A Coupled Artiﬁcial Compressibility Method for Free Surface Flows

: Modeling free surface ﬂows in a CFD context typically requires an incompressible approach along with a formulation to account for the air–water interface. Commonly, pressure-correction algorithms combined with the Volume of Fluid (VOF) method are used to describe these kinds of ﬂows. Pressure-correction algorithms are segregated solvers, which means equations are solved in sequence until convergence is accomplished. On the contrary, the artiﬁcial compressibility (AC) method solves a single coupled system of equations. Solving at each timestep a single system of equations obviates the need for segregated algorithms, since all equations converge simultaneously. The goal of the present work is to combine the AC method with VOF formulation and prove its ability to account for unsteady ﬂows of immiscible ﬂuids. The presented system of equations has a hyperbolic nature in pseudo-time, thus the arsenal of the hyperbolic discretization process can be exploited. To this end, a thorough investigation of unsteady ﬂows is presented to demonstrate the ability of the method to accurately describe unsteady ﬂows. Problems of wave propagation on constant and variable bathymetry are considered, as well as a ﬂuid structure interaction problem, where viscous effects have a signiﬁcant impact on the motion of the structure. In all cases the results obtained are compared with theoretical or experimental data. The straightforward implementation of the method, as well as its accurate predictions, shows that AC method can be regarded as a suitable choice to account for free surface ﬂows.


Introduction
Free surface flows pose significant challenges when it comes to numerical modeling. Accurate modeling of the free surface is crucial in many engineering applications ranging from ship hydrodynamics to floating platforms. Additionally, for the proper hydrodynamic and structural design of marine vessels and offshore structures, it is mandatory to take into account wave-structure interaction. To that end, several numerical methods have been developed of varying fidelity. Boundary element methods (BEM) have been widely used in the literature to model non-linear free surface waves [1]. They can provide accurate results at a small computational cost [2] and can be easily coupled with other numerical frameworks (e.g., hydro-aero-elastic solvers [3]). More sophisticated models have been also developed [4] reducing the problem size and consequently the computational cost. However, potential methods do not account for viscosity and the usual remedy is to apply viscous correction on the results [5].
When viscous effects become significant, the usual approach is to employ a computational fluid dynamic (CFD) framework. Space is discretized using structured/unstructured grids and the Reynolds Averaged Navies-Stokes (RANS) equations are solved on the underlying mesh. Two main categories of methodologies have been developed to account for free surface flows. On the one hand, the surface The paper is structured as follows. Sections 2 and 3 present Mathematical and Numerical modeling, while Section 4 includes the analysis of various method parameters. In addition, MaPFlow's predictions are validated against available experimental data. The case of wave propagation over a variable bathymetry is considered and compared to measurements available from [29]. Also, a fluid-structure interaction problem of a moonpool-type floater is considered and predictions are compared with available experimental data from [30]. Finally, in Section 5 the basic conclusions are outlined.

