Two-Phase Smoothed Particle Hydrodynamics Modelling of Hydrodynamic-Aerodynamic and Wave-Structure Interaction

: A two-phase (air and water) smoothed particle hydrodynamics (SPH) method is employed to study the hydrodynamic-aerodynamic and wave interaction with ﬁxed and ﬂoating structures in a wave basin. The method is ﬁrst veriﬁed for a classical two-phase dam-breaking. A mirror-open boundary is implemented at the top and left sides of a two-phase wave basin with a piston to generate a second-order regular wave. It is observed that, compared to the single-phase simulation, the two-phase one obtains a smoother water surface and prevents the non-physical water splash when interacting with the sloped dissipative beach. This wave basin is also used to investigate wave-structure problems such as wave interaction with a rigid cantilever beam ﬁxed to the basin bottom and downstream of the wave-maker mechanism and the dynamics of a single ﬂoating box and two ﬂoating boxes in the waves. A typical wave-structure interaction period is captured and described using pressure contours and velocity vectors at three selected instants for the wave-rigid cantilever beam case. With the increase of the structure’s height, the wave height after the structure decreases, but no evident variation is found when changing its thickness. Besides the hydrodynamics interaction, a periodical collision is observed between the two ﬂoating boxes on the wave surface.


Introduction
Typical wave-structure interaction (WSI) problems involving one or multiple solid structures interacting with surface waves and fluid flow are one of the most significant problems in oceanographic engineering with implications for the design of moored and freely floating bodies on the surfaces of open seas. For example, breakwaters and seawalls are two widely used coastal structures for enabling the safe navigation of marine vessels into harbors and maintaining calm sea-keeping conditions for portside marine operations and shore protection, respectively, as outlined in [1]. WSI also has implications for the safety and stability of floating structures such as ships or offshore platforms in the presence of unsteady wind and wave loads, as outlined in [2]. Time-varying-ocean wave loads determine the durability and resistance to fatigue; one of the key concepts in material selection, design, and manufacturing of very large floating structures (VLFS) as alternative spaces for coastal cities, as outlined in [3]. WSI is essentially a multi-physics non-linear fluid-structure interaction (FSI) problem, for which analytical solutions to the governing equations are non-existent for most practical problems, and laboratory experiments are often limited in scope. Hence, computational modeling is often a method for facilitating solutions to this problem, such as those outlined in [4,5].
Linear and weakly non-linear WSI problems can be computed accurately and economically by using the potential flow equation for inviscid irrotational flows, which has its limitations in handling highly non-linear waves and situations involving wave-breaking phenomena. Second-order wave interactions with fixed and floating structures have been investigated successfully in [6][7][8][9]. High Fidelity Computational Fluid Dynamics (CFD) based on the Eulerian Lagrangian form of Navier-Stokes equations, using mesh-based and meshless approaches, offer alternative routes for handling these highly non-linear problems. The mesh-based approaches use the finite difference (FD) or the finite volume method (FV) or the finite element (FE) spatial discretization of the flow equations for numerical solutions and are often coupled with interface-capturing/tracking techniques such as level set methods or Volume of Fluid (VOF) to capture the wave interfaces. Adaptive and moving meshes or the immersed boundary methods are also used to track the interface between the fluid and the structure. For example, Hadžić et al. [10] predicted the motion of a two-dimensional free-floating rectangular box in shallow-water waves using the FVM approach using mesh motion, mesh deformation, and sliding mesh interface approaches. Hu and Kashiwagi [11] used a constrained interpolation profile (CIP)-based Cartesian grid method to predict the strong non-linear interaction between the regular waves and a floating structure with its shape being represented by virtual particles. Chen et al. [12] used the open-source flow software OpenFOAM to compute regular wave interactions with a vertical surface piercing cylinder. Hu et al. [13] computed the extreme wave interactions with a fixed/floating truncated cylinder, and Li et al. [14] modeled the wave interactions with the marine structures.
Meshless approaches such as smoothed particle hydrodynamics (SPH) developed originally for astrophysical applications by Lucy [15] and Gingold et al. [16] have been extended to model hydrodynamic wave generation and various WSI problems. Ni and Feng [17] have used the open-source SPH-based DualSPHysics to model a two-dimensional numerical wave tank based on a source generation and absorption technology with an analytical relaxation approach using a generation zone instead of moving boundaries to model the wave based on the Stokes theory and hence limiting their model to regular waves only. Wen et al. [18] implemented an absorbing wavemaker in DualSPHysics whereby the wave generation was based on linear wave theory, thereby restricting their model to handle only regular waves. Recently, Altomare et al. [19] implemented a long-crested wave generation and absorption model in DualSPHysics to investigate second-order long-crested monochromatic and random waves by using moving boundaries to mimic the action of a piston-type wavemaker. The estimates of water surface elevation, wave orbital velocities, wave forces, and the capacity for damping the re-reflection in the flow domain have been validated against theoretical solutions and experimental results. Meringolo et al. [20] simulated wave loads using a two-dimensional diffusive weakly compressible SPH model to estimate the hydraulic characteristics of perforated breakwaters from the interaction between regular waves with fully and partially perforated breakwaters. Recently, Huang et al. [21] developed a coupled finite particle method (FPM) to combine the advantages of SPH, the δ-SPH model, and particle shifting technology (PST) to compute the unsteady hydrodynamics of both two-dimensional and three-dimensional oscillating wave surge converters.
The present work is concerned with developing and applying of a two-phase weakly compressible SPH method for various wave hydrodynamics and WSI problems. The twophase SPH is applied to the classical dam-breaking test problem to validate the approach. This is followed by modeling a second-order wave generated in a wave basin. The wave basin is then employed to investigate several WSI problems, including the wave interaction with a rigid cantilever beam fixed to the basin bottom and downstream of the wave-maker mechanism, and the dynamics of a single floating box and two floating boxes in the waves. In the framework of this two-phase SPH, our work expects to capture and understand the physics of these hydrodynamic wave-structural and dynamical interaction problems.
This paper aims to explore the difference between two-phase and single-phase flow on the fluid-structure interaction by using a smoothed particle hydrodynamics method. Our target is to study the hydrodynamic interaction between waves and fixed structures and floating boxes in the case of two-phase (air and water) flow, to and reveal the corresponding mechanism.

