In this thesis, the ship hydrodynamics during inland waterway transport and ship maneuvering are investigated using CFD (Computational Fluid Dynamics) based on OpenFoam. Validation and verification studies are carried out for the mesh convergence, time step convergence, sensitivity to turbulence models and dynamic mesh techniques. A quaternion-based 6DoF motion solver is implemented for the trim and sinkage predictions.

Environmental effects on several inland vessels (convoy 1, convoy 2, tanker) are studied using the validated numerical models. Three important aspects, the confinement effect of the waterway, head-on encounter and ship-bridge pile interaction are simulated. The testing conditions cover a wide range, including various channel dimensions, water depths, ship draughts and speeds. The ship resistance, wave pattern, Kelvin angle and wave elevation at specific positions are investigated as functions of these parameters.

Ship maneuvering is investigated using virtual captive model tests based on the MMG (Mathematical Maneuvering Group) model. An actuator disk is implemented to replace the real propeller. Open water test, rudder force test, OTT (Oblique Towing Tank test) and CMT (Circular Motion Test) of a KVLCC2 model are carried out to obtain the hydrodynamic coefficients of the propeller, rudder and ship hull. Using the obtained coefficients, system-based maneuvering simulations are carried out and validated using the free running test data.

These studies reproduce real ship tests and thus prove the validity of our numerical models. As a result, the numerical solver is promising in ship hydrodynamics and marine engineering simulations.

• Dans le cas de passage sous des piles de ponts, la résistance à l'avancement augmente au fur et à mesure que la distance "bateau-pile de ponts" diminue.

II.2. Simulations de la manoeuvrabilité par couplage CFD et systèmes dynamiques

La manoeuvre des navires a également été étudiée dans cette thèse. Elle est basée sur la résolution numérique du modèle MMG (Mathematical Maneuvering Group).

XIV

Le modèle de manoeuvrabilité actuellement implémenté dans le simulateur de navigation du CEREMA (SIMNAV), est basé sur le modèle de Norrbin, lequel ne permet pas d'étudier la manoeuvrabilité pour tout type de bateau, car les coefficients hydrodynamiques (à déterminer) de coque, de gouvernail et de propulsion (hélice) sont couplés et ne sont valables que pour un seul type de navire, en l'occurrence Esso Bernicia 190,000 DWT.

Le modèle MMG que nous proposons, dans le cadre de cette thèse, permet de surmonter ces difficultés liées aux points soulevés ci-dessus, et donc d'étudier la manoeuvrabilité de tous types de navires.

La détermination des coefficients hydrodynamiques peut être faite soit par des méthodes CFD (comme c'est le cas dans cette thèse) ou soit par des mesures expérimentales en laboratoire ou in-situ.

Ainsi, dans la présente thèse et pour la détermination par exemple des coefficients hydrodynamiques de l'hélice, il a été considéré un disque de manoeuvre numérique (actuator disk) implémenté pour remplacer l'hélice réelle.

La validation et la vérification du modèle numérique MMG avec implémentation de ce disque de manoeuvre numérique ont été faites en considérant des benchmarks internationaux issus de l'ITTC (Interbational Towing Tank Conference) validés pour des navires connus de type KVLCC2. Cela concerne les tests suivants:

• tests de force du gouvernail • tests OTT (Oblique Towing Tank)

• tests CMT (Circular Motion Test) XV

This thesis is sponsored by the China Scholarship Council (CSC) and finished in the Laboratoire Roberval, UT Compiègne-Sorbonne Université under the supervisions of Prof. Abdellatif Ouahsine and Philippe Sergent. I would like to sincerely thank Prof. Ouahsine for his encourages and helps both in my work and life. The support from Mr. Sergent is also greatly acknowledged.

Most of my simulations were done using the HPC of the Université de Strasbourg, with the help from Prof. Yannick Hoarau. We had some significant discussions on the overset method, which provided important clues for my work. The experiments in this thesis were carried out in the University of Liège. Prof. André Hage is acknowledged for providing the experimental details and discussing the maneuvering simulation methods. Dr. Shanggui Cai (Université de Technologie de Compiègne) is acknowledged for discussing the numerical methods in CFD. Dr. Shengcheng Ji (Beijing Aeronautical Science & Technology Research Institute of COMAC) is acknowledged for sharing his valuable opinions in ship simulations. Dr. Tran Khanh Toan of the Vietnam Maritime University is thanked for his previous works in the maneuvering simulations, which lays good foundations for my work. Prof. Marc Vantorre, Katrien Eloot and Dr. Evert Lataire of the Ghent University, Prof. Yasukawa of the Hiroshima University, the researchers in SIMMAN and Tokyo workshop, are thanked for sharing their experimental data and papers. Prof. Thor Fossen is thanked for sharing his experiences in maneuvering simulation. Special thanks go to Dr. Zhirong Shen (American Bureau of Shipping) for his guidance in the using and programming in OpenFoam.

I would love to dedicate my sincere acknowledgments for my wife for always accompanying me regardless of the happiness and sorrow in my life. We had our first baby girl here, which is the best thing ever. I would also thank my parents, who are the strongest support behind me.

Résumé en Français

Les travaux menés dans cette thèse concernent l'hydrodynamique et la résistance à l'avancement des bateaux en milieux confinés.

• Dans un premier temps, les vagues d'étrave et la structure de l'écoulement autour du bateau ont été simulés numériquement par CFD en tenant compte de la forme de la coque, de la vitesse du bateau et de la géométrie de la voie navigable.

• Ensuite, des scénarios de navigation en milieux confinés, qui concernent la résistance à l'avancement, le croisement et le passage sous les piles de ponts, ont été simulés et analysés également en fonction des caractéristiques cinématiques du bateau (vitesse, tirant d'eau) et géométriques (dimension de la voie navigable et des bateaux).

• Enfin, la manoeuvrabilité des bateaux en milieux confinés a été étudiée en se basant sur une procédure de couplage entre les méthodes CFD et la méthode des systèmes dynamiques. En effet la résolution directe des systèmes dynamiques requiert la détermination préalable de nombreux coefficients hydrodynamiques (∼43 coefficients), lesquels, seront déterminés par la méthode CFD dans la présente thèse. Ces coefficients hydrodynamiques sont liés à la forme de la coque, au gouvernail, à l'hélice, etc. Dans les présentes investigations nous avons adoptés le modèle MMG (Mathematical Maneuvering Group), qui, contrairement au modèle classique de Norrbin actuellement implémenté dans le simulateur SIMNAV du CEREMA-UTC, permet de simuler tout type de bateau.

I. Validation des résultats numériques

La validation des résultats obtenus est effectuée en les comparants aux mesures expérimentales, après avoir testé et vérifié le choix du modèle de turbulence et la nature du maillage (dynamique ou statique) et sa précision. La comparaison entre les résultats numériques et expérimentaux est faite:

• Tout d'abord avec un benchmark international de l'ITTC, qui consiste à comparer la hauteur des vagues et de résistance à l'avancement calculé avec des mesures expérimentales en laboratoire pour des bateaux connus KCS et KVLCC2.

• Ensuite, avec les mesures en laboratoire de l'Université de Liège (Laboratoire ANAST), qui consiste à comparer la résistance à l'avancement calculée et mesurée en fonction du tirant d'eau, de la profondeur, de la géométrie de la voie navigable et de la vitesse du bateau.

Les résultats obtenus montrent que :

• Le choix du modèle de turbulence utilisé est important pour mieux estimer la résistance à l'avancement. En effet le modèle k-ω SST reproduit avec précision les mesures expérimentales car il utilise le modèle k-ω proche de la paroi fluidecoque et le modèle k-ε loin de la paroi.

II. Présentation des résultats et discussions

II.1. Résultats des simulations CFD

Les résultats des effets du confinement sont présentés au Chapitre 5. Ils montrent que:

• La résistance à l'avancement augmente au fur à mesure que le tirant d'eau augmente, mais diminue lorsque la largeur de la voie navigable augmente.

• La résistance à l'avancement augmente lorsque la vitesse du bateau augmente.

• La résistance à l'avancement est plus importante dans le cas d'un pousseur et deux barges (convoi 2) que dans le cas d'un pousseur avec une seule barge (convoi 1).

• Dans le cas de croisement entre automoteur et un convoi 2 la résistance à l'avancement du convoi 2 augmente à cause des conditions de confinement et également à cause des vagues d'étraves générées entre les deux navires. De l'historique d'évolution de la résistance à l'avancement, on constate que, lors du processus de croisement entre les deux navires, cette résistance augmente tout d'abord, puis diminue juste après la position d'alignement des deux bateaux. Cette diminution de résistance à l'avancement est liée à l'abaissement du plan d'eau (drawdown).

Comparison between (a) AMI (Arbitrary Mesh Interface) and (b) ACMI

(Arbitrary Coupled Mesh Interface) methods. (c) shows the discretization error of the meshes at the interface. . . . . . . . . . . . . . . . .

3.4

One-dimensional domain of the overset mesh, with two-subdomains A and B. Values at nodes 6 and 7 are unknown. The arrows denote the interpolation from the donor nodes/cells to the fringe nodes/cells [START_REF] Ma | An overset mesh based multiphase flow solver for water entry problems[END_REF]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Hole cutting and DCI (Domain Connectivity Information

) of two overlapping meshes in a one-dimensional domain. A is the background mesh. B is the body-fitted mesh for a structure. The arrows denote the interpolation from the donor nodes/cells to the fringe nodes/cells. The crosses represent the hole nodes/cells [START_REF] Ma | An overset mesh based multiphase flow solver for water entry problems[END_REF]. . . . . . . . .

3.6

Flowchart of the computational procedure with overset meshes. "DCI" is the domain connectivity information. . . . . . . . . . . . . . . . . [START_REF] Fossen | Guidance and control of ocean vehicles[END_REF] Horizontal reach of bank effect y + Non-dimensional wall distance (u, v, w, p, q, r)

Linear and angular velocities of the ship (x, y, z, φ, θ, ψ)

Translational and rotational states of the ship (x g , y g , z g )

Position vector from gravity center to rotation center

(x G , y G , z G ) Coordinate of gravity center A R Profile area of rudder movable part A s , A c
Cross section areas of ship and waterway

B R Rudder span B s , B c
Widths of ship and waterway During inland shipping transport, various situations can be confronted. Ships may enter narrow and shallow waters, be influenced by wind, waves, currents, meet or overtake another ship, etc. These conditions exert additional effects on ships and change their maneuverability. Correctly establishing the corresponding numerical models based on the real conditions and reproducing the physics are quite challenging and significant to provide suggestions for inland shipping.

B wl Maximum beam of waterline C b Block coefficient C ext Extrapolated resistance coefficient C F Frictional coefficient C P Pressure coefficient C R Rudder chord length C T Total resistance coefficient C 1 , C
Ship maneuvering simulation is able to give the trajectory and motion states. Many maneuvering models exist currently. Understanding their formulations and procedures is beneficial to select the most appropriate one for applications. Obtaining the hydrodynamic coefficients is obligatory, normally using captive model tests, which are complicated and costly. Their virtual counterparts based on CFD (Computational Fluid Dynamics) have been proven feasible and efficient. The procedure from virtual captive model tests to the hydrodynamic coefficients, and finally to the maneuvering simulations needs to be clarified and understood.

CFD has been widely-used for ship-related simulations. Potential flow methods, such as the strip theory, panel method, etc. have been the state-of-the-art for such problems in the past years. RANS-based (Reynolds-Averaged Navier Stokes) solvers have thrived in recent years because of its capability of reproducing real fluid flows and acceptable computational costs. Except the fluid solver, other models, including turbulence models, dynamic mesh techniques, multiphase models, etc. are also necessary for various purposes. Many commercial and open source codes emerge nowadays. However, it is not reliable to carry out simulations without understanding the inside numerics and procedures. This is also a prerequisite if one encounters problems and fixes them during the simulations.

Objective of this thesis

This thesis aims to investigate the phenomena occurred during inland shipping transport and ship maneuvering based on CFD. The most commonly-used and advanced numerical models will be explained and validated comprehensively for ship simulations. The most important aspects include:

• The fully-confinement effect of the inland waterway will be investigated. The influences of the channel dimensions, water depths, ship draughts on the advancing resistance and ship-generated waves will be analyzed and summarized;

• The head-on encounter process between two vessels in the confined waterway will be simulated. The resistance and wave profiles during the crossing will be characterized;

• The ship-bridge pile interaction will be studied. A parameter will be defined for characterizing the resistance and wave pattern during the passing process and compared with experiments.

• Virtual captive model tests will be carried out to calculate the hydrodynamic coefficients in the MMG (Mathematical Maneuvering Group) model. An actuator disk model will be implemented to consider the propulsion. Maneuvering simulations will be conducted with these coefficients.

Outline of this thesis

This thesis is organized as follows:

-Chapter 2 presents the start-of-the-art literature review of the numerical methods and interesting topics for ship hydrodynamics;

-Chapter 3 is devoted to the framework of the numerical solver, including the multiphase models, turbulence models, finite volume discretization and solid body motion approaches;

-Chapter 4 carries out the validation and verification studies for ship simulations, specifically the boundary condition analyses, mesh and time step convergence, sensitivity to turbulence models and dynamic meshes;

-Chapter 5 simulates the ship motions in the confined waterway. The confine-ment effect, head-on encounter and ship passing bridge piles are investigated.

-Chapter 6 carries out virtual captive model tests for the MMG model to obtain the hydrodynamic coefficients. System-based maneuvering simulations are then simulated based on these results.

-Chapter 7 concludes this thesis and gives suggestions for the future work.

Chapter 2

State of the art review

Introduction

Ship simulations based on CFD (Computational Fluid Dynamics) are also referred as numerical virtual tank technique. This chapter is dedicated to a state-of-theart literature review of the relevant numerical methods and ship hydrodynamics investigations.

• The concept of the virtual tank simulations is introduced.

• The most commonly-used numerical techniques are presented.

• The interesting topics for ship hydrodynamics study are summarized, covering the most common and advanced models used nowadays.

Concept of the numerical virtual tank

To obtain the flow and force characteristics of a ship, tests have to be carried out in the basin. One commonly-used facility is the towing tank. As shown in Figure 2.1, the model is connected to a carriage to be towed (or propelled), with the allowance of predefined degrees of freedom (DoF) (e.g. heave, pitch, roll). The model has to satisfy the geometric and dynamic similarities. The geometric similarity can be satisfied as:

λ s = L f s L m s λ 2 s = S f s S m s λ 3 s = ∇ f s ∇ m s (2.1)
where λ s is the scale factor. S s and ∇ s are the surface area and volume underwater. L s is the ship length. The superscripts 'f ' and 'm' denote full-or model-scale.

Dynamic similarity means that the velocities, accelerations and forces of the model and full-scale ships have scaled magnitudes and identical directions at corresponding locations, which is physically impossible. Instead, the Froude similarity is used, i.e. the Froude numbers are kept the same:

V f L f s = V m L m s (2.2)
where V is the ship velocity. The scale effect is then partially compensated by empirical corrections. Besides, the model size should be appropriate to minimize the viscous scale effects within the capability of the towing carriage and the restrictions of the basin. Motion trackers and loadcells are installed to track the movements and to measure the forces and moments acted on the model. Probes and visualization techniques (e.g. LDV, PIV) are used to capture the wave elevations and flow fields (Du et al., 2017;[START_REF] Fossen | Handbook of marine craft hydrodynamics and motion control[END_REF]. Seakeeping studies can also be carried out by using wavemakers. Except the towing tank, the ones with special functionalities have also been built for the maneuvering study, cavitation, shallow water condition, etc. A complete list of them can be found on the website of the ITTC (International Towing Tank Conference), which is responsible of recommending procedures and guidelines for hydrodynamic experiments and simulations. While towing tank experiments have provided considerable results and suggestions [START_REF] Jagadeesh | Experimental investigation of hydrodynamic force coefficients over AUV hull form[END_REF][START_REF] Longo | Uncertainty assessment for towing tank tests with example for surface combatant DTMB model 5415[END_REF][START_REF] Mctaggart | Seakeeping of two ships in close proximity[END_REF][START_REF] Scarano | Traversing field of view and AR-PIV for mid-field wake vortex investigation in a towing tank[END_REF], the difficulties with the background noise, carriage vibration and limited running time associated with the finite length of the channel are also well-known [START_REF] Gad-El-Hak | The water towing tank as an experimental facility[END_REF]. CFD is becoming an important tool for simulating and understanding the shipgenerated flows. Compared with the experimental approaches, the numerical virtual tank has the advantages of easier extension, less noise, no scale effect, lower cost, etc. This idea is based on using the CFD theories to mimic typical towing tanks, as evidenced by the Numerical Towing Tank Symposium (NuTTS), ship hydrodynamics CFD workshops and Office of Naval Research (ONR) Symposiums on Naval Hydrodynamics (SNH). According to the purposes, virtual tanks can be classified into several types as shown in Figure 2.2, where the fluid motions are normally governed by the Navier-Stokes (NS) equations:

∂ρ ∂t + ∇ • (ρu) = 0 (2.3) ∂(ρu) ∂t + ∇ • (ρuu) = -∇(p + 2 3 µ∇ • u) + ∇ • [µ(∇u + (∇u) T )] + S u (2.4)
where u, p, ρ and µ are the fluid velocity, pressure, density and viscosity respectively. S u is the source term. They can be further simplified to inviscid or irrotational versions (Euler equations or potential theory) if necessary. Ocean waves, dynamic mesh techniques, turbulence and multiphase models can be added for various purposes. The CFD computation procedure is depicted in Figure 2.3. Modeling and numerical errors will occur during the discretization and computation process. Therefore it is obligatory to carry out grid independence study to ensure the convergence of the solution and the order of accuracy. The final results should also be validated using model test results. • Smoothed Particle Hydrodynamics

Numerical methods for the virtual tank

In Table 2.1, meshfree methods are a special class using the concept of particle collision and streaming. The connection between particles is not required, i.e. no mesh is needed for these methods. Among them, the Smoothed Particle Hydrodynamics (SPH) has been used for the sloshing tank [START_REF] Iglesias | Simulation of anti-roll tanks and sloshing type problems with smoothed particle hydrodynamics[END_REF], wave-structure interaction [START_REF] Altomare | Applicability of Smoothed Particle Hydrodynamics for estimation of sea wave impact on coastal structures[END_REF] and water entry problems [START_REF] Veen | A combined strip theory and Smoothed Particle Hydrodynamics approach for estimating slamming loads on a ship in head seas[END_REF]. Several good codes are the SPH-flow from École centrale de Nantes [START_REF] Oger | Two-dimensional SPH simulations of wedge water entries[END_REF], MLParticle-SJTU from Shanghai Jiao Tong University [START_REF] Zhang | Numerical study of interactions between waves and free rolling body by IMPS method[END_REF] and the open-source code SPHysics [START_REF] Altomare | Applicability of Smoothed Particle Hydrodynamics for estimation of sea wave impact on coastal structures[END_REF], etc. Meshfree methods are normally more time-consuming than their mesh-based counterparts but with advantages of easier representation of complex geometries and the allowance of very large deformations and displacements.

• Lattice Boltzmann Method Another special type of methods is based on the Boltzmann equation. They are generally faster and simpler than traditional CFD methods. The most representative one, the Lattice Boltzmann Method (LBM) uses a similar concept as the meshfree methods but with the collision and streaming happening on a discrete lattice mesh.

Particle distribution function is used to calculate the macroscopic variables (Du et al., 2017). This non-primitive nature also brings difficulties when implementing the velocity and pressure conditions. Since they are relatively new, there are still a lot of work to improve their performances and expand their applications. LBM has been used for free surface flow [START_REF] Janssen | A lattice Boltzmann approach for freesurface-flow simulations on non-uniform block-structured grids[END_REF] and wave [START_REF] Zhao | A lattice Boltzmann method for viscous free surface waves in two dimensions[END_REF] simulations. More complex cases are the water entry problem [START_REF] Zhang | Numerical simulation of the water-entry of body based on the Lattice Boltzmann method[END_REF], airwakes of structures [START_REF] Syms | Simulation of simplified-frigate airwakes using a lattice-Boltzmann method[END_REF] and ship-ice interactions [START_REF] Mierke | GPU-accelerated Large-Eddy Simulation of ship-ice interactions[END_REF]. However, there are still difficulties for high-Mach number flows and mesh generation of complex geometries.

• Boundary Element Method BEM is a method for solving the integral form of partial differential equations in a boundary value problem. Rather than the values throughout the whole domain, only the given boundary values are fitted into the integral equation. Thus the discretization is needed only at the boundary. The model dimension is thereby reduced by one, which dramatically increases the computational efficiency. Currently it is always used for solving the potential flow theory [START_REF] Grilli | Development of a 3D numerical wave tank for modeling tsunami generation by underwater landslides[END_REF][START_REF] Wang | An efficient numerical tank for non-linear water waves, based on the multi-subdomain approach with BEM[END_REF]. It takes the flow as inviscid and irrotational, which greatly simplifies the NS equations. However, it can only be derived for certain classes of partial differential equations. The strip theory and panel method based on 2D and 3D potential flow theory using BEM are especially useful for ship simulations because of its simplicity, high efficiency and satisfactory results. This method is currently the most widely-used in marine engineering.

(1) Strip theory

The strip theory is based on the 2D potential flow theory. It applies to slender bodies whose cross sections underwater are small compared to the overall length. It divides the submerged part into a finite number of strips. Then potential flow calculations can be done in 2D for each strip using BEM. The values are summed up to yield the 3D results. The strip theory can be used to calculate the hydrodynamic added mass, potential damping, wave loads, etc. in a seakeeping problem [START_REF] Fossen | Nonlinear time-domain strip theory formulation for low-speed manoeuvring and station-keeping[END_REF][START_REF] Xia | Non-linear wave loads and ship responses by a time-domain strip theory[END_REF]. Commonly-used commercial software are ShipX-VERES by MARINTEK and Octopus Office by Amarcon.

(2) Panel method Panel method divides the surface of marine structures or ships into discrete elements (panels). On these panels, sources and sinks (or similarly, dipole or vortex distributions) are defined which fulfill the Laplace equation. When the simple Rankine source (dipole Green function) is selected, the panel method is also called Rankine panel method or Rankine source method. This allows the application of BEM to solve a 3D potential flow problem. The velocity potentials can be calculated at last. An outstanding commercial program is WAMIT which computes the frequency-dependent added mass, potential damping coefficients, wave load transfer functions, etc. More details of the panel method can refer to [START_REF] Larsson | The principles of naval architecture series: ship resistance and flow[END_REF].

For most problems related to ship motions, potential theory is sufficient to obtain results with appropriate accuracy. Nevertheless, flows involving breaking waves and splashes can hardly be analyzed properly by potential methods. Besides, viscosity becomes important when predicting the stress and turbulent structures near the wall, such as ship maneuvering and propeller-rudder-hull interactions. It is necessary to get these details for better understanding the inner physics.

• Finite Difference Method FDM is one of the oldest numerical methods for solving partial differential equations and is intuitive for understanding. The partial derivatives in the NS equations 3.1 and 3.2 are approximated for each term based on Taylor series. To avoid the checker-board effect, staggered mesh or the Rhie-Chow interpolation on a collocated mesh should be adopted. After the discretization, a system of algebraic equations on each node will be obtained, which can be solved after introducing the boundary conditions. The final solution is an approximate one with an error between the partial derivatives and finite differences. FDM is a differential approach and we get the strong form solution at last. Many works have been done for the numerical wave tank [START_REF] Park | Fully non-linear free-surface simulations by a 3D viscous numerical wave tank[END_REF] and the sloshing phenomenon [START_REF] Chen | Time-independent finite difference analysis of fully non-linear and viscous fluid sloshing in a rectangular tank[END_REF]Kim et al., 2007) The NTNU code, called REEF3D [START_REF] Bihs | A new level set numerical wave tank with improved density interpolation for complex wave hydrodynamics[END_REF], solves the NS equations on a staggered Cartesian grid, with a ghost cell Immersed Boundary Method (IBM) for the irregular cells. The CFDShip-Iowa and SJTU code [START_REF] Wang | An efficient numerical tank for non-linear water waves, based on the multi-subdomain approach with BEM[END_REF] instead use the overset grid technique to keep the grid orthogonality.

• Finite Volume Method FVM is a natural choice for CFD problems since the equations in CFD are conservation laws, i.e. the flux entering a finite volume from one side is identical to that leaving on another side. FVM can be easily formulated to allow for unstructured meshes. As shown in Table 2.1 and Figure 2.4, there are normally two approaches for the finite volume mesh based on the domain discretization approaches: the Cell-Centered FVM (CC-FVM) and Vertex-Centered FVM (VC-FVM).

-CC-FVM means that the Control Volume (CV) is formed by the cell of the mesh (Figure 2.4(a)), with the variables storing at the cell center. The face fluxes are approximated using the values in the two adjacent cells. This approach is efficient and the resulted matrices have a low band width (the bank width is equal to the number of cell neighbors plus one), which will result in a low storage. However, the mesh topology is restricted due to the orthogonality and conjunctionality requirements. On an unstructured mesh the face center may not lie on the line joining the cell centers, which will lead to an error in the face interpolation.

-In the VC-FVM, the cells are divided and the CV is constructed within the cells rather than using the cell faces as in the CC-FVM. The variable values are stored at the vertexes. The face fluxes are calculated using the vertex values through interpolations. The VC-FVM is less efficient and has a larger matrix band width than the CC-FVM. However, the mesh topology does not have the same restrictions as the CC-FVM. Most CFD software currently uses the CC-FVM approach since its implementation in the code is straightforward. The most famous programs for marine applications are ANSYS Fluent, STAR-CCM+, FINE/Marine from École centrale de Nantes, PARNASSOS and ReFRESCO from MARIN (Maritime Research Institute Netherlands). Some widely-used open source codes are OpenFoam and SU2. The work in this thesis is based on OpenFoam. [START_REF] Moukalled | 4 Two domain discretization approaches for the FVM mesh[END_REF]. CV means the Control Volume.

Overview of ship hydrodynamics investigations 2.4.1 Ship resistance

Resistance prediction is perhaps the most general application of CFD for ships. The resistance is usually decomposed into various components, although these components may interact with each other and most of them cannot be measured individually. From a microscopic point of view, the force on the ship has a normal and tangen-tial components, which forms the pressure and frictional forces macroscopically (Figure 2.5). Since waves are generated by the ship motion, the energy will be partly taken by these waves (wave resistance), the rest will dissipate in the wake (viscous resistance). A typical resistance curve can be found in Figure 2.6. The viscous resistance is affected by the viscosity, ship velocity and wetted surface area. It increases smoothly with the ship speed. Wave resistance is caused by the shipgenerated waves. It is affected by the hull shape, beam to length ratio, ship speed, etc. The physics of wave resistance is more complicated, and small hollows and humps can be observed in the resistance curves. The position of humps is function of vessel shape and speed. The viscous component plays a more important role most times, however wave resistance can increase dramatically during motion changes (unsteady states) and may exceed the viscous resistance in extreme conditions. [START_REF] Bertram | Practical ship hydrodynamics[END_REF][START_REF] Molland | Ship resistance and propulsion[END_REF]. According to the source, the resistance can be generated by the hull, appendages and other environmental effects. Hull resistance can be predicted with high precision in the current stage. Appendages pose problems for CFD since some of them need dynamic mesh and fluid-structure interaction techniques, and introduce additional mesh refinement, which increase the numerical complexity and computational cost. Confinement effect originates from the limitations of channel side and bottom, return current will be generated and added mass effect can cause numerical instabilities (Section 2.4.3). Wave and current effects are normally known as ship seakeeping. In CFD, waves can be generated by the time-varying of velocity and phase BCs, dynamic BCs and adding source terms (Section 2.4.4). Ship-ship, shipstructure interactions etc. are all classified as ship interactions in this thesis and are introduced in Section 2.4.5.

Ship-generated waves

A moving object at the free surface will generate waves, which can be decomposed into primary and secondary waves [START_REF] Bertram | Practical ship hydrodynamics[END_REF]. The primary wave follows the Bernoulli's equation (Figure 2.7). The flows at the bow and stern decelerate and create high pressure regions. The surface elevation corresponds with the pressure distribution, yielding wave crests at the ship ends and a long wave trough along the middle. This shape is speed independent. The wave height depends quadratically on the speed. The secondary wave pattern, also known as the Kelvin wave pattern, consists of transverse and divergent waves. Early analysis is conducted for a moving pressure source as shown in Figure 2.8(a). The waves are limited to a wedge-shaped region with an angle [START_REF] Soomere | Applied wave mathematics[END_REF]:

