A Cartesian cut-cell based multiphase ﬂow model for large-eddy simulation of three-dimensional wave-structure interaction

A multiphase ﬂow numerical approach for performing large-eddy simulations of three-dimensional (3D) wave-structure interaction is presented in this study. The approach combines a volume-of-ﬂuid method to capture the air-water interface and a Cartesian cut-cell method to deal with complex geometries. The ﬁltered Navier–Stokes equations are discretised by the ﬁnite volume method with the PISO algorithm for velocity-pressure coupling and the dynamic Smagorinsky subgrid-scale model is used to compute the unresolved (subgrid) scales of turbulence. The versatility and robustness of the presented numerical approach are illustrated by applying it to solve various three-dimensional wave-structure interaction problems featuring complex geometries, such as a 3D travelling wave in a closed channel, a 3D solitary wave interacting with a vertical circular cylinder, a 3D solitary wave interacting with a horizontal thin plate, and a 3D focusing wave impacting on an FPSO-like structure. For all cases, convincing agreement between the numerical predictions and the corresponding experimental data and/or analytical or numer- ical solutions is obtained. In addition, for all cases, water surface proﬁles and turbulent vortical structures are presented and discussed.


Introduction
It is believed that extreme waves will become more common in coastal and offshore region due to climate change [1] . Detailed engineering understanding of wave-structure interaction (WSI) is a key aspect in the safe and cost-effective design of coastal and offshore structures, and marine renewable devices. Reliable numerical tools are required to predict WSI and assess the reliability and survivability of these structures due to extreme wave loads. However, not too much progress has been made to date in terms of developing and validating tools that can predict accurately three-dimensional waves interacting with complex structures, owing to nonlinearities of waves and the geometrical complexity of immersed structures. Further challenges arise in terms of the turbulence of the flow and breaking of waves before or during wave impact, with subsequent periodic air entrainment and jet-splash up [2] .
In order to roughly predict hydrodynamic loads on structures, empirical or semi-empirical methods such as the Morison equation [3] or the Froude-Krylov method [4] have been used in engineering applications. However, these methods ignore the effect of the structure and are only applicable to very simple problems.
There are a variety of mathematical models for water waves, such as the mild slope equation [5] , shallow water equation [6] , and Boussinesq equation [7,8] . In order to obtain the kinematics and dynamics of water waves, depth-resolved methods should be used. One of these is the fully nonlinear potential flow model based on Laplace's equation with inviscid and irrotational assumptions [9][10][11] . However, it is challenged when considering wave impact on structures or for wave breaking, especially when there is splash-up and/or air entrainment.
In order to immerse complex geometries for wave-structure interaction simulations, overlapping grids, boundary-fitted grids, and unstructured grids can be used. These methods provide great flexibility conforming with complex stationary or moving boundaries. However, the programming of these methods can be complicated and generating such grids is usually very cumbersome [49] . Cartesian grid methods, which can simulate flow with complex geometries on Cartesian grids, avoid these problems, such as the immersed boundary method [49] and the Cartesian cut cell method [50] . The primary advantage of Cartesian grids is that only little modification of the NSE solver is needed to account for immersed complex geometries. Examples are the partial cell treatment [51,52] , immersed boundary method with ghost cell [16] , a virtual boundary force [14] , direct forcing immersed boundaries [17] , and Cartesian cut-cell method [53,54] . Solving the NSE on fixed Cartesian grids has also the advantage of simple grid generation and allows simulating moving boundaries thereby avoiding heavily deforming grids [19,51,55,56] .
Most wave-structure interaction in engineering applications are turbulent and therefore the effect of turbulence on the mean and instantaneous flows requires to be accounted for unless all scales of turbulence are resolved, which is, for practical applications not feasible due to the exorbitant computational demands of the socalled Direct Numerical Simulation (DNS). Therefore, in many previous studies, the effects of turbulence have either been neglected or the Reynolds-averaged Navier-Stokes (RANS) equations have been solved, in which all of the unsteadiness due to turbulence is averaged out and the effects of turbulence are modelled by a so-called turbulence model. Thus, RANS models cannot provide instantaneous turbulent flow characteristics. The increase in computer power has led to the development of more powerful but more computationally demanding methods, such as the method of large-eddy simulation (LES) [57][58][59] , in which large-scale turbulence is resolved while the effect of small scale turbulence on the large scales is modelled. LES can be employed for practical problems where the Reynolds number is high and the computational domain is large. Mo et al. [60] performed single-phase LES to study wave impact on a slender cylinder. To date, two-phase flow LES for three-dimensional (3D) wave-structure interaction is rather limited [55,61] .
The objective of this paper is therefore, to refine and validate a 3D two-phase flow code (Xdolphin3D) [15,62] for large-eddy simulations of 3D wave-structure interaction using a Cartesian cut-cell method [56] , with the aim to compute accurately various quantities for wave-structure interaction, such as wave elevations, pressure fields, wave loads, water surface profiles, and turbulent vortical structures at high spatial and temporal resolution. To the best of the authors' knowledge, the combination of a volume-of-fluid method for two-phase flows, Cartesian cut-cell method to deal with complex geometries, and large-eddy simulation, has not been reported in the literature for 3D wave-structure interaction problems. Four different 3D test cases are considered in the present study and simulation results are compared with available experimental measurements and/or results of other numerical simulations.
The organisation of this paper is as follows: The description of the mathematical model for the two-phase flow and the numerical method are presented in Section 2 . The versatility, robustness and accuracy of the present two-phase flow model is demonstrated by applying it to various 3D wave-structure interaction problems involving complex geometries in Section 3 . Finally, the paper ends with conclusions in Section 4 .