Governing Equations
The VOF method uses an indicator function to describe the presence of the liquid or the gas phase. Specifically, let ρ w be the water density and ρ α the density of the air, then the volume fraction is defined as α l = (ρ m − ρ a )/ (ρ w − ρ a ). The volume fraction is equal to 0 in regions occupied by the air phase, equals to 1 when only the water phase is present and takes values between 0 and 1 near the free surface. The interface of the fluids is located in regions of rapid change of the gradient of the indicator function. The properties of the mixture fluid, density (ρ m ) and dynamic viscosity (µ m ), are described by the blending functions (1).
The free surface, as presented in Equation (2), is considered to be a material surface.
where υ = (u, v, w) is the three-dimensional velocity vector, ∇ = ∂ x , ∂ y , ∂ z is the divergence operator and t denotes the real time variable. The artificial compressibility method attempts to overcome the decoupling of the mass and momentum equations by augmenting the mass equation with a pseudo-time derivative of pressure (p), as it is shown in Equation (3). The influence of pseudo-time derivative is regulated by the artificial compressibility factor β. 1 β ∂p ∂τ + ∇ · υ = 0 (3) τ denotes the pseudo-time variable. The pseudo-time derivative corresponds to the iterative numerical procedure until convergence. In order to obtain the original equation, the pseudo-time derivative of the convergent solution should be equals to zero.
The essence of the AC method is that it assumes a relation between pressure and density during convergence (Equation (4)). This relation is similar to compressible definition of sound speed. However, in this case the parameter β is a numerical parameter that regulates convergence.
Equation (4) gives rise to a pseudo-sound speed definition. In case of one-dimensional, one-phase flows, the artificial speed of sound is given by Equation (5). The sound speed is regulated through the artificial parameter β. In Appendix A it is shown that the pseudo-sound speed regulates the convergence in a similar way that the sound speed effects the compressible equations.
Regarding the simulation of two-phase flows, the momentum equations written in a conservative form and expressed in terms of the mixture quantities are presented in Equation (6), along with the VOF equation expressed also, in conservative form.
In previous equations,σ is the stress tensor and the vector F B includes source terms and body forces.
The Equations (3) and (6) constitute a fully coupled system of equations capable of describing two-phase flows. The formed system has a hyperbolic nature and consequently numerical techniques used for such solvers can be used. By introducing the AC method the system of equation becomes hyperbolic in pseudo-time and consequently numerical techniques used for such solvers can be used. The convergent solution should satisfy the original sets of equations, by eliminating the pseudo-time derivatives. Nevertheless, the system in its original form poses several difficulties. Firstly, the density appears in the eigenvalues of the system yielding it stiff for high density ratios and secondly, the system cannot be written in a conservative form. To mitigate these perplexities, the preconditioner of Kunz [24] is employed.
By introducing, the fictitious time derivative, for the momentum equations, the artificial compressibility method can be used in time-accurate computations. Indeed, employing the dual-time stepping technique each unsteady timestep is treated as a steady state problem. Furthermore, to express the governing equations as a single coupled system of equations the time derivatives (real and fictitious) of the momentum are expressed as a sum of time derivatives for velocity and density. The governing equation can be written in the following integrated form.
The system of Equation (7) is a fully coupled system of equations. These equations express the governing equations with respect to the primitive variables Q. In order to cast the system in conservative form, the transformation matrix Γ e is used. In Equation (8) the conservative variables are given by the vector U. The three-dimension vector of velocity is denoted with υ, while p is the pressure.
The Jacobian matrix Γ e and the precondition matrix of Kunz Γ are given by (9), where ∆ρ is the density difference between the heavier and the lighter fluid and I 3×3 is the 3 by 3 identity matrix.
Finally, the inviscid and viscous fluxes can be summarized as: τ xx n x + τ xy n y + τ xz n z τ yx n x + τ yy n y + τ yz n z τ zx n x + τ zy n y + τ zz n z 0 where V n = υ · n, V g = υ vol · n, ∆V = V n − V g , while υ vol is the velocity of the control volume and n is the surface normal of the control volume. The viscous stresses τ ij are computed as where µ t is the turbulent dynamic viscosity, k is the turbulent kinetic energy and δ ij is the Kronecker delta.

Spatial Discretization
The details of the spatial discretization and the evaluation of the fluxes are given in Appendix A. In this section, the methods used for the reconstruction of the flowfield are described.
The evaluation of the inviscid fluxes in Equation (10) requires an approximation of the left and right state of a face ( Figure 1). These states are computed based on a reconstruction scheme that extrapolates the cell-centered value of the volumes in the respective face. The right and left state are designated based on the normal vector of the face, which points from the left state to the right. Due to the particular nature of the governing equations, a different reconstruction scheme is adopted for each equation. The velocity field is approximated through a piecewise linear interpolation scheme, given by Equation (12). Since the surface tension is neglected, the velocity field is continuous even across the free surface. For this reason, the gradients are retained, and no limiter is applied.
The vectors r i , r j are pointing from the cell center to the midpoint of the face, as illustrated in Figure 1.
Furthermore, the pressure field is a continuous function in space, since surface tension is neglected. However, the pressure gradient is discontinuous across the free surface, due to the density jump.
The jump condition that must hold requires ∇p ρ = 0. Many researchers [17,31] have introduced different schemes to treat this difficulty by adopting a density-based interpolation scheme. In MaPFlow, the work of Queutey et al. [17] is followed. This scheme is introduced only near the free surface, while in the rest of the computational domain a piecewise linear interpolation scheme is used, similar to Equation (12).
Finally, of great importance is the reconstruction of volume fraction field (α l ). In order to reduce the numerical diffusion, it is important to adopt a compressive reconstruction scheme. Over recent decades, numerous reconstruction schemes have been introduced which offer low numerical diffusion. The requirements should meet is boundedness and high accuracy even in large CFL numbers. Most of these reconstruction schemes are based on Leonard's Normalized Variable Diagram (NVD) [32] (e.g., CISCAM [33], HRIC [34], BICS [35]). In the present work the STACS [36] scheme is adopted as the free surface capturing scheme.