Smoothed Particle Hydrodynamics
The Lagrangian form of the conservation laws for the continuity and momentum describing the flow of a weakly compressible flow can be written as where u and ρ are the total velocity and the density field of the fluid, respectively, f is an acceleration resulting from an external force, σ is the stress tensor consisting of pressure gradient and viscous terms. In the framework of SPH, the discretized form of the continuity equation, i.e., Equation (1a) is written as where W is the smoothing kernel function, for which the cubic spline kernel outlined in Monaghan and Lattanzio [22] is used in this study. Figure 1 shows a schematic of an SPH particle in its support domain wherein the physical quantities of particle i are obtained by interpolating the quantities of the particles in its support domain using the kernel function.
The kernel function has a compact support domain Ω which is a circle of radius κ 0 h, where κ 0 is a parameter varying with the chosen kernel function shown in Figure 1, and it is equal to two for cubic spline kernel. The smoothing length h = δ s R is related to the support area over which approximations are estimated, R is the resolution, i.e., initial interparticle gap, and δ s is the smoothing factor, the value of which is typically in the range of 1.25~2 depending on the specific flow being considered. For WSI problems, convergent results are usually obtained by choosing the smoothing factor as 1.5, which has been used for most of the problems reported in this study, unless specified otherwise. The additional term δ i is introduced to improve the stability of the SPH numerical calculation [23] and to remove spurious pressure oscillations that can exist in weakly compressible SPH. The vector u ij denotes the difference between u i and u j at locations i and j, as indicated in Figure 1.
where u and ρ are the total velocity and the density field of the f acceleration resulting from an external force, σ is the stress tens gradient and viscous terms. In the framework of SPH, the discret ity equation, i.e., Equation (1a) is written as where W is the smoothing kernel function, for which the cubic Monaghan and Lattanzio [22] is used in this study. Figure 1 show particle in its support domain wherein the physical quantities of interpolating the quantities of the particles in its support doma tion. The kernel function has a compact support domain Ω whic where κ0 is a parameter varying with the chosen kernel function it is equal to two for cubic spline kernel. The smoothing length support area over which approximations are estimated, R is the terparticle gap, and δs is the smoothing factor, the value of which of 1.25~2 depending on the specific flow being considered. For W results are usually obtained by choosing the smoothing factor as for most of the problems reported in this study, unless specified o term δi is introduced to improve the stability of the SPH numeric remove spurious pressure oscillations that can exist in weakly com tor uij denotes the difference between ui and uj at locations i and 1.  The discretized form of the momentum conservation equation, i.e., Equation (1b), for the single-and two-phase flows considered in this study takes the following form as shown in [24] Du where P i and P j are the pressures at particle i and particle j, respectively, r ij is the displacement vector from particle i to particle j, g is the gravitational acceleration and µ is the dynamic viscosity of the fluid which is assumed to be constant in this study. The term Π ij denotes the artificial viscosity introduced to alleviate pressure oscillations caused by a non-physical concentration of particles and is modeled as in [25] where φ ij = hu ij r ij |rij| 2 +η 2 , ρ ij and c ij are respectively the average values of the density and the sound speed of particle i and j. η = 0.1 h is chosen to prevent the numerical divergence when particles are close to each other, α is a parameter which is set according to the flow problem with the guidance that a proper value of α should not be too large as it may lead to a large error, and should not be too small to avoid being ineffective in alleviating the oscillations in dynamic problems. For all the WSI problems considered in this study, the value of α is in the range 0.01~0.05. For the weakly compressible flow, an equation of state which is commonly used in SPH to relate the pressure to small variations of fluid density is applied where B is the reference pressure defined as B = ρ 0 c 2 γ 0 , ρ 0 is the fluid reference density (for the water and air, they are, respectively, 1000 kg/m 3 and 1.29 kg/m 3 ), c is the speed of sound in the fluid, γ 0 is a constant, and a value γ 0 = 7 is used for the problems in this study. A proper value of B is chosen to ensure that the compressibility effects are negligible [26], and for this, the speed of sound should satisfy c ≥ 10 max(||u max ||). The general form of Equations (2a) and (2b) can be written as

Boundary Conditions
For the problems of interest, both wall boundary and open boundary conditions are considered. The no-slip boundary conditions are imposed at solid walls as specified in [27]. As shown in Figure 2a, the pressure on the solid particle i is extrapolated from the neighboring fluid particle j, Energies 2022, 15, 3251 The corresponding densities are then estimated using Equation (3).
Schematic of the open boundary condition for the outflow. The translucent areas denote the supporting domains. The mirror particles are mirrored by the ghost particles concerning the outflow threshold. (a) wall boundary, the pressure on the solid particle i is extrapolated from the neighboring fluid particle j; (b) open boundary, the physical information (pressure, density, or velocity) at a ghost particle is calculated (Equation (7)) using that of the mirror particles; the information of a mirror particle is interpolated (Equation (6)) from that of the fluid particles in its supporting domain.
Open boundary conditions are usually implemented at the inflow or the outflow boundaries of the flow domain of interest by imposing zero velocity or zero pressure gradients. Recently, a versatile open boundary condition for SPH simulation outlined in [28] incorporates a buffer layer consisting of ghost particles, the width of which usually corresponds to the support domain radius. Figure 2b shows a schematic of the implementation of open boundary conditions at the outflow boundary, indicating the mirror particles mirrored from ghost particles relative to the outflow threshold for facilitating density and velocity interpolations. The new incoming particles will be added layer by layer with a uniform interval (initial interparticle gap), and the outgoing particles will be deleted. At the inflow boundaries, the velocity of the ghost particles is prescribed based on the flow problem, and the density is interpolated from the fluid or ghost particles. Both the ghost particles' velocity and density are interpolated from the fluid or ghost particles at the outflow boundaries.
The physical quantities of the mirrored particles are then projected back to the corresponding ghost particles using the first-order Taylor series expansion. The interpolation schemes used for estimating the approximate values of the physical quantities Rm and their first-order derivatives Rm,β ensures the consistency of the kernel function as pointed out in [29], and takes the form (a) wall boundary, the pressure on the solid particle i is extrapolated from the neighboring fluid particle j; (b) open boundary, the physical information (pressure, density, or velocity) at a ghost particle is calculated (Equation (7)) using that of the mirror particles; the information of a mirror particle is interpolated (Equation (6)) from that of the fluid particles in its supporting domain.
The corresponding densities are then estimated using Equation (3). Open boundary conditions are usually implemented at the inflow or the outflow boundaries of the flow domain of interest by imposing zero velocity or zero pressure gradients. Recently, a versatile open boundary condition for SPH simulation outlined in [28] incorporates a buffer layer consisting of ghost particles, the width of which usually corresponds to the support domain radius. Figure 2b shows a schematic of the implementation of open boundary conditions at the outflow boundary, indicating the mirror particles mirrored from ghost particles relative to the outflow threshold for facilitating density and velocity interpolations. The new incoming particles will be added layer by layer with a uniform interval (initial interparticle gap), and the outgoing particles will be deleted. At the inflow boundaries, the velocity of the ghost particles is prescribed based on the flow problem, and the density is interpolated from the fluid or ghost particles. Both the ghost particles' velocity and density are interpolated from the fluid or ghost particles at the outflow boundaries.
The physical quantities of the mirrored particles are then projected back to the corresponding ghost particles using the first-order Taylor series expansion. The interpolation schemes used for estimating the approximate values of the physical quantities R m and their first-order derivatives R m,β ensures the consistency of the kernel function as pointed out in [29], and takes the form where the index m denotes the mirror particles, and j represents the fluid or ghost particle inside the support domain of particle m, β and γ are the indices of dimensions repeated from one to two or three, and W mj ,γ = ∂W mj /∂x γ (∂x γ means ∂x, ∂y or ∂z) are the spatial derivatives of the kernel function for which a first-order Taylor series expansion around the location of particle m can be written as where n 0 denotes the ghost particle corresponding to the mirror particle m.