Governing equations and basic solver
The large-eddy simulation (LES) approach [59] is adopted in this study, for which the large-scale turbulence is resolved and a subgrid-scale model is employed to compute the unresolved scales of turbulence. The governing equations used for incompressible two-phase flow are based on the spatially filtered Navier-Stokes equations, given as: where the overbar · denotes the spatial filtering over the grid in Cartesian coordinates ( x, y, z ), ū = ( ū , v , w ) is the filtered velocity vector and p is the filtered pressure. t is the time, g is the gravitational acceleration vector, ρ and μ are the density and dynamic viscosity of the fluid.
The term τ sgs = ρ( ū ū − u u ) is the subgrid-scale (SGS) stress tensor and the anisotropic part of the SGS term is modelled by an eddy-viscosity model of the form [63] : where I is the unit tensor and S is the strain rate tensor of the resolved field. μ t is the turbulent eddy viscosity defined as: with the cut-off length scale ¯ = ( x y z ) 1 / 3 and the coefficient C d is calculated by the dynamic SGS model [64] in the present study. The symbol for spatial filtering ' − ' is dropped hereinafter for simplicity.

Finite volume discretisation and computational grid
The finite volume method is conservative and can deal with complex geometries, thus it is especially suitable for simulating two-phase flows and their interaction with structures (see Fig. 1 (a)) and is therefore adopted. A staggered Cartesian grid is used to discretise the governing equations, which has the advantage of strong coupling between the velocity and the pressure. Fig. 1 (b) shows a typical variable arrangement in a 3D Cartesian grid, in which the velocities are located in the centre of the face of the control volume, and the pressure, as well as all other scalar variables are stored in the cell centre. The air-water interface is captured by the volume-of-fluid method to be introduced in Section 2.2 . Complex geometries as shown for a 2D plane ( Fig. 1 (c)) and a 3D control volume ( Fig. 1 (d)) are handled with the Cartesian cut-cell method to be introduced in Section 2.3 .

Temporal discretisation
A backward difference scheme is used for the time derivative, which leads to an implicit scheme for the Navier-Stokes equations, such as the first-order backward Euler scheme: or the second-order backward Euler scheme: where t is the time step and superscripts n + 1 , n and n − 1 represent current, previous one and previous two time step, respectively. The subscript c indicates the current control volume and is the volume of the control volume of the fluid. The implicit scheme has the advantage of the solution being unconditionally stable and thus can prevent instabilities in small cut-cells as is discussed in Section 2.3 .

