In solving two-phase flows, the location of the interface between the phases is necessary to handle interface jump conditions when solving the Navier-Stokes equation. Current interface capturing and advection methods, however, suffer from various issues. The level set method uses the signed distance to the interface and the interface being the zero level set of this function allows the evolution of the level set field to be described with a simple advection equation. This means that no additional steps are required, but solving the advection equation generally does not conserve mass. On the other hand, Volume of Fluid methods utilise the local fluid fractions to represent the phase interface. For these methods, a mass conserving advection algorithm exists, but the absence of an explicit interface requires expensive reconstruction methods to be used instead. Additionally, this Volume of Fluid advection method is subject to a restrictive CFL condition on the time-step and, being dimensionally split, requires a structured grid to be used. Other Volume of Fluid or Moment of Fluid advection methods that do not have these conditions are not mass conserving. Dual interface methods that combine information from the level set and volume fractions are able to achieve higher accuracy, but these method still use Volume of Fluid advection to remain mass conserving and are thus subject to the same conditions. These methods use the level set field to create a cheaper and better justified reconstruction method, which negates one drawback of VoF methods. In this thesis, a method is formulated to allow interface advection on unstructured grids without such a strict CFL condition compared to the MCLS method. To do this, the finite volume level set advection method of the MCLS method is replaced with a nodal-modal discontinuous Galerkin method. The DG method is analogous to the Galerkin finite element methods, but the basis functions are now only valid on one element, and the solution can be multiply defined on the cell boundaries. This allows for level set advection on polygonal cells, and has the added benefit that since cells are semi-independent, level set correction now only has to be applied locally. This, like the finite volume level set advection method, is not necessarily mass conserving though, so the volume of fluid will need to solve this issue. The difference with MCLS is that no volume of fluid advection method is used, since these methods are the parts that introduce the aforementioned problems that are being avoided. Instead, a minimisation method is done on the VOF field in order to obtain a correction which can be applied to the level set and keep it mass conserving. For this minimisation method, a flux condition is used to impose additional constraints on the solution. This condition effectively attempts to match the interface intersections on a cell edge for both adjacent cells, for all cell edges that have a flux going through. The resulting method is tested for four test cases, with solid body translation and rotation, the corner flow test, and the vortex deformation test. In the translation test, there is only vertical velocity, so the flux condition does not take vertical edges into account. To stabilise this test case, level set reinitialisation is applied, which is a technique used to maintain the signed distance property of the level set function to use in the Navier-Stokes evolution method. This keeps the solution in the translation test somewhat representative of the exact interface, but makes the obtained results no longer representative of the actual method. For the other test cases, the method appears to converge when the grid is refined, except for the solid body rotation test using Euler Forward, as this result becomes unstable. However, no apparent order of convergence is present, and the effect of a more accurate time integration method is very inconsistent. To solve this, the method is altered to allow multiple level set advection steps for every optimisation routine. This makes the method better follow the order of convergence of the discontinuous Galerkin method, although this order is slightly lower for the vortex deformation test. The obtained results show that in it's current form, this optimisation based method is likely not usable for general applications, but that with additional work, this kind of method may prove to be more applicable than MCLS.

The company Deepwater Energy BV, Netherlands developed a unique water turbine called as the Oryon Watermill (OWM). It was the aim of the company to obtain an efficient numerical model for the turbine. A numerical model, which serves as a good starting point, has already been developed by Maniyara [1]. This model suffered from certain limitations and to tackle these it was realized
that a smaller and simpler problem needs to be solved. The current thesis deals with developing a numerical model for this smaller problem.
In this problem, an airfoil is hinged at its leading edge and rotating under the influence of a periodic water flow. Further, the rotation of the airfoil is restricted by a stopper. This restriction leads to a collision between the airfoil and the stopper due to which the motion of the airfoil is changed suddenly. This is a Fluid-Structure Interaction (FSI) problem. To model the structural part i.e. the airfoil motion, a non-smooth dynamical analysis is needed.
Then, the entire FSI problem is solved in a faster and efficient way using the harmonic balance method. For this purpose, the non-smooth dynamics of the airfoil have to be solved using the harmonic balance method. In the thesis, a novel approach is developed to solve the non-smooth dynamics problem using the harmonic balance method. Initially, this new approach is implemented
and verified in MATLAB on a 1D level. Later, a MATLAB-OpenFOAM coupled solver is developed to compute the solution to the airfoil problem using this newly developed approach.