Time Integration Schemes
A velocity-Verlet [30] scheme, a special form of the Newark-beta method for integrating differential equations as outlined in [31] was used for the time integration. Note that for the two-phase flows, the continuity and momentum equations of the air phase and water phase are calculated simultaneously (Equations (1a) and (1b)), hence the "fluid" is used to denote both the air and water in the time integration process. The velocity and displacement of the fluid particles are updated at the end of the first half of each time step where the superscript n denotes the time step. Then the continuity Equation (2a) is integrated using the updated velocity Equation (8a) to update the density and the displacement Equation (8b) of the fluid particles over the next half of the time step so that at the end of each time step For a floating structure governed by the rigid body motion equations, its velocity and displacement are updated after the total force F and total moment T acting on the solid body are obtained where X 0 , U, Ω, and r s are, respectively, the mass center location, the translation velocity, the rotation velocity, and the solid particle location of the structure. In Equations (8e) and (8f), the functions are used for the convenience of expressing the suffered total force and moment of a floating structure. The forces on the boundary particle i of the floating structure exerted by the neighboring fluid particle j are based on Equation (2b) where M and I are, respectively, the mass and the moment of inertia of the floating structure.
Note that for the case of two freely floating boxes with collision, an elastic collision force which will be addressed later in this section, should be added in Equations (8k) and (8l), respectively. Subsequently, the velocity of the fluid particles is updated at each half time step, Using the velocity computed at the end of the half-time step, the force is computed only once per time step to obtain the particle acceleration after the density is updated to Equation (8c). To obtain a proper time step, several criteria outlined in [25] must be satisfied for the stability reasons of the numerical scheme. Namely, the CFL-condition related to the smoothing length h, the maximum sound speed c max (c max must be estimated according to the specific flow problems before calculating, and the actual value may be smaller for obtaining better numerical results), the maximum flow speed u max , the viscosity v as well as the body force g, are adopted to estimate the global time step. It requires, (9c)

Collision Force Model
When the two boxes are very close to each other, a collision occurs and this is modeled based on the discrete element method (DEM) [32]. As indicated in Figure 3, the contact force F c acting on box A resulting from a collision with Box B can be decomposed into the normal and tangential components F n and F t , respectively. Both the forces can be further decomposed into a repulsion force F r and a damping force F d , which are respectively caused by the deformation of the material and the dissipation of energy during the deformation. The normal forces can be expressed as [33] where k n,AB is the stiffness constant of the boxes; δ AB = R-r AB is the box overlap (see Figure 3, R is the initial interparticle gap or the resolution); e AB is the unit vector between particles of the two boxes with the minimum distance (see Figure 3, the distance between point A and point B); . δ AB = v AB · e AB is the rate of normal deformation (see Figure 3, v AB is the velocity difference between point A and point B); γ n,AB is the damping constant. The stillness and damping constants are, k n,AB = 4 3 E * √ R * and γ n,AB =C n 6M * E * √ R * , respectively, where C n is in order of 10 −5 [34]; and M* are, respectively, half of the height and mass of the boxes. The tangential forces are given similarly to [35] δ t AB e t AB with k t,AB = 2 7 k n,AB and γ t,AB = 2 7 γ n,AB ; the superscript or subscript t denotes the tangential direction; µ f,AB is the kinetic friction coefficient. In the present study, the boxes are assumed to be made of hardwood, corresponding to Young modulus E b = 2 × 10 8 N/m 2 , Poisson ratio v b = 0.2, and the kinetic friction coefficient µ f,AB = 0.7.

Numerical Results and Discussion
In this section, the SPH simulations of several two-phase WSI examples are discussed. First, a classical two-phase dam-breaking problem is addressed to show the convergence issues of the SPH method. Next, a two-phase SPH simulation of a second-order wave generated by a piston in a wave basin is validated to capture the wavy air-water interface. This is followed by the two-phase SPH simulation of the interaction between the wave generated in the wave basin and a rigid cantilever beam fixed to the bottom of the wave basin and downstream of the wavemaker mechanism. Subsequently, the two-phase SPH simulation of the interaction between the wave and a floating box to capture the wave-induced motion of the box is conducted. The two-phase SPH simulation of the interaction between the wave and two identical floating boxes with the possibility of elastic collision between them is considered finally.

SPH Simulation of Dam-Breaking
The two-phase dam-breaking problem is a standard benchmark test case that has several hydrodynamic phenomena such as a surface break-up, water jet, and impact pressure. This test was chosen for validating the accuracy and stability of the present twophase SPH code given the large difference in the densities of the two phases, namely air and water. The computational domain is shown in Figure 4, the height and length of the computational domain are, respectively, 0.5 m and 1.6 m. The water is initially contained in a rectangular domain (dam) with a length L of 0.6 m and a height H of 0.3 m (shaded in light green) implying a volume of water per unit width of 0.18 m 3 . The remaining area of the computational domain (shaded in light blue) is assumed to be filled with air. Both air and water are assumed to be at rest initially and t = 0 is defined as the instant at which the water dam is allowed to break free and fall. The no-slip condition is imposed on the solid walls. The total number of particles N used for this computation is 158916 (5316 for solid walls, 34560 for water, and 119040 for air). A smoothing length h = 1.5 R with a resolution R = 0.0025 m is used. Following the work of Gong et al. [36], the laminar viscosity is ig-