Temporal Discretization
The present section focuses on the time discretization process. Let Q * be the vector of unknowns. An implicit formulation of the problem can be defined as where R * is the unsteady residual (14), which includes the spatial terms of Equation (7), here denoted as R D i Q * , and the unsteady term.
If n is the index of the known time solution and k is the internal iterator of the steady problem, then at the beginning of every steady iteration k = 0, the vector of unknowns is initialized as Q * = Q n and until convergence they do not satisfy the original unsteady problem. The convergent solution of the problem is obtained when R * → 0 and hence Q * → Q n+1 .
The unsteady term is discretized as a series expansion of successive levels backwards in time (BDF schemes) [37].
When the control volume changes in time, the Geometric Conservation Law (GCL) [38] should be satisfied (Equation (16) Using a similar backwards differentiation in time as in Equation (15) the GCL can be written as where the residual of the GCL is defined as To ensure that the GCL is satisfied, Equation (17) is applied directly to the discretization of the unsteady term and thus Equation (15) becomes In MaPFlow two successive levels of solution are retained, yielding a second order accurate scheme in time. The fictitious time derivative of equation is discretized using a first-order backward difference scheme To facilitate convergence the local time stepping technique is used. The local pseudo-timestep is determined by whereΛ c,i is the convective spectral radii and it is defined bŷ

Turbulence Modeling
Regarding turbulence modeling, the k − ω SST model of Menter [39] is used in this work. A well-known downside of RANS models when employed for free surface flows is the exponential growth of the turbulent kinetic energy and eddy viscosity [40,41]. The overproduction of turbulence viscosity, especially around the interface between water and air, is triggered by the shear layer in the free surface and results to an artificial damping of the propagated waves. Devolder et al. [40] proposed a buoyancy term (23) to the turbulent kinetic energy equation.
Although this technique improves drastically the results, when free surface waves are considered in more complex flows (e.g., with surface piercing geometries) the problem of the overproduction persists [41]. For this reason, the Kato-Launder production limiter [42] is also used.

Wave Generation and Absorption
When considering a numerical wave tank, the generation of the desired wave profile as well as the effective radiation of the outgoing waves outside the computational domain are two non-trivial tasks. The generation of steadily progressive waves is performed by forcing the numerical solution, in a specific part of the computational domain, to follow a wave solution provided by a wave theory. Furthermore, an artificial damping of the waves is required near the boundaries of the domain. The boundary conditions, assume a uniform field and thus any physical disturbance created inside the domain should not reach the farfield boundary. Various techniques have been considered for practical implementations. Among them, is the coarsening of the computational mesh near the boundaries. Although this technique is able to smear the solution, the mesh coarsening is case-dependent, and it does not guarantee zero reflection. In most cases, the boundaries are supplied with damping zones where the farfield conditions are imposed.
In this work, the numerical generation and absorption of the free surface waves is performed through source terms applied in the momentum equations, in specific zones of the computational domain near the farfield boundaries. Typically, these zones extend for a few wavelengths. The general form of the source terms is given by Equation (24). This source term drives the solution to the imposed υ tar velocity vector.
The effect of the source terms is regulated through the C nwt function, its form is given by Equation (25), while in Figure 2 it is plotted for different values of α and n. The maximum value of the function is regulated through parameter α, while its spatial distribution through n parameter. The function is zero away from the boundaries of the computational domain, scales exponentially inside the specified zone and it reaches its maximum value α at the boundary. This function was originally proposed in [43] where the desired numerical solution was imposed explicitly. In the present work, the desired numerical solution (e.g., a specific wave profile in the generation zone) is obtained implicitly, meaning that the solver needs to converge to that solution through the numerical procedure. The source term (25), should be able to drive the solution to the desired one, while at the same time maintain the convergence properties of the method. Large values of C nwt may lead to numerical instabilities, while small values of α may not be able to drive the solution. Typical values for the exponent n is between 2 and 5 and usually α multiplier is not greater than 200. In the first section of the numerical results, it will be shown that there is an acceptable range of values for the coefficient C nwt where no numerical instabilities appear.
The source term of Equation (25) is a function of the non-dimensional space variable x r , which depends on the starting position x s of the specified zone and the end point of the zone x e .
The wave generation is performed by imposing the velocity υ tar as given by an appropriate wave theory. In this work, the wave solution is obtained from the semi-analytical method of Fenton's Stream Function Theory [44].
In case of the wave damping zones, the target is to minimize the transverse velocity components. In a two-dimensional context, where the wave propagation is along the x-axis, the velocity vector υ tar takes the form υ damp = 0, v in f . This method has been assessed extensively as a wave damping technique [45]. In the present work, wave generation is accomplished in a similar manner.

Volume Fraction Boundedness
The volume fraction equation is an advection equation. The values of α l should always remain inside the definition bounds [0 : 1]. However, due to numerical instabilities its value may surpass its domain of definition. This can be triggered by two factors. First, for the derivation of the conservative form of the VOF equation a divergence free field was assumed. This may not hold when distorted meshes are used, near the sharp edges of the wall boundary or when dominant sources are applied (such as in damping and generation zones). Secondly, another reason for violating the domain of definition is the spurious oscillations caused by the reconstruction scheme of the VOF field. High order reconstruction scheme may introduce unphysical oscillations near sharp discontinuities [46]. This may be treated by employing a first-order approximation which would ensure boundedness, but at the cost of excessive diffusion. Much work has been done to remedy this behavior, especially in the frameworks of the Total Variation Diminishing (TVD) schemes [47] and of the Normalized Variable Formulation [32]. Unfortunately, in practical implementations the requirements for boundedness of the reconstruction scheme may not be satisfied. Finally, it should be noticed that both arguments are correlated with the timestep resolution and thus small Courant number are usually adopted.
In MaPFlow a relaxed limit (of order 1e −3 ) to the volume fraction is applied. Through numerical investigation it was seen that the VOF clipping was large during the initial stages of simulations; however, when a steady or periodic state was reached the clipping was reduced significantly. In all cases, the introduced mass error was negligible.

Numerical Results
In this section, numerical results based on the artificial compressibility method are presented. First, the ability of the solver to accurately generate and propagate waves on constant bathymetry is demonstrated. Following a grid and timestep sensitivity study, a parametric study regarding the influence of the artificial compressibility parameter is presented. Afterwards, the interaction of waves with variable bathymetry is considered and the results are compared with experimental data. Finally, a two-dimensional fluid-structure interaction problem is presented. In the latter case, the artificial compressibility methodology is coupled with a dynamic solver and the interaction of moonpool-type floater with incident waves is examined. The amplitudes of the motions are compared with experimental data.

Numerical Wave Tank
The propagation of a cnoidal wave at constant bathymetry is considered. The wave period is equal to T = 5 s, the wave height is H = 0.05 m, while the depth of the numerical wave tank is d = 0.5 m. Using the stream function theory of Fenton, the wavelength is found equals to l = 11.08 m.

Wave Generation and Absorption
First, the influence of the generation parameters α, n, as defined in Equation (25) is investigated. In Table 1 six test cases are defined. The same grid configuration is used in all cases. In the direction of the wave propagation, the mesh is uniform almost everywhere with 150 cells per wave period, except in the damping zone where the mesh coarsens based on a geometric rule. In the vertical direction, 20 cells per wave height are employed. Regarding the timestep discretization, a constant timestep is chosen corresponding to 800 timesteps per wave period (dt = 0.0625 s). The surface elevation, for all cases, inside the generation zone are shown in Figure 3a and is compared with the analytical solution provided by the Stream Function Theory. It can be seen that the waveform is distorted for large values of α. Furthermore, discrepancies are noted near the farfield boundary. This occurs due to a lagging between the boundary conditions and the time window that solver need to converge to the imposed wave solution. However, in all cases the solution converges to the desired wave profile as depicted in Figure 3b, where the free surface elevation is plotted against time and it is compared with its analytical form.  Tables 2 and 3 present the error of the amplitude of the first four harmonics and the error of the phase angle in degrees between the analytical and numerical solution. Overall, the error of the amplitude at the dominant frequency is below 0.5% for all cases, while for large values of α an amplification in higher frequencies is evident. The exponent n does not seem to have any significant influence on the results. The general recommendation is that the source terms should have small influence in the system of equations, in order to enhance the stability of the solver. For this reason, large values of exponent n and small values of the parameter α are desired. On these grounds, an accurate representation of the wave profile is obtained, while keeping the influence of the source terms small, by choosing α = 60 and n = 3.5. Next the damping performance of the solver is evaluated. A parametric study is conducted by varying the values of the parameters α, n and the length of the damping zone. The case studies are defined in Table 4. The asterisk, in the last two cases, indicates that the damping zone is coarser. In Figure 4a the surface elevation inside the damping zone is shown for the six test cases, while in Figure 4b the free surface elevation during a wave period at the beginning of the damping zone is presented. As illustrated in Figure 4a, the elevation of the free surface near the boundary is completely damped in case 4, contrary to cases 1 and 3. It is also worth noticing that for the same parameters α, n (cases 2 and 5) the grid coarsening does not affect the damping of the wave. However, all cases produce the same wave profile at the beginning of the damping zone.  Furthermore, in order to examine in detail the effect of the damping zone in the wave flume, the amplitude of the 1st harmonic of the free surface elevation is plotted in Figure 4c, at various stations across the computational domain. As expected, in all cases the amplitude decays due to numerical diffusion. However, in case 3 an oscillatory behavior of the harmonics amplitudes is noticed. This is caused by reflections at the boundary of the domain. Also, in case 1 the amplitude of the 1st harmonic is not monotonically decreasing in space. This can be attributed to reflections at the farfield boundary, as in the previous case. In all other cases, the differences are regarded negligible. A damping zone extending six wavelengths is excessive for practical implementations while large values of α may cause numerical instabilities. For this reason, α = 120, n = 3.5 is chosen (case 2) for the rest of the work.
Most of the results showed that the generation and damping technique can produce accurate results. The differences noted for different values of α and n can be considered small. In case of wave generation, the numerical procedure was able to provide acceptable results in all cases however for large values of α (>300) numerical instabilities were noticed. Regarding wave damping, the method can absorb the outgoing waves, provided that the length of the zones is sufficient. Finally, the damping sensitivity study indicates that the damping performance is more sensitive to the chosen value of α. For this reason, a larger value of α, compared to wave generation, is chosen.

Grid and Timestep Independence
Once α and n are decided the next step is to perform a grid and a timestep sensitivity study. The grid parameterization is based on the number of cells per wavelength (λ/dx) and the number of cells per wave height near the free surface (H/dy FS ). Predictions from four grids are compared in Figure 5 with respect to free surface elevation in Figure 5a and the elevation at a specific station during a wave period in Figure 5b.
In Figure 5c the amplitude of the first harmonic of the propagating wave is depicted. The figure shows that in case of the fine grid (blue line) the amplitude of the wave elevation remains almost constant across the domain. The coarse grid (red line) produces significant numerical diffusion and the amplitude of the waves is decreasing. The yellow line corresponds to a mesh with highly skewed cells near the free surface. It should be noticed that although the green line corresponds to a coarser grid, the numerical diffusion introduced is smaller. Due to very thin cells near the free surface, the CFL number is large and consequently the numerical diffusion introduced by the temporal discretization is significant [48]. The computational mesh with 150 cells per wavelength and 20 cells in the wave height (green line) exhibits similar characteristics with the finer grid (blue line) and consequently, it is kept for the rest of the simulations. In Figure 5d the amplitude of the first harmonic of the cnoidal wave can be found. Three different timesteps are considered, based on the wave period. Predictions for 400, 800 and 1600 timesteps per period are shown. It is evident that that even when 400 timesteps per period are chosen the relative error after 8 wavelengths is ≈ 2.5%. Predictions with 800 and 1600 timesteps per wave period are in very close agreement and thus the former is selected in the present work.

Influence of the Artificial Compressibility β Parameter
The next case is to study the impact of the artificial compressibility parameter β on wave propagation. In general, β regulates the coupling between pressure and density during convergence, i.e., β regulates the pseudo-sound speed (c). On the one hand, small values of β lead to small values of c thus the system is strongly coupled, on the other hand large values of β lead to a loose coupling of the equations. Because the parameter β affects the convergence of the method, it is imperative to quantify how the choice of β affects the predictions.
Several guidelines for choosing β can be found in the literature. For one-phase incompressible flows, Turkel [22] proposed that β should scale with the square of the local velocity, and hence β varies across the domain. Special care should be taken near stagnation points where its value should be safeguarded. In case of wave propagation, the velocities under the crest are small and so the β parameter could not depend on the local characteristics of the flow. Nevertheless, the flows considered in his work did not involve gravity. Kunz et al. [24] applied the AC method for cavitation problems and suggested values of β ranging from 5 − 15 · υ 2 ∞ . In order to understand the effect of β on the predictions, we conduct a parametric study, using five different constant values of β = 1, 10, 50, 100 and 1000.
The numerical setup is based on the previous remarks. The results are summarized in Figure 6. Figure 6a shows the free surface elevation 5 wavelengths downstream of the generation zone. It is evident that the predictions yield large differences for the two extreme values of β (β = 1 and β = 1000). However, for values ranging from 10-100 the predictions yield negligible differences. This is also exhibited in Figure 6b where the amplitude of the first harmonic for the various β values is shown. For values β = 1 and β = 1000 elongation of the wavelength is observed. This is an artifact that implies undesired compressibility effects. To further investigate the effect of β on the solution, the 2nd and 3rd harmonics are provided in Figure 6. It is apparent that within the range β = 10 − 100 the solution is insensitive to the chosen value for both the 2nd (Figure 6c) and 3rd (Figure 6d) harmonic. Nonetheless, for β = 1 the differences are larger but still within an acceptable limit. For β = 1000 the solutions is clearly affected.
From the preceding analysis, it is evident that for the present formulation predictions are insensitive of value of β within an acceptable range. With these considerations and taking into account the guidelines of the literature, a value of β = 10 is chosen for the rest of this work.

Wave Interaction with Variable Bathymetry
We shall now consider the interaction of a wave with a trapezoidal bathymetry. This is a standard benchmark to evaluate the solver's capability to accurately simulate dispersive phenomena. The original experiment was setup and conducted by Beji and Battjes [49], where seven wave gauges were placed along the wave tank. Later, Dingemans [29] repeated the same experiment using four more probes. The location of the wave probes are depicted in Figure 7. The generated wave profile has height h = 2.0 cm and period T = 2.02 s.  The timestep resolution was based on the incident wave period, and the parametric study included the values dt = T/400, T/800 and T/1600. The spatial mesh should have enough resolution to account for the smallest ripples of the irregular wave caused by the interaction with the bathymetry. The sensitivity study concluded for the timestep resolution to dt = 2.5 e −3 s which corresponds to ≈800 steps per period and to a grid of total size of 3.2 · 10 3 cells. The grid in the x direction is uniform with dx = 5 mm, in the y direction the mesh is dense near the free surface (dy = 1 mm) and coarsens gradually towards the bottom.
Regarding the numerical setup, the artificial compressibility factor is equal to β = 10. The pseudo-CFL number is equal to 50 and at each time true iteration 10 dual steps were executed for the convergence of the pseudo-steady state problem. At the bottom, no slip boundary condition is applied.
The predicted free surface elevation is also presented in Figure 7 (continuous red line). Due to shoaling, the wavelength is compressed, and energy is transferred to higher harmonics. The deepening after the bar releases the energy from the higher harmonics and an irregular wave profile is obtained.
The signals recorded at the wave probes are compared against the experimental data in Figure 8, at eleven stations.
At the first probe (Figure 8a), the two signal appear to have a small difference in amplitude. However, as the comparison at the other probes suggests it does not affect the overall comparison. At probe 7 (Figure 8g) a small phase shift is observed between the simulation and the measurements. Nevertheless, the comparison of free surface displacement is very good at all stations even when higher harmonic are excited after the shoaling (Figure 8f-k).

Moonpool-Type Floater
In this section, the response of a moonpool-type floater in incident wave is investigated and compared with experimental data from [30]. Under wave excitation the floater is displaced and the free surface inside the moonpool is subjected to a vertical motion, commonly called piston motion, while additional sloshing modes may appear. The piston motion inside the moonpool reduces the motion of the structure induced by the incident waves. This problem poses additional challenges since the motion of the floater is induced by the underlying flowfield. It is noted that viscous effect are significant in this case since the flow separation at the sharp edges directly affects the free surface elevation inside the moonpool. Employing potential methods in this case results in overshoots of the piston mode and thus viscous corrections must be included [50,51].
The motion of the floater is described by Newton's second law. At each unknown timestep n + 1, Equation (26) should be satisfied.
The vector x contains the unknown displacements and rotations and the double dots imply twice differentiation with respect to time. The matrix M is the mass matrix and vector F tot includes the sum of the forces applied to the body (e.g., hydrodynamic forces, external forces, etc.). The above system of equations constitutes a non-linear system of equations. A linearization is imposed based on the assumption of small perturbations in each pseudo-timestep. As a result a linearized δ-formulation is obtained (Equation (27)).
In the above equation F k tot is the force vector calculated at the known timestep k and for the correction This linearization implies an iterative process during pseudo-timesteps. At each k-iteration the correction components δ x are computed based on the hydrodynamics forces provided by the solver. The iterative procedure ends when both solvers converge, implying a strong coupling between them.
The second order ordinary differential equations (ODE) (27) are solved with the Newmark-β method with coefficient values β = 0.25 and γ = 0.5. This method uses a second order Taylor expansion variant, in which a predictor step is intervened based on the coefficients β and γ.
Since the geometry, is not stationary an appropriate grid technique needs to be adopted to account for the moving boundaries. In cases of small amplitude motion, the deforming grid technique could be used. In MaPFlow the nodes of the mesh are moved based on a damping term as described in [52].
In the experiments conducted by [30] two barges of rectangular cross section of 20 cm were mounted together to form a moonpool of 10 cm. The breadth of the barges was 14 mm smaller than the width of the experimental wave tank and thus 3D effects are limited.
Following the notation of the experiments, the test case is defined on the y-z plane. In a two-dimensional context, a body has three degrees of freedom, namely the sway motion (motion along the y-axis, η 2 ), the heave motion (motion along the z-axis, η 3 ) and the roll motion (motion around the x-axis, η 4 ).
The setup of the experiment is depicted in Figure 9. The structure is mounted with horizontal mooring lines to prevent it from drifting. At the end of these lines springs were connected to allow the structure to move. Pre-tension was applied to the springs, which was much larger than the tension due to motion, thus the restoring effect of the mooring system was always active. A coupling between sway and roll motion is assumed because the center of gravity of the model and the mounting point of the mooring lines were not in the same height. Considering the abovementioned remarks, the rigid body motion is described by Equations (28).
In the above equations m is the mass of the structure, I is the inertia of momentum about its center of gravity, n is the normal unit vector defined on the boundary S of the body, n y , n z the components of n, p is the fluid pressure andτ is the viscous stress tensor. Lastly, K 22 and K 44 are the spring constants for the heave and roll motion, respectively and K 24 = K 42 are the spring constants for coupling of the heave and roll motion due to the level arm formed between the axis of the restoring forces and the center of gravity. A more detailed report for the experimental setup and the physical modeling of the rigid body motion can be found in [53]. The experiment considered three types of waves based on the wave height to wavelength ratio for several wave periods. In the present study, only results for waves with ratio H/λ = 1/30 for wave periods between the range [0.6 : 1.2 s] are presented.
The flow is considered fully turbulent for all wave periods. This is valid for a large amplitude of motion, where wave-breaking occurs. In smaller amplitudes, the production of turbulent viscosity is lower resulting in an almost laminar approximation. The grid and the timestep is the same for all cases. In the stream-wise, direction the grid is almost uniform with size of dy = 0.008 m and in the transverse direction in the near free surface regime dz = 0.001 m. On the damping zone region, the mesh coarsens with a geometric rule to effectively damp the outgoing waves. The grid around the structure has an O-grid topology which is dense to account for the vortex formation near the sharp edges. The total size of the mesh is 3 · 10 5 cells. Based on a sensitivity study, the timestep is selected equals to 0.001 s .
In Figure 10 the results are compared with the experimental data. In Figure 10a-c the Response Amplitude Operators (RAOs) of the motions are presented, while Figure 10d presents the amplitude of the space averaged free surface elevation inside the moonpool. Regarding the sway motion, the amplitudes scale linearly with respect to the wave period and the results coincide with the experimental ones. The heave amplitudes appear to have a resonant behavior for wave periods between T = 0.65 s and T = 0.8 s. This behavior is depicted in the amplitudes of the free surface as well. Small discrepancies are observed near this resonant region. Furthermore, the Figure 10c shows the amplitudes of the roll motion, except for large wave periods the present method produces similar results with the experiments. However, larger amplitudes of roll motion are obtained for wave periods between T = 1 s and T = 1.2 s. This could be caused due to the assumption of two-dimensional flow. Additional drag forces could be act on the structure caused by the liquid between the side of the structure and the wave tank.
In Figure 11 vorticity near the edges of the structure is depicted where the flow separation at the corners of the structure is evident. The flow separation at the lower entrance of the moonpool regulates the resonant behavior of the piston mode which is apparent for wave periods between T = 0.75 s and T = 0.85 s. Finally, the pressure field upstream of the structure is presented in Figure 12. The figure illustrates the pressure distribution at the incident wave, the free surface elevation inside the moonpool, as well as the pressure drop due to vortex formation near the sharp edges of the structure.

Conclusions
Concluding, in the present work a two-phase flow solver was formulated based on the artificial compressibility approach and the volume of fluid method. In the first part, the derivation of the hyperbolic system of equations was presented. A detailed analysis of the time discretization process was illustrated, while emphasis was given to account for the moving boundaries. Regarding the spatial discretization the finite volume method was used, while an approximate Riemann solver was used for the evaluation of the inviscid fluxes. Finally, a simple and efficient implicit methodology to generate and absorb free surface waves was illustrated.
A series of numerical test cases was performed to demonstrate the efficiency and consistency of the solver. A simple wave propagation problem on constant bathymetry was chosen to illustrate the ability of the solver to generate, propagate and damp free surface waves. Also, the same configuration was used to investigate the effect of the artificial compressibility parameter-β. The results indicate that there is a wide range of values where no significant differences were observed. Regarding the last two validation examples, the predictions of the solver were compared to available experimental data. In case of the variable bathymetry the solver was able to perform well and account for the higher harmonics developed due to shoaling.
In the final test case, it was shown that the present methodology was able to produce adequate results for a surface piercing problem with wave-structure interaction. The solver was able to handle the fluid-structure interaction problem accurately being in close agreement with the measurements.
As a last remark, the proposed methodology compared to the alternatives provided by the literature solves a single system of equations in a coupled manner. The formulation enables simulation of free surface flows using a hyperbolic perspective. The method performed well in a variety of cases, providing a viable alternative for modeling free surface flows within the CFD framework. Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Discretization of the Hyperbolic System of Equations
The equations of flow (7) are discretized-based on the finite volume method. Control volumes are defined in every cell of the mesh with its center being the geometric center of the mesh elements. The unknown variables Q are expressed over a control volume D i as The surface terms of Equation (7) are considered constant in each face of the control volume. Thus, the surface integrals are approximated through a sum of surface terms evaluated at the midpoint of every face. Furthermore, the volume terms are considered constant in each D i . As a result, the spatial terms are computed based on the following equation.
The inviscid fluxes are computed through the approximate Riemann solver of Roe [54]. As mentioned in Section 2, the eigenvalues of the inviscid Jacobian depend on the density field. For this reason, the preconditioned matrix Γ is used in the evaluation of fluxes. The preconditioned Jacobian A c is defined as where Γ is the preconditioned matrix (9) and Γ −1 is the inverse matrix. The matrix A c is given by Equation (A4).
0 n x n y n z 0 n x ρ m (n x u + ∆V) ρ m n y u ρ m n z u u∆V∆ρ n y ρ m n x v ρ m n y v + ∆V ρ m n z v v∆V∆ρ n z ρ m n x w ρ m n y w ρ m (n z w + ∆V) w∆V∆ρ 0 α l n x α l n y α l n z ∆V The convective fluxes are computed, at a face f , using the Roe approximate Riemann solver as The indices · R , · L denote the right and the left state of the face, as defined by the normal vector. The Jacobian matrix is computed based on the absolute values of the eigenvalues of the system. Specifically, let R −1 be the matrix constructed by arranging the left eigenvectors of A c in a column wise order, R be the matrix constructed by arranging the right eigenvectors in a row wise order and Λ be the diagonal matrix constituted by the absolute values of the eigenvalues, then the Jacobian matrix can be written as The inviscid Jacobian of Equation (A6) is computed based on the Roe averaged quantities at face f , given by (A7). v = The eigenvalues of the Jacobian matrix A c are given by equation (A8). The eigenvalues have similar form as in the case of the compressible equations. The difference is the pseudo-sound speed does not depend on the flow characteristics, but on the choice of β parameter and the hyperbolic definition of the equations is valid only during pseudo-time iterations.
In previous equations, c is the artificial sound speed, where in case of two-phase flows, is defined as Finally, the eigenvectors of the Jacobian (A3) are given by Equations (A10) and (A11). 1 c g β n x y 1 − n y x 1 + ∆V (vx 1 − uy 1 )] 0 − 1 ρ m 2c c g (β + λ 4 V n ) 1 2c c g βc p n x 1 2c c g βc p n y 1 2c c g βc p n z 0 1 ρ m 2c c g (β + λ 3 V n ) 1 2c c g βc m n x 1 2c c g βc m n y 1 2c c g βc m n z 0 In the preceding equations c m , c p and c g are expressed as while for the unit vectors x 1 = (x 1 , y 1 , z 1 ) and x 2 = (x 2 , y 2 , z 2 ) should hold x 1 · n = x 2 · n = x 3 · n = 0.