sin α k = 1 + 2k w h w sinh -1 (2k w h w ) 3 -2k w h w sinh -1 (2k w h w ) (2.5)
where α k is the Kelvin wave angle, k w is the wave number, h w is the water depth. In deep water (h w → ∞), sin α k → 1/3, then α k ≈ 19.47 • . This angle is independent of the hull shape. The actual ship waves is formed by the interference of these simple elements. The wave system interference of a typical hull form is shown in Figure 2.8(b). Two waves starting with a crest are formed at the bow and stern, while other waves starting with a trough are formed at the middle parts. The sum of these wave profiles, together with the primary wave pattern (Figure 2.7) forms the actual shipgenerated waves. In confined areas, the waves will be modified because the water depth h w comes into effect in 2.5, which will be discussed in Section 2.4.3. [START_REF] Molland | Ship resistance and propulsion[END_REF]. V and λ denote the velocity and the wavelength.

Confinement effect

When a vessel is navigated in waters with limited width and depth, it experiences resistance increase and drops of speed and propulsion efficiency compared with open water, leading to higher propulsion power and fuel consumption. The maneuvering also becomes harder, leading to larger turning radius. The distance between the keel and the channel bottom decreases known as "squat" (Figure 2.9), a combination of trim and sinkage. The flow around the vessel is altered, causing a return current and depression of water level. Important thresholds of the confinement effects are shown in Table 2.2. The condition when h w /T d < 1.5 is called shallow water. Three flow regimes can be identified referring the depth-based Froude number (F r h ), subcritical, critical (F r h = 1) and supercritical (Figure 2.10). The phenomena of the three regimes are different. In the subcritical condition, the ship advancement pushes away the water at the bow and leaves a void at the stern, creating a high and a low pressure region respectively. The return current under the ship decreases the pressure, resulting in a downward force on the hull. And the so-called squat appears. The generated ship waves are able to overtake the downstream flow and travel in an upstream direction. On the contrary, in the supercritical condition, the vessel will rise and the surface waves cannot travel upstream. However, it has been confirmed that there exists a limiting speed that self-propelled vessels cannot exceed in restricted waterways regardless of the power applied, known as the Schijf limiting speed [START_REF] Schijf | Influence on the form and dimensions of the cross-section of the canal, of the form, of the speed and the propulsion system of vessels[END_REF].

Resistancce

The resistance in open, shallow and confined waters are different, whose typical curves can be observed in Figure 2.10. For the case of shallow water, no lateral restriction exists, only the undulotary effect caused by the bottom takes place. The key parameter is h w /T d (Table 2.2). There exists a critical speed (F r h = 1) where the ship resistance increases dramatically. Around the critical speed the ship motions are unsteady. The range 0.90 < F r h < 1.1 is also defined as the transcritical region.

When in confined water where the width is also restricted, the hydraulic effect appears besides the undulotary effect. The ratio A c /A s (Table 2.2) dominates the resistance change with a subcritical (2.6) and a supercritical speed (2.7) [START_REF] Pompée | About modelling inland vessels resistance and propulsion and interaction vessel-waterway key parameters driving restricted/shallow water effects[END_REF][START_REF] Constantine | On the movement of ships in restricted waterways[END_REF]:

Subcritical F r sub h = 2sin
Arcsin(1 -m b ) 3

1.5 (2.6) Supercritical F r super h = 2sin π -Arcsin(1 -m b ) 3
1.5

(2.7) where m b = A s /A c is the blockage ratio. The weight distribution of the confinement effect can be represented by a decreasing exponential function (Figure 2.11) [START_REF] Lataire | Longitudinally directed bank effects[END_REF]:

ω w = e - ξy |y| y inf l + ξz |z| T d (2.8)
where ξ y and ξ z are determined from experiments. y inf l = B s (5F r h + 5) is the horizontal reach of the bank effect [START_REF] Lataire | Ship-bank interaction induced by irregular bank geometries[END_REF]. T d and B s are the ship draught and beam. As shown in Figure 2.11, a water particle close to the hull will have a high value, meaning high confinement effect. Ships can experience a speed drop of 30% in shallow water and up to 60% in fully-confined area (Tezdogan et al., 2016;[START_REF] Barrass | Ship stability for masters and mates[END_REF].

Squat

Squat is caused by the reduction in Under Keel Clearance (UKC) between a vessel at rest and underway due to the accelerated flow around the moving body. It greatly influences the ship resistance and maneuverability (Figure 2.9). [START_REF] Gourlay | A brief history of mathematical ship-squat prediction, focussing on the contributions of EO Tuck[END_REF][START_REF] Gourlay | Mathematical and computational techniques for predicting the squat of ships[END_REF]; [START_REF] Millward | A review of the prediction of squat in shallow water[END_REF]; [START_REF] Varyani | Squat effects on high speed craft in restricted waterways[END_REF]. However, they can only give estimates for limited conditions. Comparably, numerical simulations are more accurate. Potential flow simulations can be found in [START_REF] Gourlay | Application of potential flow methods to ship squat in different canal widths[END_REF]; [START_REF] Gourlay | Ship squat in water of varying depth[END_REF][START_REF] Gourlay | Slender-body methods for predicting ship squat[END_REF][START_REF] Gourlay | Mathematical and computational techniques for predicting the squat of ships[END_REF]; Yuan & Incecik (2016) (2016). However, numerical simulations have difficulties for very confined conditions because of the added mass effect (Förster et al., 2007;[START_REF] Devolder | A review of the implicit motion solver algorithm in OpenFOAM to simulate a heaving buoy[END_REF][START_REF] Dunbar | Development and validation of a tightly coupled CFD/6-DOF solver for simulating floating offshore wind turbine platforms[END_REF], which leads to the divergence of the simulation. Detailed discussions on the added mass effect can be found in Section 3.4.3.

Ship waves

Ship waves in shallow and fully-confined water areas are different, but tightly connected. The influence of water depth dominates the generated wave pattern as function of the depth-Froude number (Figure 2.12). In shallow water, traditional Kelvin wave patterns can be observed at the subcritical region. As the ship speed approaches the critical speed, the wave angle approaches 90 • , i.e. the ship waves propagate parallel to the sailing line. At the supercritical region, the transverse waves are left behind and only divergent waves are present under this condition [START_REF] Pethiyagoda | Time-frequency analysis of ship wave patterns in shallow water: modelling and experiments[END_REF]. The bow waves folds backwards to the ship with respect to the Froude number as:

α k = arcsin(1/F r h ) (2.9)
The wave angle changing with the Froude number is shown in Figure 2.12. The blue line follows 2.5 and the brown line follows 2.9. In a fully-confined waterway, the wave angle follows that in the shallow water. However, the waves will be reflected at the side wall (Figure 2.13).

Figure 2.12: Dependence of the wave angle (α k ) on the depth-Froude number (F r h ) [START_REF] Soomere | Applied wave mathematics[END_REF][START_REF] Pethiyagoda | Time-frequency analysis of ship wave patterns in shallow water: modelling and experiments[END_REF].

Ship seakeeping

Seakeeping performance of a ship refers to the ability to remain and carry out its duty in all conditions, such as in wind, waves and current. Experimental method is still the most reliable. The usage of CFD keeps growing and has become a powerful alternative, where creating these conditions (wind, wave, current) accurately is a prerequisite. For maneuvering analysis, flow details are not necessary. Only the forces and moments on the ships are formulated empirically, theoretically or numerically.

Wind

Early studies directly estimate the wind forces and moments (F wind , M wind ) as functions of the air density (ρ a ), reference speed (V r ) and projected area (A T , A L ) as follows:

F wind = 1 2 C F w ρ a V 2 r A T M wind = 1 2 C M w ρ a V 2 r A L L s (2.10)
with different coefficients (C F w , C M w ) obtained from model tests for each component. Several models for merchant ships, larger tankers, medium-sized ships, VLCCs (Very Large Crude Carriers), moored ships and floating structures can be found in [START_REF] Fossen | Guidance and control of ocean vehicles[END_REF][START_REF] Fossen | Handbook of marine craft hydrodynamics and motion control[END_REF]. Note that a different convention from 2.10 is used by the OCIMF (1994) (Oil Companies International Marine Forum) but the overall idea is similar. These empirical models are normally used for maneuvering analysis. As a more accurate approach, CFD has been principally used to investigate the wind loads [START_REF] Koop | Predicting wind loads for FPSO tandem offloading using CFD[END_REF][START_REF] Wnęk | CFD assessment of the wind loads on an LNG carrier and floating platform models[END_REF][START_REF] Mousaviraad | CFD prediction of ship response to extreme winds and/or waves[END_REF] and air wakes of the superstructure [START_REF] Snyder | USNA ship air wake program overview. 29th AIAA Applied Aerodynamics Conference[END_REF][START_REF] Thornber | Implicit large eddy simulation of ship airwakes[END_REF][START_REF] Reddy | Numerical simulation of ship airwake[END_REF]Yesilel & Edis, 207). LBM is also used to simulate the airwakes of simplified ship models in recent years [START_REF] Syms | Simulation of simplified-frigate airwakes using a lattice-Boltzmann method[END_REF]. Wind is normally generated at the inlet uniformly [START_REF] Wnęk | CFD assessment of the wind loads on an LNG carrier and floating platform models[END_REF] or with atmospheric boundary layer [START_REF] Koop | Predicting wind loads for FPSO tandem offloading using CFD[END_REF]. Experimental works of the wind effect can be found in [START_REF] Mousaviraad | CFD prediction of ship response to extreme winds and/or waves[END_REF], [START_REF] Snyder | USNA ship air wake program overview. 29th AIAA Applied Aerodynamics Conference[END_REF] and [START_REF] Wnęk | CFD assessment of the wind loads on an LNG carrier and floating platform models[END_REF]. Since the air density is too small compared with water, the wind effect is always not a hot spot for ships.

Wave

Modern wave-makers (piston-or flap-type) are able to create harmonic waves. The superposition of sinusoidal waves of different frequency can create irregular waves similar to natural wind seas. Several experimental data can be found for the DTMB 5415/5512, KVLCC2, Wigley, DTC (Duisburg Test Case) and S175 in regular head waves [START_REF] Hino | CFD workshop Tokyo 2005 workshop on numerical ship hydrodynamics[END_REF][START_REF] Gui | Forces, moment, and wave pattern for surface combatant in regular head waves Part I. Measurement systems and uncertainty assessment[END_REF][START_REF] Gui | Forces, moment and wave pattern for surface combatant in regular head waves Part II. Measurement results and discussions[END_REF][START_REF] Lee | EFD and CFD for KCS heaving and pitching in regular head waves[END_REF][START_REF] Journée | Experiments and calculations on 4 Wigley hull forms in head waves[END_REF][START_REF] Sigmund | Numerical and experimental investigation of added resistance of different ship types in short and long waves[END_REF][START_REF] Fujii | Experimental study on the resistance increase of a ship in regular oblique waves[END_REF][START_REF] Fonseca | Experimental investigation of the nonlinear effects on the vertical motions and loads of a containership in regular waves[END_REF] , KCS (KRISO Container Ship) and ONRT (ONR Tumblehome Ship) in regular head and oblique waves [START_REF] Larsson | The principles of naval architecture series: ship resistance and flow[END_REF][START_REF] Carrica | Computations of self-propulsion free to sink and trim and of motions in head waves of the KRISO Container Ship (KCS) model[END_REF][START_REF] Simonsen | EFD and CFD for KCS heaving and pitching in regular head waves[END_REF]. Experiments in irregular waves can refer to [START_REF] Nakamura | Propulsive performance of a container ship in waves[END_REF], [START_REF] Umeda | Model experiments of ship capsize in astern seas[END_REF], [START_REF] Bulian | On the nonlinear modeling of parametric rolling in regular and irregular waves[END_REF] and [START_REF] Lee | Theoretical and experimental study on dynamic behavior of a damaged ship in waves[END_REF].

For system-based maneuvering analysis, formulas for calculating the wave forces should be defined, using either the two-time scale model or the hybrid approach (ITTC, 2011;[START_REF] Tello Ruiz | Ship Manoeuvring in Waves: A literature review[END_REF]. Both methods are based on the general assumption of a possible separation of fluid effects into a contribution of the conventional maneuvering in calm water and the wave-induce motions. However, their evaluations are completely different. The two-time scale methods separate the basic ship motion into two parts: the one for high frequency wave-induced motion and the other for the low frequency maneuvering motion, based on the fact that the wave-induced motion is generally much faster than the maneuvering motion. The motion separation brings a favorable merit that the second-order wave forces can be evaluated in a more accurate manner. Several works falling into this category can be found in [START_REF] Skejic | A unified seakeeping and maneuvering analysis of ships in regular waves[END_REF], [START_REF] Yasukawa | Simulation of ship maneuvering in waves (1st report: turning motion)[END_REF], [START_REF] Seo | Numerical analysis on ship maneuvering coupled with ship motion in waves[END_REF][START_REF] Fossen | Guidance and control of ocean vehicles[END_REF] and Zhang et al. (2017). The hybrid approach, on the other hand, integrates the low frequency maneuvering motion and the high frequency wave induced motion into a generic set of equations to describe the maneuvering in waves. Several works falling into this classification are [START_REF] Bailey | A unified mathematical model describing the manoeuvring of a ship travelling in a seaway[END_REF], [START_REF] Fang | A nonlinear mathematical model for ship turning circle simulation in waves[END_REF], [START_REF] Sutulo | A unified nonlinear mathematical model for simulating ship manoeuvring and seakeeping in regular waves[END_REF], [START_REF] Lin | Numerical simulations of ship maneuvering in waves[END_REF] and [START_REF] Subramanian | A time-domain strip theory approach to maneuvering in a seaway[END_REF]. Among these works, [START_REF] Yasukawa | Simulation of ship maneuvering in waves (1st report: turning motion)[END_REF] and [START_REF] Fossen | Guidance and control of ocean vehicles[END_REF] were dedicated to irregular waves; [START_REF] Seo | Numerical analysis on ship maneuvering coupled with ship motion in waves[END_REF], [START_REF] Skejic | A unified seakeeping and maneuvering analysis of ships in regular waves[END_REF], [START_REF] Subramanian | A time-domain strip theory approach to maneuvering in a seaway[END_REF], [START_REF] Sutulo | A unified nonlinear mathematical model for simulating ship manoeuvring and seakeeping in regular waves[END_REF] and Zhang et al. (2017) included the second-order wave loads. The wave forces can be divided into viscous force and inertial force, where the latter consists of the diffraction, radiation and Froude-Krylov forces arising from the potential flow theory. They can also be expressed as a sum of first, second and even higher order terms. The linear first order terms are important for predicting the motion responses, while the second order ones are responsible for the drift.

In CFD, three approaches can be adopted for the wave generation, static boundary condition, moving boundary and internal generation method. ( 1) The first one simply implements the free surface level and velocity as boundary conditions. This approach is the most widely-used and can be found in many works [START_REF] Higuera | Realistic wave generation and active wave absorption for Navier-Stokes models: Application to OpenFOAM[END_REF][START_REF] Jacobsen | A wave generation toolbox for the open-source CFD library: OpenFoam[END_REF]. The applicable range of different wave theories can refer to the classic Le [START_REF] Méhauté | An introduction to hydrodynamics and water waves[END_REF] graph. Irregular waves can normally be treated as a superposition of a series of linear waves based on the wave spectrum.

Commonly-used wave spectrum is the Pierson-Moskowitz spectrum, JONSWAP spectrum, Torsethaugen spectrum and white noise [START_REF] Fossen | Handbook of marine craft hydrodynamics and motion control[END_REF][START_REF] Shen | An irregular wave generating approach based on naoe-FOAM-SJTU solver[END_REF]. ( 2) The second method using the moving boundary techniques for wave generation is similar with the experimental facilities [START_REF] Higuera | Three-dimensional numerical wave generation with moving boundaries[END_REF]. The generated wave type can be controlled by the number and frequency of the bound-ary motions.

(3) The internal wave generation is to incorporate the mass source functions into the NS equations. By using different source functions, all kinds of wave trains can be generated [START_REF] Lin | Internal wave-maker for Navier-Stokes equations models[END_REF].

Along with the wave generation, the wave absorption becomes a problem, where the wave reflections and backflows can disturb the internal flow. Several methods for this purpose are as follows: (1) increase the domain size and grid damping, thus increase the numerical discretization and iteration errors to diffuse the waves;

(2) numerical beaches, where a slope in the domain bottom is used which leads to wave breaking and energy dissipation; (3) sponge layer, where momentum sinks are included in the governing equations to damp the waves propagating through the predefined zones. These techniques are also called passive methods. They are effective but trade with the computational cost of the increased domain. The counterpart of them is the active absorption methods. The wave velocity is first identified at the outlet. Two options are then possible to absorb the waves: (1) an equal velocity with opposite direction is imposed; (2) the boundary is moved [START_REF] Perić | Reliable damping of free-surface waves in numerical simulations[END_REF][START_REF] Higuera | Realistic wave generation and active wave absorption for Navier-Stokes models: Application to OpenFOAM[END_REF]. The main advantage of this kind of absorption is that it prevents the energetic and mean water levels to increase unbounded, while adding no significant computational cost.

Current

Currents are horizontal and vertical circulation of waters produced by gravity, wind friction and water density variation. They have the characteristics of low frequency and long term on ships compared with waves. The forces and moments can be included in two ways for the maneuvering analysis: (1) modify the ship velocity using the current velocity [START_REF] Fossen | Guidance and control of ocean vehicles[END_REF][START_REF] Fossen | Handbook of marine craft hydrodynamics and motion control[END_REF], where no explicit force equations are needed for this approach; (2) propose empirical equations similar with 2.10 based on model tests [START_REF] Ocimf | Prediction of Wind and Current Loads on VLCCs[END_REF]. These models are rather simplified and the interactions with waves are neglected. Some works studying the wave-current interactions can be found in [START_REF] Umeyama | Dynamic-pressure distributions under Stokes waves with and without a current[END_REF][START_REF] Umeyama | Coupled PIV and PTV measurements of particle velocities and trajectories for surface waves following a steady current[END_REF] and [START_REF] Singh | Combined effect of wave and current in rough bed free surface flow[END_REF]. The influence of the wave-current interaction on floating structures has been studied by [START_REF] Lin | Wave-current interaction with a vertical square cylinder[END_REF], [START_REF] Kang | Wave-current interaction with a vertical square cylinder at different Reynolds numbers[END_REF] and [START_REF] Park | Three-dimensional numerical wave tank simulations on fully nonlinear wave-current-body interactions[END_REF] using CFD. A THOBEM (Time-domain Higher-Order Boundary Element Method) is also developed by [START_REF] Liu | Wave-current interactions with three-dimensional floating bodies[END_REF] for this purpose. Some theoretical works can be found in [START_REF] Zhao | Interaction between waves and current on a two-dimensional body in the free surface[END_REF] and [START_REF] Li | Ship waves on uniform shear current at finite depth: wave resistance and critical velocity[END_REF].

Ship interactions

All the interactions between the ship and other structures, such as channel bank, bottom or another ship are summarized in this section. They can alter the maneuverability and controllability, and introduce additional hydrodynamic effects for ships.

Ship-bank interaction

When the ship moves near the bank, the water between ship and bank will be accelerated. The high velocity in this area results in a pressure drop (Figure 2.14), which is governed by the Bernoulli principle. The pressure difference between the port and starboard sides causes suction towards the closest bank, an effect known as the bank suction. But the lateral force on the aft part is always larger than that on the fore part of the ship. A yaw moment pushing the bow away is generated, known as the bank cushion. When the ship is attracted towards the bank, the fluids in the confined space will be squeezed, which generate a force directed away from the bank. The smaller the confinement the greater this effect will be, higher speeds and smaller under keel clearances might also aggravate the effect. When this effect becomes larger than the Bernoulli effect, suction may turn into repulsion. Together with the propulsion and hydrodynamic effects, the ship will be pushed away from the bank. During this process, the stern is prone to collide with the bank, which can be counteracted by the rudder action (Du et al., 2017;[START_REF] Torben | The development of a method to predict ship-bank interaction in real-time[END_REF]. Several empirical models for predicting the bank-induced forces can be found in [START_REF] Lee | Investigation of ship maneuvering with hydrodynamic effects between ship and bank[END_REF], [START_REF] Vantorre | Experimental investigation of ship-bank interaction forces[END_REF], Du et al. (2017Du et al. ( , 2018)), [START_REF] Duffy | Modelling of ship-bank interaction and ship squat for shiphandling simulation[END_REF] and [START_REF] Yasukawa | Bank effect on ship maneuverability in a channel with varying width[END_REF]. Ghent University has conducted perhaps the most comprehensive experiments and analyses for the bank effects [START_REF] Vantorre | Experimental investigation of ship-bank interaction forces[END_REF][START_REF] Lataire | Ship-bank interaction induced by irregular bank geometries[END_REF][START_REF] Lataire | Longitudinally directed bank effects[END_REF][START_REF] Lataire | The Influence of the ship's speed and distance to an arbitrarily shaped bank on bank effects[END_REF][START_REF] Lataire | Experiment based mathematical modelling of ship-bank interaction[END_REF][START_REF] Lataire | Impact of banks on ship squat. 4th MASHCON-International Conference on Ship Manoeuvring in Shallow and Confined Water with Special Focus on Ship Bottom Interaction[END_REF]. CFD studies were also carried out by researchers, where RANS-based works can be found in [START_REF] Mucha | Ship-bank interaction of a large tanker and related control problems[END_REF]

.2 Ship overtaking

To characterize the overtaking process, a parameter denoting the relative position between two ships is normally defined:

ξ ss = x ss /L m (2.11)
where x ss is the longitudinal distance between the two vessels.

L m = (L s1 + L s2 )/2
is the average of the ship lengths. The overall overtaking process can be divided into 5 phases (Figure 2.15) [START_REF] Kokarakis | Hydrodynamic interaction analysis in marine accidents[END_REF][START_REF] Vantorre | Model test based formulations of ship-ship interaction forces[END_REF]Du et al., 2018). When the faster ship first catches up with the slower ship (Phase 1: ξ ss = -1), the longitudinal force drags the slower ship and speeds up the faster ship. The two vessels are laterally repelled and a yaw moment induces a 'bow-out' effect for the slower ship and a 'bow-in' effect for the faster ship. The repulsion may be caused by the 'pushing' of the faster ship transferred by the fluid. At the same time, the pressure drop between the two ships accelerates the water, imposing an attraction force for both ships. In Phase 2, the attraction of this Bernoulli effect becomes dominant. When passing abeam (Phase 3: ξ ss = 0), both ships will experience strong turning forces. Their bows will be repelled and sterns will be attracted. In Phase 4, the faster ship overtakes the slower one. The longitudinal forces reverse, showing a drag for the faster ship and a thrust for the slower ship. A 'bow-in' moment is induced for the slower ship and a 'bow-out' one for the faster ship. In the final phase (ξ ss = 1), the two vessels are repelled laterally. Phases 2-4 present the highest risk of collision since both yaw moment and suction force appear, drawing the ships closer to each other. Counter helm may be necessary for safety. The relative strength of the interaction varies with the type of hulls involved and their drafts. Accidents are more likely to occur on a smaller vessel like a tugboat. Models for predicting the overtaking forces and moments of [START_REF] Brix | Manoeuvring technical manual[END_REF], [START_REF] Varyani | New generic mathematical model to predict hydrodynamic interaction effects for overtaking maneuvers in simulators[END_REF][START_REF] Varyani | New empirical and generic models to predict interaction forces for several ships in encounter and overtaking manoeuvres in a channel[END_REF] 2015) etc. Slender ship theory [START_REF] Yeung | On the interactions of slender ships in shallow water[END_REF] and potential flow methods [START_REF] Xiang | Maneuvering of two interacting ships in calm water[END_REF], such as the Rankine source method [START_REF] Söding | Analysis of overtaking manoeuvres in a narrow waterway[END_REF]Yuan et al., 2015), Hess-Smith panel method [START_REF] Sutulo | Simulation of the hydrodynamic interaction forces in close-proximity manoeuvring[END_REF], three-dimensional translating-pulsating source method [START_REF] Xu | Numerical study on wave loads and motions of two ships advancing in waves by using three-dimensional translating-pulsating source[END_REF], were also adopted. A mathematical collision model for overtaking can be found in [START_REF] Curtis | A ship collision model for overtaking[END_REF]. Special cases of overtaking like lightering operations [START_REF] De Decker | Ship-ship interaction during lightering operations[END_REF][START_REF] Lataire | Ship to ship interaction forces during lightering operations. 2nd International Conference on Ship Manoeuvring in Shallow and Confined Water: Ship to ship interaction[END_REF][START_REF] Zou | Numerical predictions of ship-to-ship interaction in shallow water[END_REF], a ship passing a stationary ship [START_REF] Varyani | Prediction of load on mooring ropes of a container ship due to the forces induced by a passing bulk carrier[END_REF]Wang & Zou, 2014) and two ships parallel traveling (Yuan et al., 2015) are also of interest. For the head-on encounter between two ships, a parameter similar with 2.11 can be defined to characterize the relative positions [START_REF] Kokarakis | Hydrodynamic interaction analysis in marine accidents[END_REF][START_REF] Vantorre | Model test based formulations of ship-ship interaction forces[END_REF]. In Phase 1 (ξ ss = -1), additional resistances will be induced. Both ships are laterally repelled and "bow-out" in this stage. In Phase 2, the lateral forces and yaw moments turn to the opposite directions because of the Bernoulli effect. Longitudinal resistances and yaw moments disappear temporarily in Phase 3 (ξ ss = 0). In Phase 4, the longitudinal forces change directions and exert resistances again on both ships. Bow-out phenomenon appears since the sterns are attracted to each other. In the final Phase (ξ ss = 1), lateral forces change to repulsions and yaw moments have the bow-in effect. It can be observed that the longitudinal force always drags the ship motions.

Head-on encounter

Several famous models for prediction the encounter forces and moments are [START_REF] Kijima | The interaction effects between two ships in the proximity of bank wall[END_REF], [START_REF] Varyani | New empirical and generic models to predict interaction forces for several ships in encounter and overtaking manoeuvres in a channel[END_REF] and [START_REF] Vantorre | Model test based formulations of ship-ship interaction forces[END_REF]. Their applications can refer to [START_REF] Lee | Assessment of safe navigation including the effect of ship-ship interaction in restricted waterways[END_REF][START_REF] Lee | A numerical study on ship-ship interaction in shallow and restricted waterway[END_REF] and [START_REF] Lima | Modelling of close proximity manoeuvres in shallow water channels[END_REF]; [START_REF] Lima | Study of ship-to-ship interaction in shallow water with account for squat phenomenon[END_REF]. CFD simulations of the encounter maneuver can be found in [START_REF] Lo | Numerical simulation of hydrodynamic interaction produced during the overtaking and the head-on encounter process of two ships[END_REF] and [START_REF] Wan | Overset-RANS computations of two surface ships moving in viscous fluids[END_REF], and potential flow theory in [START_REF] Zhou | Computation of ship hydrodynamic interaction forces in restricted waters using potential theory[END_REF].

Ship maneuvering