Numerical Results and Discussion
In this section, the SPH simulations of several two-phase WSI examples are discussed. First, a classical two-phase dam-breaking problem is addressed to show the convergence issues of the SPH method. Next, a two-phase SPH simulation of a second-order wave generated by a piston in a wave basin is validated to capture the wavy air-water interface. This is followed by the two-phase SPH simulation of the interaction between the wave generated in the wave basin and a rigid cantilever beam fixed to the bottom of the wave basin and downstream of the wavemaker mechanism. Subsequently, the two-phase SPH simulation of the interaction between the wave and a floating box to capture the waveinduced motion of the box is conducted. The two-phase SPH simulation of the interaction between the wave and two identical floating boxes with the possibility of elastic collision between them is considered finally.

SPH Simulation of Dam-Breaking
The two-phase dam-breaking problem is a standard benchmark test case that has several hydrodynamic phenomena such as a surface break-up, water jet, and impact pressure. This test was chosen for validating the accuracy and stability of the present two-phase SPH code given the large difference in the densities of the two phases, namely air and water. The computational domain is shown in Figure 4, the height and length of the computational domain are, respectively, 0.5 m and 1.6 m. The water is initially contained in a rectangular domain (dam) with a length L of 0.6 m and a height H of 0.3 m (shaded in light green) implying a volume of water per unit width of 0.18 m 3 . The remaining area of the computational domain (shaded in light blue) is assumed to be filled with air. Both air and water are assumed to be at rest initially and t = 0 is defined as the instant at which the water dam is allowed to break free and fall. The no-slip condition is imposed on the solid walls. The total number of particles N used for this computation is 158916 (5316 for solid walls, 34560 for water, and 119040 for air). A smoothing length h = 1.5 R with a resolution R = 0.0025 m is used. Following the work of Gong et al. [36], the laminar viscosity is ignored. The artificial viscosity term α of 0.05 and a fixed computational time step of 5 × 10 −5 s have been used for this problem. The evolution of the solution in the computational domain is then tracked in terms of a non-dimensional time step t defined as t(g/H) 0.5 . Note that the theory of the two-and single-phase simulation is integrated into Equation (4), hence we only need to remove the air particles when considering the single-phase flow problem.   [37] at the same instant. It can be seen that after impacting the vertical wall w the pressure sensor is located (see Figure 4), the water height reaches a relatively higher level than that predicted by a single-phase SPH and also with that predicted in [37]. An area of air entrapment can also be seen here. Apart from the area of the air bubble entrapped by the backward plunging breaker, it can be seen that the outline of the water height from the single-phase flow is comparable with that from the two-phase flow. It is observed that the two-phase SPH simulation forms a smoother interface of entrapped air due to the weak compressibility of the air phase there. Figure 5. Comparison of single-and two-phase SPH results with that of Chen et al. [37] for dambreaking flow at dimensionless time t(g/H) 0.5 = 6.4. The hollow "triangles" and "diamonds" denote the surfaces of the water obtained from the present single-phase simulation and reference [37], respectively. The blue solid "dots" denote the water particles, in which the air phase has been hidden for a distinct water-air interface. Figure 6 compares the computed spatial pressure contours from the present twophase SPH model with that from reference [37] at the same instant, and while there are differences in the contour levels, the spatial distribution is similar. This is due to the dif-   [37] at the same instant. It can be seen that after impacting the vertical wall w the pressure sensor is located (see Figure 4), the water height reaches a relatively higher level than that predicted by a single-phase SPH and also with that predicted in [37]. An area of air entrapment can also be seen here. Apart from the area of the air bubble entrapped by the backward plunging breaker, it can be seen that the outline of the water height from the single-phase flow is comparable with that from the two-phase flow. It is observed that the two-phase SPH simulation forms a smoother interface of entrapped air due to the weak compressibility of the air phase there.   [37] at the same instant. It can be seen that after impacting the vertical wall w the pressure sensor is located (see Figure 4), the water height reaches a relatively higher level than that predicted by a single-phase SPH and also with that predicted in [37]. An area of air entrapment can also be seen here. Apart from the area of the air bubble entrapped by the backward plunging breaker, it can be seen that the outline of the water height from the single-phase flow is comparable with that from the two-phase flow. It is observed that the two-phase SPH simulation forms a smoother interface of entrapped air due to the weak compressibility of the air phase there. Figure 5. Comparison of single-and two-phase SPH results with that of Chen et al. [37] for dambreaking flow at dimensionless time t(g/H) 0.5 = 6.4. The hollow "triangles" and "diamonds" denote the surfaces of the water obtained from the present single-phase simulation and reference [37], respectively. The blue solid "dots" denote the water particles, in which the air phase has been hidden for a distinct water-air interface. Figure 6 compares the computed spatial pressure contours from the present twophase SPH model with that from reference [37] at the same instant, and while there are differences in the contour levels, the spatial distribution is similar. This is due to the different wall boundary conditions, in which a comparatively higher pressure is obtained in Figure 5. Comparison of single-and two-phase SPH results with that of Chen et al. [37] for dambreaking flow at dimensionless time t(g/H) 0.5 = 6.4. The hollow "triangles" and "diamonds" denote the surfaces of the water obtained from the present single-phase simulation and reference [37], respectively. The blue solid "dots" denote the water particles, in which the air phase has been hidden for a distinct water-air interface. Figure 6 compares the computed spatial pressure contours from the present twophase SPH model with that from reference [37] at the same instant, and while there are differences in the contour levels, the spatial distribution is similar. This is due to the different wall boundary conditions, in which a comparatively higher pressure is obtained in our simulations (see Figure 7). Note that the air and the water are presented in Figure 6. The splashed water particles (see Figure 6a) with uneven distribution may be affected by the higher pressure (compared to that of the reference [37]). Note that the pressure below the cavity is comparatively larger than that acting on the cavity top. This pressure difference leads to the collapse of the cavity for the single-phase flow without the air (see Figures 5 and 6a).     Figure 4, corresponding to single-and two-phase SPH simulations from the present study and also with that from Chen et al. [37] and Buchner [38]. For the two-phase SPH simulations, two resolutions have been considered, i.e., R = 0.005 m and 0.0025 m, to gauge the convergence of the pressure temporal variation. It can be seen that the choice of the resolution of 0.0025 results in a smoothening of the pressure variation and improves the predicted temporal pressure evolution agreement with that of Chen et al. [37]. Compared to Buchner's experimental data [38], it appears that the present SPH models give a reliable prediction of the pressure evolution. The pressure of the two-phase simulation lags a little from the single-phase result in time scale, which is due to the compression of the surrounding air. The simulation results also lag from the experimental data, which may be due to the present SPH's weakly compressibility. The pressure peak of the sensor in our simulation appears at t = 6.4 when the backward plunging breaker entraps the air and creates the cavity (see Figures 5 and 6). Subsequently, it gradually reaches the expected hydrostatic pressure.
(a) (b) Figure 6. Comparison of pressure for two-phase SPH simulation of dam-breaking flow at dimensionless time t(g/H) 0.5 = 6.4. (a) present study; (b) study of Chen et al. [37]. The measurement unit is m in the x-and y-axis.   Figure 4, corresponding to single-and two-phase SPH simulations from the present study and also with that from Chen et al. [37] and Buchner [38]. For the two-phase SPH simulations, two resolutions have been considered, i.e., R = 0.005 m and 0.0025 m, to gauge the convergence of the pressure temporal variation. It can be seen that the choice of the resolution of 0.0025 results in a smoothening of the pressure variation and improves the predicted temporal pressure evolution agreement with that of Chen et al. [37]. Compared to Buchner's