Spatial discretisation
The finite volume discretisation of the advection term is obtained as: where the subscript f denotes the corresponding face of the control volume, A is the area of the face available for the fluid and m = ρu · n A is the mass flux through the face.
It is worth mentioning that it is attractive to use a consistent and conservative form for the momentum advection method [62,65,66] . In this study, the m f is obtained from the interpolation of the mass fluxes, which is already available at the faces of the continuity control volumes. u f is approximated with the highresolution scheme with flux limiting and more details can be found in [62] .
The finite volume discretisation of the diffusion term is obtained as: where the viscosity on the face is obtained by the harmonic mean, ∂ u ∂n is calculated by the finite difference approach, τ w is the shear stress of the wall on the face of the control volume, and A solid is the area occupied by the solid. The finite volume discretisation of the source term is obtained as: where the value in the centre of the control volume ( ρ c ) is obtained by the volume averaging of two values on the face of the control volume and the pressure gradient is calculated from two neighbouring control volumes.
Substituting all the approximated derivatives into Eq. (2) and subtracting the continuity equation ∂ρ/ ∂t + ∇ · (ρu ) = 0 multiplied by u n +1 P , leads to where a u is the coefficient, b u P is the source term, the subscripts P and nb denote the variables in the present and neighbouring control volume, respectively.

Pressure-velocity coupling
The PISO algorithm [67] is employed for pressure-velocity coupling which calculates the pressure correction term twice. This term is used to update pressure and velocities at the end of each time step to satisfy the continuity equation. A brief summary is shown below.
For a guessed pressure p * , the discretised momentum equations can be solved to produce fluid velocities u * , which satisfy To obtain the pressure correction, the updated fluid velocities are substituted into the discretised continuity equation Eq. (1) and the resulting pressure correction equation has the following form: where the term b P , is the mass residual from the predicted velocities, p P is the pressure correction and a p P , nb are the matrix coefficients for this equation.
In the PISO algorithm [67] , a second correction step is introduced as: where the coefficients are the same as in the first pressure correction equation shown in Eq. (12) , however the mass residual is based on the value of first velocity correction u . After solving the first and second pressure correction equations, the quantities of the current time step are updated as: where

Interface capturing method
The volume-of-fluid method is employed here to capture the air-water interface in the two-phase flow solver during the simulation. F is the volume fraction defined as: F = 1 , if a control volume contains water only; 0 , if a control volume contains air only.
The air-water interface is then in cells where 0 < F < 1. A particle on the surface stays on the surface and the volume fraction F has a zero material derivative: After interface capturing for the volume fraction field, the momentum equation ( Eq. (2) ) is closed with the constitutive relations for the density and dynamic viscosity of the fluid as given by: where the superscripts 'w' and 'a' denote water or air, respectively.
A key requirement for simulating two-phase flows is to track/capture the shape of the interface. Numerous methods have been proposed and used for the simulation of free-surface/twophase flows. However, the VOF method for capturing the interface is one of the most popular approaches due to its advantages: mass conservation, computational efficiency and easy implementation. From a general point of view, there are two classes of algorithms to solve the F transport equation ( Eq. (17) ): algebraic and geometric computation [68] . Considering the efficiency in algebraic computation, the high resolution VOF scheme Compressive Interface Capturing Scheme for Arbitrary Meshes (CICSAM) is employed to capture the air-water interface for two-phase flows. CICSAM is a high resolution scheme based on the normalised variable diagram used by Leonard [69] . It contains two high resolution schemes to calculate the volume fraction F on the face of the control volume and the weighting factor is based on the angle between the interface and the direction of motion. Refer to [15,70] for more details.

Cartesian cut-cell method
The cut-cell method proposed by Xie and Stoesser [56] is employed here. The geometry of the solid is described by a level set function LS ( x, y, z ), in which the boundary of the solid is represented by LS(x, y, z) = 0 while the fluid domain has the value of LS ( x, y, z ) > 0 and the solid domain is LS ( x, y, z ) < 0. The cut-cell interface between the fluid and solid is calculated by a piecewise linear interface, which is a straight-line in 2D and a plane in 3D (as shown in Fig. 1 ). For each 3D Cartesian grid cell, the area of each side surface A and the total volume available for the fluid is calculated. A θ function is defined here, the value of which is 1 for a point accessible to fluid and 0 for a point inside a solid. The average of this function over a control volume or cell face is the fraction of the volume or area available to the fluid in a standard control volume.
Regarding numerical implementation, no modification of the equations in Section 2.1 are required for full fluid cells. In cut cells around the structures, the mass flux has to be modified by the θ function on the boundary. If A f = 0 (such as the cells in the left of the cylinder in Fig. 1 (c)), there is no mass flux through the east face and the advective flux is obtained as m f = 0 . In addition, the diffusive flux as well as the cell volume have to be modified in cut cells. The details for the Cartesian cut-cell method can be found in [56] . It is worth mentioning that an implicit scheme for time integration is employed to prevent numerical instabilities in small cut cells. The cut-cell method has also been validated for some LES studies for turbulent open-channel flows [71][72][73] .