Ship maneuvering test is an important topic for evaluating the maneuverability of a vessel, which refers to the ability of a vessel to keep or change its motion state under control actions, including the straight line stability, course-keeping ability, course-changing ability, yaw checking ability, turning ability, stopping ability, etc. Correspondingly, turning test, zig-zag test, stopping test, spiral test and pull-out test can be carried out to determine these abilities. Detailed standards have been defined by the International Maritime Organization [START_REF] Imo | Standards for ship manoeuvrability[END_REF]. ITTC (2017,2008,2017,2014,2017) has comprehensively concluded the methods and procedures for the maneuvering tests. These methods are classified as in Figure 2.17 in this thesis. The maneuvering analysis can be conducted either by the free running test or mathematical modeling (in full-or model-scale). Free running test is carried out with functional propeller and rudder. The motion states are recorded and analyzed to characterize the maneuverability directly. Experimental data of ship maneuvering are very limited, several exist for the tankers [START_REF] Crane | Maneuvering trials of 278000-DWT tanker in shallow and deep waters[END_REF][START_REF] Broglia | Simulation of turning circle by CFD: Analysis of different propeller models and their effect on manoeuvring prediction[END_REF][START_REF] Broglia | Turning ability analysis of a fully appended twin screw vessel by CFD. Part I: Single rudder configuration[END_REF][START_REF] Dubbioso | Turning ability analysis of a fully appended twin screw vessel by CFD. Part II: Single vs. twin rudder configuration[END_REF], DTMB 5415, KVLCC1, KVLCC2, KCS [START_REF] Stern | Experience from SIM-MAN 2008-the first workshop on verification and validation of ship maneuvering simulation methods[END_REF][START_REF] Steinwand | Manövrierversuche mit einem Modell des Kriso Container Ship (KCS) Modell[END_REF][START_REF] Carrica | Direct simulation and experimental study of zigzag maneuver of KCS in shallow water[END_REF], ONRT [START_REF] Araki | Estimating maneuvering coefficients using system identification methods with experimental, system-based, and CFD free-running trial data[END_REF][START_REF] Bottiglieri | Uncertainty assessment for free-running model cases at the IIHR wave basin[END_REF][START_REF] Elshiekh | Maneuvering characteristics in calm water and regular waves for ONR Tumblehome[END_REF][START_REF] Sadat-Hosseini | CFD, system-based and EFD study of ship dynamic instability events: surf-riding, periodic motion, and broaching[END_REF][START_REF] Sanada | Trajectories for ONR Tumblehome maneuvering in calm water and waves[END_REF], SR108 [START_REF] Yasukawa | Simulation of ship maneuvering in waves (1st report: turning motion)[END_REF][START_REF] Yasukawa | Simulation of ship maneuvering in waves (2nd report: zigzag and stopping maneuvers)[END_REF] [START_REF] Bhushan | Model-and full-scale URANS simulations of Athena resistance, powering, seakeeping, and 5415 maneuvering[END_REF][START_REF] Carrica | CFD analysis of broaching for a model surface combatant with explicit simulation of moving rudders and rotating propellers[END_REF][START_REF] Carrica | Direct simulation and experimental study of zigzag maneuver of KCS in shallow water[END_REF][START_REF] Carrica | Turn and zigzag maneuvers of a surface combatant using a URANS approach with dynamic overset grids[END_REF][START_REF] Mofidi | Simulations of zigzag maneuvers for a container ship with direct moving rudder and propeller[END_REF][START_REF] Dubbioso | Analysis of propeller bearing loads by CFD. Part I: straight ahead and steady turning maneuvers[END_REF][START_REF] Dubbioso | Turning ability analysis of a fully appended twin screw vessel by CFD. Part II: Single vs. twin rudder configuration[END_REF][START_REF] Broglia | Simulation of turning circle by CFD: Analysis of different propeller models and their effect on manoeuvring prediction[END_REF][START_REF] Broglia | Turning ability analysis of a fully appended twin screw vessel by CFD. Part I: Single rudder configuration[END_REF][START_REF] Wang | CFD simulations of free running ship under course keeping control[END_REF][START_REF] Wang | Numerical simulations of zigzag maneuver of free running ship in waves by RANS-Overset grid method[END_REF][START_REF] Shen | Dynamic overset grids in OpenFOAM with application to KCS self-propulsion and maneuvering[END_REF]. Modeled propellers were sometimes used to lighten the computational load [START_REF] Jacquin | Simulation of unsteady ship maneuvering using free-surface RANS solver[END_REF][START_REF] Muscari | Trajectory prediction of a selfpropelled hull by unsteady RANS computations[END_REF][START_REF] Muscari | Analysis of propeller bearing loads by CFD. Part II: Transient maneuvers[END_REF][START_REF] Carrica | Turn and zigzag maneuvers of a surface combatant using a URANS approach with dynamic overset grids[END_REF][START_REF] Dubbioso | Analysis of propeller bearing loads by CFD. Part I: straight ahead and steady turning maneuvers[END_REF][START_REF] Jasak | CFD validation and grid sensitivity studies of full scale ship self propulsion[END_REF], which is commonly-used in commercial codes (FINE-MARINE, STAR-CCM+, etc.). Apart from the free running test, the mathematical modeling approach, also called the system-based method, is widely-used, where calculating the hydrodynamic coefficients is critical. This method is based on the solving of the rigid body equations:

m[ u -vr + wq -x g (q 2 + r 2 ) + y g (pq -ṙ) + z g (pr + q)] = X m[ v -wp + ur -y g (r 2 + p 2 ) + z g (qr -ṗ) + x g (qp + ṙ)] = Y m[ ẇ -uq + vp -z g (p 2 + q 2 ) + x g (rp -q) + y g (rq + ṗ)] = Z I x ṗ + (I z -I y )qr -( ṙ + pq)I xz + (r 2 -q 2 )I yz + (pr -q)I xy +m[y g ( ẇ -uq + vp) -z g ( v -wp + ur)] = K I y q + (I x -I z )rp -( ṗ + qr)I xy + (p 2 -r 2 )I zx + (qp -ṙ)I yz +m[z g ( u -vr + wq) -x g ( ẇ -uq + vp)] = M I z ṙ + (I y -I x )pq -( q + rp)I yz + (q 2 -p 2 )I xy + (rq -ṗ)I zx +m[x g ( v -wp + ur) -y g ( u -vr + wq)] = N (2.12)
where the first three equations are for the translational motions, the last three for the rotational motions. Thr forces and moments (X, Y, Z, K, M, N ) on the RHS are expanded using Taylor series. (u, v, w, p, q, r) are the linear and angular velocities.

I x , I y , I z , ... are the moments of inertia. (x g , y g , z g ) is the vector from the gravity center to the rotation center. Their derivatives regarding the variables (velocity, yaw rate, rudder angle, etc.) bring about the hydrodynamic coefficients. Other forces generated by wind, waves, current, ship-ship interactions, etc. can also be included.

Their calculations can refer to the former parts. Several models exist for calculating the hydrodynamic coefficients. The linearized models [START_REF] Davidson | Turning and course keeping qualities[END_REF][START_REF] Nomoto | On the steering qualities of ships[END_REF] were used in early studies. Nonlinear models [START_REF] Bibliography Abkowitz | Lectures on ship hydrodynamics-steering and manoeuvrability[END_REF][START_REF] Norrbin | Theory and observation on the use of a mathematical model for ship maneuvering in deep and confined waters[END_REF][START_REF] Blanke | Ship propulsion losses related to automatic steering and prime mover control[END_REF][START_REF] Yasukawa | Introduction of MMG standard method for ship maneuvering predictions[END_REF][START_REF] Ross | A novel manoeuvering model based on low-aspect-ratio lift theory and Lagrangian mechanics[END_REF] are of the trend currently. Among them, the models of [START_REF] Bibliography Abkowitz | Lectures on ship hydrodynamics-steering and manoeuvrability[END_REF] and MMG [START_REF] Yasukawa | Introduction of MMG standard method for ship maneuvering predictions[END_REF] are the most commonly-used. For the MMG model, the forces of the hull, rudder, propeller are separated, and so their hydrodynamic coefficients, while for the [START_REF] Bibliography Abkowitz | Lectures on ship hydrodynamics-steering and manoeuvrability[END_REF] model, they are treated as a whole. Models coupling the ship steering and rolling can be found in Vam Ameronge (1981); [START_REF] Son | On the coupled motion of steering and rolling of a high speed container ship[END_REF]; [START_REF] Christensen | A linearized state-space model of steering and roll of a high-speed container ship[END_REF]. The methods for determining the hydrodynamic coefficients are distinguished in Figure 2.17. (1) The empirical formula method is the most convenient but less accurate. Good reviews of these methods can be found in ITTC (2002); [START_REF] Vantorre | Manoeuvring coefficients for a container carrier in shallow water: an evaluation of semi-empirical formulae[END_REF].

(2) Captive model test refers to the forced ship motions using PMM (Planar Motion Mechanism). Standard procedures have been defined by ITTC (2017ITTC ( , 2014)).

Experimental data can be found for the KVLCC1, KVLCC2, KCS, DTMB 5415, ONRT on the website of the Tokyo 2015 CFD Workshop in ship hydrodynamics and SIMMAN 2008, 2014, 2019[START_REF] Stern | Experience from SIM-MAN 2008-the first workshop on verification and validation of ship maneuvering simulation methods[END_REF][START_REF] Yasukawa | Introduction of MMG standard method for ship maneuvering predictions[END_REF][START_REF] Yoon | Phase-averaged stereo-PIV flow field and force/moment/motion measurements for surface combatant in PMM maneuvers[END_REF][START_REF] Simonsen | Maneuvering predictions in the early design phase using CFD generated PMM data[END_REF][START_REF] Sakamoto | URANS simulations of static and dynamic maneuvering for surface combatant: Part 1. Verification and validation for forces, moment, and hydrodynamic derivatives[END_REF] (3) System identification can be used to calculate or correct the hydrodynamic coefficients based on the free running test results [START_REF] Stern | Experience from SIM-MAN 2008-the first workshop on verification and validation of ship maneuvering simulation methods[END_REF]ITTC, 2014;[START_REF] Luo | Parametric identification of ship maneuvering models by using support vector machines[END_REF][START_REF] Luo | Parameter identification of ship maneuvering model based on support vector machines and particle swarm optimization[END_REF][START_REF] Rajesh | System identification for nonlinear maneuvering of large tankers using artificial neural network[END_REF][START_REF] Sutulo | An algorithm for offline identification of ship manoeuvring mathematical models from free-running tests[END_REF]. To apply this method, standard maneuvers should be conducted to obtain system inputs (rudder angle, propeller revolution, etc.) and outputs (ship velocities, heading, yaw rate, etc.). The advantage of system identification is that all the maneuvering coefficients can be estimated by one or a few free running trials as opposed to numerous captive model tests. Classical system identification techniques for ship maneuvering include least squares method [START_REF] Nonaka | Estimating maneuvering coefficients from free-running trials[END_REF][START_REF] Lasdon | Design and testing of a generalized reduced gradient code for nonlinear programming[END_REF][START_REF] Rhee | Estimation of manoeuvring coefficients from PMM test by genetic algorithm[END_REF], model reference method [START_REF] Hayes | Parametric identification of nonlinear stochastic systems applied to ocean vehicle dynamics[END_REF][START_REF] Van Amerongen | Adaptive steering of ships-A model reference approach[END_REF], extended Kalman Filter method [START_REF] Shi | Identification of ship maneuvering model using extended Kalman Filters[END_REF][START_REF] Abkowitz | Measurement of hydrodynamic characteristics from ship maneuvering trials by system identification[END_REF][START_REF] Hwang | Application of system identification to ship maneuvering[END_REF], maximum likelihood method [START_REF] Åström | Identification of ship steering dynamics[END_REF], recursive prediction error method [START_REF] Zhou | Identification of nonlinear marine systems[END_REF][START_REF] Källström | Experiences of system identification applied to ship steering[END_REF], annealing search method [START_REF] Tiano | Multivariable identification of ship steering and roll motions[END_REF]. However, some problems exist in these methods, such as sensitivity to the initial values, ill-conditioned solutions, simultaneous drift of parameters and multicollinearity of regressive parameters. Several new approaches have also been proposed for the modeling of ship maneuvering, including the Markov process theories, statistical linearization techniques [START_REF] Roberts | Stochastic estimation methods for non-linear ship roll motion[END_REF], frequency spectrum analysis [START_REF] Selvam | A frequency domain system identification method for linear ship maneuvering[END_REF][START_REF] Bhattacharyya | Parametric identification for nonlinear ship maneuvering[END_REF], optimization techniques [START_REF] Di Mascio | Investigation of twin-screw naval ships maneuverability behavior[END_REF][START_REF] Khanh | Assessment of ship manoeuvrability by using a coupling between a nonlinear transient manoeuvring model and mathematical programming techniques[END_REF]Du et al., 2017), genetic algorithm [START_REF] Sutulo | An algorithm for offline identification of ship manoeuvring mathematical models from free-running tests[END_REF], Particle Swarm Optimization (PSO) [START_REF] Chen | Parameters identification for ship motion model based on particle swarm optimization[END_REF], Artificial Neural Networks (ANN) [START_REF] Haddara | Parametric identification of manoeuvring models for ships[END_REF][START_REF] Mahfouz | Effect of the damping and excitation on the identification of the hydrodynamic parameters for an underwater robotic vehicle[END_REF][START_REF] Hess | Simulation of ship maneuvers using recursive neural networks[END_REF][START_REF] Moreira | Dynamic model of manoeuvrability using recursive neural networks[END_REF][START_REF] Rajesh | System identification for nonlinear maneuvering of large tankers using artificial neural network[END_REF][START_REF] Yin | On-line prediction of ship roll motion during maneuvering using sequential learning RBF neural networks[END_REF], and Support Vector Machine (SVM) [START_REF] Drezet | Support vector machines for system identification[END_REF][START_REF] Goethals | Identification of MIMO Hammerstein models using least squares support vector machines[END_REF][START_REF] Gretton | Support vector regression for black-box system identification[END_REF][START_REF] Luo | Parameter identification of ship maneuvering model based on support vector machines and particle swarm optimization[END_REF][START_REF] Luo | Parametric identification of ship maneuvering models by using support vector machines[END_REF][START_REF] Wang | System identification modelling of ship manoeuvring motion based on ε-support vector regression[END_REF]. Among them, ANN and SVM are applicable both in the parameter identification and blind prediction of ship maneuvering. Differently, for the parameter identification, the hydrodynamic coefficients can be obtained directly by SVM whereas another regressor (e.g., the least square method) is required by ANN. For the blind prediction of ship maneuvering, it can be confirmed that SVM gains better generalization ability than ANN. Sensitivity analysis can also be conducted to determine the most influential variables for identification [START_REF] Blanke | A sensitivity approach to identification of ship dynamics from sea trial data[END_REF]Wang et al., 2014;[START_REF] Khanh | Assessment of ship manoeuvrability by using a coupling between a nonlinear transient manoeuvring model and mathematical programming techniques[END_REF]Du et al., 2017).

Concluding remarks

In this chapter, the numerical methods for the virtual tank and the principal topics of ship hydrodynamics study were presented and summarized.

For the numerical methods, the most important ones in CFD were introduced, especially those commonly-used for maritime applications, such as BEM, FDM and FVM. The work in this thesis is based on FVM.

For the ship hydrodynamics investigations, the ship resistance, ship-generated waves, confinement effect, seakeeping, ship interactions and maneuvering were reviewed, which are the most relevant topics with this thesis. The most general and advanced methods for each topic were presented, providing solid theoretical and engineering support for this work. Other topics such as the hull-propeller-rudder interaction, ship motion control, hull form design and optimization, etc. are also of interest but out of the range of this thesis.

Chapter 3

Framework of the numerical solver

Introduction

In this chapter, the details of the numerical solver are presented. The following aspects are emphasized:

• The RANS equations with multiphase and turbulence models are presented.

• The finite volume discretization procedure including the deduction of the pressure Poisson equation and the pressure-velocity coupling are explained.

• The ship is viewed as rigid body, and the commonly-used mesh motion approaches and procedures are concluded, for either prescribed or 6DoF motions.

Governing equations

For marine simulations, the incompressible NS equations (conservative form) are normally used:

∇ • u = 0 (3.1) ∂(ρu) ∂t + ∇ • (ρuu) = -∇p + ∇ • (µ∇u) + S u (3.2)
The following analyses and discretizations are based on these equations.

RANS

Since DNS (Direct Numerical Simulation) and LES (Large-Eddy Simulation) are computationally expensive, Reynolds decomposition is always used for the NS equations 3.1 and 3.2 to separate the velocity fluctuations from the mean flow velocity. The famous Reynolds-Averaged Navier-Stokes (RANS) equations can then be obtained [START_REF] Jasak | Error analysis and estimation for the finite volume method with applications to fluid flows[END_REF]:

∇ • u = 0 (3.3) ∂(ρu) ∂t + ∇ • (ρuu) = -∇p + ∇ • (µ∇u) + S u + R(u , u ) (3.4)
where u is the time average of the velocity u. R(u , u ) is the so-called Reynolds stress term. The Boussinesq assumption is then adopted to relate the Reynolds stress tensor to the mean velocity gradient:

R(u , u ) = ∇ • (µ t [∇u + (∇u) T ] - 2 3 ρk t I) = ∇ • (µ t ∇u) -∇ • ( 2 3 ρk t I) (3.5)
The RHS of 3.4 can then be rearranged:

-∇p + ∇ • (µ∇u) + S u + R(u , u ) = -∇ • (pI + 2 3 ρk t I) + ∇ • [(µ + µ t )∇u] + S u = -∇p * + ∇ • (µ ef f ∇u) + S u (3.6)
µ ef f is the effective viscosity, k t is the turbulent kinetic energy and µ t is the eddy viscosity obtained from a turbulence model. Gravity (ρg) and surface tension (f σ ) are normally included as the source terms for multiphase problems. Omitting the overline, 3.3 and 3.4 can be reformulated as:

∇ • u = 0 (3.7) ∂(ρu) ∂t + ∇ • (ρuu) = -∇p * + ∇ • (µ ef f ∇u) + ρg + f σ (3.8)

Multiphase model

The VOF (Volume-of-fluid) approach is used for multiphase flow simulations, together with an artificial compression term. The transport equation reads [START_REF] Rusche | Computational fluid dynamics of dispersed two-phase flows at high phase fractions[END_REF]:

∂α p ∂t + ∇ • [α p (u -u g )] + ∇ • [α p (1 -α p )u r ] = 0 (3.9)
where α p is the phase fraction, which takes values within the range 0 ≤ α p ≤ 1. α p = 0 and α p = 1 correspond with air and water respectively. u g is the grid velocity taking into account the mesh motion. It can be seen that the compression term (the last term on the LHS) only takes effect within the interface. This term is able to compress the free surface towards a sharper one. u = α p u w + (1 -α p )u a is the effective velocity, and u r = u w -u a is the relative velocity between the two phases, where the subscripts 'w' and 'a' denote water and air respectively. The density and dynamic viscosity are calculated according to:

ρ = α p ρ w + (1 -α p )ρ a µ = α p µ w + (1 -α p )µ a (3.10)
The surface tension term f σ is calculated as:

f σ = σ w κ w ∇α p (3.11)
where σ w is the surface tension coefficient (0.07 kg/s 2 in water). κ w is the curvature of the free surface interface, determined from the volume of fraction by κ w = -∇ • (∇α p /|∇α p |).

Turbulence models

Various turbulence models have been implemented in OpenFoam, including Reynoldsaveraged stress (RAS), LES or detached-eddy simulation (DES) models. The common two-equation turbulence models, such as k-ε, k-ω, SST k-ω, etc., are directly used for the simulations in this thesis. The initial values of the turbulent kinetic energy k t , the turbulent dissipation ε t and ω t are calculated as follows (Du et al., 2018):

k t = 3 2 (UI t ) 2 ω t = k t ν ν t ν -1 ε t = C µ ρk 2 t µ µ t µ -1 (3.12)
where the turbulence intensity I t , the eddy viscosity ratio νt ν , µt µ and the turbulence constant C µ are taken to be 1%, 1 and 0.09. Sensitivity analyses will be conducted in the next chapter to choose an appropriate turbulence model. Wall functions are used in the near-wall region, so the non-dimensional wall distance of the first layer mesh needs to be controlled within the range of 30≤ y + ≤200.

Finite volume discretization 3.3.1 Equation discretization

The equation discretization is to transform the governing equations 3.1 and 3.2 into a set of algebraic equations, which is solved to get the approximate solutions at pre-defined locations and times. The continuity equation 3.1 and the momentum equation 3.2 are integrated over the CV and time: [START_REF] Jasak | Error analysis and estimation for the finite volume method with applications to fluid flows[END_REF].

t+∆t t V P ∇ • udV dt = 0 (3.13) t+∆t t V P ∂(ρu) ∂t dV + V P ∇ • (ρuu)dV dt = t+∆t t V P ∇ • (µ∇u)dV + V P S u dV dt -∇p (3.14)
The pressure gradient term is not discretized currently in the spirit of the Rhie-Chow procedure. The discretization of each term is shown below (Figure 3.1):

• Continuity equation:

V P ∇ • udV = S dS • u ≈ f S • u f = 0 (3.15)
• Temporal term:

V P ∂(ρu) ∂t dV ≈ ρ P u n P -u o P ∆t V P (3.16)
where the superscripts 'n' and 'o' denote the new and old time values, i.e. u n = u(t + ∆t), u o = u(t).

• Convection term:

V P ∇ • (ρuu)dV = S dS • (ρuu) ≈ f S • (ρu) f u f = f F f u f (3.17)
where F f = S•(ρu) f is the face flux. u f can be evaluated using face interpolation from cell center values.

• Diffusion term:

V P ∇ • (µ∇u)dV = S dS • (µ∇u) ≈ f µS • (∇u) f (3.18)
where (∇u) f denotes the gradient at the face f obtained from interpolation.

Corrections should be made here for non-orthogonal meshes.

• Source term:

V P S u dV = S P V P u P + S C V P (3.19)
It can be seen that the source term is linearized. S P and S C are the coefficients obtained from the linearization. Therefore, a semi-discretized momentum equation is obtained:

t+∆t t ρ P u n P -u o P ∆t V P + f F f u f dt = t+∆t t f µS • (∇u) f + S P V P u P + S C V P dt -∇p (3.20)
The time integral can select the explicit, implicit or the Crank-Nicholson scheme, which depends on the utilization of the old, new or the blended values. By adopting the Euler implicit scheme, the final discretized continuity and momentum equations can be written as:

f S • u n f = 0 (3.21) ρ P u n P -u o P ∆t V P + f F f u n f = f µS • (∇u) n f + S P V P u n P + S C V P -∇p (3.22)

Derivation of the pressure equation

In 3.21 and 3.22, there is no equation for the pressure, which will be deduced here. Equation 3.22 is simplified as:

A P u n P = H u -∇p (3.23)
The H u term consists of two parts: the 'transport part' related with all the neighbors and the 'source part' including the contributions of the transient term and all other source terms apart from the pressure gradient. Then:

u n P = H u A P - ∇p A P (3.24)
The velocity on the cell face u f is calculated as the face interpolation of this equation:

u n f = ( H u A P ) f -( 1 A P ) f (∇p) n f (3.25)
A pressure Poisson equation can be obtained by substituting 3.24 into the continuity equation 3.1:

∇ • 1 A P ∇p = ∇ • H u A P (3.26)
The discretization of this Poisson equation can be realized by substituting 3.25 into 3.21 and the final form of the discretized incompressible NS equations is:

f S • ( 1 A P ) f (∇p) n f = f S • ( H u A P ) f (3.27) A P u n P = H u - f Sp n f (3.28)

Pressure-velocity coupling

From equations 3.27 and 3.28, it can be seen that the velocity and the pressure depend on each other. As a result, a decoupling procedure is needed for them.

Commonly-used methods are PISO (Pressure Implicit with Split Operator) and SIMPLE (Semi-Implicit Method for Pressure Linked Equations). A new approach called PIMPLE (merged PISO-SIMPLE) is also implemented in OpenFoam. The procedure of the PIMPLE algorithm are explained here (Figure 3.2):

(1) The pressure and velocity fields (p o , u o ) are built based on the old values from initialization or last time step.

(2) An under-relaxed form of the momentum equation 3.28 is solved:

A P α r u u * P = H u - f Sp n f + 1 -α r u α r u A P u o P (3.29)
where α r u is the velocity under-relaxation factor. This procedure is called the momentum predictor. A new velocity field u * can be obtained.

(3) The face flux F * f can then be approximated as:

F * f = S • ( H u A P ) f (3.30)
which is needed in the RHS of the pressure Poisson equation 3.27:

f S • ( 1 A P ) f (∇p) * f = f F * f (3.31)
This procedure is called the pressure solution. A new pressure field p * can be obtained.

(4) With the new pressure field p * , the face flux can then be updated according to 3.25:

F n f = F * f -S • ( 1 A P ) f (∇p) * f (3.32) (5)
The pressure solution is under-relaxed explicitly:

p n = p * + α r p (p * -p o ) (3.33)
where α r p is the pressure under-relaxation factor. This step can be repeated to account for the mesh non-orthogonality (nNonOrthogonalCorrectors). [START_REF] Du | Hydrodynamics prediction of a ship in static and dynamic states[END_REF] The velocity is corrected using 3.27 and the new pressure p n :

u n P = H u A P - 1 A P f Sp n f (3.34)
This step is also called the momentum corrector.

Step 3-6 can be repeated several times (nCorrectors). This is the characteristic of the PISO algorithm, which consists of an implicit momentum predictor and a series of pressure solutions and explicit velocity corrections.

(7) All the parameters with the superscript 'n' (F n f , p n , u n ) are the final ones at the end of the PIMPLE process. The whole process from step 2 can also be repeated (nOuterCorrectors). Therefore, the PIMPLE algorithm is an iterative method merging the PISO and SIMPLE methods. The SIMPLE is the PIMPLE without nCorrectors, and the PISO is the one without under-relaxations and nOuterCorrectors. PIMPLE allows the usage of large time steps, which breaks the limit of the CFL (Courant-Friedrichs-Lewy) condition and greatly reduces the simulation time.

Solid body motion techniques

To account for the mesh motion, the momentum equation 3.8 is reformulated using ALE (Arbitrary Lagrangian-Eulerian) as:

∂(ρu) ∂t + ∇ • (ρ(u -u g )u) = -∇p * + ∇ • (µ ef f ∇u) + ρg + f σ (3.35)
where u g is the grid velocity, which can be determined by the space conservation law:

∂ ∂t V dV - S u g • dS = 0 (3.36)
Normally, solid body motions are either prescribed or solved based on the rigid body motion equations, where the former is simpler since the state of the next time step is given. To realize these motions, dynamic mesh techniques should be adopted.

Dynamic mesh methods

Dynamic mesh methods are necessary to realize the solid body motions. Common techniques include sliding mesh, deforming mesh, overset, layer addition/removal, Immersed Boundary Method (IBM), remeshing, etc. The first three techniques have no topological changes of the mesh. Remeshing method [START_REF] Menon | A numerical study of droplet formation and behavior using interface tracking methods[END_REF] is available for OpenFoam but not mature enough for engineering applications. The MRF (Multiple Reference Frame) is a good approach for steady-state simulations with no need of mesh motions, but it is not capable of reproducing unsteady effects.

Sliding mesh

Sliding mesh enables simulation across disconnected, adjacent mesh domains. Two techniques, AMI (Arbitrary Mesh Interface) for OpenFoam and GGI (General Grid Interface) for foam-extend, are available respectively. Depending on the meshes are fully or partially overlapped, AMI is further divided into AMI (Arbitrary Mesh Interface) and ACMI (Arbitrary Coupled Mesh Interface) as in Figure 3.3. Since the meshes have to match with each other at the interface, sliding meshes are mostly used for prescribed solid body motions.

During the calculation, each face at the sliding interface accepts contributions from overlapping faces of the neighbor patch, with the weights defined as a fraction of the intersecting areas. For each face, the sum of the weights should be 1. Interpolation weights at the interface are constructed such that the interpolation results in conservation. Conservation errors are introduced as the sum of weights deviates from 1 where the patch geometries are not well matched (Figure 3.3). Consistent and conservative discretization across the interface is achieved using weighted interpolation.

For the flow variables from the master patch to the slave patch [START_REF] Beaudoin | Development of a generalized grid interface for turbomachinery simulations with OpenFOAM[END_REF][START_REF] Farrell | Conservative interpolation between volume meshes by local Galerkin projection[END_REF]:

φ S i = n W Mn-S i φ Mn (3.37)
For the ones from the slave patch to the master patch:

φ M j = m W Sm-M j φ Sm (3.38)
In order for the interface discretization to remain conservative, the following three constraints should be obeyed:

W Mn-S i = 1.0 (3.39) W Sm-M j = 1.0 (3.40) W Mn-S i |S Mn | = W Sm-M j |S M j | = |S M ∩S | (3.41)
with the additional symmetry constraint:

if W Mn-S i > 0 ⇒ W S i -Mn > 0 (3.42)
where φ S and φ M are the variables of the slave and master patch. The value of the weighting factors is basically the percentage of the surface intersection between two overlapping faces deduced from 3.41:

W M -S i = |S M ∩S i | |S Mn | ∈ [0.0, 1.0] (3.43) W S-M i = |S S∩M j | |S Sn | ∈ [0.0, 1.0] (3.44)
Because of the domain discretization, the two meshes may not overlap with each other, making the weighting factors underestimated (smaller than 1.0), i.e. the AMI will not be conservative. It is therefore suggested to use enough mesh resolution at the interface and use similar mesh size between the master and slave patches for the sake of accuracy. If this cannot be satisfied, remedies can be adopted. For AMI, corrections for low weight ("lowWeightCorrection") can be used in the AMI boundary conditions. GGI rescales the face weighting factors so they will sum up to 1.0.

Deforming mesh

Deforming mesh techniques are used to diffuse the movements of the solid body into surrounding meshes, thus keeping the mesh quality. In OpenFoam, the mesh deformation techniques by solving the displacement-based or velocity-based Laplace equations are adopted [START_REF] Jasak | Dynamic mesh handling in OpenFOAM applied to fluid-structure interaction simulations[END_REF]:

∇ • (γ d ∇d g ) = 0 ∇ • (γ d ∇u g ) = 0 (3.45)
where γ d is the diffusivity coefficient based on the mesh quality or the cell center distance to the selected boundary. d g and u g are the grid displacement and grid velocity, which can be obtained by the prescribed motion or the solving of rigid body equations. Then the new mesh positions will be updated according to:

x new = x old + d g x new = x old + u g ∆t (3.46)
The velocity-based approach is comparably more robust since it needs one more step of integration, which reduces the numerical oscillations.