Two-Phase Wave Generation in a Wave Basin
This section is concerned with the generation of waves in a two-phase flow in a wave basin. Based on the assumptions of irrotational and incompressible fluid and constant pressure at the free surface, the Biesel transfer function outlined in [39] and also in [40] can be used to describe the relation between wave amplitude and wavemaker displacement. The far-field solution for the free-surface elevation is expressed as where H 0 is the wave height, ω = 2π/T is the angular frequency, k = 2π/λ is the wavenumber, T is the wave period, λ is the wavelength and the initial phase δ 0 is given by a random number between 0 and 2π. For a piston-type wavemaker, the Biesel transfer function under the hypothesis of monochromatic sinusoidal waves is expressed as where S 0 is the piston stroke and d is the water depth. The time series of the piston motion is then given by Madsen [41] has developed a simple second-order wavemaker theory to generate relatively long second-order Stokes waves that preserve shape as they propagate. To generate second-order regular waves, an additional term is added to Equation (13) to account for the piston displacement so that In the present study, the second-order Stokes waves are generated using a piston-type wavemaker as shown in Figure 8, which shows the wave basin configuration and the computational domain for the two-dimensional two-phase wave basin. In the following content, based on the above wave making theory, we first generate the second-order Stokes wave using the SPH method, and then validate our results against the theoretical solution. The selection of the wave parameters is based on the existing literature research and the specific problems we face in our project. This step is to lay the foundation for the later study of the interaction between the waves and fixed structures, and between the waves and the freely floating boxes. particles for single-phase SPH, and 53161 particles (31261 particles for the air and 28584 particles for the water) for the two-phase SPH simulations. The artificial viscosity parameter α in Equation (2b) and the computational time step are set to 0.01 and 1 × 10 −4 s, respectively. Note that in this case, the cantilever beam does not exist. Figure 9 compares the computed temporal variation of the free surface wave elevations recorded at the wave gauges WG1 and WG3 placed at x = 2 m and 4 m, respectively. It is observed that the results from the present single-phase and two-phase SPH and Du-alSPHysics match the theoretical solution to a good extent, although the wave troughs are slightly overestimated (less than 3%) by SPH, indicating that the waves are properly generated, and the effect of the air on the wave elevations appears to be negligible. Note that the results of DualSPHysics are obtained by us, rather than from the literature. Figure 10 compares the computed velocity components from theory and DualSPHysics corresponding to the same case shown in Figure 9. The horizontal and vertical velocity components are monitored at a height d2 of 0.12 m above the basin in the WG3 profile. Apart from the horizontal velocity, which is slightly underestimated at the peak, the results could provide a satisfactory accuracy. Note that in this case, the cantilever beam does not exist. Figure 9 compares the computed temporal variation of the free surface wave elevations recorded at the wave gauges WG1 and WG3 placed at x = 2 m and 4 m, respectively. It is observed that the results from the present single-phase and two-phase SPH and DualSPHysics match the theoretical solution to a good extent, although the wave troughs are slightly overestimated (less than 3%) by SPH, indicating that the waves are properly generated, and the effect of the air on the wave elevations appears to be negligible. Note that the results of DualSPHysics are obtained by us, rather than from the literature. Figure 10 compares the computed velocity components from theory and DualSPHysics corresponding to the same case shown in Figure 9. The horizontal and vertical velocity components are monitored at a height d 2 of 0.12 m above the basin in the WG3 profile. Apart from the horizontal velocity, which is slightly underestimated at the peak, the results could provide a satisfactory accuracy. the results of DualSPHysics are obtained by us, rather than from the literature. Figure 10 compares the computed velocity components from theory and DualSPHysics corresponding to the same case shown in Figure 9. The horizontal and vertical velocity components are monitored at a height d2 of 0.12 m above the basin in the WG3 profile. Apart from the horizontal velocity, which is slightly underestimated at the peak, the results could provide a satisfactory accuracy.  Figure 11a,b shows the contours of the horizontal velocity contours for single-phase and two-phase wave generation at the instant t = 16 s, respectively. In Figure 11b, the air particles are not shown deliberately to understand the physics predicted by the singleand two-phase SPH. It can be observed that the same velocity patterns and free-surface profiles are predicted by the two approaches from a global viewpoint. However, if one focuses on the circled regions in the vicinity of the sloped dissipative beach area labeled and shown in Figure 11c,d, respectively, it is observed that the water particles of the single-phase simulation splash more severely than those of the two-phase simulation. This indicates that air in the vicinity could smoothen the air-water interface by complementing the particles for the calculation in the vicinity of the water surface area where the particle supporting domain is truncated for the single-phase problem.  Figure 11a,b shows the contours of the horizontal velocity contours for single-phase and two-phase wave generation at the instant t = 16 s, respectively. In Figure 11b, the air particles are not shown deliberately to understand the physics predicted by the single-and two-phase SPH. It can be observed that the same velocity patterns and free-surface profiles are predicted by the two approaches from a global viewpoint. However, if one focuses on the circled regions in the vicinity of the sloped dissipative beach area labeled and shown in Figure 11c,d, respectively, it is observed that the water particles of the single-phase simulation splash more severely than those of the two-phase simulation. This indicates that air in the vicinity could smoothen the air-water interface by complementing the particles for the calculation in the vicinity of the water surface area where the particle supporting domain is truncated for the single-phase problem. Figure 12a shows the contours of the density predicted using the two-phase SPH for the same case at the instant t = 16 s. Both air and water phases can be distinguished because of the large density difference between the two phases. In the vicinity of the wave crest, as indicated in Figure 12c, the entraining of air into the water creates a circulatory flow, indicating the continuity of velocity at the air-water interface. and two-phase SPH. It can be observed that the same velocity patterns and free-surface profiles are predicted by the two approaches from a global viewpoint. However, if one focuses on the circled regions in the vicinity of the sloped dissipative beach area labeled and shown in Figure 11c,d, respectively, it is observed that the water particles of the single-phase simulation splash more severely than those of the two-phase simulation. This indicates that air in the vicinity could smoothen the air-water interface by complementing the particles for the calculation in the vicinity of the water surface area where the particle supporting domain is truncated for the single-phase problem.  Figure 12a shows the contours of the density predicted using the two-phase SPH for the same case at the instant t = 16 s. Both air and water phases can be distinguished because of the large density difference between the two phases. In the vicinity of the wave crest, as indicated in Figure 12c, the entraining of air into the water creates a circulatory flow, indicating the continuity of velocity at the air-water interface.   Figure 12a shows the contours of the density predicted using the two-phase SPH for the same case at the instant t = 16 s. Both air and water phases can be distinguished because of the large density difference between the two phases. In the vicinity of the wave crest, as indicated in Figure 12c, the entraining of air into the water creates a circulatory flow, indicating the continuity of velocity at the air-water interface.