Initial and boundary conditions
Initial and boundary conditions are as follows. Water waves are generated at the inlet based on wave theories [74] and velocities over time are prescribed as Dirichlet boundary condition. The time history of the water surface elevation is prescribed in terms of the volume fraction at the inlet, i.e. in cells above the interface F = 0 and in cells below the interface F = 1 whilst cells containing the interface the volume fraction is calculated from the water surface elevation. On the solid boundary, the wall model [75] is used for the near-wall treatment in the LES. At the outlet, a radiative boundary condition is applied and a damping zone can also be employed. As both fluids, i.e. air and water, are solved simultaneously, the kinematic and dynamic free surface boundary conditions are inherently implemented at the air-water interface. In the computation, the initial flow field as well as the hydrostatic pressure at t = 0 are prescribed depending on the flow to be solved.

Results and discussion
A convergence study of the two-phase flow approach is carried first, by solving the benchmark problem of a 3D travelling solitary wave without a structure with the goal to validate the Navier-Stokes solver together with the interface capturing method at various spatial and temporal resolutions. The validation of the cutcell method is then demonstrated by studying several 3D wavestructure interaction problems, such as a solitary wave traveling over a vertical circular cylinder, a solitary wave traveling over a horizontal thin plate, and a focusing wave impacting on an FPSOlike structure. The numerical results are compared with available experimental data and detailed 3D water surface profiles and turbulent vortical structures are visualised.

Convergence study -solitary wave run-up
A solitary wave propagating at constant water depth in a canal [76] is simulated, considered to be a classical benchmark problem for wave simulations with the aim to compute viscous damping and wave run-up on a vertical wall. For this case, the turbulence model is switched off and only the molecular viscosity is included to quantify the amount of artificial viscous damping of the approach.
The computational setup of this problem is shown in Fig. 2 , where h is the still water depth and the domain size L x × L y × L z is 20 h × 2 h × h in the streamwise, vertical, and spanwise direction, respectively. The same parameters as used in [76] are considered, where the theoretical wave speed C w = gh = 1 . 0 m / s and the Reynolds number is Re = 5 × 10 4 . The no-slip boundary condition is applied at the left, right, bottom and top walls of the domain while the slip boundary condition is applied for the walls in the spanwise direction. In the simulations, the quantities A c and A runup , which represent the wave amplitude at the centre of the domain and the wave run-up amplitude at the right wall, are used to assess the performance of the method. Initially, the water surface and zero velocity in the entire domain are prescribed with a Boussinesq profile at the left wall as: where A 0 is the maximum elevation at the left wall. The initial wave moves towards the right due to gravity and can be considered as a solitary wave after 6.0 s and this time is set as t = 0 for the remaining propagation of the solitary wave. The computational domain is discretised by a uniform grid of 400 × 200 × 20 in the streamwise, vertical, and spanwise directions, respectively. A small The results in the finest mesh are used as reference.
time step t = 0 . 001 s is used with the first-order backward Euler scheme and the simulation is run for a total of 20 s. Fig. 3 (a) shows snapshots of the solitary wave propagating along the channel for the case A 0 /h = 0 . 4 at various instants in time. There is slight wave damping due to the viscous effects and the predicted wave speed is 1.05 m/s, which is within 5% of the theoretical value and agrees with the previous study by [76] . The energy transformation and dissipation compared to the initial state are plotted in Fig. 3 (b). After release of the water column near the left wall, some of the initial potential energy is converted into kinetic energy until the solitary wave is fully established. After that, both potential and kinetic energy remain almost constant until kinetic energy transforms into potential energy during the wave runup near the right wall until the wave reaches its maximum height. The total energy decreases due to viscous damping and friction.
In order to assess grid convergence, additional simulations are performed using the same time step t = 0 . 001 s for three additional different grid sizes, namely 100 × 50 × 5, 200 × 100 × 10, and 800 × 400 × 40, respectively. Fig. 4 (a) presents water surface profiles in the centre of the channel at t = 4 s as computed on four different grids. It is shown that apart from the coarse mesh result, there are only minor differences for the other three simulation results and the profiles collapse onto each other for the finer grids. In order to quantify the difference between simulation results, the numerical result from simulation on the finest grid are used as a reference solution (as there is no analytical solution), and the root mean square (RMS) errors are calculated and shown in Table 1 .
In order to assess the sensitivity of the results due to the time step, simulations are carried out on the grid mentioned above (400 × 200 × 20), but employing three different time steps for both first-order and second-order time integration schemes,   Table 2 .
In order to quantify the viscous damping of the wave, three additional waves with different A 0 / h are simulated and the results are compared with the analytical solutions proposed by [77] as where A max is the amplitude of the solitary wave and A 0max is the initial elevation. Fig. 5 plots the simulated waves together with the analytical solution demonstrating that the simulation results agree convincingly well with the analytical solutions for smaller wave and this is because the analytical solution obtained from the perturbation theory is valid only for A 0max / h ≤ 0.1 and it is similar to the results obtained in [76] .  The results in the smallest time step are used as reference. Another quantity to examine is the wave run-up height at the right wall. Fig. 6 shows the predicted run-up as a function of wave height for nine different initial elevations together with numerical simulations from other researchers. Good agreement is obtained demonstrating the capability of the present approach to simulate wave propagation and interaction with solid boundaries without introducing marked numerical dissipation.
Finally, mass conservation of the simulations is calculated and it is found that the error of the mass is less than 0.001% for A 0 /h = 0 . 1 , 0.0015% for A 0 /h = 0 . 2 , and 0.002% for A 0 /h = 0 . 4 at the end of each simulation.