Using the Laplace equation, the mesh deformation adjacent to the moving boundary is the largest, which may deteriorate the mesh quality potentially. This problem can be reduced by assigning different diffusivity models (quality-or distance-based) in the Laplace equations 3.45 [START_REF] Jasak | Automatic mesh motion for the unstructured finite volume method[END_REF]. Other mesh deformation techniques based on the algebraic mesh motion and Radial Basis Function (RBF) are also available in OpenFoam [START_REF] Casadei | Dynamic-mesh techniques for unsteady multiphase surfaceship hydrodynamics[END_REF][START_REF] Jasak | Dynamic mesh handling in OpenFOAM applied to fluid-structure interaction simulations[END_REF].

Overset

Overset, overlapping grid, or Chimera is a technique to construct meshes for different parts of the domain separately. These part meshes overlap with each other and communicate using interpolation at the interfaces. Best meshing strategies can be adopted for each part and mesh quality can be sustained since no mesh morphing occurs. Figure 3.4 illustrates the essential idea of the overset method using a one-dimensional domain with Dirichlet boundary conditions at the left and right sides. This domain is split into two sub-domains A and B. However, the solutions (φ) at the node 6 and 7 for the two sub-domains are unknown. They are usually named fringe nodes/cells. Their values are calculated using interpolation from the surrounding values. For node 6, several donor nodes/cells (7, 8, 9, 10 here) are selected, and the explicit interpolation formula can be given as:

φ n+1 6 = (β 7 φ n 7 + β 8 φ n 8 + β 9 φ n 9 + β 10 φ n 10 )/β 6 (3.47)
where β i are interpolation coefficients. The superscripts n and n + 1 represent the old and new time steps respectively. Implicit interpolation can also be used. This formula can be calculated multiple times along with the prediction and correction steps in one time step during the simulation. Same interpolation procedure can be adopted for the node 7. A more complex case with a structure is presented in Figure 3.5. Two domains are generated, the background mesh A and the body-fitted mesh B. Another procedure called hole cutting is needed for this case. Before the real computation, the mesh A is scanned to mark and block the nodes/cells (8 and 9 here) at the same position as the structure. They are marked as the hole nodes/cells. Their neighboring nodes 7 and 10 could also be blanked out to reduce the number of points taking part in the computation. The valid points adjacent to the hole cells (6 and 11) can be selected as fringe nodes/cells. For mesh B, point 17 and 20 can be selected as fringe cells.

Of course the points at the two ends of the mesh B can also be marked as hole nodes/cells to reduce the overall mesh numbers. To ensure a smooth transition of flow variables between the sub-domains, more fringe nodes/cells can be selected.

Then the procedure for selecting the interpolation nodes, i.e. the Domain Connectivity Information (DCI) can be executed. These important information should be calculated and stored before solving the flow problem. The subsequent procedures are interpolation and computation as in Figure 3.4. For stationary problems, hole cutting and DCI can be processed only once. For moving boundaries, these steps need to be repeated to update these information. The computational procedure can be found in Figure 3.6. Native overset functionality in OpenFoam only executes the hole cutting for the background mesh as in Figure 3.5. The commercial overset package Suggar++ [START_REF] Noack | SUGGAR: a general capability for moving body overset grid assembly[END_REF]Noack et al., 2009) has been connected with OpenFoam [START_REF] Shen | Dynamic overset grids in OpenFOAM with application to KCS self-propulsion and maneuvering[END_REF]Noack & Boger, 2009;[START_REF] Miller | An overset grid method for fluid-structure interaction[END_REF]. Several in-house overset codes are also available [START_REF] Chandar | Assessment of interpolation strategies and conservative discretizations on unstructured overset grids in OpenFOAM[END_REF][START_REF] Lu | An overlapping domain decomposition based near-far field coupling method for wave structure interaction simulations[END_REF][START_REF] Ma | An overset mesh based multiphase flow solver for water entry problems[END_REF][START_REF] Gross | Development and Application of the overset grid library Bellerophon[END_REF].

Layer addition/removal

Layer addition/removal, or dynamic cell layering is able to add or remove layers of cells in accordance with the boundary motion. The thickness of the mesh layer is monitored during the simulation. Cells will be added when the thickness of the stretched layer rises above the maximum threshold and removed when the thickness of the squeezed layer falls below the minimum threshold. This strategy keeps the mesh quality since the majority of the grid remains fixed. It is especially useful for translating boundary motions like piston engines. Currently OpenFoam only supports layer addition/removal at planar surfaces. It can be used with ACMI for translational motion simulations.

Immersed boundary method

IBM (Immersed Boundary Method) discretizes the solid boundary into Lagrangian points, and represents the influence of the solid on the fluid using a boundary force. This force is added into the momentum equation as a source term to take effects. The treatment of the fluid and solid meshes can then be separated, namely the Eulerian discretization for the fluid and the Lagrangian discretization for the solid. The fluid mesh does not need to conform to the solid geometry. The solid as such seems "immersed" in the fluid. Simple structured mesh can be used for the fluid. This approach eliminates the mesh quality problem caused by the body-conforming mesh method and the mesh distortion and interpolation errors of the remeshing and mesh deforming process. The mathematical formulation of IBM is (Figure 3.7):

∂(ρu) ∂t + ∇ • (ρuu) = -∇p + ∇ • (µ∇u) + f IBM (3.48) Ω f u(x, t)δ(x -X(s, t))dx = U b (3.49) f IBM (x, t) = Γs F(s, t)δ(x -X(s, t))ds (3.50)
in the domain Ω f , where U b is the solid boundary velocity. δ is the Dirac delta function. X(s, t) represents the solid boundary position. F(s, t) and f IBM (x, t) are the boundary force defined on the Lagrangian position X(s, t) and the corresponding force on the Eulerian frame. It can be seen that the boundary force acts as a Lagrange multiplier to satisfy the no-slip wall condition 3.49 at the immersed interface. The evaluation of this force differs one IBM from another. Recent reviews can be found in [START_REF] Iaccarino | Immersed boundary technique for turbulent flow simulations[END_REF]; [START_REF] Mittal | Immersed boundary methods[END_REF]; [START_REF] Sotiropoulos | Immersed boundary methods for simulating fluid-structure interaction[END_REF].

Communication between the two frames of reference is performed via two steps:

• Interpolation, where the physical quantities in the Eulerian mesh are interpolated onto the Lagrangian points to estimate the boundary force;

• Spreading, where the boundary force previously calculated on the Lagrangian points is spread back on the Eulerian mesh.

A discrete forcing IBM is implemented in foam-extend and OpenFoam [START_REF] Riahi | A pressurecorrected Immersed Boundary Method for the numerical simulation of compressible flows[END_REF][START_REF] Constant | An immersed boundary method in OpenFOAM: Verification and validation[END_REF]. However, the loss of accuracy in the vicinity of immersed boundary has been observed. Particularly, IBM solver exhibits strong instability next to the immersed boundary for multiphase flows. The crux of the problem is the inability of IBM to match gradients of different variables at the immersed boundary, and polynomial fitting of the VOF field which contains a stepfunction. A new Immersed Boundary Surface method based on the cut-cell approach has been implemented and validated [START_REF] Jasak | Immersed boundary surface method in foam-extend[END_REF]. Several user-implemented ver-sions can also be found in

Prescribed solid body motion

A ship has six degrees of freedom (Figure 3.8), with three translations (surge, sway, heave) and three rotations (roll, pitch, yaw). Prescribed motions can be realized by their combinations, which are useful for captive model tests. Since the trajectory of the rigid body is already known, sliding mesh, layer addition/removal, overset and IBM techniques can all be adopted for this case.

6DoF motion solver

To calculate the 6DoF motions of the ship, a set of rigid body equations should be solved [START_REF] Shen | Dynamic overset grids in OpenFOAM with application to KCS self-propulsion and maneuvering[END_REF][START_REF] Carrica | Ship motions using single-phase level set with dynamic overset grids[END_REF][START_REF] Fossen | Handbook of marine craft hydrodynamics and motion control[END_REF]:

m[ u -vr + wq -x g (q 2 + r 2 ) + y g (pq -ṙ) + z g (pr + q)] = X m[ v -wp + ur -y g (r 2 + p 2 ) + z g (qr -ṗ) + x g (qp + ṙ)] = Y m[ ẇ -uq + vp -z g (p 2 + q 2 ) + x g (rp -q) + y g (rq + ṗ)] = Z I x ṗ + (I z -I y )qr -( ṙ + pq)I xz + (r 2 -q 2 )I yz + (pr -q)I xy +m[y g ( ẇ -uq + vp) -z g ( v -wp + ur)] = K I y q + (I x -I z )rp -( ṗ + qr)I xy + (p 2 -r 2 )I zx + (qp -ṙ)I yz +m[z g ( u -vr + wq) -x g ( ẇ -uq + vp)] = M I z ṙ + (I y -I x )pq -( q + rp)I yz + (q 2 -p 2 )I xy + (rq -ṗ)I zx +m[x g ( v -wp + ur) -y g ( u -vr + wq)] = N (3.51)
The 6DoFs of a ship (Figure 3.8) are denoted as (x 1 , x 2 ) = (x, y, z, φ, θ, ψ), where x 1 and x 2 are the translational and rotational motions respectively. The linear and angular velocities are represented as (v 1 , v 2 ) = (u, v, w, p, q, r). (x g , y g , z g ) is the vector from the center of rotation to the center of gravity. I x , I y and I z are the moments of inertia around the center of rotation. The forces and moments can be integrated from the pressure and shear stress on the solid surface. The rigid body accelerations can then be calculated, which are further integrated to obtain the velocities and displacements. Deforming mesh technique can be used for small displacements, while overset and IBM has no such limitations. Important concepts when solving the 6DoF equations are the coordinate transformation between the earth and body fixed systems (Figure 3.8) and the numerical integration techniques.

Coordinate transformation

• Euler angle approach This is the traditional approach based on the three basic rotation matrices:

R x (φ) =    1 0 0 0 cosφ -sinφ 0 sinφ cosφ    (3.52) R y (θ) =    cosθ 0 sinθ 0 1 0 -sinθ 0 cosθ    (3.53) R z (ψ) =    cosψ -sinψ 0 sinψ cosψ 0 0 0 1    (3.54)
The linear and angular velocities in the body-fixed frame can be transformed to the earth-fixed frame by:

ẋ1 =    ẋ ẏ ż   = R x (φ)R y (θ)R z (ψ)    u v w    =    cosθcosψ sinφsinθcosψ -cosφsinψ cosφsinθcosψ + sinφsinψ cosθsinψ sinφsinθsinψ + cosφcosψ cosφsinθsinψ -sinφcosψ -sinθ sinφcosθ cosφcosθ       u v w    = R e 1 v 1 (3.55) ẋ2 =    φ θ ψ   =    1 sinφtanθ cosφtanθ 0 cosφ -sinφ 0 sinφ/cosθ cosφ/cosθ       p q r    = R e 2 v 2 (3.56)
Notice that this equation fails when θ = 90 • , which is known as the gimbal lock phenomenon. To avoid the singularity problem caused by this phenomenon, another approach based on the unit quaternions can be adopted.

• Unit quaternion approach

A unit quaternion is defined as a complex number with one real part η q and three imaginary parts ε q = [ε 1 , ε 2 , ε 3 ], which satisfies:

η 2 q + ε 2 1 + ε 2 2 + ε 2 3 = 1 (3.57)
The transformation for linear and angular velocities can be realized by:

ẋ1 =    ẋ ẏ ż   =    1 -2(ε 2 2 + ε 2 3 ) 2(ε 1 ε 2 -ε 3 η q ) 2(ε 1 ε 3 + ε 2 η q ) 2(ε 1 ε 2 + ε 3 η q ) 1 -2(ε 2 1 + ε 2 3 ) 2(ε 2 ε 3 -ε 1 η q ) 2(ε 1 ε 3 -ε 2 η q ) 2(ε 2 ε 3 + ε 1 η q ) 1 -2(ε 2 1 + ε 2 2 )       u v w    = R q 1 v 1 (3.58) ẋ2 =    φ θ ψ   = 1 2       -ε 1 -ε 2 -ε 3 η q -ε 3 ε 2 ε 3 η q -ε 1 -ε 2 ε 1 η q          p q r    = R q 2 v 2 (3.59)
Quaternions and Euler angles can be translated to each other [START_REF] Fossen | Handbook of marine craft hydrodynamics and motion control[END_REF].

• Splitting approach

This approach is commonly-used in Molecular Dynamics (MD). The basic rotations 3.52-3.54 are modified to consecutive rotations R 1 -R 5 as follows:

R 1 = R x ( 1 2 φ) = R x ( 1 2 ∆t π 1 I 1 ) R 2 = R y ( 1 2 θ) = R y ( 1 2 ∆t π 2 I 2 ) R 3 = R z (ψ) = R z ( 1 2 ∆t π 3 I 3 ) R 4 = R y ( 1 2 θ) = R y ( 1 2 ∆t π 2 I 2 ) R 5 = R x ( 1 2 φ) = R x ( 1 2 ∆t π 1 I 1 ) (3.60)
where I 1 , I 2 and I 3 are the moment of inertia; π 1 , π 2 and π 3 are the corresponding components of the angular momentum π s . R x , R y and R z are the rotation matrices in 3.52-3.54. Quaternions can also be used instead of Euler angles for efficiency.

This approach avoids the gimbal lock phenomenon by splitting the angle increment into two more steps. The symmetric decomposition insures time-reversibility, and conserves total angular momentum [START_REF] Dullweber | Symplectic splitting methods for rigid body molecular dynamics[END_REF]. They are used to update the angular momentum and orientation matrices in OpenFoam:

π s = R 1 π s π s = R 2 π s π s = R 3 π s π s = R 4 π s π s = R 5 π s (3.61) Q s = Q s R T 1 Q s = Q s R T 2 Q s = Q s R T 3 Q s = Q s R T 4 Q s = Q s R T 5 (3.62)

Numerical integration techniques

The rigid body equations should be integrated to obtain the motion parameters of the ship. Traditional integration techniques such as Runge-Kutta methods, predictorcorrector schemes, and Gragg-Bulirsch-Stoer extrapolation are commonly-used. The Fifth-order Cash-Karp embedded Runge-Kutta scheme has been used in foam-extend [START_REF] Gatin | Enhanced coupling of solid body motion and fluid flow in finite volume framework[END_REF][START_REF] Jasak | Monolithic coupling of the pressure and rigid body motion equations in computational marine hydrodynamics[END_REF]. The symplectic-type methods are adopted in OpenFoam, which is proven more stable and enables the use of larger time steps [START_REF] Dullweber | Symplectic splitting methods for rigid body molecular dynamics[END_REF]. They are based on the splitting approach of the rotation matrices. Three numerical integration methods are implemented, the Verlet/leapfrog method, Newmark method and Crank-Nilcolson method, where the combination of the verlet/leapfrog with the splitting method is the standard symplectic method.

• Verlet/leapfrog method This method is normally used in Molecular Dynamics (MD), where a molecule can be described by a center of mass position vector q s , an orientation matrix Q s , and vectors for its linear and angular momentum (p s and π s ).

In the first step we need to calculate linear force vectors f n and torques τ n at time t n . Then we can start the integration and propagate the angular and linear momenta from t n to t n+1 :

π n+1/2 s = π n s + 1 2 ∆tτ n (3.63) p n+1/2 s = p n s + 1 2 ∆tf n (3.64)
Then we move the center of mass position a full time step:

q n+1 s = q n s + ∆tp n+1/2 s m (3.65)
Now we apply consecutive rotations R 1 to R 5 to all angular momenta and update all orientation matrices according to 3.61 and 3.62. After obtaining q n+1 s and Q n+1 s we can now calculate the forces f n+1 and torques τ n+1 at time t n+1 and propagate the angular and linear momenta another half time step:

π n+1 s = π n+1/2 s + 1 2 ∆tτ n+1 (3.66) p n+1 s = p n+1/2 s + 1 2 ∆tf n+1 (3.67)
This is the end of one integration step. This procedure is able to maintain the symplecticity of the integration.

• Newmark method

The Newmark method is widely used in numerical evaluation of the dynamic response of structures and finite element analysis. After the force vectors f n and torques τ n at time t n are calculated, the linear and angular momentum are calculated from t n to t n+1

π n+1 s = π n s + γ N τ n+1 + (1 -γ N )∆tτ n (3.68) p n+1 s = p n s + γ N f n+1 + (1 -γ N )∆tf n (3.69)
Both the old and new forces and torques are considered in this procedure. It is proven that γ N should take the value 0.5 to avoid the spurious damping effect [START_REF] Newmark | A method of computation for structural dynamics[END_REF]. The position vector is calculated according to:

q n+1 s = q n s + ∆tp n s + ∆t 2 (0.5 -β N )∆tf n + ∆t 2 β N ∆tf n+1 (3.70)
The value of β N normally equals 0.25. The orientation is corrected at last using 3.60.

• Crank-Nilcolson method

The Crank-Nilcolson procedure is similar with the Newmark method. The linear and angular momentum are first calculated:

π n+1 s = π n s + ∆t[α 1 τ n+1 + (1 -α 1 )τ n ] (3.71) p n+1 s = p n s + ∆t[α 1 f n+1 + (1 -α 1 )f n ] (3.72)
Then the position is calculated as:

q n+1 s = q n s + ∆t[α 2 p n+1 + (1 -α 2 )p n ] (3.73)
The orientation is corrected at last using 3.60. The values of α 1 and α 2 are taken to be 0.5.

6DoF solution procedures

Based on the two coordinate frames in Figure 3.8, the overall 6DoF computation procedure in OpenFoam can be summarized as:

1. Calculate the hydrodynamic forces and moments in the earth fixed system from the flow solutions;

2. Transform the forces and moments, including gravity forces, into body fixed system. Equations 3.51 are solved to calculate the linear and angular accelerations. Add constraints and restraints if necessary;

3. Integrate the accelerations to linear and angular momentum, new position and orientation using the numerical integration techniques;

4. Transform the orientation back to the earth system;

5. Calculate new velocities/positions of the mesh points on the solid boundary;

6. Update the mesh positions using the dynamic mesh solver in Section 3.4.1;

7. Current iteration is complete. New flow solutions will be calculated and the next iteration starts.

Enhanced coupling between solid and fluid

The 6DoF motion simulations is within the area of Fluid Structure Interaction (FSI). Two distinct FSI approaches are possible, the monolithic and partitioned approach. The former brings the fluid and solid physics together into a single set of equations to be solved, leading to a single, fully dependent solver that cannot be implemented directly. The latter separates the two solvers through the use of a coupling scheme. This scheme is based on successive solutions of the fluid and structure solvers: the coupling is loose/weak if only one solution per time step is required for each field, and strong if an iterative procedure is used to ensure the convergence of the coupled solution. Its best advantage is that it is a modular approach that greatly facilitates the implementation of different solvers and allows distributed computations [START_REF] Song | A partitioned coupling scheme extended to structures interacting with high-density fluid flows[END_REF].

The original 6DoF motion solver in OpenFoam uses the loose coupling approach.

The fluid solution is passed on to the motion solver libraries only once at the end of every time step to determine its position and orientation [START_REF] Chow | Strongly coupled partitioned six degreeof-freedom rigid body motion solver with Aitken's dynamic under-relaxation[END_REF].

The major drawback of the partitioned approach is that when high density fluids are involved, convergence is no longer guaranteed, and divergence is generally observed, regardless of the chosen time step and discretization scheme [START_REF] Van Brummelen | Added mass effects of compressible and incompressible flows in fluid-structure interaction[END_REF][START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF]Förster et al., 2007). It has been reported that this divergence is dependent on the size of the fluid domain [START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF], and the density ratios between the fluid and the structure (Förster et al., 2007). [START_REF] Song | A partitioned coupling scheme extended to structures interacting with high-density fluid flows[END_REF] stated that the divergence is best explained in terms of the ratio between the inertia and the integration of the pressure force on the structure, denoted as the added mass effect. This "artificial" added mass effect originates from the artifact of the loose coupling between the CFD and 6DoF solvers, which will reduce the stability of the solver and results in divergence.

To solve this problem, several approaches have been proposed [START_REF] Hughes | Lagrangian-Eulerian finite element formulation for incompressible viscous flows[END_REF]; [START_REF] Fernández | A projection algorithm for fluid-structure interaction problems with strong added-mass effect[END_REF]; [START_REF] Song | A partitioned coupling scheme extended to structures interacting with high-density fluid flows[END_REF]), including the semi-implicit discretization [START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF], Aitken's under-relaxation [START_REF] Küttler | Fixed-point fluid-structure interaction solvers with dynamic relaxation[END_REF], new partitioned procedures based on Robin-type transmission conditions [START_REF] Badia | Fluid-structure partitioned procedures based on Robin transmission conditions[END_REF], interface Newton-Krylov solver [START_REF] Michler | An interface Newton-Krylov solver for fluid-structure interaction[END_REF], multigrid technique [START_REF] Van Brummelen | Space/time multigrid for a fluid-structure-interaction problem[END_REF], and modified mass approach [START_REF] Song | A partitioned coupling scheme extended to structures interacting with high-density fluid flows[END_REF][START_REF] Tezduyar | Finite element methods for flow problems with moving boundaries and interfaces[END_REF][START_REF] Tezduyar | Space-time finite element techniques for computation of fluid-structure interactions[END_REF][START_REF] Connell | Flapping dynamics of a flag in a uniform stream[END_REF][START_REF] Van Brummelen | Partitioned iterative solution methods for fluidstructure interaction[END_REF]. [START_REF] Chow | Strongly coupled partitioned six degreeof-freedom rigid body motion solver with Aitken's dynamic under-relaxation[END_REF] implemented the second order Adams-Bashforth-Moulten method instead of the original leapfrog scheme in OpenFoam along with the Aitken's underrelaxation. [START_REF] Dunbar | Development and validation of a tightly coupled CFD/6-DOF solver for simulating floating offshore wind turbine platforms[END_REF] implemented a solver which tightly coupled the fluid and 6DoF equations of motion using sub-iterations in the SIMPLE loop (nOuter-Correctors in Figure 3.2) and Aitken's relaxation to eliminate the artificial added mass instability. [START_REF] Gatin | Enhanced coupling of solid body motion and fluid flow in finite volume framework[END_REF] calculated the 6DoF motion equations after each pressure correction step further within the PISO loop (nCorrectors in Figure 3.2). [START_REF] Devolder | A review of the implicit motion solver algorithm in OpenFOAM to simulate a heaving buoy[END_REF] included the added mass into the leapfrog scheme using the modified-mass approach [START_REF] Song | A partitioned coupling scheme extended to structures interacting with high-density fluid flows[END_REF], but the added mass is constant and should be adjusted with the cases. [START_REF] Jasak | Monolithic coupling of the pressure and rigid body motion equations in computational marine hydrodynamics[END_REF] introduced a monolithic approach in OpenFoam where the updated rigid body boundary velocity is added directly into the pressure equation. Some of these schemes will be implemented to overcome the artificial added mass effect.

Chapter 4

Validation and verification studies

Introduction

In this chapter, numerical simulations of ships will be conducted for validations and verifications. These results will give helpful suggestions for the following works in this thesis.

• The boundary conditions in OpenFoam are tested to find the most appropriate conditions of their usage.

• Mesh and time step convergence studies are carried out to find a good meshing strategy and an appropriate time step. Different turbulence models are tested to characterize their influences on the simulation results.

• Prescribed ship motions and 6DoF motions are validated to test the influences of the dynamic mesh techniques.

Boundary condition tests

Boundary conditions (BCs) should be physically realistic to guarantee the stability and convergence of the solution. A simple interface-piercing cylinder case (Figure 4.1) is designed to select the appropriate ones. Many types of BCs have been implemented in OpenFoam, the possible candidates will be tested one by one. The basic BCs used can be found in Table 4.1. The SST k-ω turbulence model is used for testing in this case. The time derivative term is discretized using the implicit Euler scheme. The convection term in the momentum equation is discretized with second order linear upwind scheme, and the convection term in the VOF equation is discretized with second order TVD (Total Variation Diminishing) scheme with van Leer's flux limiter. Second order central differencing scheme with non-orthogonal correction is used for the discretization of the diffusion term. One outer correctors are used per time step, and two pressure corrections are employed in this case. All equations are solved to the tolerance of O(-8).

Velocity boundary conditions

Possible outlet BCs of the velocity can be found in Table 4.2. The "outletPhase-MeanVelocity" is designed specially for towing tank ship simulations in OpenFoam to maintain the outlet water level as the inlet. However, this BC is found to be instable under certain conditions. It can cause weird free surface shapes and lead to divergence at last (Figure 4.2). The "pressureInletOutletVelocity" BC is normally used with the "fixedValue" inlet pressure BC, which makes its convergence curve completely different from that of the "zeroGradient" type BCs. The other two BCs have similar convergence histories. The "inletOutlet" BC is normally the same as "zeroGradient", but it switches to "fixedValue" if the velocity vector next to the boundary aims inside the domain (backward flow), where that "fixedValue" is set as the inlet value. Therefore, as shown in Figure 4.3, "inletOutlet" BC doesnot allow backflow, which seems plausible for some cases. However, it is found to cause problems when the boundary moves using dynamic mesh techniques. "zeroGradient" BC is comparably more stable and suitable for most conditions.

Pressure boundary conditions

The pressure BCs should be in accordance with the velocity BCs. The candidates can be found in Table 4.3. "fixedFluxPressure" sets the pressure gradient such that the flux on the boundary is specified by the velocity BC. When it is used with a "zeroGradient" type velocity BC, large fluctuations of the solution will be observed (Figure 4.4 and 4.5). Although the convergence history seems alright, "zeroGradient" at the inlet and outlet should be avoided, which result in an ill-posed problem. The combination of "zeroGradient" at inlet and "fixedValue" at the outlet should be robust and is always suggested by some researchers. However, it is found to cause the drop of the water level during the ship simulations. Therefore, "fixedFluxPressure" and "zeroGradient" are used since the pressure has already been set by the "totalPressure" BC at the top.

Other boundary conditions

The candidates of alpha.water, k and nut are shown in Tables 4.4∼4.6. It can be seen in Figure 4.6 that the solution is not sensitive to them. For the wall (Table 4.7), "movingWallVelocity" and "fixedValue" have equivalent convergence curves as long as the velocity values are the same. For the sides (Table 4.7), "symmetry", "symmetryPlane" and "slip" are the same. However, the "symmetryPlane" BC is not suggested because it is sensitive to the smoothness of the boundary, which may cause problems for post-processing. A wall BC is also designed for the sides in Table 4.8. It will not influence the solution as long as the domain is large enough.

Convergence study 4.3.1 Mesh convergence

The mesh convergence study is conducted for a scaled KCS (KRISO container ship) model, whose geometrical parameters are shown in Table 4.9. Figure 4.7 presents the computational domain. The regions of the free surface, bow and stern are refined to resolve the water-air interface and complex geometries. Only a half domain in the y direction is simulated by using the symmetry boundary condition at the midPlane.

Other BCs can refer to the analyses in Section 4.2. The turbulence model and schemes adopted in this case are the same as in Section 4.2. Three outer correctors are used per time step, and one pressure correction are employed in this case. Three types of grids are designed in Table 4.10. The mesh is generated using the snappyhexMesh utility. A background mesh is first generated and then the meshes at the target geometry are split and snapped. Boundary layer meshes are added at last, which are controlled to make the non-dimensional wall distance in the range 30 < y + < 200. A refinement ratio √ 2 is used for the background mesh as suggested by the ITTC (2017). Since the mesh near the hull is refined to a certain level based on the background mesh, the final mesh ratio is larger than √ 2 in Table 4. The convergence ratio (R C ) and order of accuracy (p C ) can be calculated according to:

R C = ε 21 C ε 32 C (4.1) p C = ln(|ε 32 C /ε 21 C |) ln(r f ) (4.2)
= C 2 T -C 1 T and ε 32 C = C 3 T -C 2
T are the change of resistance coefficients (C T ). The subscripts 1, 2 and 3 denote the the fine, medium and coarse meshes.

The experimental result of the total resistance coefficient (C exp T ) is also listed in Table 4.10 for comparison, which is close to the simulation results. The frictional and pressure resistance coefficients are also calculated and compared with the ITTC 1957 friction line in Table 4.11. Three possible conditions can be diagnosed: Hence the simulations in this study are monotonically convergent. The Richardson extrapolation is adopted to define the extrapolated resistance coefficient (C 21 ext ), approximate relative error (e 21 a ), extrapolated relative error (e 21 ext ) and fine-grid convergence index (GCI 21 f ine ). The numerical uncertainty for the total resistance in this work is 1.348%.

• 0 < R C < 1 Monotonic convergence; • R C < 0 Oscillatory convergence; • R C > 1 Divergence.
C 21 ext = r p C f C 1 T -C 2 T r p C f -1 (4.
3)