Two-Phase SPH Simulation of the Wave-Rigid Beam Interaction
The interaction of the waves created in the wave basin with a rigid cantilever beam fixed perpendicularly to the base of the wave basin and downstream of the wavemaker is addressed in this section. As shown in Figure 8, a rigid structure is set in the water, which is activated here to consider the WSI problem. The rigid cantilever beam of 0.2 m height and 0.03 m thickness is located at 2.8 m from the initial position of the piston wavemaker. The height of the beam can be changed to allow it to be submerged below the wavy interface or to protrude the wavy interface. A pressure sensor is placed at the height of 0.1 m (i.e., at half of the beam height). A second-order regular wave is generated using the same parameters and SPH settings and number of particles as outlined in Section 3.2. Figure 13 compares the temporal variation of the non-dimensional pressure P* = (P/(ρ 0 gd)) recorded by this pressure sensor computed using the single-phase and two-phase SPH simulations using the resolution R = 0.01 m developed in this study and also with that predicted using the single-phase mode in DualSPHysics using resolutions R = 0.01 m and 0.005 m.

Two-Phase SPH Simulation of the Wave-Rigid Beam Interaction
The interaction of the waves created in the wave basin with a rigid cantilever beam fixed perpendicularly to the base of the wave basin and downstream of the wavemaker is addressed in this section. As shown in Figure 8, a rigid structure is set in the water, which is activated here to consider the WSI problem. The rigid cantilever beam of 0.2 m height and 0.03 m thickness is located at 2.8 m from the initial position of the piston wavemaker. The height of the beam can be changed to allow it to be submerged below the wavy interface or to protrude the wavy interface. A pressure sensor is placed at the height of 0.1 m (i.e., at half of the beam height). A second-order regular wave is generated using the same parameters and SPH settings and number of particles as outlined in Section 3.2. Figure 13 compares the temporal variation of the non-dimensional pressure P* = (P/(ρ0 gd)) recorded by this pressure sensor computed using the single-phase and two-phase SPH simulations using the resolution R = 0.01 m developed in this study and also with that predicted using the single-phase mode in DualSPHysics using resolutions R = 0.01 m and 0.005 m. During the initial stages of the wave generation, the pressure recorded at the sensor appears erratic. After 3 s when the waves are fully developed, there appears to be some temporal quasi-periodicity of the non-dimensional pressure and with a period of about 1.3 s. There appears to be no significant differences between the pressure variations computed using the present single-phase and two-phase SPH simulations, indicating that the effect of the air phase on the submerged structure may be negligible.
Three representative time instants t1, t2, and t3 of 12.2 s, 12.6 s, and 13 s, respectively, (as shown in Figure 13) were selected to plot and compare the instantaneous contours of the non-dimensional pressure in the vicinity of the rigid beam structure computed using the single-and two-phase SPH codes from the present study which are shown in Figure  14. A pressure peak is achieved at t1 = 12.2 s when the wave crest passes over the structure's top, as shown in Figure 14a,b. After this, it can be seen that the pressure drops gradually as the wave crest passes over the structure at t2 = 12.6 s, and this can be seen in Figure  14c,d. Subsequently, as shown in Figure 14e,f, a temporary backflow can be seen to occur at t3 = 13 s, thereby resulting in a decay of the non-dimensional pressure to the trough. Since the wave elevation after passing the structure is significant for assessing the equivalence of the effectiveness of a barrier for say coastal protection, the temporal evolution of the free wavy surface at wave gauge WG3 computed using the single-and two-phase SPH from the present study are compared with the theoretical result (without the rigid beam structure) in Figure 15 with the same resolution. During the initial stages of the wave generation, the pressure recorded at the sensor appears erratic. After 3 s when the waves are fully developed, there appears to be some temporal quasi-periodicity of the non-dimensional pressure and with a period of about 1.3 s. There appears to be no significant differences between the pressure variations computed using the present single-phase and two-phase SPH simulations, indicating that the effect of the air phase on the submerged structure may be negligible.
Three representative time instants t 1 , t 2, and t 3 of 12.2 s, 12.6 s, and 13 s, respectively, (as shown in Figure 13) were selected to plot and compare the instantaneous contours of the non-dimensional pressure in the vicinity of the rigid beam structure computed using the single-and two-phase SPH codes from the present study which are shown in Figure 14. A pressure peak is achieved at t 1 = 12.2 s when the wave crest passes over the structure's top, as shown in Figure 14a,b. After this, it can be seen that the pressure drops gradually as the wave crest passes over the structure at t 2 = 12.6 s, and this can be seen in Figure 14c,d. Subsequently, as shown in Figure 14e,f, a temporary backflow can be seen to occur at t 3 = 13 s, thereby resulting in a decay of the non-dimensional pressure to the trough. Since the wave elevation after passing the structure is significant for assessing the equivalence of the effectiveness of a barrier for say coastal protection, the temporal evolution of the free wavy surface at wave gauge WG3 computed using the single-and two-phase SPH from the present study are compared with the theoretical result (without the rigid beam structure) in Figure 15 with the same resolution.