3D solitary wave interacting with a vertical circular cylinder
A solitary wave interacting with a vertical circular cylinder is considered next. This case has been studied previously both experimentally [78] and numerically [18,22] . In the physical experiment, the wave tank was 7.  The simulation is set up to replicate the laboratory measurement and the schematic view of this problem is shown in Fig. 7 . In the experiment detailed water surface elevations were measured in the vicinity of the cylinder as well as horizontal wave loadings on the cylinder. The locations of the measuring points in terms of polar coordinates are shown on the top in Fig. 7 . in which four gauging points are along the radial direction with θ measured in the clockwise direction and θ = 0 o is the upstream wave direction.
The computational domain L x × L y × L z is 32 R × 2 h × 12 R with the cylinder in the centre and it is discretised by a uniform mesh 1280 × 64 × 480 in the streamwise, vertical, and spanwise directions, respectively. In the simulation, the solitary wave is generated from the inlet based on analytical solutions and the case with a wave height of H/h = 0 . 4 is considered here.
Figs. 8 , 9 , 10 plot the time history of the water surface elevation as predicted numerically together with experimental data of [78] , where the water surface elevation η is normalised by the wave height H and the time t is normalised by h/g . All numerical results are shifted at an instant in the time axis so that the predicted peak wave coincides with the experimental measurement away from the cylinder at r/R = 2 . 92 and θ = 100 o . Reasonably good agreement is obtained between the LES results and the  experimental measurements. However, the wave heights at gauge points r/R = 1 . 03 , i.e. close to the cylinder, are over-predicted in the present simulation. However, we found that this is also observed in other numerical simulations [18,22,78] . Compared to previous simulations, the predicted wave heights appear closer to the experimental measurements than previous studies due to the high resolution of the simulation and an improved representation of viscous effects and hence improved energy dissipation.  Fig. 11 plots the time history of the dimensionless horizontal wave force on the cylinder as predicted by the simulation and as measured in the experiment of [78] . Overall good agreement is obtained, especially in terms of the magnitude and phase until t(g/h ) 1 / 2 = 15 . After that, the force is overestimated due to a higher wave height predicted in the simulation. Fig. 12 presents snapshots of the water surface profiles together with turbulent vortical structures, here identified via isosurfaces of λ 2 , i.e. the second invariant of the velocity gradient tensor [79] , to identify vortex cores. The solitary wave is nearly two-dimensional before approaching the cylinder and most of the vortical struc- tures are in the vicinity of the wave crest and the cylinder. During the solitary wave's impact on the cylinder, significant wave-runup is observed and more turbulent vortical structures are generated around the cylinder. When the wave passes the cylinder, it is split into two parts wrapping around the cylinder and a reflected wave is generated that propagates radially outwards, and this can be noticed from the strong vortical structures generated upstream. After that, oscillatory waves are generated followed by a negative wave in the upstream. Downstream of the cylinder, the two parts of the wave meet the main wave passing the cylinder and develop into a three-dimensional wave that propagates towards the outlet. Turbulent structures are significantly decreased in the vicinity of the cylinder while they are enhanced around the solitary wave.