The accurate approximation of the surface tension force is paramount for continuum surface models in the field of computational fluid dynamics for multiphase flow where surface tension is relevant. This involves being able to accurately calculate the curvature at the interface. This study focuses on the use of convolution in smoothing the VOF colour field in order to obtain better approximations of the curvature. Given the sudden jump in values of the VOF colour field, the calculation of its derivative for the curvature is sensitive to errors, given the large values of high order terms that determine the truncation error. To deal with this problem, convolution of this abruptly varying field can be used to create a smoother transition. The curvature approximation of a circular interface improved as the support of the convolution was increased.
It was proven analytically that, for these interfaces, the original curvature is retrieved from the convoluted field. Interfaces along which the curvature varies were also considered, and it was found that there is a critical convolution support that minimizes the error in the curvature, given that the choice of the support length can modify the curvature that is estimated.
An algorithm was implemented in OpenFOAM that calculates the convolution of the VOF colour field. The resulting smoothed field was then used to calculate the curvature, which is needed for the surface tension force of the system. The simulations of a two-dimensional rising bubble resulted in more accurate results for the circularity and the rising velocity, when compared to the original OpenFOAM implementation with no smoothing. With the convolution algorithm, the terminal velocity deviated only 0.01% from a well-accepted benchmark case, a great improvement when compared to the 4.2% difference when no smoothing was used. However, simulations of a static bubble in zero-gravity rapidly resulted in unphysical flow, manifested as a wavy interface, when a convolution support larger than 2 cells was chosen. An improvement of the estimation of the surface tension force direction may be needed for this behaviour to disappear.

For many types of floating structures, roll motion is the most important wave induced motion. More specifically, roll motion is of high significance for barges, in order to correctly predict the acceleration of the structure caused by roll motion and determine the procedure of sea-fastening and off-loading. The correct forecast of roll motion requires accurate estimation of roll damping. At the same time, the main sources of roll damping are the creation of waves, skin friction and creation of vortices (eddies), due to roll motion. For this reason, the physics of the problem cannot be fully described by a flow model based on potential theory, as it is highly dependent on vorticity and viscous effects. As a consequence, potential theory algorithms underestimate roll damping, which results in over prediction of roll motion and thus in a conservative estimate of workability. It is also important to mention that the vorticity and the viscous effects, lead to nonlinear damping characteristics. More specifically, the roll damping moment is controlled by its odd numbered harmonics, the first of which is dominant. Consequently, it is common to express the damping moment in an equivalent linearized form, equal to the first harmonic.
During the last years simulating the flow around the rolling vessel using viscous flow algorithms is becoming more and more popular. Using CFD (Computational Fluid Dynamics) is a cheap and fast way to create a “numerical” wave tank and perform numerical decay tests, forced roll simulations, and roll response simulations in regular waves. Of course the underlying algorithm should be validated, before any commercial or scientific use. In this thesis, roll damping will be estimated by performing virtual forced oscillation tests, and the computed roll damping coefficients will be presented as a function of roll amplitude. The virtual forced oscillation tests have been performed with the open-source CFD software OpenFOAM. Additionally, numerous numerical experiments are performed in order to choose the optimum discretization schemes, turbulence model, mesh and time configurations. Finally, Ikeda’s experimental data is used in order to validate the viscous flow algorithm and the numerical model. Two methodologies have been used in order to determine the nonlinear roll damping. In the first case the free surface is included in the viscous flow model, using the Volume of Fluid method. The second approach disregards any free surface effect and the total damping is calculated as a superposition of viscous (by viscous flow algorithm) and wave (by potential theory algorithm) damping. Both approaches are compared with Ikeda’s experimental data and it is concluded that both methods are able to capture accurately the linearized roll damping coefficients, for various amplitudes. Finally, viscous flow simulations of the full scale barge, in order to calculate the roll damping coefficients, are notably time consuming. For this reason, in order to reduce the computational time, the methodology which neglects the free surface effect, is chosen, as it is a good compromise between time efficiency and accuracy.