Figure 15.
Single-and two-phase SPH computed temporal evolution of the free surface at WG3 immediately after the wave has passed the structure. Figure 16 shows the impact of structural shape parameters such as thickness and height on the wave profile. After passing the structure, the wave appears to have a period of 1.3 s, consistent with that of the wave generated without the structure in the basin. The variation in the choice of structural parameters such as the height (d3) and thickness (S1) (different structural sizes) can be considered to assess the effectiveness of the structure for coastal protection. As shown in Figure 16, the wave height appears to decrease as the structure height (d3) is increased. Increasing the thickness, i.e., S1 = 0.12 m appears to have no impact on the temporal variation of the wave heights, probably because the beam is assumed rigid for this study, which may not be the case if the beam is assumed to be elastic. Figure 15. Single-and two-phase SPH computed temporal evolution of the free surface at WG3 immediately after the wave has passed the structure. Figure 16 shows the impact of structural shape parameters such as thickness and height on the wave profile. After passing the structure, the wave appears to have a period of 1.3 s, consistent with that of the wave generated without the structure in the basin. The variation in the choice of structural parameters such as the height (d 3 ) and thickness (S 1 ) (different structural sizes) can be considered to assess the effectiveness of the structure for coastal protection. As shown in Figure 16, the wave height appears to decrease as the structure height (d 3 ) is increased. Increasing the thickness, i.e., S 1 = 0.12 m appears to have no impact on the temporal variation of the wave heights, probably because the beam is assumed rigid for this study, which may not be the case if the beam is assumed to be elastic.
Energies 2022, 14, x FOR PEER REVIEW 19 of 28 Figure 16. Comparison of the present single-and two-phase SPH computed temporal evolution of the free surface at WG3 immediately after the wave has passed the structure with different geometrical sizes computed using SPH code and theory without the structure in the basin.

Two-Phase SPH Simulation of Wave Interaction with a Freely Floating Box in a Wave Basin
This section considers the hydrodynamic wave interaction with a freely floating rectangular box in the context of the present single-and two-phase SPH. This problem has already been experimented and simulated using single-phase SPH by Ren et al. [42] and Figure 16. Comparison of the present single-and two-phase SPH computed temporal evolution of the free surface at WG3 immediately after the wave has passed the structure with different geometrical sizes computed using SPH code and theory without the structure in the basin.

Two-Phase SPH Simulation of Wave Interaction with a Freely Floating Box in a Wave Basin
This section considers the hydrodynamic wave interaction with a freely floating rectangular box in the context of the present single-and two-phase SPH. This problem has already been experimented and simulated using single-phase SPH by Ren et al. [42] and Domínguez et al. [43] using DualSPHysics. Figure 17 shows a schematic of the test problem showing a freely floating box (Box A) with a length l b of 0.3 m, a height h b of 0.2 m, and unit width in the spanwise direction in the wave basin (Box B will be activated in the next section). The box with a density of 500 kg/m 3 is initially half-submerged in the water and then activated by a wave generated from the wavemaker with a period T of 1.2 s and a wave height H 0 of 0.1 m. The effect of this will be to cause the freely floating box to move in the direction (surge) and perpendicular to the direction (heave) of wave propagation and the rotation (pitch) concerning the mass center (counter-clockwise is positive). The wavemaker generates a second-order Stokes wave in the basin, and the waves are absorbed in the 1:4-sloped dissipative beach at the end of the basin. No-slip boundary conditions are implemented on the bottom wall and the dissipative beach. Buffers 1 and 2 denote the open boundaries treated as discussed in Sections 3.2 and 3.3. Other numerical parameters remain the same as those mentioned in Section 3.3 if these are not specified. The total number of particles used for single-phase SPH is 70023 and for the two-phase SPH is 32770. A resolution of 0.01 m suggested in [42] is used in the present study. The computed temporal variation of heave, surge, and pitch for the floating box using single-and two-phase SPH from the present study are shown in Figure 18, respectively, and are compared with the experimental data for the same problems from [42]. Comparing the results between single-phase and two-phase simulations, the observed slight differences indicate that the single-phase model is better than the two-phase model to reflect the two-phase wave-floating structure interaction problem when the air is still in the infinite field. Still, we find that the single-phase model induces a slightly greater fluctuation of the box than the two-phase one. Figure 19a,d shows the computed contours of the density and the non-dimensional pressure P = P/(ρ 0 gd)) using the two-phase SPH of the present study at two selected times instants of t 1 and t 2 of 6.1 s and 6.8, respectively, as indicated in Figure 18a. It can be seen that the air-water interface in the vicinity of the box is captured clearly, and the box bobs up as shown in Figure 19c with higher pressure at the bottom of the wave basin than when it bobs down as shown in Figure 19d.
OR PEER REVIEW 20 of 28 6.1 s and 6.8, respectively, as indicated in Figure 18a. It can be seen that the air-water interface in the vicinity of the box is captured clearly, and the box bobs up as shown in Figure 19c with higher pressure at the bottom of the wave basin than when it bobs down as shown in Figure 19d.

Wave Interaction with Two Freely Floating Boxes in a Basin
In this section, two freely floating boxes in the wave basin are investigated. This problem has not been fully resolved numerically due to the complexity of combining twophase flow, WSI, and the collisions between the two boxes. Here, we release a box (Box A) as in Section 3.4 and add an extra box (Box B) behind it (initially half-submerged in the water, see Figure 17), with an initial distance of d5 from the mass center of Box B to that of