e 21 a = C 1 T -C 2 T C 1 T (4.4) e 21 ext = C 21 ext -C 1 T C 21 ext (4.5) GCI 21 f ine = 1.25e 21 a r p C f -1 (4.6)
Figure 4.8 shows the wave elevation contours of the three cases. The waves are resolved better with a finer mesh. Considering both the accuracy and the simulation time, the mesh generation strategy 2 is used for the following studies. The wave profiles of the case 2 agree well with experiments in Figure 4.9.

Time step convergence 4.3.2.1 Steady case

The time convergence study is carried out for the KVLCC2 model (Table 4.12). The computational domain and the meshes near the bow and stern can be found in Figure 4.10, with the total cell number 818,372. The turbulence model and schemes in this case are the same as in Section 4.2. Five outer corrections and two pressure corrections are used. As in Figure 4.11, four time steps are designed. Because of the use of the PIMPLE scheme, large time steps can be adopted. The maximum Courant number in this study reaches 60 and the case is still convergent. 20 processors are used for this case. It can be seen that the simulation time is greatly reduced. The results of the four cases seem all reasonable. The frictional and pressure resistances are also calculated and compared with the ITTC 1957 friction lines in Table 4.13. However, the accuracy of the case 4 (∆t = 1 × 10 -1 s) is clearly not enough. Although the trend is correct, large fluctuations can be found in the convergence history and the wave profiles. A very small time step (∆t = 1 × 10 -4 s) can capture small changes during the convergence, but the simulation time is enormously increased. As a compromise, time steps 1 × 10 -2 and 10 -3 s are preferred in this work. Steady solutions are sought in this case. Another transient case with different time steps are also studied in the following section.

Unsteady case

The water entry phenomenon of a free falling sphere is simulated for the time step convergence of a transient case. The computational domain can be found in Figure 4.12, with the sphere placing 1 m above the free surface. The density of the sphere is 500 kg/m 3 (half of the water density). The total grid number is 665,409. Six outer correctors and four pressure correctors are used since transient behaviors are important in this case, with five different time steps 1 × 10 -2 , 5 × 10 -3 , 1 × 10 -3 , 5 × 10 -4 , 1 × 10 -4 s. No turbulence model is adopted in this case. Since no wave absorption techniques are used, the boundaries are placed far from the sphere, and the mesh size are incremented from the sphere so that waves can be diffused.

Figure 4.12: Computational domain of the water entry case of a free falling sphere. R s is the radius of the sphere.

The displacement and vertical force histories are validated (Figure 4.13). The overall results agree with the work of [START_REF] Shen | Numerical simulation of sphere water entry problem based on VOF and dynamic mesh methods[END_REF]. The converge histories of the displacement and vertical force for the time steps 1 × 10 -3 , 5 × 10 -4 and 1 × 10 -4 s are close. The impact forces when the sphere first touches the free surface are predicted close to the work of [START_REF] Shen | Numerical simulation of sphere water entry problem based on VOF and dynamic mesh methods[END_REF]. Discrepancies can be seen after about 4 s, where our results diffuse more. This may be caused by the fact that Shen & Wan (2011) used sponge layers around the sphere, which would have less influences of wave reflections. Although the mesh used here is designed specifically to reduce this phenomenon, the waves cannot be fully diffused by this method. However, both works converge to the same values at last. Small time steps (≤ 1 × 10 -3 s in this case) are better to be used for transient simulations.

Sensitivity analyses of turbulence models

Three two-equation turbulence models, k-ε, k-ω and SST k-ω, are selected in this part. The mesh and BC setups are the same as in Figure 4.10. One outer correction and two pressure corrections are used. Their convergence histories in Figure 4.14 are similar. The SST k-ω model produces the most accurate resistance coefficient among them, which will be adopted in the following studies. The frictional and pressure resistances are also compared in Figure 4.14. It can be seen that the viscous force contributes more to the total resistance. The pressure parts of the three models are nearly the same. The difference mainly lies in the prediction of the frictional force.

Dynamic mesh validations 4.4.1 Prescribed ship motion

To validate the prescribed dynamic mesh techniques, two cases are designed in this part (Figure 4.15). In the first case, the ship is fixed and a velocity is added for the fluids at the inlet as the former studies. In the second case, the ship together with the computational domain is moved. Three outer corrections and two pressure corrections are used. Their convergence histories are similar. Both cases can correctly predict the ship resistance. The difference between them mainly lies in the pressure resistance, which can be caused by the added mass effect when ship moves. A 6DoF motion solver based on the quaternion approach (Section 3.4.3) is implemented in this section with the tightly coupling scheme in [START_REF] Dunbar | Development and validation of a tightly coupled CFD/6-DOF solver for simulating floating offshore wind turbine platforms[END_REF]. A KCS model (Table 4.9) with free sinkage and trim is simulated in this case (Du et al., 2018). The case setups are shown in Table 4.14, where six different ship velocities are designed. The SST k-ω model is used in this case. The grid number is 1,074,514. The numerical schemes, corrector numbers and convergence criterion are the same as the water entry case (Section 4.3.2.2). The resistance coefficient (C T ), sinkage (σ s ) and trim (τ t ) are calculated and compared with experiments [START_REF] Larsson | Numerical ship hydrodynamics: an assessment of the Gothenburg 2010 workshop[END_REF]. The resistance coefficient is calculated as:

6DoF motion

C T = F x 1 2 ρU 2 in S w (4.7)
where F x , U in and S w are the total resistance, ship speed and wetted surface area at rest. In Figure 4.17, the resistance coefficients are underestimated for F r ≤ 0.260. The sinkage is underestimated for F r ≤ 0.195 and overestimated for F r ≥ 0.260. The trim prediction is fairly well, with only the case F r = 0.108 underestimated and F r = 0.282 overestimated. The largest error of the three parameters are about 3.0%, 6.4% and 7.1% respectively. Overall, the implemented 6DoF model can predict the ship resistance and motion states with satisfactory accuracy.

Figure 4.17: Comparison of the resistance coefficient (C T ), sinkage (σ s ) and trim (τ t ) with experimental data [START_REF] Larsson | Numerical ship hydrodynamics: an assessment of the Gothenburg 2010 workshop[END_REF]. F r is the Froude number based on the ship length.

Concluding remarks

In this chapter, validation and verification studies were carried out to test the validity and accuracy of the numerical solvers.

Tests of various boundary conditions in OpenFoam were carried out. Suggestions for the conditions of their usage were given to ensure the convergence and validity of the simulations. Velocity and pressure boundaries were found to influence significantly on the computations, while others have small influences.

Two open ship models, KCS and KVLCC2 were used. The mesh convergence was tested based on the suggested procedures of ITTC using three mesh resolutions. The second one is as a result the most appropriate considering the accuracy and the computational cost.

The time step convergence was conducted for the steady state of the KCS and the unsteady case of a free falling sphere. Time steps 1 × 10 -2 and 10 -3 s were suggested for the former and small time steps (≤ 1 × 10 -3 s) for the latter. The SST k-ω turbulence model is preferred based on the sensitivity analyses.

Two tests were carried out for the dynamic motion solvers, a prescribed one and a 6DoF one. These results can prove the validity and accuracy of these models and give suggestions for the following studies.

Chapter 5

Influences of environmental effects

Introduction

In this chapter, the environmental effects on ships are simulated using the numerical solvers introduced and validated in the former chapters. The results will be compared with the experimental data. The following three effects are investigated:

• The confinement effect induced by the channel bank and bottom;

• Head-on encounter between two inland vessels in the confined waterway;

• Interaction between the ship and bridge piles during the crossing.

Confinement effect

Experimental details

The experiments were carried out in the towing tank of the University of Liège (Figure 5.1). The basin has been adapted with an aluminum bottom and PVC banks. Two convoys were tested in this part, which consists of a pusher and one/two barges (Figure 5.2). They are promising in increasing the payload by at least 20% and decreasing the costs by at least 30% for the European inland waterway transport system. The body lines of the pusher and the barge are shown in Figures 5.3 and 5.4. Their geometrical parameters are presented in Table 5.1. Reduced models with the scale ratio 1/25 are used during the experiments. The barge and pusher are connected rigidly as shown in Figure 5.1(b).

The experimental conditions are all shown in Figure 5.5 and listed in Table 5.2. The bank slope is 2:1 to resemble the real waterways. Three water depths, bank positions and two draughts (empty and charged) are designed to account for various channel and transport conditions. The weights of the pusher with two different draughts are 54.77 and 144.11 kg; while those of the barge are 102.81 and 265.22 kg. The testing speed is in the range of 0∼1 m/s, which means that all our tests are within the subcritical and critical regions (Figure 2.10). The Schijf limited speed is calculated for reference:

F r Schijf h = 8 cos 3 π 3 + arccos(1 -A s /A c ) 3 (5.1)
where A s and A c are the cross section areas of the ship and channel. The Schijf limiting speed is for self-propelled ships. Since the ship models in our experiments are towed, some cases exceed this limit. The parameter h w /T d is within 3∼6, A c /A s 4.263∼44.211 and B c /B s 2.305∼7.368, which show that our tests cover a wide range of confinement conditions.

W b T d h w V max F r max h V s V sub c h w /T d A c /A s B c /B

Computational setups

The computational domain is shown in Figure 5.6, with two ship lengths upstream and three downstream. Overall BC setups are the same as in Section 4.2, except the channel bottom is set to wall. 'OutletPhaseMeanVelocity' BC is used to maintain the water level. The mesh information for water depth h w =0.24 m is shown in Table 5.

Results and analyses

The resistance results of the convoy 1 and 2 are shown in Figures 5.7 and 5.8 respectively. The numerical simulations agree well with experiments. The cases with smaller water depth, narrower channel width and larger ship draught have higher advancing resistance since they are more confined in the waterway. Return flows appear around the hull and added mass effect induces higher resistance. The ship resistance is plotted against the blockage ratio (m b = A s /A c ) in Figure 5.9.

For the same velocity, ship resistance shows linear relationship with the blockage ratio. And the curve seems to have a larger slope with a higher ship speed. This means that a ship with higher speed is more sensitive to the confinement change of the waterway. Even for the same ship, its performance changes when the draught is different. In Figure 5.10, the resistance changes more intensively with a larger draught value. Therefore, when the ship is navigated in a confined waterway, its performance should be reevaluated with different loadings. Special attention should be paid to high loading and confinement conditions since the ship experiences large resistance and becomes hard to maneuver. The ship waves of the convoy 2 are more complicated since another barge creates additional waves. Because the vessel is longer, the interfering region of the ship waves is also larger, creating more complex wave patterns. The wave angles are measured in Figure 5.13, which seem irrelevant with the channel width. This agrees with the traditional viewpoint. According to Section 2.4.3, the wave angle in shallow water is only influenced by the depth-based Froude number. However, it is clearly affected by the draught here. The vessel with a larger draught has a larger wave angle. In Figures 5.11 and 5.12, the wave reflections appear earlier in the channel and create more interferences in the downstream. To account for the draught, the original depth-Froude number is suggested to modify using the ship draught h T = h w -T d :

F r T h = V √ gh T (5.2)
With this modification, the equations for the curve in Figure 5.13(b) are changed to:

sin α k = 1 + 2k w h T sinh -1 (2k w h T ) 3 -2k w h T sinh -1 (2k w h T ) (5.3)
α k = arcsin(1/F r T h ) (5.4)
k w h T relates with F r T h as:

tanh k w h T k w h T [3 -2k w h T sinh -1 (2k w h T )] = F r T 2 h
(5.5) The three dotted lines in Figure 5.13(a) denote the theoretical wave angles. Line 1 is the original wave angle which is irrelevant with the draught. Line 2 and 3 are new values calculated based on the modified water depth, which agree better with the simulation results. Water level increases at the bow and maintains a long wave along the hull. A maximum wave elevation can be observed at the connection part between the two vessels. This is caused by the superposition of the stern wave of the front part and the bow wave of the back part. The cases with larger draught seem to have more intense wave changes, and the channel confinement further amplifies this effect. When the channel width is smaller, the ship-generated waves become higher at the bow and lower along the hull. The ship waves are more complicated for the convoy 5.3 Head-on encounter between two ships

Experimental details

The head-on encounter between the convoy 2 and an inland tanker is tested in this part (Figure 5.15). During the experiment, the convoy 2 is driven by the trolley and the tanker by a rope wrapped around a drum at the end of the basin. The hull form of the tanker can be seen in Figure 5.16. A motor providing constant power rotates the drum, allowing the tanker to move at a constant speed of 0.50 m/s with a draught 0.07 m and weight 148.52 kg. A system consisting of pulleys attached to both ends of the channel bottom and a tensioned rope connecting the front to the back of the tanker is used to keep a straight path. As illustrated in Figure 5.17, the two vessels cross at a distance of 0.37 m with the bank width 1.44 and 2.88 m. The center of the canal is adjusted to the crossover axis of the models (0.185 m from the edge of each model). The water depth is 0.18 m and the tests were carried out for the convoy 2 under light and loaded conditions.

Detailed case configurations can be found in Table 5.4.

Computational setups

The computational domain is designed as the experiments (Figure 5.18). ACMI technique is used to move the ships. Two moving parts are created separately for the convoy 2 and the tanker. The meshes can be observed in Figure 5.19. Detailed mesh numbers for each case are shown in Table 5.4. Five outer corrections and two pressure corrections are used. 90 processors are used in this part and the overall simulation durations are about 28 cpu hours for a1 and a2, 47 cpu hours for b1 and b2.

Results and analyses

To characterize the passing process, a parameter ξ representing the relative position between the two ships is defined as in Section 2.4.5.2:

ξ ss = x ss /L m (5.6)
where x ss is the longitudinal distance between the ships. L m = (L convoy2 +L tanker )/2 is the averaged length of the convoy 2 and the tanker. Figure 5.20 demonstrates the three most important positions. When ξ ss = -1, the head of the two vessels meet.

When ξ ss = 0, the centers of the two vessel coincides at the same longitudinal (x) position. When ξ ss = 1, the head-on encounter ends. The positions within -1 < ξ ss < 1 during the crossing will be emphasized. Only the average resistance values of the convoy 2 were calculated during the experiments. In Figure 5.21, the simulation values agree well with experiments. The resistance histories of the simulation are plotted in , whose values are modified using the average resistance to concentrate on the crossing process. It can be observed that the longitudinal force caused by the ship-ship interaction is positive in the first stage (ξ ss < 0) and negative in the second stage (ξ ss > 0), which means the longitudinal force always drags the vessels. Large ship speed, draught and small channel dimension will enhance the interaction during the crossing. The ship waves at the position ξ ss = 0 are specially shown in Figures 5.25 and 5.26, which are clearly modified by the appearance of another ship. The superposition of these waves creates a very complicated wave pattern in the confined channel. In Figure 5.25, the ship waves of the tanker is disturbed and nearly not observable any more since its velocity is small. The drawdown on the channel banks can also be observed. It is responsible for the erosion of the soil and the failure of the bank structure. As in Section 5.2.3, higher velocity, greater draught and smaller confinement area will induce greater wave elevation changes of the free surface, and also larger drawdown on the banks. Since the velocity of the convoy 2 is higher, the drawdown on its side is greater than that of the tanker. The wave profiles at the channel center are extracted for comparison. The sampling line is depicted in Figure 5.27. The wave elevation changes are clearer with the ship speed, draught and channel width by the wave profiles ). Higher vessel speeds make the wave profiles higher in front of the vessel and lower in the middle. The longer and faster vessel dominates the flow and wave fields. During inland shipping, vessels will inevitably pass bridges, where the bridge piles exert stronger confinement apart from the channel banks. The ship dynamics and maneuverability change during this process. The increase of the squat may lead to grounding and the asymmetrical moment on the hull may result in collision.

In this part, the convoy 2 passing the bridge piles in the confined waterway is tested and simulated. The experimental setups can be found in Figure 5.29 and Table 5.5.

All the tests are done with a channel width 1.44 m. Two pile distances, two draughts and three water depths are tested within the subcritical conditions (F r h < 1). The piles are composed of two concrete blocks of dimensions 40 cm × 16 cm whose edges are cut.

Computational setups

The computational domain is shown in Figure 5.30. Symmetry BC is used to reduce the domain size. ACMI technique is adopted here, whose positions are marked as red lines. The mesh is shown in Figure 5.31, where the free surface, ship hull and bridge regions are refined to capture the water level and flow fields. The mesh numbers are listed in Table 5.5. Five outer corrections and two pressure corrections are employed. Each simulation uses 52 processors and takes about 12 cpu hours to finish.

Figure 5.32: Relative position (ξ sb ) between the convoy 2 and bridge piles.

Results and analyses

To characterize the passing process, a parameter ξ representing the relative position between the ship and the bridge is defined: When ξ sb = 0, the centers of the vessel and the pile coincide at the same x position.

ξ sb = x sb /L m (5.7)
When ξ sb = 1, the passing ends. The positions within -1 < ξ sb < 1 will be emphasized.

The resistance during the passing process is averaged. Figure 5.33 shows the resistance comparison between experiments and simulations, which agrees well. With a smaller water depth, the vessel is more confined, leading to a higher resistance, which is also the case for the draught.

The detailed resistance change with pile distances (D b ), vessel speed (V ) and water depth (h w ) are further analyzed for the simulations. In Figure 5.34, with two different pile distances, the overall trends of the resistance are similar. The resistance changes intensively and has a minimum value near the position ξ sb = 0. This means that the ship is the most unstable at this position, which should be paid attention to during the real maneuvering. A window with shorter period of time are also zoomed to demonstrate the convergence of the resistance within several time steps. The resistance change with different vessel speeds is shown in Figure 5.35. Clearly a higher speed results in a larger resistance on the vessel. The minimum resistance appears at the same position (ξ sb = 0) regardless of the vessel speeds. The influence of the water depth is similar. A smaller water depth has a larger resistance, but the positions of the minimum value seem unchanged. The most dramatic change always lies at the position (ξ sb = 0) where the center of the vessel and the piles coincide during the passing process. Two distances between the piles 0.7 m and 0.8 m are designed. Since their values are small, no significant differences can be observed in Figure 5.41. However, they are clear enough in the extracted profiles (Figure 5.42). The pile will create a more confined condition besides the channel and the water level will be elevated. A smaller pile distance will create a more confined condition, which will make the water level higher. Nevertheless, the pile distance will not influence the overall wave pattern because the wave crests and troughs are at the same positions. This is also the reason why no significant differences can be seen in Figure 5.41 except the change of the isoline densities.

Concluding remarks

In this chapter, the influences of the confinement, ship-ship interaction and shipstructure interaction were simulated and compared with experimental data.

For the confinement effect, two convoy models in a fully-confined channel were investigated with three water depths, three channel widths, two draughts and several vessel speeds. A wide range of real conditions were covered by these combinations. The simulations were conducted without sinkage and trim. The advancing resistance agrees well with experiments. The relationship between the ship resistance and the blockage ratio is found to be linear. Higher speeds and larger draughts make the slope of the linear relation higher, which means that the resistance changes more intensively under these conditions. The theory of the Kelvin angle was modified to account for the draught, which turns out to coincide with the real conditions. Since the convoys consist of the barge and pusher, these wave patterns superimpose with each other and create a very complex waves in the confined channel. The head-on encounter between the convoy 2 and the inland tanker was simulated using the ACMI technique under the conditions of two channel widths, two draughts and two vessel speeds. The resistance change was analyzed, which is found to always drag the ships during the crossing. The ship wave analysis shows that larger vessel speed, draught and smaller channel width contribute more to the wave elevation and drawdown on the banks.

The ship passing bridge piles was simulated to investigate the influence of the appearance of structures in the confined channel. The case design included one channel width, two bridge pile distances, two draughts, three water depths and several ship speeds. A parameter ξ sb similar with the ship-ship interaction was defined to characterize the passing process. The resistance shows a minimum value at the position ξ sb = 0, regardless of the conditions. The ship waves near the bridge pile become complicated since the waves are reflected by the wall of the pile and superimpose with the original wave patterns. Two wave crests are observed there.

Overall, the advancing resistance, ship wave pattern, Kelvin angle were emphasized for these studies. The findings provide important clues for the ship motions in the confined waterway and ship interactions.

Chapter 6

Ship maneuvering analyses

Introduction

In this chapter, ship maneuvering analysis is carried out based on the virtual captive model tests, the following aspects are emphasized:

• Commonly-used maneuvering models are analyzed. Particularly, the formulations and procedures of the MMG (Mathematical Maneuvering Group) model is presented;

• Virtual captive model tests are conducted to calculate the hydrodynamic coefficients of the MMG model;

• System-based maneuvering simulations are carried out.

Maneuvering models

System-based maneuvering simulations can predict the ship trajectories, yaw angles, etc. within an acceptable time range, which are important for ship motion prediction and evaluation. As introduced in Section 2.4.6, it solves the rigid body equations directly, where the hydrodynamic forces on the ship are calculated from the hydrodynamic derivatives obtained from either empirical formula, captive model tests or system identification methods. Therefore, the hydrodynamic derivatives play an important role for ship maneuvering simulations. A 3DoF rigid body equation is introduced in this chapter for simplification:

m( u -vr -x G r 2 ) = X(u, v, r, u, v, ṙ, δ, ...) (6.1) m( v + ur + x G ṙ) = Y (u, v, r, u, v, ṙ, δ, ...) (6.2)
I z ṙ + mx G ( v + ur) = N (u, v, r, u, v, ṙ, δ, ...) (6.3)
The RHS can be expanded using Taylor series:

X =X 0 + (u -U ) ∂X ∂u + v ∂X ∂v + r ∂X ∂r + u∂X ∂ u + v ∂X ∂ v + ṙ ∂X ∂ ṙ + δ ∂X ∂δ + 1 2! (u -U ) ∂ ∂u + v ∂ ∂v + r ∂ ∂r + u ∂ ∂ u + v ∂ ∂ v + ṙ ∂ ∂ ṙ + δ ∂ ∂δ (2)
X + . . .

+ 1 n! (u -U ) ∂ ∂u + v ∂ ∂v + r ∂ ∂r + u ∂ ∂ u + v ∂ ∂ v + ṙ ∂ ∂ ṙ + δ ∂ ∂δ (n) X + . . . (6.4) Y =Y 0 + (u -U ) ∂Y ∂u + v ∂Y ∂v + r ∂Y ∂r + u∂Y ∂ u + v ∂Y ∂ v + ṙ ∂Y ∂ ṙ + δ ∂Y ∂δ + 1 2! (u -U ) ∂ ∂u + v ∂ ∂v + r ∂ ∂r + u ∂ ∂ u + v ∂ ∂ v + ṙ ∂ ∂ ṙ + δ ∂ ∂δ (2)
Y + . . .

+ 1 n! (u -U ) ∂ ∂u + v ∂ ∂v + r ∂ ∂r + u ∂ ∂ u + v ∂ ∂ v + ṙ ∂ ∂ ṙ + δ ∂ ∂δ (n) Y + . . . (6.5
)

N =N 0 + (u -U ) ∂N ∂u + v ∂N ∂v + r ∂N ∂r + u∂N ∂ u + v ∂N ∂ v + ṙ ∂N ∂ ṙ + δ ∂N ∂δ + 1 2! (u -U ) ∂ ∂u + v ∂ ∂v + r ∂ ∂r + u ∂ ∂ u + v ∂ ∂ v + ṙ ∂ ∂ ṙ + δ ∂ ∂δ (2)
N + . . .

+ 1 n! (u -U ) ∂ ∂u + v ∂ ∂v + r ∂ ∂r + u ∂ ∂ u + v ∂ ∂ v + ṙ ∂ ∂ ṙ + δ ∂ ∂δ (n) N + . . . (6.
6) The derivatives are usually expressed as:

∂X ∂u = X u , ∂X ∂v = X v , ∂X ∂r = X r , ∂X ∂ u = X u, • • • (6.7)
Normally, only some terms in Equation 6.4-6.6 are retained. Early models only keep the first order terms, which are referred as linear models. Details of various models can be found in Section 2.4.6. Among them, the models of MMG [START_REF] Yasukawa | Introduction of MMG standard method for ship maneuvering predictions[END_REF] and [START_REF] Bibliography Abkowitz | Lectures on ship hydrodynamics-steering and manoeuvrability[END_REF] are the most commonly-used. They are emphasized in the following parts.

Abkowitz model

The Abkowitz mathematical model directly expands the forces and moment X, Y and N in equations 6. 1-6.3 as (Strøm-Tejsen & Chislett, 1966):

X =X * + X u u + X u ∆u + X uu ∆u 2 + X uuu ∆u 3 + X vv v 2 + X vvu v 2 ∆u + X rr r 2 + X rru r 2 ∆u + X vr vr + X vru vr∆u + X δδ δ 2 + X δδu δ 2 ∆u + X vδ vδ + X vδu vδ∆u + X rδ rδ + X rδu rδ∆u (6.8) Y =Y * + Y ṙ ṙ + Y v v + Y u ∆u + Y uu ∆u 2 + Y uuu ∆u 3 + Y v v + Y vvv v 3 + Y vrr vr 2 + Y vδδ vδ 2 + Y vu v∆u + Y vuu v∆u 2 + Y r r + Y rrr r 3 + Y vvr v 2 r + Y rδδ rδ 2 + Y ru r∆u + Y ruu r∆u 2 + Y δ δ + Y uδ ∆uδ + Y uuδ δ∆u 2 + Y δδδ δ 3 + Y vvδ v 2 δ + Y rrδ r 2 δ (6.9) N =N * + N ṙ ṙ + N v v + N u ∆u + N uu ∆u 2 + N uuu ∆u 3 + N v v + N vvv v 3 + N vrr vr 2 + N vδδ vδ 2 + N vu v∆u + N vuu v∆u 2 + N r r + N rrr r 3 + N vvr v 2 r + N rδδ rδ 2 + N ru r∆u + N ruu r∆u 2 + N δ δ + N uδ ∆uδ + N uuδ δ∆u 2 + N δδδ δ 3 + N vvδ v 2 δ + N rrδ r 2 δ ( 6

MMG model

The coordinate system of MMG model is shown in Figure 6.1. Similar with 6.1-6.3, the 3DoF maneuvering equation of MMG model with the added masses (m x , m y ) and added moment of inertia (J z ) is written as [START_REF] Yasukawa | Introduction of MMG standard method for ship maneuvering predictions[END_REF]:

(m + m x ) u -(m + m y )v m r -x G mr 2 = X (6.11) (m + m y ) vm + (m + m x )ur + x G m ṙ = Y (6.12) (I zG + x 2 G m + J z ) ṙ + x G m( vm + ur) = N m (6.13)
Differently, the forces and moments of the MMG model on the RHS are separated based on the hull, propeller and rudder as follows: 6.16) This separation allows for more accurate evaluation of each forces. Other forces induced by wind, waves, current, ship interactions, etc. can also be added in a similar way. It should be noted that this separation also brings errors since some phenomena are essentially coupled, such as the hull-ruder-propeller interaction, wave-current interaction, wave-wave interaction, etc. Corrections can be considered for such situations. The hydrodynamic coefficients in the MMG model are deduced from captive model tests in this chapter.

X = X H + X R + X P (6.14) Y = Y H + Y R (6.15) N m = N H + N R (

Hydrodynamic forces on the hull

The hydrodynamic forces on the ship hull are non-dimensionalized as:

X H = 1 2 ρL pp T d U 2 X H (v m , r ) (6.17) Y H = 1 2 ρL pp T d U 2 Y H (v m , r ) (6.18) N H = 1 2 ρL 2 pp T d U 2 N H (v m , r ) (6.19)
where the dimensionless hydrodynamic forces are calculated as the 1st and 3rd order polynomial functions of v m and r :

X H (v m , r ) = -R 0 + X vv v 2 m + X vr v m r + X rr r 2 + X vvvv v 4 m (6.20) Y H (v m , r ) =Y v v m + Y R r + Y vvv v 3 m + Y vvr v 2 m r + Y vrr v m r 2 + Y rrr r 3 (6.21) N H (v m , r ) =N v v m + N R r + N vvv v 3 m + N vvr v 2 m r + N vrr v m r 2 + N rrr r 3 (6.22)
The lateral velocity v m and yaw rate r are non-dimensionalized as v m = v m /U and r = rL pp /U respectively. Other forces and moments are non-domensionalized the same way as 6.17-6.19. The non-dimensionalization is based on the Prime system II in Table A.1 of the Appendix.

Propeller induced forces

Only the longitudinal force induced by the propeller is considered: 6.23) where the thrust deduction factor t P is assumed to be constant at a given propeller load for simplicity. The propeller thrust T P is written as:

X P = (1 -t P )T P (
T P = ρn 2 P D 4 P K T (J P ) (6.24) K T (J P ) = k 2 J 2 P + k 1 J P + k 0 (6.25)
where k 0 , k 1 , k 2 are calculated from the propeller open water tests.

During the maneuvering simulations, the propeller advance ratio J P is calculated as:

J P = u(1 -w P ) n P D P (6.26) 1 -w P 1 -w P 0 = 1 + {1 -exp(-C 1 |β P |)}(C 2 -1) (6.27)
where w P is the wake coefficient at propeller position, w P 0 is that in straight moving. β P = β -x P r is the geometrical inflow angle to the propeller. β is the hull drift angle at midship. x P is the longitudinal coordinate of the propeller position. β = tan -1 (-v m /u) is the hull drift angle at midship. C 1 and C 2 are the constants of wake characteristics determined from experiments.

Rudder induced forces

The rudder induced forces are written as: 6.30) where t R and a H are the steering resistance deduction factor and rudder force increase factor. x R and x H are the longitudinal coordinates of the rudder position and the acting point of the additional lateral force. They represent the interactions between hull and rudder. The rudder normal force F N is expressed as: 6.31) where ρ is the water density. A R is the profile area of the rudder movable part. f α is the rudder lift gradient coefficient. The rudder inflow velocity U R and the angle α R are calculated as: 6.33) where δ is the rudder angle. u R and v R are the longitudinal and lateral inflow velocity components to the rudder:

X R = -(1 -t R )F N sinδ (6.28) Y R = -(1 + a H )F N cosδ (6.29) N R = -(x R + a H x H )F N cosδ (
F N = 1 2 ρA R U 2 R f α sinα R (
U R = u 2 R + v 2 R (6.32) α R ≈ δ - v R u R (
u R = ε R u(1 -w P ) η R 1 + κ R 1 + 8K T πJ 2 P -1 2 + (1 + η R ) (6.34) v R = U γ R β R (6.35) β R = β -l R r
is the effective inflow angle to the rudder. η R = D P /H R is the ratio of propeller diameter to rudder span. ε R = (1 -w R )/(1 -w P ) is the wake fraction ratio at rudder position to that at propeller position. κ is an experimental constant. An actuator disk model is implemented to account for the propulsion [START_REF] Stern | A viscous-flow approach to the computation of propeller-hull interaction[END_REF]. Instead of using a real propeller, the actuator disk accelerates and rotates the fluids passing the defined region in the axial and tangential directions, which can greatly reduce the computational costs. The volume forces in the two directions are expressed as: 6.36) where r d is the radius of the disk, r * d = (r d -r h )/(1 -r h ), r d = r d /R P . A x and A θ are computed as:

Virtual captive model tests

f bx = A x r * d 1 -r * d , f bθ = A θ r * d 1 -r * d r * d (1 -r h ) + r h (
A x = 105 8 T P π∆(3R H + 4R P )(R P -R H ) (6.37) A θ = 105 8 Q P π∆(R P -R H )(3R P + 4R H ) (6.38)
R H and R P are the radius range of the disk. The thrust T P is obtained from the self propulsion tests [START_REF] Yasukawa | Introduction of MMG standard method for ship maneuvering predictions[END_REF] and the torque Q P can be calculated from the open water test (Figure 6.18). When the ship moves, the position of the actuator disk should also be updated based on the two coordinate frames in Figure 6.1:

x = (x 0 -x 0G ) cos ψ + (y 0 -y 0G ) sin ψ y = -(x 0 -x 0G ) sin ψ + (y 0 -y 0G ) cos ψ z = z 0 (6.39)
It can be seen that the upper equations only update in the horizontal plane as the captive model tests without sinkage and trim.

Computational setups

The hydrodynamic derivatives of the hull are determined using OTT (Oblique Towing Test) and CMT (Circular Motion Test). PMM (Planar Motion Mechanism) can also be used, but the obtained derivative values change remarkably due to the influence of motion frequency and amplitude in the tests [START_REF] Yasukawa | Introduction of MMG standard method for ship maneuvering predictions[END_REF]. It should be noted that the simulations here were carried out without sinkage and trim. The SST k-ω turbulence model is adopted in this case. The rudder angle δ = 0 • in these tests.

• Oblique Towing Tank test

In OTT, the hull drift angle β is changed in the range of 0 ∼ 20 • with 4 • interval. The computational domain is designed as in Figure 6.2. During the OTT simulations, the same mesh is used while the direction of the inlet velocities is changed to realize different oblique conditions. This avoids remeshing when changing the ship drift angle. The actuator disk is located as the real propeller (Table 6.5). The case setups are shown in Table 6. • Circular Motion Test In CMT, the non-dimensional yaw rate r is changed in the range of 0 ∼ 0.8 with 0.2 interval at a certain drift angle. 6.2). The position of the actuator disk is updated each time step according to 6.39. The non-dimensional forces and moment are calculated in . Only the simulations with positive drift angles were carried out. The negative drift angle part is projected directly from the positive ones. The ship has larger lateral force and yaw moment with a larger drift angle. The lateral force and yaw moment in the OTT agree the best with the experimental results [START_REF] Yasukawa | Introduction of MMG standard method for ship maneuvering predictions[END_REF].

For the longitudinal force, the values with negative drift angles are underestimated since they are projected directly from the positive ones.

• Circular Motion Test

The ship-generated waves, velocity contours generated by the actuator disk and the force and moment histories when the yaw rate r = 0 are shown in Figures 6.9-6.11. The wave elevation is higher with a smaller drift angle. The wave elevation is clearly influenced by the volume force of the actuator disk. The cross section in z direction is extracted to show the velocity contours in Figure 6.10. A similar effect to a propeller can be observed. The velocity turns to the inside of the turning circle in CMT. The convergences of the surge and lateral forces are steady, while that of the yaw moment is periodic. Only the first 30 s is shown for clarity. The amplitudes of the forces and moment increase with the drift angle, but the period remains the same. The non-dimensionalized forces and moment are shown in .

• Hydrodynamic coefficient calculation Equations 6.11-6.13 should be modified to express the actual forces on the hull:

(m + m x ) u = X * H + X R + X P (6.40) 6.42) where the approximation u ≈ 1 is employed. The hull forces are written as:

(m + m y ) v m + x G m ṙ = Y * H + Y R (6.41) (I zG + x 2 G m + J z ) ṙ + x G m v m = N * H + N R (
X * H = X H + (m + m y )v m r + x G m r 2 (6.43) 6.44) 6.45) By taking 6.20∼6.22 into the last equations, the following expressions can be obtained to calculate the hydrodynamic derivatives: Least squares method is adopted to fit these equations using the interpolated simulation results in . The fitting surfaces are shown in Figure 6.12. The coefficient values can be found in Table 6.3, which agree with experiments except several coefficients (X rr , Y rrr , N rrr ). This is because the experiments covered a wider range comparably. In terms of (X vr +m +m y ), (X rr +x G m ), (Y R -m -m x ) and (N R -x G m ), the added mass components are calculated as follows:

Y * H = Y H -(m + m x )r (
N * H = N H -x G m r (
X * H = -R 0 + X vv v 2 m + (X vr + m + m y )v m r + (X rr + x G m )r 2 + X vvvv v 4 m (6.46) Y * H =Y v v m + (Y R -m -m x )r + Y vvv v 3 m + Y vvr v 2 m r + Y vrr v m r 2 + Y rrr r 3
• m is given from the displacement volume of the ship.

• m x , m y and J z are estimated by [START_REF] Motora | On the measurement of added mass and added moment of inertia for ship motions[END_REF][START_REF] Motora | On the measurement of added mass and added moment of inertia for ship motions[END_REF] chart.

The final hydrodynamic derivatives and added mass values of the hull are shown in Table 6.3. The computational meshes can be seen in Figure 6.13. AMI technique is used for the propeller rotation. The rudder, propeller and its downwind are refined to capture the flow characteristics. The velocity contours at the propeller position with different rudder angles are shown in Figure 6.14. The flows are apparently accelerated by the propeller and directed by the rudder, which is also the mechanism of the advancing and turning of the ship. In Figure 6.15, the longitudinal force, lateral force, yaw moment and rudder normal force are calculated and compared with experiments [START_REF] Yasukawa | Introduction of MMG standard method for ship maneuvering predictions[END_REF]. Except the longitudinal force, others are approximately linear with the rudder angle. The propeller thrust remains the same at different rudder angles. Based on 6.49-6.51, the relations between the forces and moment on the hull and the rudder forces are calculated in Figure 6.16. The hydrodynamic coefficients t R , a H , x H can then be obtained from the slopes. The calculated rudder coefficients can be found in Table 6.4. Deviations can be clearly observed since only five cases were simulated in our work. More simulations can be done to reduce the computational and fitting errors. Other coefficients of the rudder are determined as follows:

Results and analyses

• The Longitudinal coordinate of the rudder position x R is approximately -0.5L pp , i.e. x R = -0.5.

• Rudder lift gradient coefficient f α = 6.13Λ/(Λ+ 2.25) is estimated using [START_REF] Fujii | Experimental researches on rudder performance (2)[END_REF] formula with a practical treatment for mariner rudder by [START_REF] Hirano | Open water performance of semi-balanced rudder[END_REF]. Λ = B R /C R is the aspect ratio of the rudder including the horn part, where B R , C R are the rudder span and chord length respectively.

• κ and ε are determined from the rudder force test results in straight moving under various propeller loads. Their values are directly taken from [START_REF] Yasukawa | Introduction of MMG standard method for ship maneuvering predictions[END_REF] here.

• γ R and l R are calculated by fitting the data from oblique towing tank tests and flow straightening coefficient tests. Since these tests are easy to introduce uncertainties, the coefficient values are taken directly from [START_REF] Yasukawa | Introduction of MMG standard method for ship maneuvering predictions[END_REF]. The computational domain is designed as in Figure 6.17. The domain extends 5 propeller diameters (D P ) upstream, 10 diameters downstream and 10 diameters in the lateral direction. The region near the propeller is refined to capture the flow details. The SST k-ω model is adopted in this case. The grid number and time step are 1,089,926 and 1 × 10 -3 s respectively. The outer correctors and pressure correctors are 18 and 3. All the domain rotates with a prescribed rotation velocity.

The inlet velocity (U in ) is varied to obtain different advance coefficient (J P ):

J P = U in n P D P (6.52) Two methods are used during the open water tests, the single-run approach and multi-run approach [START_REF] Xing | Computational towing tank procedures for single run curves of resistance and propulsion[END_REF]. The multi-run approach is the normal one where a constant inlet velocity is used. The single-run approach ramps up the inlet velocity with a small acceleration such that the time derivatives become negligible, thus the solution is in a quasi-steady state. The obtained result covers a wide range in a single run, which is time-saving. After the simulations, the thrust (T P ) and torque (Q P ) can be measured, the corresponding coefficients can be calculated as:

K T = T P ρn 2 P D 4 P K Q = Q P ρn 2
P D 5 P (6.53)

Results and analyses

The simulation results can be found in Figure 6.18. Both approaches agree with experiments, except for the beginning of the single run approach. This is because the initial parameters may not be realistic and the flows should evolve with the simulation. The values out of the experimental range (J P >0.8) are also predicted well. The vortical structures of Q v =300 s -2 , the second invariant of the velocity gradient tensor, with different advance coefficients are shown in Figure 6.19. The tip vortexes are resolved well. However, because of RANS, detailed vortical structures

• The wake characteristic constants C 1 and C 2 are taken from [START_REF] Yasukawa | Introduction of MMG standard method for ship maneuvering predictions[END_REF] directly.

• The Longitudinal coordinate of propeller position x P is approximately -0.5L pp , i.e. x P = -0.5.

• Effective wake in straight moving w P 0 is assumed to be 0.40 for L7-model and 0.35 for full-scale as in [START_REF] Yasukawa | Introduction of MMG standard method for ship maneuvering predictions[END_REF].

Maneuvering simulations

Using the hydrodynamic coefficients obtained in Tables 6.3,6.4 and 6.6, the maneuvering simulations are conducted for the KVLCC2 L7 model (Table 6 and yaw rate (r) are analyzed. For the turning circle tests , the ship velocity and heading angle are predicted accurately. But the maximum value of the yaw rate is underestimated, which leads to the deviation of the trajectory. The zigzag tests are more difficult since the rudder angle changes with time . Phase shift can be clearly observed. Generally, the maneuvering simulations agree with the experimental results, which proves the validity of our methods. Chapter 7

Conclusions and future work

Conclusions

Inland shipping transport and ship maneuvering were simulated based on CFD. Numerical solvers, including the turbulence models, boundary conditions and dynamic mesh techniques, meshing strategy and time step selection were validated using the surface-piercing cylinder, static KVLCC2, KCS with free sinkage and trim and free falling sphere case. These results provide valuable suggestions and reliability for the simulations in this thesis.

The confinement effect of the inland waterway was simulated using the convoy 1 and convoy 2 with various combinations of channel dimensions, water depths, ship draughts and speeds. Narrower channel width, smaller water depth, larger draught and higher speed were found to increase the advancing resistance and the wave elevations. The ship resistance and the blockage ratio show linear relations, whose slope increases with higher speed and draught. This means that the ship resistance changes more intensively under these conditions. The theoretical Kelvin angle is modified with the ship draught, which corresponds well with the simulation results.

The head-on encounter between the convoy 2 and the inland tanker was simulated using the ACMI technique. The ship resistance and wave pattern were analyzed based on the relative position parameter ξ ss . The average resistance agrees with experiments. The advancing resistance was analyzed, which was found to always drag the involving ships during the head-on encounter. The wave elevations at the channel center and the drawdown on the bank were analyzed as functions of the ship speed, channel width and draught, whose results agree with the confinement cases.

The ship passing bridge piles was simulated. A parameter ξ sb similar with the headon encounter was defined to characterize the passing process. The resistance change during the crossing was found to have a minimum value when the centers of the ship and pile coincide (ξ sb = 0). The original wave pattern is disturbed by the bridge piles. Two crests appear at the pile position because of the wave reflections and superpositions.

Virtual captive model tests were simulated, including the propeller open water test, rudder force test, OTT and CMT. An actuator disk model was implemented and validated to account for the propulsion. The hydrodynamic coefficients of the propeller, rudder and hull were calculated based on the MMG model. System-based maneuvering simulations were carried out and validated using experimental data.

The results presented how ships are influenced by the channel bank and bottom, passing ships and bridge piles. They provide useful suggestions for navigating under these conditions. These models and procedures allow us to simulate many shiprelated phenomena and investigate their physics based on CFD.

Future work

Besides the numerical and engineering contributions in this thesis, many improvements are also necessary:

-The inland shipping simulations and virtual captive model tests were carried out without sinkage and trim, which can be achieved by combining the dynamic mesh solvers in OpenFoam. This should increase the accuracy of the results.

-Real propellers instead of the actuator disk can be used to consider the thrust and torque influences on ships more realistically in the captive model tests.

-Other dynamic mesh techniques, such as Overset, IBM, etc. are good alternatives besides the ACMI method. They have no such constraints as ACMI for ship motions, which are possible to simulate real ship motions and maneuvering directly without the hydrodynamic coefficients.

-Other environmental influences, such as wind, waves and current can be added for seakeeping simulations.

With these improvements, this numerical solver will be able to simulate more conditions for ships. And the simulation results can be greatly improved for better predictions of their physics.

  to dedicate this thesis to my loving family.

3. 7

 7 Comparison between the (a) body-fitted mesh and (b) IBM (Immersed Boundary Method) mesh. Ω f denotes the fluid domain. Γ s is the solid boundary. The points are the discretized Lagrangian points of the solid boundary. The box is the supporting box at one Lagrangian point. The arrows represent the spreading and interpolation procedures of IBM. . 3.8 Six degrees of freedom and two coordinate frames of a ship. . . . . . 4.1 Computational domain of the interface-piercing cylinder case. . . . . 4.2 Convergence histories of the drag (F x ) with different velocity boundary conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Comparison of (a) zeroGradient and (b) inletOutlet boundary conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Convergence histories of the drag (F x ) with different pressure boundary conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Separated zoom view of the drag convergence history (F x ) with different pressure boundary conditions. . . . . . . . . . . . . . . . . . . 4.6 Convergence histories of the drag (F x ) with different phase fraction (alpha.water), turbulent kinetic energy (k), eddy viscosity (nut), wall and side boundary conditions. . . . . . . . . . . . . . . . . . . . . . . 4.7 Computational domain (a) of the KCS model and the meshes near (b) the bow and (c) the stern of the medium mesh (grid 2). . . . . . . . . 4.8 Wave elevation contours of (a) grid 1; (b) grid 2; (c) grid 3 in the mesh convergence tests. . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi 4.9 Wave profiles of the medium mesh (grid 2) at three lateral positions (y/L pp ) (a) 0.0741; (b) 0.1509; (c) 0.4224 and (d) on the hull surface. Experimental data can be found in Hino (2005). . . . . . . . . . . . . 4.10 Computational domain (a) of the KVLCC2 model and the meshes near (b) the bow and (c) the stern. . . . . . . . . . . . . . . . . . . . 4.11 (a) Convergence histories of the drag (F x ); (b) total resistance coefficients (C T ); (c) wave profiles on the hull and (d) simulation times (t) with different time steps (∆t) using 20 processors. Experimental data can be found in Pereira et al. (2007); Kim et al. (2001). L pp is the ship length between perpendiculars. . . . . . . . . . . . . . . . . . . . . .4.12 Computational domain of the water entry case of a free falling sphere. R s is the radius of the sphere. . . . . . . . . . . . . . . . . . . . . . . 4.13 Comparison of the (a) displacement (D m ) and (b) vertical force (F z ) histories of the free falling sphere (Shen & Wan, 2011). ∆t is the time step. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.14 (a) Convergence histories of the drag (F x ) and (b) comparison of resistance coefficients with different turbulence models. C T , C F , C P are the total, frictional and pressure resistance coefficients respectively. The experimental data can be found in Pereira et al. (2007). . . . . . . 4.15 Sketch of the (a) static and (b) prescribed motion cases. . . . . . . . . 4.16 (a) Convergence histories of the drag (F x ) and (b) comparison of the resistance coefficients between the static and prescribed motion cases. C T , C F , C P are the total, frictional and pressure resistance coefficients respectively. The experimental data can be found in Pereira et al. (2007). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.17 Comparison of the resistance coefficient (C T ), sinkage (σ s ) and trim (τ t ) with experimental data (Larsson et al., 2013). F r is the Froude number based on the ship length. . . . . . . . . . . . . . . . . . . . .5.1 Photos of the (a) confined channel and (b) convoy test of the experiments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Geometries of the (a) convoy 1 (1 pusher + 1 barge) and (b) convoy 2 (1 pusher + 2 barges). . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Body lines of the barge. . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Body lines of the pusher. . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Channel dimension and water depth setups for the convoy tests. . . . . 5.6 (a) Computational mesh of the convoy 1 case. (b) and (c) show the meshes near the stern and the bow. . . . . . . . . . . . . . . . . . . . vii 5.7 Experimental and numerical resistances of the convoy 1. T d , W b , R x , h w and V are the draught, channel bottom width, resistance, water depth and ship speed respectively. . . . . . . . . . . . . . . . . . . . 5.8 Experimental and numerical resistances of the convoy 2. T d , W b , R x , h w and V are the draught, channel bottom width, resistance, water depth and speed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.9 Ship resistance (R x ) in function of the blockage ratio (m b = A s /A c ) of the (a) convoy 1; (b) convoy 2 with different vessel speeds (V ). The draught T d = 0.04 m. The lines are fitting curves of the data. . . . . . 5.10 Ship resistance (R x ) in function of the blockage ratio (m b = A s /A c ) with different ship draughts (T d ). (a) and (b) belong to the convoy 1. (c) and (d) belong to the convoy 2. (a) and (c) correspond with a velocity 0.45 m/s. (b) and (d) correspond with a velocity 0.80 m/s. The lines are fitting curves of the data. . . . . . . . . . . . . . . . . . 5.11 Ship-generated waves of the convoy 1 (vessel speed V = 0.91 m/s, water depth h w = 0.24 m). The draught T d = 0.04 m for (a), (c), (e); 0.1 m for (b), (d), (f). . . . . . . . . . . . . . . . . . . . . . . . . . . 5.12 Ship-generated waves of the convoy 2 (vessel speed V = 0.91 m/s, water depth h w = 0.24 m). The draught T d = 0.04 m for (a), (c), (e); 0.1 m for (b), (d), (f). . . . . . . . . . . . . . . . . . . . . . . . . . . 5.13 (a) Wave angles (α k ) of the convoy 1 and convoy 2 with different draughts (T d ) and channel bottom widths (W b ). The vessel speed V = 0.91 m/s. (b) Theoretical prediction of the wave angle depending on the depth-Froude number (F r h ). Line/point 1 is the original theoretical wave angle which is irrelevant with the draught. Lines/points 2 and 3 are the new values calculated based on the modified equation 5.3. 5.14 Wave elevation at the lateral position (y = 0.3 m) of the convoy 1 and convoy 2. The draughts of (a), (b), (c) are 0.04 m; (d), (e), (f) are 0.1 m. The channel widths of (a), (d) are 0.72 m; (b), (e) are 1.44 m; (c), (f) are 2.88 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.15 Photos of the (a) head-on encounter between the convoy 2 and an inland tanker. (b) The tanker is towed by a rope connecting with (c) a rotating drum. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.16 Body lines of the inland tanker. . . . . . . . . . . . . . . . . . . . . . 5.17 Channel dimension and water depth setups for the head-on encounter. 5.18 Computational domain and ACMI positions for the head-on encounter case between two inland vessels. . . . . . . . . . . . . . . . . . . . . viii 5.19 Computational mesh and ACMI positions. . . . . . . . . . . . . . . . 5.20 Relative position (ξ ss ) between the two inland vessels during the headon encounter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.21 Comparison between experimental and numerical resistances of the convoy 2 during the head-on encounter. R x , T d , W b , V are the ship resistance, draught, channel bottom width and speed respectively. The draught of the tanker is 0.07 m. . . . . . . . . . . . . . . . . . . . . . 5.22 Resistance (R x ) of the convoy 2 with the relative position (ξ ss ) under different vessel speeds (V convoy2 ). The velocity of the tanker is 0.5 m/s. The channel bottom width is 1.44 m. The draughts of the convoy 2 and the tanker are 0.04 and 0.07 m. . . . . . . . . . . . . . . . . . . 5.23 Resistance (R x ) of the convoy 2 with the relative position (ξ ss ) under different channel widths (W b ). The velocities of the tanker and the convoy 2 are 0.5 and 0.91 m/s respectively. The draughts of the convoy 2 and the tanker are 0.04 and 0.07 m. . . . . . . . . . . . . . 5.24 Resistance (R x ) of the convoy 2 with the relative position (ξ ss ) under different ship draughts (T d ). The velocities of the tanker and the convoy 2 are 0.5 and 0.91 m/s respectively. The channel bottom width and the draught of the tanker are 1.44 and 0.07 m. . . . . . . . . . . . 5.25 Ship-generated waves with different ship speeds (ξ ss = 0). The velocity of the convoy 2 is (a) 0.91; (b) 0.68 m/s. The velocity of the tanker is 0.5 m/s. The channel bottom width is 1.44 m. The draughts of the convoy 2 and the tanker are 0.04 and 0.07 m. . . . . . . . . . . 5.26 Ship-generated waves with different draughts and channel widths (ξ ss = 0). The velocities of the convoy 2 and the tanker are 0.91 and 0.5 m/s. The channel width is 1.44 m for (a) and (b), 2.88 m for (c) and (d). The draugtht of the convoy 2 is 0.04 m for (a) and (c), 0.1 m for (b) and (d). The draught of the tanker is 0.07 m. . . . . . . . . . . . . . . 5.27 Ship wave profiles with different speeds (V convoy2 ) of the convoy 2. The speed of the tanker is 0.5 m/s. The channel bottom width is 1.44 m. The draughts of the convoy 2 and the tanker are 0.04 and 0.07 m. . 5.28 Ship wave profiles with different draughts (T d ) and channel widths (W b ). The wave elevations are modified using the initial draught for comparison. The velocities of the convoy 2 and the tanker are 0.91 and 0.5 m/s respectively. The draught of the tanker is 0.07 m. . . . . 5.29 Photo of the convoy 2 passing the bridge piles. . . . . . . . . . . . . . 5.30 Computational domain of the convoy 2 passing the bridge piles. . . . ix 5.31 Computational mesh of the convoy 2 passing the bridge piles. . . . . . 5.32 Relative position (ξ sb ) between the convoy 2 and bridge piles. . . . . . 5.33 Experimental and numerical resistances of the convoy 2 passing the bridge piles. R x , h w , T d , V are the resistance, water depth, draught and vessel speed respectively. . . . . . . . . . . . . . . . . . . . . . . 5.34 Ship resistance (R x ) with the relative position (ξ sb ) under different pile distances (D b ). The velocity, draught and water depth of the convoy 2 are 0.91 m/s, 0.04 m and 0.24 m respectively. . . . . . . . . . 5.35 Ship resistance (R x ) with the relative position (ξ sb ) under different vessel speeds (V ). The draught and water depth of the convoy 2 are 0.04 m and 0.24 m respectively. The pile distance is 0.8 m. . . . . . . 5.36 Ship resistance (R x ) with the relative position (ξ sb ) under different water depths (h w ). The velocity and draught of the convoy 2 are 0.91 m/s and 0.04 m respectively. The pile distance is 0.8 m. . . . . . . . 5.37 Ship-generated waves with different speeds. ξ sb = 0 for (a), (b), (c). ξ sb = 1 for (d), (e), (f). V = 0.907 m/s for (a), (d). V = 0.793 m/s for (b), (e). V = 0.684 m/s for (c), (f). The draught, water depth and pile distance are 0.04, 0.24 and 0.8 m respectively. . . . . . . . . . . 5.38 Ship wave profiles at the lateral position y = 0.26 m with different speeds. The red line denotes the sampling position of the wave. (b) is the zoomed view of the box in (a). The draught, water depth and pile distance are 0.04, 0.24 and 0.8 m respectively. The dotted lines in (b) represent the position of the pile. . . . . . . . . . . . . . . . . . . . . 5.39 Ship-generated waves with different water depths. ξ sb = 0 for (a), (b). ξ sb = 1 for (c), (d). h w = 0.24 m for (a), (c). h w = 0.18 m for (b), (d). The vessel speed, draught and pile distance are 0.91 m/s, 0.04 m and 0.8 m respectively. . . . . . . . . . . . . . . . . . . . . . . . . . 5.40 Ship wave profiles at the lateral position y = 0.26 m with different water depths. The water levels are adjusted using the initial value for comparison. The vessel speed, draught and pile distance are 0.91 m/s, 0.04 m and 0.8 m respectively. The dotted lines denote the positions of the vessel and pile. . . . . . . . . . . . . . . . . . . . . . . . . . . 5.41 Ship-generated waves with different distances of the bridge piles. ξ sb = 0 for (a), (b). ξ sb = 1 for (c), (d). D b = 0.8 m for (a), (c). D b = 0.7 m for (b), (d). The vessel speed, draught and water depth are 0.91 m/s, 0.04 m and 0.24 m respectively. . . . . . . . . . . . . . . . . . . . . x 5.42 Ship wave profiles at the lateral position y = 0.26 m with different bridge pile distances. The vessel speed, draught and water depth are 0.91 m/s, 0.04 m and 0.24 m respectively. The dotted lines denote the vessel and pile positions. . . . . . . . . . . . . . . . . . . . . . . 6.1 Coordinate systems of the MMG model. . . . . . . . . . . . . . . . . 6.2 Computational domain and boundary conditions for the OTT (Oblique Towing Tank) and CMT (Circular Motion Test). . . . . . . . . . . . . 6.3 Sketch of the CMT (Circular Motion Test). r, R, U are the yaw rate, turning radius and ship speed. . . . . . . . . . . . . . . . . . . . . . . 6.4 Velocity contours (a, c, d) and volume force (b) at the actuator disk position (Drift angle β = 20 • ). U in is the inlet velocity. . . . . . . . . 6.5 Ship-generated waves with different drift angles (β) in the Oblique Towing Tank (OTT) tests. The drift angle (a) β = 4 • ; (b) β = 8 • ; (c) β = 12 • ; (d) β = 16 • ; (e) β = 20 • . . . . . . . . . . . . . . . . . . . . 6.6 Non-dimensional surge force (X * H ) in OTT and CMT. . . . . . . . . . 6.7 Non-dimensional lateral force (Y * H ) in OTT and CMT. . . . . . . . . 6.8 Non-dimensional yaw moment (N * H ) in OTT and CMT. . . . . . . . . 6.9 Ship-generated waves in CMT (Circular Motion Tests). The drift angle (a) β = 0 • ; (b) β = 4 • ; (c) β = 12 • ; (b) β = 20 • . The yaw rate r = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.10 Velocity contours at the actuator disk position in CMT (Circular Motion Tests). The drift angle (a) β = 0 • ; (b) β = 4 • ; (c) β = 12 • ; (b) β = 20 • . The yaw rate r = 0. . . . . . . . . . . . . . . . . . . . . . 6.11 Convergence histories of the surge force (X * H ), lateral force (Y * H ) and yaw moment (N * H ) in CMT (Circular Motion Tests) (yaw rate r = 0). 6.12 Fitting surfaces of the OTT (Oblique Towing Tank) and CMT (Circular Motion Tests) results. v m , r are the dimensionless lateral velocity and yaw rate. X * H , Y * H , N * H are the dimensionless forces and moment on the hull. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.13 Computational domain and mesh for the rudder force tests. (a) shows the positions of the hull, rudder and propeller; (b) zooms the AMI (Arbitrary Mesh Interface) position; (c) and (d) show the mesh refinements at two cross sections of the propeller. . . . . . . . . . . . . . . 6.14 Velocity contours at the propeller position in the rudder force tests. The ruder angle (a) δ = 5 • ; (b) δ = 10 • ; (c) δ = 15 • ; (d) δ = 20

58 4 . 4 5 . 1

 4451 Possible outlet boundary conditions for the phase fraction (alpha.water). 60 4.5 Possible inlet and outlet boundary conditions for the turbulent kinetic energy (k). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.6 Possible inlet and outlet boundary conditions for the eddy viscosity (nut). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.7 Possible velocity (U ) and pressure (p_rgh) boundary conditions for the wall (Figure 4.1). . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.8 Possible boundary conditions of sides (Figure 4.1). . . . . . . . . . . 60 4.9 Geometrical parameters of the KCS model. . . . . . . . . . . . . . . 63 4.10 Grid convergence results of the total resistance coefficient. Experimental data can be found in Hino (2005). . . . . . . . . . . . . . . . . 63 4.11 Comparison of the frictional (C F ) and pressure (C P ) coefficients in the mesh convergence tests. The frictional coefficient is estimated by the ITTC 1957 friction line C F = 0.075/(log 10 Re-2) 2 . C P = C T -C F . 63 4.12 Geometrical parameters of the KVLCC2 model. . . . . . . . . . . . . 66 4.13 Comparison of the frictional (C F ) and pressure (C P ) coefficients in the time step convergence tests. The frictional coefficient is estimated by the ITTC 1957 friction line C F = 0.075/(log 10 Re -2) 2 . C P = C T -C F . ∆t is the time step. . . . . . . . . . . . . . . . . . . . . . . 68 4.14 Case setups of the trim and sinkage predictions for the KCS model. . . 71 xiii Geometrical parameters of the full-and model-scale ships. . . . . . . 5.2 Test configurations and related channel parameters. The units of the channel bottom width W b , ship draught T d and water depth h w are [m]. The units of the maximum speed V max , Schijf limiting speed V s and subcritical speed V sub c are [m/s]. A c and A s are the cross section areas of the channel and the ship. B c = A c /h w and B s are the average channel width and the ship breadth. . . . . . . . . . . . . . . . . . . . 5.3 Mesh numbers of the convoy tests when the water depth h w =0.24 m. T d , W b and N mesh are the ship draught, channel bottom width and mesh number respectively. . . . . . . . . . . . . . . . . . . . . . . . 5.4 Test configurations and mesh numbers of the head-on encounter case. The units of the channel width W b , draught T d and water depth h w are [m]. The unit of the maximum velocity V max is [m/s]. T d is the draught of the convoy 2. N mesh is the mesh number. . . . . . . . . . . 5.5 Test configurations and mesh numbers of the convoy 2 passing the bridge piles. The units of the channel bottom width W b , distance between bridge piles D b , ship draught T d and water depth h w are [m]. The unit of the maximum velocity V max is [m/s]. N mesh is the mesh number. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1 Principle parameters of the KVLCC2 tanker. . . . . . . . . . . . . . . 6.2 Setups of the OTT (Oblique Towing Tank) and CMT (Circular Motion Test). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Hydrodynamic derivatives and added masses of the hull. . . . . . . . 6.4 Hydrodynamic coefficients obtained from the rudder force tests. . . . 6.5 Geometrical parameters of the KP458 propeller. . . . . . . . . . . . . 6.6 Hydrodynamic coefficients of the propeller. . . . . . . . . . . . . . . A.1 Standard normalization systems

Figure 2 . 1 :

 21 Figure 2.1: Sketch of the towing tank.

Figure 2 . 2 :

 22 Figure 2.2: Numerical virtual tanks for various purposes.

Figure 2 . 3 :

 23 Figure 2.3: Flowchart of the CFD computation process.

Figure 2 . 4 :

 24 Figure 2.4: Two domain discretization approaches for the FVM mesh[START_REF] Moukalled | 4 Two domain discretization approaches for the FVM mesh[END_REF]. CV means the Control Volume.

Figure 2 . 5 :

 25 Figure 2.5: Components of the ship resistance[START_REF] Bertram | Practical ship hydrodynamics[END_REF][START_REF] Molland | Ship resistance and propulsion[END_REF].

Figure 2 . 6 :

 26 Figure 2.6: Typical curves of the ship resistance.

Figure 2 .

 2 Figure 2.7: Velocity, pressure and surface elevation of the primary wave pattern.

Figure 2 .

 2 Figure 2.8: (a) Kelvin wave pattern of a moving pressure source; (b) wave system interference of a real ship hull[START_REF] Molland | Ship resistance and propulsion[END_REF]. V and λ denote the velocity and the wavelength.

Figure 2 .

 2 Figure 2.9: Schematic of the ship squat, wave change and return current. (a) Cross section; (b) side view. "UKC" is the Under Keel Clearance.

Figure 2 .

 2 Figure 2.10: Typical curves of the ship resistance in open, shallow and confined waters[START_REF] Pompée | About modelling inland vessels resistance and propulsion and interaction vessel-waterway key parameters driving restricted/shallow water effects[END_REF].

Figure 2 .

 2 Figure 2.11: Confinement effect distribution in the cross section[START_REF] Lataire | Longitudinally directed bank effects[END_REF].

  ; Duffy (2008); Yao & Zou (2010), and CFD simulations in He et al. (2016); Jachowski (2008); Mucha et al. (2016, 2014); Sergent et al. (2015); Terziev et al. (2018); Tezdogan et al.

Figure 2 .

 2 Figure 2.13: Wave pattern in a fully-confined waterway(Caplier, 2016).

Figure 2 .

 2 Figure 2.14: Sketch of the bank-induced effect.

  ,[START_REF] Sian | Assessment of Ship-bank interactions on LNG tanker in shallow water[END_REF],[START_REF] Kaidi | Numerical estimation of bankpropeller-hull interaction effect on ship manoeuvring using CFD method[END_REF],[START_REF] Zou | CFD prediction and validation of ship-bank interaction in a canal. 2nd International conference on Ship Manoeuvring in Shallow and Confined Water: Ship to ship interaction[END_REF],[START_REF] Ma | Hydrodynamic interaction among hull, rudder and bank for a ship sailing along a bank in restricted waters[END_REF],[START_REF] Torben | The development of a method to predict ship-bank interaction in real-time[END_REF], and potential flow methods inYuan 2.4.5

Figure 2 .

 2 Figure 2.15: Sketch of the ship overtaking process. ξ ss represents the relative position between two ships. F x , F y , M z denote the longitudinal force, lateral force and yawing moment on the ship respectively.

  [START_REF] Vantorre | Model test based formulations of ship-ship interaction forces[END_REF],[START_REF] Lima | Modelling of close proximity manoeuvres in shallow water channels[END_REF];[START_REF] Lima | Study of ship-to-ship interaction in shallow water with account for squat phenomenon[END_REF] and potential flow theory[START_REF] Kijima | The interaction effects between two ships in the proximity of bank wall[END_REF] are the most commonly-used. Several review and application of these models can be found in ITTC(2002),[START_REF] Lima | Modelling of close proximity manoeuvres in shallow water channels[END_REF][START_REF] De Decker | Ship-ship interaction during lightering operations[END_REF] and[START_REF] Lee | Assessment of safe navigation including the effect of ship-ship interaction in restricted waterways[END_REF];[START_REF] Lee | The investigation of ship maneuvering with hydrodynamic effects between ships in curved narrow channel[END_REF];[START_REF] Lee | A numerical study on ship-ship interaction in shallow and restricted waterway[END_REF]. CFD simulations for overtaking were carried out by[START_REF] Lo | Numerical simulation of hydrodynamic interaction produced during the overtaking and the head-on encounter process of two ships[END_REF][START_REF] Fonfach | Numerical study of shipto-ship interaction forces on the basis of various flow models[END_REF][START_REF] Mousaviraad | Ship-Ship interactions in calm water and waves. Part 2: URANS validation in replenishment and overtaking conditions[END_REF],Lee (

Figure 2 .

 2 Figure 2.16: Sketch of the head-on encounter process. ξ ss represents the relative position between two ships. F x , F y , M z denote the longitudinal force, lateral force and yawing moment on the ship respectively.

Figure 2 .

 2 Figure 2.17: Classifications of ship maneuvering test methods.

Figure 3 . 1 :

 31 Figure 3.1: Parameter definitions in the finite volume discretization[START_REF] Jasak | Error analysis and estimation for the finite volume method with applications to fluid flows[END_REF].

Figure 3 . 2 :

 32 Figure 3.2: Pressure-velocity coupling procedure (PIMPLE).

  i and j denote the ith and jth facets at the slave and master patches. n is the number of master facets neighboring the slave facet i, and m is the number of the slave facet neighboring the master facet j. W M -S and W S-M are the master to slave and slave to master weighting factor. |S M |, |S S | and |S M ∩S | are the surface area of the master facet, slave facet and the intersection between the master and slave facets.

Figure 3 . 3 :

 33 Figure 3.3: Comparison between (a) AMI (Arbitrary Mesh Interface) and (b) ACMI (Arbitrary Coupled Mesh Interface) methods. (c) shows the discretization error of the meshes at the interface.

Figure 3 . 4 :

 34 Figure 3.4: One-dimensional domain of the overset mesh, with two-subdomains A and B. Values at nodes 6 and 7 are unknown. The arrows denote the interpolation from the donor nodes/cells to the fringe nodes/cells (Ma et al., 2018).

Figure 3 . 5 :

 35 Figure 3.5: Hole cutting and DCI (Domain Connectivity Information) of two overlapping meshes in a one-dimensional domain. A is the background mesh. B is the body-fitted mesh for a structure. The arrows denote the interpolation from the donor nodes/cells to the fringe nodes/cells. The crosses represent the hole nodes/cells (Ma et al., 2018).

Figure 3 . 6 :

 36 Figure 3.6: Flowchart of the computational procedure with overset meshes. "DCI" is the domain connectivity information.

Figure 3 .

 3 Figure 3.7: Comparison between the (a) body-fitted mesh and (b) IBM (Immersed Boundary Method) mesh. Ω f denotes the fluid domain. Γ s is the solid boundary. The points are the discretized Lagrangian points of the solid boundary. The box is the supporting box at one Lagrangian point. The arrows represent the spreading and interpolation procedures of IBM.

Figure 3 . 8 :

 38 Figure 3.8: Six degrees of freedom and two coordinate frames of a ship.

Figure 4 . 1 :

 41 Figure 4.1: Computational domain of the interface-piercing cylinder case.

Figure 4 . 2 :

 42 Figure 4.2: Convergence histories of the drag (F x ) with different velocity boundary conditions.

Figure 4 . 3 :

 43 Figure 4.3: Comparison of (a) zeroGradient and (b) inletOutlet boundary conditions.

Figure 4 . 4 :

 44 Figure 4.4: Convergence histories of the drag (F x ) with different pressure boundary conditions.

Figure 4 . 5 :

 45 Figure 4.5: Separated zoom view of the drag convergence history (F x ) with different pressure boundary conditions.

Figure 4 . 6 :

 46 Figure 4.6: Convergence histories of the drag (F x ) with different phase fraction (alpha.water), turbulent kinetic energy (k), eddy viscosity (nut), wall and side boundary conditions.

  10. The ship is fixed in the domain and a constant velocity (2.196 m/s) is added at the inlet for the fluids. The corresponding Froude (F r) and Reynolds numbers (Re) are 0.26 and 1.4×10 7 . The grid convergence of the total resistance coefficient is analyzed.

Figure 4 . 7 :

 47 Figure 4.7: Computational domain (a) of the KCS model and the meshes near (b) the bow and (c) the stern of the medium mesh (grid 2).

Figure 4 .

 4 Figure 4.8: Wave elevation contours of (a) grid 1; (b) grid 2; (c) grid 3 in the mesh convergence tests.

Figure 4 .

 4 Figure 4.9: Wave profiles of the medium mesh (grid 2) at three lateral positions (y/L pp ) (a) 0.0741; (b) 0.1509; (c) 0.4224 and (d) on the hull surface. Experimental data can be found in Hino (2005).

Figure 4 .

 4 Figure 4.10: Computational domain (a) of the KVLCC2 model and the meshes near (b) the bow and (c) the stern.

Figure 4 .

 4 Figure 4.11: (a) Convergence histories of the drag (F x ); (b) total resistance coefficients (C T ); (c) wave profiles on the hull and (d) simulation times (t) with different time steps (∆t) using 20 processors. Experimental data can be found in Pereira et al. (2007); Kim et al. (2001). L pp is the ship length between perpendiculars.

Figure 4 .

 4 Figure 4.13: Comparison of the (a) displacement (D m ) and (b) vertical force (F z ) histories of the free falling sphere[START_REF] Shen | Numerical simulation of sphere water entry problem based on VOF and dynamic mesh methods[END_REF]. ∆t is the time step.

Figure 4 .

 4 Figure 4.14: (a) Convergence histories of the drag (F x ) and (b) comparison of resistance coefficients with different turbulence models. C T , C F , C P are the total, frictional and pressure resistance coefficients respectively. The experimental data can be found in Pereira et al. (2007).

Figure 4 .

 4 Figure 4.15: Sketch of the (a) static and (b) prescribed motion cases.

Figure 4 .

 4 Figure 4.16: (a) Convergence histories of the drag (F x ) and (b) comparison of the resistance coefficients between the static and prescribed motion cases. C T , C F , C P are the total, frictional and pressure resistance coefficients respectively. The experimental data can be found in Pereira et al. (2007).

Figure 5 . 1 :

 51 Figure 5.1: Photos of the (a) confined channel and (b) convoy test of the experiments.

Figure 5 . 2 :

 52 Figure 5.2: Geometries of the (a) convoy 1 (1 pusher + 1 barge) and (b) convoy 2 (1 pusher + 2 barges).

Figure 5 . 3 :

 53 Figure 5.3: Body lines of the barge.

Figure 5 . 4 : 1 (Figure 5 . 5 :

 54155 Figure 5.4: Body lines of the pusher.

Table 5 . 2 :

 52 Test configurations and related channel parameters. The units of the channel bottom width W b , ship draught T d and water depth h w are [m]. The units of the maximum speed V max , Schijf limiting speed V s and subcritical speed V sub c are [m/s]. A c and A s are the cross section areas of the channel and the ship. B c = A c /h w and B s are the average channel width and the ship breadth.

Case

3 .

 3 The turbulence model and schemes in this case are the same as in Section 4.2. Five outer corrections and two pressure corrections are adopted. The simulations use 50 processors and take at most 42 cpu hours to finish. It should be noted that the simulations were carried out without trim and sinkage. The enhanced coupling approaches of Dunbar et al. (2015), Gatin et al. (2017) and Devolder et al. (2015) were implemented but with no success. This is probably due to the large added mass effect induced by the fully-confined condition as introduced in Section 3.4.3. Better FSI alternatives are Song et al. (2013); Jasak et al. (2017); Tezduyar (2001); Tezduyar et al. (2006); Connell & Yue (2007); Van Brummelen (2011) but not tested in this thesis.

Figure 5 . 6 :

 56 Figure 5.6: (a) Computational mesh of the convoy 1 case. (b) and (c) show the meshes near the stern and the bow.

Figure 5 . 7 :

 57 Figure 5.7: Experimental and numerical resistances of the convoy 1. T d , W b , R x , h w and V are the draught, channel bottom width, resistance, water depth and ship speed respectively.

Figure 5 . 8 :

 58 Figure 5.8: Experimental and numerical resistances of the convoy 2. T d , W b , R x , h w and V are the draught, channel bottom width, resistance, water depth and speed.

Figure 5 .

 5 Figure 5.9: Ship resistance (R x ) in function of the blockage ratio (m b = A s /A c ) of the (a) convoy 1; (b) convoy 2 with different vessel speeds (V ). The draught T d = 0.04 m. The lines are fitting curves of the data.

Figure 5 .

 5 Figure 5.10: Ship resistance (R x ) in function of the blockage ratio (m b = A s /A c ) with different ship draughts (T d ). (a) and (b) belong to the convoy 1. (c) and (d) belong to the convoy 2. (a) and (c) correspond with a velocity 0.45 m/s. (b) and (d) correspond with a velocity 0.80 m/s. The lines are fitting curves of the data.

Figure 5 .

 5 Figure 5.14 extracts the wave elevations at the lateral position y = 0.3 m near the convoys. Regardless of the draught value, the overall wave profiles are similar.Water level increases at the bow and maintains a long wave along the hull. A maximum wave elevation can be observed at the connection part between the two vessels. This is caused by the superposition of the stern wave of the front part and the bow wave of the back part. The cases with larger draught seem to have more intense wave changes, and the channel confinement further amplifies this effect. When the channel width is smaller, the ship-generated waves become higher at the bow and lower along the hull. The ship waves are more complicated for the convoy

Figure 5 .

 5 Figure 5.11: Ship-generated waves of the convoy 1 (vessel speed V = 0.91 m/s, water depth h w = 0.24 m). The draught T d = 0.04 m for (a), (c), (e); 0.1 m for (b), (d), (f).

Figure 5 .

 5 Figure 5.12: Ship-generated waves of the convoy 2 (vessel speed V = 0.91 m/s, water depth h w = 0.24 m). The draught T d = 0.04 m for (a), (c), (e); 0.1 m for (b), (d), (f).

Figure 5 .

 5 Figure 5.13: (a) Wave angles (α k ) of the convoy 1 and convoy 2 with different draughts (T d ) and channel bottom widths (W b ). The vessel speed V = 0.91 m/s. (b) Theoretical prediction of the wave angle depending on the depth-Froude number (F r h ). Line/point 1 is the original theoretical wave angle which is irrelevant with the draught. Lines/points 2 and 3 are the new values calculated based on the modified equation 5.3.

Figure 5 .

 5 Figure 5.14: Wave elevation at the lateral position (y = 0.3 m) of the convoy 1 and convoy 2. The draughts of (a), (b), (c) are 0.04 m; (d), (e), (f) are 0.1 m. The channel widths of (a), (d) are 0.72 m; (b), (e) are 1.44 m; (c), (f) are 2.88 m.

Figure 5 .

 5 Figure 5.15: Photos of the (a) head-on encounter between the convoy 2 and an inland tanker. (b) The tanker is towed by a rope connecting with (c) a rotating drum.

Figure 5 .

 5 Figure 5.16: Body lines of the inland tanker.

Figure 5 .

 5 Figure 5.17: Channel dimension and water depth setups for the head-on encounter.

Figure 5 .

 5 Figure 5.18: Computational domain and ACMI positions for the head-on encounter case between two inland vessels.

Figure 5 .

 5 Figure 5.19: Computational mesh and ACMI positions.

Figure 5 .

 5 Figure 5.20: Relative position (ξ ss ) between the two inland vessels during the headon encounter.

Figure 5 .

 5 Figure 5.21: Comparison between experimental and numerical resistances of the convoy 2 during the head-on encounter. R x , T d , W b , V are the ship resistance, draught, channel bottom width and speed respectively. The draught of the tanker is 0.07 m.

Figure 5 .

 5 Figure 5.22: Resistance (R x ) of the convoy 2 with the relative position (ξ ss ) under different vessel speeds (V convoy2 ). The velocity of the tanker is 0.5 m/s. The channel bottom width is 1.44 m. The draughts of the convoy 2 and the tanker are 0.04 and 0.07 m.

Figure 5 .

 5 Figure 5.23: Resistance (R x ) of the convoy 2 with the relative position (ξ ss ) under different channel widths (W b ). The velocities of the tanker and the convoy 2 are 0.5 and 0.91 m/s respectively. The draughts of the convoy 2 and the tanker are 0.04 and 0.07 m.

Figure 5 .

 5 Figure 5.24: Resistance (R x ) of the convoy 2 with the relative position (ξ ss ) under different ship draughts (T d ). The velocities of the tanker and the convoy 2 are 0.5 and 0.91 m/s respectively. The channel bottom width and the draught of the tanker are 1.44 and 0.07 m.

Figure 5 .

 5 Figure 5.25: Ship-generated waves with different ship speeds (ξ ss = 0). The velocity of the convoy 2 is (a) 0.91; (b) 0.68 m/s. The velocity of the tanker is 0.5 m/s. The channel bottom width is 1.44 m. The draughts of the convoy 2 and the tanker are 0.04 and 0.07 m.

Figure 5 .

 5 Figure 5.26: Ship-generated waves with different draughts and channel widths (ξ ss = 0). The velocities of the convoy 2 and the tanker are 0.91 and 0.5 m/s. The channel width is 1.44 m for (a) and (b), 2.88 m for (c) and (d). The draugtht of the convoy 2 is 0.04 m for (a) and (c), 0.1 m for (b) and (d). The draught of the tanker is 0.07 m.

Figure 5 .

 5 Figure 5.27: Ship wave profiles with different speeds (V convoy2 ) of the convoy 2. The speed of the tanker is 0.5 m/s. The channel bottom width is 1.44 m. The draughts of the convoy 2 and the tanker are 0.04 and 0.07 m.

Figure 5 .

 5 Figure 5.28: Ship wave profiles with different draughts (T d ) and channel widths (W b ). The wave elevations are modified using the initial draught for comparison. The velocities of the convoy 2 and the tanker are 0.91 and 0.5 m/s respectively. The draught of the tanker is 0.07 m.

Figure 5 .

 5 Figure 5.29: Photo of the convoy 2 passing the bridge piles.

Figure 5 .

 5 Figure 5.30: Computational domain of the convoy 2 passing the bridge piles.

Figure 5 .

 5 Figure 5.31: Computational mesh of the convoy 2 passing the bridge piles.

Figure 5 .

 5 Figure 5.33: Experimental and numerical resistances of the convoy 2 passing the bridge piles. R x , h w , T d , V are the resistance, water depth, draught and vessel speed respectively.

Figure 5 .

 5 Figure 5.34: Ship resistance (R x ) with the relative position (ξ sb ) under different pile distances (D b ). The velocity, draught and water depth of the convoy 2 are 0.91 m/s, 0.04 m and 0.24 m respectively.

Figure 5 .

 5 32 demonstrates the three most important positions. When ξ sb = -1, the vessel meets the piles.

Figure 5 .

 5 Figure 5.35: Ship resistance (R x ) with the relative position (ξ sb ) under different vessel speeds (V ). The draught and water depth of the convoy 2 are 0.04 m and 0.24 m respectively. The pile distance is 0.8 m.

Figure 5 .

 5 Figure 5.36: Ship resistance (R x ) with the relative position (ξ sb ) under different water depths (h w ). The velocity and draught of the convoy 2 are 0.91 m/s and 0.04 m respectively. The pile distance is 0.8 m.

Figure 5 .

 5 Figure 5.37: Ship-generated waves with different speeds. ξ sb = 0 for (a), (b), (c). ξ sb = 1 for (d), (e), (f). V = 0.907 m/s for (a), (d). V = 0.793 m/s for (b), (e). V = 0.684 m/s for (c), (f). The draught, water depth and pile distance are 0.04, 0.24 and 0.8 m respectively.

Figure 5 .

 5 Figure 5.38: Ship wave profiles at the lateral position y = 0.26 m with different speeds. The red line denotes the sampling position of the wave. (b) is the zoomed view of the box in (a). The draught, water depth and pile distance are 0.04, 0.24 and 0.8 m respectively. The dotted lines in (b) represent the position of the pile.

Figure 5 .

 5 Figure 5.39: Ship-generated waves with different water depths. ξ sb = 0 for (a), (b). ξ sb = 1 for (c), (d). h w = 0.24 m for (a), (c). h w = 0.18 m for (b), (d). The vessel speed, draught and pile distance are 0.91 m/s, 0.04 m and 0.8 m respectively.

Figure 5 .

 5 Figure 5.40: Ship wave profiles at the lateral position y = 0.26 m with different water depths. The water levels are adjusted using the initial value for comparison. The vessel speed, draught and pile distance are 0.91 m/s, 0.04 m and 0.8 m respectively. The dotted lines denote the positions of the vessel and pile.

Figure 5 .

 5 Figure 5.41: Ship-generated waves with different distances of the bridge piles. ξ sb = 0 for (a), (b). ξ sb = 1 for (c), (d). D b = 0.8 m for (a), (c). D b = 0.7 m for (b), (d). The vessel speed, draught and water depth are 0.91 m/s, 0.04 m and 0.24 m respectively.

Figure 5 .

 5 Figure 5.42: Ship wave profiles at the lateral position y = 0.26 m with different bridge pile distances. The vessel speed, draught and water depth are 0.91 m/s, 0.04 m and 0.24 m respectively. The dotted lines denote the vessel and pile positions.

  .10) where ∆u = u -u 0 is the disturbance in surge velocity. u 0 is the surge velocity in the initial state of forward motion. Depending on the ship model and maneuvering conditions, several terms can be neglected and the equations are thus simplified. Studies use PMM (Planar Motion Mechanism) to determine the hydrodynamic coefficients in the Abkowitz model can be found in Sakamoto et al. (2012); Hajivand & Mousavizadegan (2015); Sadat-Hosseini et al. (2011); Kim et al. (2015); Liu et al. (2018); El Moctar et al. (2014); Shenoi et al. (2014); Yoon (2009).

Figure 6 . 1 :

 61 Figure 6.1: Coordinate systems of the MMG model.

2 .

 2 Five drift angles are designed in this study. The inlet velocity and revolution rate of the KP458 propeller are 0.76 m/s and 17.95 rps. The thrust value is obtained from the experiment in[START_REF] Yasukawa | Introduction of MMG standard method for ship maneuvering predictions[END_REF]. The torque is obtained from the open water test, which will be introduced later, according to the thrust.

Figure 6 . 2 :

 62 Figure 6.2: Computational domain and boundary conditions for the OTT (Oblique Towing Tank) and CMT (Circular Motion Test).

Figure 6 .

 6 3 presents the typical movement of the CMT test. The ship model is towed on a circular arc line. Since the yaw rate r = U/R, changing the turning radius R is equivalent to changing r under the same ship speed U = 0.76 m/s. The surge force, lateral force and yaw moment are measured during the tests. The thrust and torque are determined the same way as the former part (Table

Figure 6 . 3 :

 63 Figure 6.3: Sketch of the CMT (Circular Motion Test). r, R, U are the yaw rate, turning radius and ship speed.

Figure 6 . 4 :

 64 Figure 6.4: Velocity contours (a, c, d) and volume force (b) at the actuator disk position (Drift angle β = 20 • ). U in is the inlet velocity.

Figure 6 . 5 :

 65 Figure 6.5: Ship-generated waves with different drift angles (β) in the Oblique Towing Tank (OTT) tests. The drift angle (a) β = 4 • ; (b) β = 8 • ; (c) β = 12 • ; (d) β = 16 • ; (e) β = 20 • .

Figure 6 . 6 :

 66 Figure 6.6: Non-dimensional surge force (X * H ) in OTT and CMT.

Figure 6 .

 6 Figure 6.7: Non-dimensional lateral force (Y * H ) in OTT and CMT.

Figure 6 .

 6 Figure 6.8: Non-dimensional yaw moment (N * H ) in OTT and CMT.

Figure 6 .

 6 Figure 6.9: Ship-generated waves in CMT (Circular Motion Tests). The drift angle (a) β = 0 • ; (b) β = 4 • ; (c) β = 12 • ; (b) β = 20 • . The yaw rate r = 0.

Figure 6 .

 6 Figure 6.10: Velocity contours at the actuator disk position in CMT (Circular Motion Tests). The drift angle (a) β = 0 • ; (b) β = 4 • ; (c) β = 12 • ; (b) β = 20 • . The yaw rate r = 0.

Figure 6 .

 6 Figure 6.11: Convergence histories of the surge force (X * H ), lateral force (Y * H ) and yaw moment (N * H ) in CMT (Circular Motion Tests) (yaw rate r = 0).

  =N v v m + (N R -x G m )r + N vvv v 3 m + N vvr v 2 m r + N vrr v m r 2 + N rrr r 3(6.48)

Figure 6 .

 6 Figure 6.12: Fitting surfaces of the OTT (Oblique Towing Tank) and CMT (Circular Motion Tests) results. v m , r are the dimensionless lateral velocity and yaw rate. X * H , Y * H , N * H are the dimensionless forces and moment on the hull.

Figure 6 .

 6 Figure 6.13: Computational domain and mesh for the rudder force tests. (a) shows the positions of the hull, rudder and propeller; (b) zooms the AMI (Arbitrary Mesh Interface) position; (c) and (d) show the mesh refinements at two cross sections of the propeller.

Figure 6 .

 6 Figure 6.14: Velocity contours at the propeller position in the rudder force tests. The ruder angle (a) δ = 5 • ; (b) δ = 10 • ; (c) δ = 15 • ; (d) δ = 20 • .

Figure 6 .

 6 Figure 6.15: Simulation results of the rudder force tests. X , Y , N m , F N , T , δ are the dimensionless longitudinal force, lateral force, yaw moment around midship, rudder normal force, propeller thrust and rudder angle respectively.

Figure 6 .

 6 Figure 6.16: Analysis results of the rudder force tests. X , Y , N m , F N , δ are the dimensionless longitudinal force, lateral force, yaw moment around midship, rudder normal force and rudder angle respectively.

Figure 6 .

 6 Figure 6.20: Turning circle maneuver (δ = 35 • ) of the KVLCC2 L7 model. U, ψ, r, L pp are the absolute velocity, heading angle, yaw rate and ship length respectively.

Figure 6 .

 6 Figure 6.21: Turning circle maneuver (δ = -35 • ) of the KVLCC2 L7 model. U, ψ, r, L pp are the absolute velocity, heading angle, yaw rate and ship length respectively.

Figure 6 .

 6 Figure 6.22: Zigzag maneuver (δ = 10 • / -10 • ) of the KVLCC2 L7 model. The dashed line represents the rudder angle. U, ψ, r, L pp are the absolute velocity, heading angle, yaw rate and ship length respectively.

Figure 6 .

 6 Figure 6.23: Zigzag maneuver (δ = -10 • /10 • ) of the KVLCC2 L7 model. The dashed line represents the rudder angle (δ). U, ψ, r, L pp are the absolute velocity, heading angle, yaw rate and ship length respectively.

Figure 6 .

 6 Figure 6.24: Zigzag maneuver (δ = 20 • / -20 • ) of the KVLCC2 L7 model. The dashed line represents the rudder angle (δ). U, ψ, r, L pp are the absolute velocity, heading angle, yaw rate and ship length respectively.

Figure 6 .

 6 Figure 6.25: Zigzag maneuver (δ = -20 • /20 • ) of the KVLCC2 L7 model. The dashed line represents the rudder angle (δ). U, ψ, r, L pp are the absolute velocity, heading angle, yaw rate and ship length respectively.

  Simulation results of the rudder force tests. X , Y , N m , F N , T , δ are the dimensionless longitudinal force, lateral force, yaw moment around midship, rudder normal force, propeller thrust and rudder angle respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.16 Analysis results of the rudder force tests. X , Y , N m , F N , δ are the dimensionless longitudinal force, lateral force, yaw moment around midship, rudder normal force and rudder angle respectively. . . . . . . 6.17 Computational domain and mesh for the open water test. . . . . . . . 6.18 Open water test results for the KP458 propeller. K T , K Q are the thrust and torque coefficients. J P is the advance ratio. . . . . . . . . . . . . 6.19 Vortexes (Q v =300 s -2 ) generated by the propeller with different advance coefficients: (a) J P =0.8; (b) J P =0.85; (c) J P =0.9. . . . . . . . . 6.20 Turning circle maneuver (δ = 35 • ) of the KVLCC2 L7 model. U, ψ, r, L pp are the absolute velocity, heading angle, yaw rate and ship length respectively. . . 6.25 Zigzag maneuver (δ = -20 • /20 • ) of the KVLCC2 L7 model. The dashed line represents the rudder angle (δ). U, ψ, r, L pp are the absolute velocity, heading angle, yaw rate and ship length respectively. . . Classifications of the most representative numerical methods. NS equations here include the potential flow theory, Euler equations and all other versions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2 Thresholds of the appearance of confinement effects. A s and A c are the cross sections of the vessel and waterway. h w and T d are the water depth and ship draught. B c and B s are the channel width and ship breadth. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 4.1 Basic boundary conditions. . . . . . . . . . . . . . . . . . . . . . . . 56 4.2 Possible outlet boundary conditions for the velocity (U ). . . . . . . . 57 4.3 Possible inlet and outlet boundary conditions for the pressure (p_rgh).

• . . . . xi 6.15 pp are the absolute velocity, heading angle, yaw rate and ship length respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.21 Turning circle maneuver (δ = -35 • ) of the KVLCC2 L7 model. U, ψ, r, L pp are the absolute velocity, heading angle, yaw rate and ship length respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.22 Zigzag maneuver (δ = 10 • / -10 • ) of the KVLCC2 L7 model. The dashed line represents the rudder angle. U, ψ, r, L pp are the absolute velocity, heading angle, yaw rate and ship length respectively. . . . . . 6.23 Zigzag maneuver (δ = -10 • /10 • ) of the KVLCC2 L7 model. The dashed line represents the rudder angle (δ). U, ψ, r, L pp are the absolute velocity, heading angle, yaw rate and ship length respectively. . . 6.24 Zigzag maneuver (δ = 20 • / -20 • ) of the KVLCC2 L7 model. The dashed line represents the rudder angle (δ). U, ψ, r, L xii List of tables 2.1

  . . . . . . . . . . . . .

	IBM SPH r	Immersed Boundary Method Smoothed Particle Hydrodynamics Yaw rate
	IHI SVAP r d	Hokkaido University, Japan Schiffbau-Versuchsanstalt Potsdam GmbH, Germany Radius of the actuator disk
	IIHR SVM r f IMO THOBEM Time-domain Higher-Order Boundary Element Method Iowa Institute of Hydraulic Research, USA Support Vector Machine Refinement ratio Nomenclature International Maritime Organization t Time
	INSEAN TVD t P	Instituto Nazionale per Studi ed Esperienze di Architettura Navale, Italy Total Variation Diminishing scheme Thrust deduction factor
	ITTC UKC t R	International Towing Tank Conference Under Keel Clearance Steering resistance deduction factor
	KCS VC-FVM Vertex Centered-Finite Volume Method KRISO Container Ship u Fluid velocity
	KRISO VLCC u g List of abbreviations Korea Research Institute of Ships and Ocean Engineering Very Large Crude Carrier Grid velocity LBM Lattice Boltzmann Method VOF Volume-Of-Fluid u R , v R Longitudinal and lateral inflow velocity components to rudder
	LDV XFEM v m	Laser Doppler Velocimetry Extended Finite Element Method Lateral velocity at midship in the MMG model
	ACMI LES w P	Arbitrary Coupled Mesh Interface Large-Eddy Simulation Wake coefficient at propeller position
	ALE LHS w P 0	Arbitrary Lagrangian-Eulerian Left-Hand Side Wake coefficient at propeller position in straight moving
	AMI LGA List of symbols Arbitrary Mesh Interface Lattice Gas Automata w R Wake coefficient at rudder position
	ANN LSMM x H	Artificial Neural Network Least-Squares Meshfree Method Longitudinal coordinate of the acting point of additional lateral force
	BC BEM MARIN x P a H BSHC MD x R d g MLPG x ss e a CC-FVM Cell Centered-Finite Volume Method Boundary Conditions Maritime Research Institute Netherlands Rudder force increase factor Longitudinal coordinate of propeller position Boundary Element Method Molecular Dynamics Grid displacement Longitudinal coordinate of rudder position Bulgarian Ship Hydrodynamic Center, Bulgaria Meshfree Local Petrov-Galerkin Approximate relative error Longitudinal distance between two ships CFD MMG Mathematical Maneuvering Group Extrapolated relative error x sb Longitudinal distance between ship and bridge pile e ext Computational Fluid Dynamics CFL CMT CTO NuTTS Numerical Towing Tank Symposium g Gravity acceleration Centrum Techniki Okretowej, Poland NTNU Norwegian University of Science and Technology f σ Surface tension term Circular Motion Test NS Navier-Stokes f α Rudder lift gradient coefficient Courant-Friedrichs-Lewy MRF Multiple Reference Frame f bx , f bθ Body force in axial and tangential directions of actuator disk model y inf l
	CV DCI OCIMF h T DES ONR h w DNS ONRT k t DoF OTT k w DTC PIMPLE k 0 , k 1 , k 2 DUGKS PISO l R PIV m	Control Volume Oil Companies International Marine Forum Modified water depth based on ship draught Domain Connectivity Information Office of Naval Research Water depth Detached-Eddy Simulation ONR Tumblehome Ship Turbulent kinetic energy Direct Numerical Simulation Oblique Towing Tank test Wave number Degrees of Freedom Merged PISO-SIMPLE Coefficients of the propeller Duisburg Test Case Pressure Implicit with Split Operator Effective longitudinal coordinate of rudder position Discrete Unified Gas Kinetic Scheme Particle Image Velocimetry Mass
	EFG FDM PMM m b FEM PSO m x , m y FSI RANS n b FV-LBM RAS n P RHS p	Element Free Galerkin Planar Motion Mechanism Blockage ratio Finite Difference Method Particle Swarm Optimization Added mass Finite Element Method Reynolds-Averaged Navier Stokes Number of propeller blades Fluid Structure Interaction Reynolds-Averaged Stress Propeller revolution Finite Volume-Lattice Boltzmann Method Right-Hand Side Fluid pressure
	FVM GFEM SIMPLE p C GGI HSVA SJTU SNH p s , π s q s	Finite Volume Method Semi-Implicit Method for Pressure Linked Equations Order of accuracy Generalized Finite Element Method General Grid Interface Hamburgische Schiffbau-Versuchsanstalt GmbH, Germany Shanghai Jiao Tong University Linear and angular momentum Symposiums on Naval Hydrodynamics Position vector of a rigid body
		xiv xv xvi xvii

Table 2 .

 2 1: Classifications of the most representative numerical methods. NS equations here include the potential flow theory, Euler equations and all other versions.

	Theories	Methods Classifications
			Smoothed Particle Hydrodynamics (SPH)
			Element Free Galerkin (EFG)
		Meshfree	Least-squares Meshfree Method (LSMM) Meshfree Local Petrov-Galerkin (MLPG)
			Local Radial Point Interpolation Method (LRPIM)
			• • •
		FDM	Standard Finite Difference Method General Finite Difference Method
	NS equations		Standard Finite Element Method
		FEM	Extended Finite Element Method (XFEM) Generalized Finite Element Method (GFEM)
			• • •
		FVM	Cell Centered Finite Volume Method (CC-FVM) Vertex Centered Finite Volume Method (VC-FVM)
			Dual Boundary Element Method (Dual BEM)
		BEM	Galerkin BEM
			Dual Reciprocity BEM
		LGA	Lattice Gas Automata
	Boltzmann		Standard Lattice Boltzmann Method
	equation	LBM	Finite Volume Lattice Boltzmann Method (FV-LBM)
			Discrete Unified Gas Kinetic Scheme (DUGKS)
	Many computational methods have been proposed for simulating fluid flows. The
	most representative ones are concluded in Table 2.1. The classifications here are
	based on the flow theories and spatial discretization techniques. Among them, Finite
	Element Method (FEM) is mainly used for solids. There are some flow solvers
	based on FEM, such as ADINA and ANSYS CFX. But its utilization in fluids is still
	restricted. Boundary Element Method (BEM), Finite Difference Method (FDM) and

Table 2 .

 2 2: Thresholds of the appearance of confinement effects. A s and A c are the cross sections of the vessel and waterway. h w and T d are the water depth and ship draught. B c and B s are the channel width and ship breadth.

	Parameters Start of confinement effect Important confinement Highly confined
	A c /A s	50	7-8	4
	h w /T d	15	3-4	1.5 (Shallow water)
	B c /B s	50-200	10-15	4

  Iowa Institute of Hydraulic Research, USA) and INSEAN (Instituto Nazionale per Studi ed Esperienze di Architettura Navale, Italy), where most were tested for the zig-zag and turning maneuvers, several for the spiral and stopping maneuvers. Maneuvering in shallow water was also tested in SIMMAN 2014. Detailed experimental data can be found on the website of the Tokyo 2015 CFD Workshop in ship hydrodynamics,SIMMAN 2008SIMMAN , 2014SIMMAN , 2019. CFD simulations nowadays are able to accurately reproduce the experimental results. Because of the numerical complexity and computational cost, very limited number of CFD simulations is documented in the literature, where RANS-based simulations are dominant

	(Hamburgische Schiffbau-Versuchsanstalt GmbH), CTO (Ship Hydrodynamics Di-
	vision, Centrum Techniki Okretowej, Poland), MARIN (Maritime Research In-
	stitute Nederlands), SVAP (Schiffbau-Versuchsanstalt Potsdam GmbH, Germany),
	BSHC (Bulgarian Ship Hydrodynamic Center, Bulgaria), IHI (Hokkaido University,
	Japan), IIHR (
	tested by HSVA

  Blais et al. (2016); Yi et al. (2016); Municchi & Radl (2017); McIntyre (2011); Albadawi et al. (2016); Singh et al. (2015).

Table 4 . 1 :

 41 Basic boundary conditions.

		U	p_rgh	alpha.water
	inlet	fixedValue	fixedFluxPressure fixedValue
	outlet	outletPhaseMeanVelocity	zeroGradient	variableHeightFlowRate
	top	pressureInletOutletVelocity totalPressure	inletOutlet
	cylinder movingWallVelocity	fixedFluxPressure zeroGradient
	sides	symmetry	-	-
		k	nut	omega
	inlet	fixedValue	fixedValue	fixedValue
	outlet	inletOutlet	zeroGradient	inletOutlet
	top	inletOutlet	zeroGradient	inletOutlet
	cylinder kqRWallFunction	nutkWallFunction omegaWallFunction
	sides	symmetry	-	-

Table 4 .

 4 2: Possible outlet boundary conditions for the velocity (U ).

	U	No. outlet
		1	outletPhaseMeanVelocity
	BCs	2 3	pressureInletOutletVelocity zeroGradient
		4	inletOutlet

Table 4 .

 4 3: Possible inlet and outlet boundary conditions for the pressure (p_rgh).

	p_rgh No. inlet	outlet
		1	zeroGradient	fixedValue
		2	zeroGradient	fixedFluxPressure
	BCs	3 4	zeroGradient fixedFluxPressure fixedValue zeroGradient
		5	fixedFluxPressure fixedFluxPressure
		6	fixedFluxPressure zeroGradient

Table 4 .

 4 4: Possible outlet boundary conditions for the phase fraction (alpha.water).

	alpha.water No. outlet
		1	variableHeightFlowRate
	BCs	2	zeroGradient
		3	inletOutlet

Table 4 .

 4 5: Possible inlet and outlet boundary conditions for the turbulent kinetic energy (k).

	k	No. inlet	outlet
		1	fixedValue	zeroGradient
	BCs	2 3	fixedValue turbulentIntensityKineticEnergyInlet zeroGradient inletOutlet
		4	turbulentIntensityKineticEnergyInlet inletOutlet

Table 4 .

 4 6: Possible inlet and outlet boundary conditions for the eddy viscosity (nut).

		nut No. inlet	outlet
		BCs	1 2	fixedValue zeroGradient calculated calculated
	Table 4.7: Possible velocity (U ) and pressure (p_rgh) boundary conditions for the
	wall (Figure 4.1).
	wall No. U	p_rgh
		1	movingWallVelocity fixedFluxPressure
	BCs	2 3	movingWallVelocity zeroGradient fixedValue fixedFluxPressure
		4	fixedValue	zeroGradient

Table 4 .

 4 8: Possible boundary conditions of sides (Figure 4.1).

	sides No. U , p_rgh, alpha.water, k, nut, omega
		1	symmetry
	BCs	2 3	symmetryPlane slip
		4	wall

Table 4 .

 4 9: Geometrical parameters of the KCS model.

	Parameters	Symbols	Full-scale Model-scale
	Scale factor	λ s	1	31.599
	Length between perpendiculars	L pp [m]	230	7.2786
	Length of waterline	L wl [m]	232.5	7.3570
	Maximum beam of waterline	B wl [m]	32.2	1.0190
	Draft	T d [m]	10.8	0.3418
	Displacement volume	∇ s [m 3 ]	52,030	1.6490
	Block coefficient	C b	0.6505	0.6505
	Longitudinal center of buoyancy, fwd+ L cb (%L pp ) -1.48	-1.48
	Wetted surface area without rudder	S w [m 2 ]	9530	9.5441
		K xx /B wl	0.40	0.40
	Moment of Inertia	K yy /B wl	0.25	0.25
		K zz /B wl	0.25	0.25

Table 4 .

 4 10: Grid convergence results of the total resistance coefficient. Experimental data can be found in[START_REF] Hino | CFD workshop Tokyo 2005 workshop on numerical ship hydrodynamics[END_REF].

	Parameters	Values
	Grid 1; 2; 3 r f C exp (×10 -3 ) 3.56 4,032,156; 1,706,703; 712,650 √ 2 T C 1 T (×10 -3 ) 3.554
	C 2 T (×10 -3 )	3.539
	C 3 T (×10 -3 )	3.513
	R C	0.600
	p C	1.473
	C 21 ext (×10 -3 ) 3.579
	e 21 a (%)	0.447
	e 21 ext (%)	0.666
	GCI 21 f ine (%) 1.348

Table 4 .

 4 11: Comparison of the frictional (C F ) and pressure (C P ) coefficients in the mesh convergence tests. The frictional coefficient is estimated by the ITTC 1957 friction line C F = 0.075/(log 10 Re -2) 2 . C P = C T -C F .

	Parameters	ITTC 1957 Grid 1 Grid 2 Grid 3
	C F (×10 -3 ) 2.832	2.822 2.810 2.789
	C P (×10 -4 ) 7.28	7.32	7.29	7.24
	where ε 21 C

Table 4 .

 4 12: Geometrical parameters of the KVLCC2 model.

	Parameters	Symbols	Full-scale Model-scale
	Scale factor	λ s	1	58
	Length between perpendiculars	L pp [m]	320	5.5172
	Length of waterline	L wl [m]	325.5	5.6121
	Maximum beam of waterline	B wl [m]	58	1.0000
	Draft	T d [m]	20.8	0.3586
	Displacement volume	∇ s [m 3 ]	312,622	1.6023
	Block coefficient	C b	0.8098	0.8098
	Longitudinal center of buoyancy, fwd+ L cb (%L pp ) 3.48	3.48
	Wetted surface area without rudder	S w [m 2 ]	27194	8.0838
		K xx /B wl	0.40	0.40
	Moment of Inertia	K yy /B wl	0.25	0.25
		K zz /B wl	0.25	0.25

Table 4 .

 4 13: Comparison of the frictional (C F ) and pressure (C P ) coefficients in the time step convergence tests. The frictional coefficient is estimated by the ITTC 1957 friction line C F = 0.075/(log 10 Re -2) 2 . C P = C T -C F . ∆t is the time step.

	Parameters	ITTC 1957 ∆t 1	∆t 2	∆t 3	∆t 4
	C F (×10 -3 ) 3.453	3.315 3.313 3.340 2.101
	C P (×10 -3 ) 0.657	0.973 0.935 0.951 1.550

Table 4 .

 4 14: Case setups of the trim and sinkage predictions for the KCS model.

	Case No.	1	2	3	4	5	6
	U in [m/s] 0.914	1.281	1.647	1.922	2.196	2.380
	F r	0.108	0.152	0.195	0.227	0.260	0.282
	Re	6.66×10 6 9.32×10 6 1.20×10 7 1.40×10 7 1.60×10 7 1.73×10 7

Table 5 .

 5 3: Mesh numbers of the convoy tests when the water depth h w =0.24 m. T d , W b and N mesh are the ship draught, channel bottom width and mesh number respectively.

	s

Table 5 .

 5 4: Test configurations and mesh numbers of the head-on encounter case. The units of the channel width W b , draught T d and water depth h w are [m]. The unit of the maximum velocity V max is [m/s]. T d is the draught of the convoy 2. N mesh is the mesh number.

	Case W b	T d	h w	V max F r max h	N mesh
	a1 a2	1.44	0.04 0.18 0.91 0.1 0.18 0.793 0.597 10,258,709 0.685 10,134,445
	b1 b2	2.88	0.04 0.18 0.91 0.1 0.18 0.91	0.685 15,221,105 0.685 15,345,369

Table 5 . 5 :

 55 Test configurations and mesh numbers of the convoy 2 passing the bridge piles. The units of the channel bottom width W b , distance between bridge piles D b , ship draught T d and water depth h w are [m]. The unit of the maximum velocity V max is [m/s]. N mesh is the mesh number.

	Case W b D b	T d	h w V max F r max h	N mesh
	a1				0.12 0.80	0.738 1,284,531
	a2			0.04	0.18 0.91	0.685 1,423,626
	a3 a4 a5 b1	1.44	0.7	0.1	0.24 0.91 0.18 0.80 0.24 0.91 0.12 0.80	0.593 1,426,806 0.602 1,340,223 0.593 1,437,711 0.738 1,278,281
	b2			0.04	0.18 0.91	0.685 1,419,086
	b3 b4 b5		0.8	0.1	0.24 0.91 0.18 0.80 0.24 0.91	0.593 1,534,598 0.602 1,331,889 0.593 1,471,276

Table 6 .

 6 1: Principle parameters of the KVLCC2 tanker. of the captive model tests. The parameters of this ship can be found in Table6.1. The virtual captive model tests are carried out at propeller condition with the rudder. Ship speed U 0 is set to 0.76 m/s (15.5 kn in full-scale). Dimensionless variables are used in some cases and the Prime system II is adopted for their nondimensionalizations (TableA.1 in Appendices).

	Parameter	Symbol	L3 model L7 model Full-scale
	Scale ratio	λ s	1/110	1/45.7	1.00
	Length between perpendiculars	L pp [m] 2.902	7.00	320.0
	Ship breadth	B wl [m] 0.527	1.27	58.0
	Ship draught	T d [m]	0.189	0.46	20.8
	Displacement volume	∇ s [m 3 ] 0.235	3.27	312,600
	Longitudinal gravity center position x G [m]	0.102	0.25	11.2
	Block coefficient	C b	0.810	0.810	0.810
	Propeller diameter	D P [m] 0.090	0.216	9.86
	Rudder span length	H R [m] 0.144	0.345	15.80
	Profile area of rudder movable part	A R [m 2 ] 0.00928	0.0539	112.5

Table 6 .

 6 2: Setups of the OTT (Oblique Towing Tank) and CMT (Circular Motion Test).

	Parameters	Symbols	Values
	Inlet velocity	U in [m/s] 0.76
	Ship drift angle β [ • ]	0, 4, 8, 12, 16, 20
	Thrust	T [N ]	4.350, 4.321, 3.810, 4.054, 3.735, 3.784
	Torque	Q [N m]	0.0884, 0.0872, 0.0822, 0.0845, 0.0815, 0.0820

Table 6 . 3 :

 63 Hydrodynamic derivatives and added masses of the hull.

	Coefficient	Exp	Sim	Coefficient	Exp	Sim
	R 0	0.022 0.019 Y vrr	-0.391 -0.548
	X vv	-0.040 -0.061 Y rrr	0.008 0.059
	X vr + m + m y 0.518 0.525 N v	-0.137 -0.139
	X vr	0.002 0.007 N R -x G m -0.059 -0.070
	X rr + x G m	0.021 0.009 N R	-0.049 -0.060
	X rr	0.011 -0.001 N vvv	-0.030 -0.027
	X vvvv	0.771 0.730 N vvr	-0.294 -0.143
	Y v	-0.315 -0.355 N vrr	0.055 0.099
	Y R -m -m x	-0.233 -0.218 N rrr	-0.013 -0.002
	Y R	0.083 0.084 m x	0.022 -
	Y vvv	-1.607 -1.161 m y	0.223 -
	Y vvr	0.379 0.288 J z	0.011 -

Table 6 .

 6 4: Hydrodynamic coefficients obtained from the rudder force tests. R (β R < 0) γ R (β R > 0) l R Propeller open water tests 6.3.3.1 Computational setupsPropeller open water test is used to determine the coefficients k 0 , k 1 and k 2 . The open water test of a KP458 propeller is carried out to obtain the curves of the thrust and torque coefficients, which are useful to determine the thrust and torque values for the actuator disk model. The geometrical parameters of the propeller can be found inTable 6.5.

	Coefficients t R	a H	x H
	Exp	0.387	0.312	-0.464
	Sim	0.123	0.284	-0.232
	Coefficients f α	x R	κ	ε
	Exp	2.747	-0.5	0.50	1.09
	Coefficients γ Exp 0.395	0.640	-0.710

Figure 6.17: Computational domain and mesh for the open water test.

6.3.3

Table 6 .

 6 5: Geometrical parameters of the KP458 propeller.

	Parameters	Symbols Values
	No. of blades n b	4
	Diameter	D P [m] 0.090
	Pitch ratio	P P /D P	0.721
	Area ratio	A e /A 0	0.431
	Hub ratio	d h /D P	0.155
	Rotation		Right hand

Table 6 .

 6 6: Hydrodynamic coefficients of the propeller.

	Coefficient	k 0	k 1	k 2
	Exp	0.2931 -0.2753	-0.1385
	Multi-run (Sim) 0.2754 -0.2595	-0.1334
	Single-run (Sim) 0.2720 -0.2513	-0.1347
	Coefficient	t P
	Exp	0.220
	Sim	0.209
	Coefficient	C 1	C 2 (β P < 0)	C 2 (β P > 0)
	Value	2.0	1.1	1.6
	Coefficient	x P	w P 0 (L7 model) w P 0 (full-scale)
	Value	-0.50	0.40	0.35

  .1). The initial velocity U 0 is 0.76 m/s. Propeller revolution is assumed to keep the U 0 constant. Turning circle tests (δ = ±35 • ) and zigzag tests (δ = ±10 • / ± 20 • ) are simulated using the MMG model in Section 6.2.2. The free running model test experiments were carried out by MARIN (Maritime Research Institute Netherlands) and used for comparison here). The dimensionless trajectory, absolute velocity (U = √ u 2 + v 2 ), heading angle (ψ)

& Incecik (2016) and[START_REF] Torben | The development of a method to predict ship-bank interaction in real-time[END_REF].

Acknowledgments

Concluding remarks

In this chapter, the details of the numerical solver were explained. Traditionallyused RANS solver and VOF multiphase model were presented. The calculations of kinetic energy and turbulent dissipation in turbulence models in this thesis were shown.

The equation discretization procedure of FVM was explained. The deduction of the pressure Poisson equation and the pressure-velocity coupling were included because of their extraordinary importance in CFD.

Dynamic mesh techniques including the sliding mesh, deforming mesh, overset, layer addition/removal and IBM were summarized. The background theories of these methods and their usage for prescribed or 6DoF motions were explained in detail. The 6DoF simulation procedures were presented, where the most important parts, coordinate transformation and numerical integration methods were emphasized. The implementation procedure in OpenFoam was listed at last. as in [START_REF] Carrica | Self-propulsion computations using a speed controller and a discretized propeller with dynamic overset grids[END_REF] are not captured. This is a trade-off between the flow details, accuracy and computational costs. The values of k 0 , k 1 , k 2 are calculated according to 6.25 and shown in Table 6.6. They are captured accurately by both the multi-and single-run approaches. The single-run case seems less accurate because the initial evolving process influences the results at the beginning. The determination of other propeller coefficients are introduced below:

• The calculation of the thrust deduction factor t P is calculated according to equation 6.49. The R 0 value has already calculated. The torque (T ) can be obtained from Figure 6.15. It can be seen that it is not affected by the rudder angle. The thrust deduction coefficient is calculated in Table 6.6.

Concluding remarks

In this chapter, the maneuvering analysis was carried out based on the MMG model, which separates the hull, rudder and propeller forces. The formulations and procedures were presented in detail.

Virtual captive model tests were conducted using a L3 KVLCC2 model. The actuator disk model was implemented into OpenFoam to replace a real propeller, which is beneficial to reduce the computational cost. It was validated using a case with real propeller. Open water test, rudder force test, OTT and CMT were carried out to obtain the hydrodynamic coefficient of the propeller, rudder and hull.

Based on the obtained hydrodynamic coefficients, maneuvering simulations were conducted and compared with free running test results. The agreement proves the validity and accuracy of our model and procedure.

Appendix A

Standard normalization systems

The most commonly-used normalization systems are shown in Table A.1. They are used for the non-dimensionalization of variables. The Prime system II is adopted in this thesis. [START_REF] Fossen | Guidance and control of ocean vehicles[END_REF]).

Unit

Prime system I Prime system II Bis system Length L L L Mass