The Interface Capturing method, which is a finite volume method as formulated by Queutey and Vissoneau for free surface immiscible, incompressible multiphase flows employs a collocated arrangement of unknowns
and achieves a discrete force balance for the case when the interface coincides with the faces of the control volumes in the computational domain. This constraint limits the applicability of the method. Furthermore,
the authors do not provide the exact formulation of the operators involved in the pressure velocity coupling. In the present research, a balanced-force numerical method is formulated, applicable for an interface that
neither has to coincide nor be aligned with the faces of the control volumes. The approach consists of the reconstruction of the values of the flow variables at the interface based on the interface jump conditions, with which the limit values of the normal derivatives at the interface are calculated. Furthermore, the construction of the operators of the discrete system is delineated to achieve a discrete force balance, by incorporating the reconstructed flow variables and employing a discretization which complies with the interface jump conditions. It is sufficient for a stationary discrete formulation to comply with the differential equation and the interface jump conditions. However, to apply this approach to solve unsteady flow problems the influence of the reformulated operators on the stability properties of the system should also be investigated.
The properties of the individual operators are analyzed as well as their behaviour when they are embedded in the complete solver algorithm. Results are shown for both steady and unsteady test cases and compared with numerical results obtained with OpenFOAM. The resulting framework avoids the occurrence of
spurious velocities as it discretely complies with the interface conditions.

The current research aims to numerically predict the performance characteristics of the Oryon Water Mill. The simulation mechanism to be developed, begins at the structural solver part of a partitioned FSI solver, upon which the unique functional characteristics of the turbine are embedded. The modelling is carried out to account for the turbine’s functional characteristics like, multiple nested rigid body FSI dynamics, non-smooth dynamics and lamella partial torque contributions. The solver is built to resolve the turbine dynamics in two dimensions and re-scale the final result with respect to the turbine core height.
The results of the simulation show the predicted torque to be of similar trend as that of the experimental data collected by MARIN, while being an over-prediction by a factor of two. The optimal performance range of the turbine is predicted within an error bound of 0.1RPS. The overprediction of the torque is considered to be a result of the 2D solution not accounting for the flow leakage in the third and unaccounted for direction. The time variant fluid solution is found to be unstable due to the spikes observed in the torque time signal and is a consequence of the perfect momentum sink assumption employed in the modelling of the lamella’s non-smooth dynamics. The lamella dynamics are observed to not convergence to fully periodic behaviour and the sensitivity of the lamella’s dynamics is considered to be the root cause. The velocity Verlet based structural solver does not provide fluid-structure coupling iterations, the inclusion of which, could improve the stability of the fluid solution.
This endeavour establishes a foundation for the predictive numerical simulation of the OWM in the form of baseline numerical simulation results and a developed simulation framework.

This research entails the numerical simulation of physical flow instabilities that can occur in two-phase pipe flows with a new efficient algorithm. The fluids are assumed to be immiscible, and the flow is incompressible and isothermal in a straight circular pipe section with a certain inclination. The numerical algorithm that was developed consists of a flow solver and a sharp interface model that solve the Navier-Stokes equations in cylindrical coordinates. The simulation results obtained with the new method are validated through comparison with other models and with experiments. The focus lies on obtaining an accurate and efficient algorithm that can ultimately be used for Direct Numerical Simulations or Large Eddy Simulations of turbulent two-phase pipe flows.