Wave Interaction with Two Freely Floating Boxes in a Basin
In this section, two freely floating boxes in the wave basin are investigated. This problem has not been fully resolved numerically due to the complexity of combining twophase flow, WSI, and the collisions between the two boxes. Here, we release a box (Box A) as in Section 3.4 and add an extra box (Box B) behind it (initially half-submerged in the water, see Figure 17), with an initial distance of d 5 from the mass center of Box B to that of Box A. The two boxes maintain the length of 0.3 m, the height of 0.2 m, and the density of 500 kg/m 3 . The resolution of 0.01 m is adopted for our two-phase simulations. The other simulation setup remains the same, as indicated in Figure 17. By varying the location of Box B, two relative distances, i.e., d 5 = 0.5 m and 1.1 m, are initially given, corresponding to the motion of the two boxes as plotted in Figures 20 and 21. Note that the experimental data of a single box (Box A) in the reference [42], and single-phase simulation (two boxes) results obtained using DulSPHysics with a resolution of 0.007 m are plotted here for a better comparison. For Box A with a greater initial distance (1.1 m) to Box B, the effect of Box B on the motion of Box B can be neglected (compare to single box simulation results). On a smaller initial distance (0.5 m), Box A is hindered by the hydrodynamic interaction with Box B in the direction (surge) of wave propagation (see Figure 20 b). Regarding Box B, it is observed that the heaving amplitude is somewhat smaller than that of a single box (compare to the experiment data), as shown in Figure 21a, indicating that some wave energy is dissipated by Box A. It is observed that Box B with a smaller initial distance (0.5 m) to Box A presents a more significant surging displacement than that of the case with a greater one. This can be explained as the smaller initial distance contributes to a more vital hydrodynamic interaction in pushing Box B (Figure 21b) and hindering Box A (see Figure 20b) in the surging direction. Figure 22 presents the time evolution for the relative distances between the two boxes with different initial locations in which the final distances tend to remain stable periods. The differences between our two-phase simulation results and those single-phase results obtained using DualSPHysics may result from the existence of the air phase. greater one. This can be explained as the smaller initial distance contributes to a more vital hydrodynamic interaction in pushing Box B (Figure 21b) and hindering Box A (see Figure  20b) in the surging direction. Figure 22 presents the time evolution for the relative distances between the two boxes with different initial locations in which the final distances tend to remain stable periods. The differences between our two-phase simulation results and those single-phase results obtained using DualSPHysics may result from the existence of the air phase.      Since the initial distances aforementioned are too large to collide, a significantly smaller initial distance, i.e., d 5 = 0.32 m, is adopted for considering the collision between the two boxes. The collision force model adopted here has been addressed as outlined in Section 2.4. After the collision for the first time, the two boxes remain attached, corresponding to the coincidence of peaks in the surge direction as shown in Figure 23b. The variation of the relative distance and the motion for the two boxes reflects a typical interaction process in which the two boxes collide with each other, periodically taking the contact points as the fulcrum. This process is clearly described through the two snapshots at t 1 = 5.3 s and t 2 = 5.95 s (see Figure 24), as shown in Figure 25. Since the initial distances aforementioned are too large to collide, a significantly smaller initial distance, i.e., d5 = 0.32 m, is adopted for considering the collision between the two boxes. The collision force model adopted here has been addressed as outlined in Section 2.4. After the collision for the first time, the two boxes remain attached, corresponding to the coincidence of peaks in the surge direction as shown in Figure 23b. The variation of the relative distance and the motion for the two boxes reflects a typical interaction process in which the two boxes collide with each other, periodically taking the contact points as the fulcrum. This process is clearly described through the two snapshots at t1 = 5.3 s and t2 = 5.95 s (see Figure 24), as shown in Figure 25.

Conclusions
In this paper, an SPH method was used to simulate several two-phase WSI examples. We validate a two-phase dam-breaking and a second-order wave generated by a piston in a wave basin to capture the wavy air-water interface, and the interaction between the wave and a floating box. Subsequently, we first studied the interaction between the wave and a rigid cantilever beam, and the interaction between the wave and two identical floating boxes with the possibility of elastic collision between them.
For the dam-breaking case, a smoother cavity interface was observed due to the weak compressibility of the local air phase for the two-phase simulation. Without the air phase, the cavity collapses due to the pressure difference between its top and bottom. The pressure evolution of a sensor on the right wall using single-and two-phase SPH simulations shows agreement with the experimental data and the available numerical results, indicating that our numerical model can provide a reliable solution to this problem. Based on the

Conclusions
In this paper, an SPH method was used to simulate several two-phase WSI examples. We validate a two-phase dam-breaking and a second-order wave generated by a piston in a wave basin to capture the wavy air-water interface, and the interaction between the wave and a floating box. Subsequently, we first studied the interaction between the wave and a rigid cantilever beam, and the interaction between the wave and two identical floating boxes with the possibility of elastic collision between them.
For the dam-breaking case, a smoother cavity interface was observed due to the weak compressibility of the local air phase for the two-phase simulation. Without the air phase, the cavity collapses due to the pressure difference between its top and bottom. The pressure evolution of a sensor on the right wall using single-and two-phase SPH simulations shows agreement with the experimental data and the available numerical results, indicating that our numerical model can provide a reliable solution to this problem. Based on the

Conclusions
In this paper, an SPH method was used to simulate several two-phase WSI examples. We validate a two-phase dam-breaking and a second-order wave generated by a piston in a wave basin to capture the wavy air-water interface, and the interaction between the wave and a floating box. Subsequently, we first studied the interaction between the wave and a rigid cantilever beam, and the interaction between the wave and two identical floating boxes with the possibility of elastic collision between them.
For the dam-breaking case, a smoother cavity interface was observed due to the weak compressibility of the local air phase for the two-phase simulation. Without the air phase, the cavity collapses due to the pressure difference between its top and bottom. The pressure evolution of a sensor on the right wall using single-and two-phase SPH simulations shows agreement with the experimental data and the available numerical results, indicating that our numerical model can provide a reliable solution to this problem. Based on the Biesel transfer function theory, the two-phase and second-order regular wave is generated by moving a piston in a wave basin and implementing the open boundary condition at the top and left of the domain of the air phase. The obtained wave elevation evolution and the velocity components using the single-and two-phase SPH at the wave gauges match the theoretical solution and those obtained by the single-phase DualSPHysics to a reasonable extent. Even though the same velocity patterns and free-surface profiles were observed for the one-and two-phase simulation in a global sight, the air phase can smooth the water-air interface by complementing the particles for the calculation in the vicinity of the water surface area where the particle supporting domain is truncated for the single-phase problem. The effect of different structural sizes on the wave elevation was investigated, and it was found that the increase of structure height decreases the wave elevation. On the motion of a floating box in the presence of the wave, our results show an agreement with the available experimental data, and only slight differences have been observed between the single-and two-phase simulations. On the motion of two floating boxes in the wave, it is observed that Box B with a smaller initial distance (0.5 m) to Box A presents a more significant surging displacement than that of the case with a greater one, which can be explained as the smaller initial distance contributes to a stronger hydrodynamic interaction in pushing Box B and hindering Box A in the surging direction. For the case of collision between the two boxes, we find that the two boxes collide periodically with taking the contact points as the fulcrum. Funding: This manuscript received no external funding.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.