3D solitary wave propagating over a thin horizontal plate
A solitary wave traveling over a thin horizontal plate is investigated next, a case that has been studied previously by both ex-perimentally [80] and numerically [22,80] and to the best of our knowledge these previous numerical studies [22,80] have all employed two-dimensional domains.
The computation is set up to replicate the laboratory measurement and the schematic view of this problem is shown in Fig. 13 . A thin plate of length L = 1 . 156 m, width L z = 0 . 6 m and thickness δ = 0 . 01 m is submerged in a wave tank. The submergence is d and the water depth is h = 0 . 2 m. The corresponding Reynolds number based on the plate thickness and wave speed is approximately 1.57 × 10 4 . In the experiment four wave gauges (WG1-WG4) and ten pressure gauges (P1-P10) along the plate were placed to measure the wave elevation and pressure on the plate, respectively, and their locations are marked in Fig. 13 . The computational domain L x × L y × L z is 11 (m) × 0.32 (m) × 0.6 (m) with the leading edge of the plate located at x = 4 m, and it is discretised by a uniform mesh 1600 × 160 × 96 in the streamwise, vertical, and spanwise directions, respectively. In the simulation, the solitary wave is generated from the inlet based on analytical solutions and the case Fig. 14 plots water surface elevations at four gauging points as computed together with experimental data, and the 2D numerical results of [80] , in which the water surface elevation is normalised by the wave height H and the time is normalised by the wave period T 0 and t = 0 is when the peak wave passes the gauge point WG2. It can be seen from WG1 that a correct incident wave is generated upstream and that the reflecting wave is well captured by the present simulation. There is a wave transformation above the plate during wave propagation and the oscillatory wave is reduced at WG4 when the wave has passed over the plate. Overall, very good agreement is obtained between the present simulations and the experimental data in terms of wave amplitude and phase. The LES results are similar to the numerical results obtained by the 2D COBRAS model [80] . However it is worth noting that when compared to the experimental data, the maximum wave elevation is slightly better predicted by the LES than the 2D simulation without turbulence model in [80] , demonstrating that the LES model is able to capture the energy dissipation during wave-structure interaction.
The dynamic characteristics during wave impact is quantified in Fig. 15 presenting normalised pressure along the centreline of the plate as computed here, as experimentally measured, as well as predicted by the 2D numerical model COBRAS of [80] . The LES results are similar to the COBRAS results, and both of them agree well with the experimental data. Better agreement is obtained for the pressure above the plate than the points below the plate. Although the magnitude of the pressure is overpredicted by the present simulation towards the lower end of the plate, the phase characteristics are well captured similar to the COBRAS model predictions.
The net vertical force acting on the plate and the moment about the plate's centre is plotted in Fig. 16 . The predicted vertical force agrees with the experimental data, again similar to the predictions of the COBRAS model. However, there is some discrepancy between simulation and experiment in terms of the moment, but similar to the result obtained from the COBRAS model. Lo and Liu [80] mentioned that this might be due to the different spatial resolution used in experiment and simulation, i.e. there are fewer values available to compute the moment in the experiment than in the simulation. Fig. 17 shows snapshots for the water surface profile before, during and after the solitary wave impacting on the horizontal plate at four different time. It can be seen that the plate has a little effect on the solitary wave and there is no significant change of the wave form when passing the structure. During wave-structure interaction, the vortical structures identified by the λ 2 method [79] in the flow are also shown. Vortical structures are normally generated around the two edges of the plate and in the vicinity of the wave crest and trough. Large vortices are generated by the wave approaching the plate and being advected when the wave passes the plate.

3D focusing wave impact on a FPSO-like structure
A 3D focusing wave, a form of an extreme wave, impacting on a complex FPSO-like structure is considered. This case was a blind test during the Collaborative Computational Project in Wave Structure Interaction (CCP-WSI), in which ten different numerical models were employed and their results were compared with experimental data. The computational model is set up to replicate the wave and geometric properties of the laboratory experiment [30] , in which detailed free surface profiles and pressure on the structure were measured. The schematic view of the computational domain is shown in Fig. 18 , where the origin is located at the still water level of the inlet along the central plane. A smaller section was selected from the test Fig. 15. Normalised pressure along the centreline of the plate as computed, experimentally measured, and as predicted by the 2D numerical model of [80] . Fig. 16. Normalised force and moment as computed, experimentally measured, and as predicted by the 2D numerical model of [80] .   section in the experiment in order to avoid wall effects and also to save computational effort. The computational domain L x × L y × L z is 12 (m) × 3.3 (m) × 3 (m) with the FPSO-like structure (1.2 m long, 0.3 m wide, and 0.3 m deep) located at x = 2 . 37 m away from the inlet, and it is discretised by a uniform mesh 320 × 320 × 160 in the streamwise, vertical, and spanwise directions, respectively. The corresponding Reynolds number based on the structure height and wave speed is approximately 1.8 × 10 6 .
In order to generate the focusing wave, linear wave theory is used here, where both the time history of wave elevations and velocities for each component (244 components in total) is specified for the water corresponding to the experiment, while the velocities are set to zero for the air. In order to check whether the correct focusing wave is generated at the inlet, the time history of the wave elevation is compared with the experimental measurements without the structure at wave gauge WG1 in Fig. 19 . It can be seen that a reasonable agreement is obtained from the linear wave the-ory and better wave generation can only be obtained by using the self-correction wavemaker technique [61] . The present LES model has been validated against the experimental measurements as well as other numerical models in terms of water surface elevation and pressure at various locations around the FPSO as well as the efficiency in [30,61] , here the focus is on the kinematics and dynamics of the wave-structure interaction and the resulting turbulence field. Fig. 20 shows snapshots of the water surface profile and turbulent vortical structures identified by the λ 2 method [79] . It can be seen that some vortical structures are generated around the structures at the beginning when small waves pass the FPSO. The wave amplitude gradually increases until the focusing wave impacts on the structure. During impact, more turbulent vortical structures are generated in the vicinity of the crest and trough of the focusing wave. When the focusing wave passes the structure, there is a small reflected wave generated upstream the bow, which can be also observed from the turbulent structures around the FPSO.

Conclusions
In this paper, a two-phase flow approach has been introduced to be employed for investigating 3D wave-structure interaction problems. The large-eddy simulation approach is adopted, in which the space-filtered Navier-Stokes equations are solved on Cartesian grids. The dynamic Smagorinsky subgrid-scale model is employed to compute sub-grid scale stresses and to account for the effects of the unresolved small scales of turbulence on the large scales. The finite volume method is utilised to discretise spatial deriva-tives and the PISO algorithm for the pressure-velocity coupling. An implicit backward difference discretisation is used for the time derivative. The air-water interface is being captured using the high resolution VOF scheme CICSAM, and the Cartesian cut-cell method is implemented to deal with complex geometries immersed in the Cartesian fluid grid.
With the goal to thoroughly validate the in-house code Xdol-phin3D, several complex multi-dimensional flows have been studied. A travelling solitary wave in a closed channel is considered first, and run-up height and viscous damping characteristics have been compared with experimental data and an analytical solution. Good agreement has been found, spatial and temporal convergence has been demonstrated and energy dissipation and mass conservation has been assessed. The capability of the present LES-based two-phase flow approach to predict accurately and reliably wavestructure interaction problems for geometrically-complex structures (such as a circular cylinder, a horizontal thin plate, and an FPSO-like structure) has been demonstrated. Numerically predicted water surface elevations, pressure at selected locations on the structure, and acting wave forces and moments on the structure have been compared with available experimental data or the results of other numerical simulations, respectively, an overall convincing agreement has been obtained. Furthermore, snapshots of wave surfaces and turbulent vortical structures, as identified by the λ 2 method, have been presented and discussed, visualising the complex turbulent flow fields around the structures during WSI.
This study demonstrates the capability of the present multiphase flow model to predict 3D wave-structure interaction. The model can act as a complementary approach to experimental investigations to gain further insight into the kinematics and dynamics of three-dimensional wave-structure interaction problems.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.