Two-Dimensional Free-Surface Flow Modeling for Wave-Structure Interactions and Induced Motions of Floating Bodies

: In this study, the level set (LS) and immersed boundary (IB) methods were integrated into a Navier–Stokes equation two-phase ﬂow solver, to investigate wave-structure interactions and induced motions of ﬂoating bodies in two dimensions. The movement of an interfacial boundary of two ﬂuids, even with severe free-surface deformation, is tracked by using the level set method, while an immersed object inside a ﬂuid domain is treated by the IB method. Both approaches can be implemented by solving the Navier–Stokes equations for viscous laminar ﬂows with embedded objects in ﬂuids. For accurate treatment of the solid–liquid phase, appropriate source terms as forcing functions to take into account the hydrodynamic e ﬀ ects on the body boundaries are added into the governing equations. The integrated compact interfacial tracking techniques between the interfaces of gas–liquid phase and the solid–liquid phase allow the use of a combined Eulerian Cartesian and Lagrangian grid system. Problems related to the ﬂuid-structure interactions and induced motions of a ﬂoating body, such as (a) a dam-break wave over a dry bed; (b) a dam-break wave over either a submerged semicircular or rectangular cylinder; (c) wave decomposition process over a trapezoid breakwater; (d) a free-falling wedge into a water body; and (e) wave packet interacting with a ﬂoating body are selected to test the model performance. For all cases, the computed results are found to agree reasonably well with published experimental data and numerical solutions. For the case of modeling wave decomposition process, improved solutions are obtained. The detailed features of ﬂow phenomena described by the physical variables of velocity, pressure and vorticity are presented and discussed. The present two-phase ﬂow model is proved to have the advantage of simulating the cases with induced severe interfacial oscillations and coupled gas (or air) motions where the single-phase model may miss the contributions of the air motions on the interfaces. Additionally, the proposed method with uses of the LS and IB methods is demonstrated to be able to achieve the reliable predictions of complex ﬂow ﬁelds. the instantaneous vorticity contours of two moving cylinders for Re = 100 at four selected instants. The color-coded vorticity contours range from − 5 to 5. In this transient simulation study, the results of hydrodynamic e ﬀ ect and vortex shedding phenomena under the condition of two cylinders moving against each other in viscous ﬂuid are found to be completely di ﬀ erent when compared to those for a single moving body at the same Reynolds number. Due to the viscous e ﬀ ect within the boundary layers of the solid bodies, the moving rigid bodies are noticed to play an important role in generating the vortex ﬂow in the immersed viscous ﬂuid domain.


Introduction
Obtaining solutions of a hydrodynamic model for the studies of complex two-phase flows, such as the propagation of surge fronts, movement of free-surface by a moving body, or waves pass over either completely or partially immersed structures, is an important subject in coastal and ocean engineering applications [1,2]. The key issue of modeling free-surface flows is the treatment of moving grids assigned in the domain. Figure 1 depicts a schematic diagram showing a fixed Cartesian grid system that is used for solving the Navier-Stokes equations. Under the assumption of incompressible fluid medium and considering the effects of surface tension between two fluids and an immersed body, the governing equations for two-phase flows can be formulated as follows: (2) detailed descriptions of the physical variables given in Equation (2) can be found in Sethian and Smereka [37]. For the level set function, ϕ , t , typically, when ϕ 0, it denotes the gas region, while ϕ 0 is for the liquid region. The fluid properties, such as the density (ρ and viscosity μ can be interpolated according to the following: ρ ϕ ρ H ϕ ρ 1 H ϕ and μ ϕ μ H ϕ μ 1 H ϕ (3) where H ϕ is the Heaviside function [38].
H ϕ 0 1 2 1 ϕ ε 1 π sin πϕ ε 1 ϕ ε |ϕ| ε ϕ ε (4) The interfacial thickness of the Heaviside function is approximately 2ε, where ε is a small parameter that is related to the order of the size of the mesh cells close to the interface and can be selected through the numerical tests. The density (ρ and viscosity (μ with the subscripts l and g represent, respectively, the corresponding properties of liquid and gas. Here, the computational domain contains two different fluids, and δ is the Dirac delta function. σ is the coefficient of surface tension, k is the curvature of an interface, f denotes hydrodynamic forces induced by an immersed body, A b as defined later represents the wave absorbing coefficient, g is the gravitational acceleration vector, t is time, and φ denotes the level set function. The 2D velocity vector, u, describes the velocity components in x and z directions, and p is the pressure. More detailed descriptions of the physical variables given in Equation (2) can be found in Sethian and Smereka [37].
For the level set function, φ(x, t), typically, when φ < 0, it denotes the gas region, while φ > 0 is for the liquid region. The fluid properties, such as the density (ρ) and viscosity (µ) can be interpolated according to the following: where H(φ) is the Heaviside function [38].
The interfacial thickness of the Heaviside function is approximately 2ε, where ε is a small parameter that is related to the order of the size of the mesh cells close to the interface and can be selected through the numerical tests. The density (ρ) and viscosity (µ) with the subscripts l and g represent, respectively, the corresponding properties of liquid and gas. The level set function (φ) describing the interface that moves with the fluid particles is governed by the transport equation: ∂φ ∂t + u·∇φ = 0 (5) For the numerical accuracy, a fourth-order Runge-Kutta integration in time and WENO [32] scheme in space are adopted to discretize Equation (5).
In Equation (2), the smoothed delta function can be written as follows: Moreover, the curvature k(φ) can be computed as follows: With the consideration of a wave-absorbing region, the proposed absorbing coefficient (A b ) by Lin and Liu [39] is adopted: (8) where A b denotes the absorbing coefficient, and x xt and x ab are the starting position and length of the absorbing region. C α and n c are the empirical damping related coefficients. Here, the coefficients of C α = 200 and n c = 10, as recommended by [39], are used in numerical simulations.

Numerical Formulations
The governing equations (Equations (1) and (2)) are solved by using the projection method [40], where the predicted intermediate velocity field with the known values obtained at the n (previous) time level are firstly calculated. According to the two-step calculation procedure and Adams-Bashforth time-stepping approach, the solutions of velocity vector (u * ) at each intermediate (between n and n + 1) time level are determined by solving Equation (2) with the expressions of temporal difference as follows:û − u n ∆t = RHS (9) u * − u n ∆t = RHS + f û (10) where ∆t is the time step, and u n is the velocity vector at the n time level. The terms with the effects of convection, diffusion, surface tension, wave absorbing, and gravity are combined with a notation of RHS = 1.5A n − 0.5A n−1 : In Equation (10), f is an added term describing force density induced by the immersed bodies. Its calculation procedure, based on the IB method, is provided in the following section. With a flow passing over a rigid body, the calculation domain is separated by two phases, where one is the ordinary fluid flow surrounding the rigid body, and the other is the body phase with an inclusion of a body forcing. As can be seen in Equation (10), the velocity vectors,û, are first computed without considering the effect of immersed boundary. Then, the forcing term, f, is computed by a continuous forcing IB method with the computed velocity vectors,û. By imposing a body force, f, in those cells, which has the property of a signed distance function near the immersed boundary, the intermediate velocity field is updated from Equation (10).
Taking the divergence of Equations (1) and (2), the pressure field can be obtained by solving the Poisson equation, which is derived by satisfying the continuity equation: The pressure Neumann boundary conditions should be satisfied. We have the following: Once the velocity vectors at the intermediate time level are determined, Equation (12) can be solved, using a spline alternating direction implicit method (SADI) with a standard tri-diagonal matrix solver applied for the calculations of the pressure field. Here, the SADI method is used in the z direction. The iteration consists of a sweep through the i index for each (i, k) pair, solving implicitly for the pressures in the k-index direction, where i and k represent, respectively, the grid indices along x and z directions. Then, the velocity vectors at the new (n + 1) time level can be obtained by solving the following equation: In terms of the evaluation of the level set function (φ) for capturing the positions of the two-phase interfaces, Equation (5), when expressed in Lagrangian description, is given as follows: where L(φ) in Equation (15) is an approximation (e.g., WENO approximation in these lecture notes). Using a uniform grid in one space dimension as an example, Equation (15) results in an ordinary differential equation (ODE) with a discrete conservative finite difference scheme. Then, the fourth-order Runge-Kutta scheme, as shown below, is used to advance the interfaces: where φ 0 is the coordinates of the interfacial points from the previous time step, and k 1 , k 2 , k 3 , and k 4 are obtained as follows: Here, the computed u n+1 are used as the advection needed velocity field for the calculation of the level set function.

Immersed Boundary (IB) Method
The IB method basically follows the concept that the effect of an immersed structure on its surrounding flow is modeled through a force density vector (f), which is included in the Navier-Stokes equations. Figure 2 depicts an example plot showing the arrangement of the Eulerian and Lagrangian grid points used for transferring the values of variables between them through an interpolation procedure with a weighted sum of regularized Dirac delta function, δ h . Here, the body boundary (fluid-solid interface) is distributed with discrete surface markers (Lagrangian grid points). The force density vector (f) on the right-hand side of Equation (2) can be formulated as follows: where F(X l , t) is the force density vector at the Lagrangian points with coordinates X l along the body boundary. The F(X l , t) values at Lagrangian points are computed beforehand and then transferred into the Eulerian meshes for f(x, t), using Equation (21). The discrete Dirac delta function δ h developed by Peskin [23] is adopted in this study. The expression of δ h for a 2D case is given as follows: where δ 1D h (x 1 − X 1 ) = 1 h ϕ(r 1 ) is the regularized one-dimensional delta function. The formulations of ϕ(r), proposed both by Roma et al. [25] and Yang et al. [29], are as follows, respectively:  where , t is the force density vector at the Lagrangian points with coordinates along the body boundary. The , t values at Lagrangian points are computed beforehand and then transferred into the Eulerian meshes for , , using Equation (21). The discrete Dirac delta function δ developed by Peskin [23] is adopted in this study. The expression of δ for a 2D case is given as follows: where δ x X φ r is the regularized one-dimensional delta function. The formulations of φ r , proposed both by Roma et al. [25] and Yang et al. [29], are as follows, respectively: These equations are used together with the selection of the influential radius, r, to be either 1.5∆h (Equation (23)) or 2∆h (Equation (24)) for the numerical sensitivity study. The force density vector, F(X l , t), at the interfacial Lagrangian marker points is calculated based on the time rate change of velocity between the body velocity (U b ) and the Lagrangian point velocity (U Γ ): where the velocity vector (U b ) of an immersed body is determined from the equations of motion. Accordingly, the force density vector (f) at the Eulerian grid points can be calculated by using a spreading scheme involving the Lagrangian marker points as follows: where ∆V l is the control volume defined about the m-th Lagrangian marker. The resulting Eulerian force (f) from Equation (27) is then substituted into Equation (10) for solution computation.

Rigid Body Dynamics
The movements of a body subject to the external forces and moments are governed by the equations of motion for linear and angular momentum of a rigid body. In this vertical two-dimensional modeling study, the equations of motion for a body subject to external forcing can be formulated by considering the effects of inertial motion and internal mass. According to [41], the equations describing the translational and rotational motions of a moving body are expressed as follows: Here, V b , I b , and ρ b are volume, moment of inertia, and density of the body, respectively; ρ is the fluid density; r represents the position vector relative to the body centroid, g is the vector of gravitational acceleration, and the symbol × denotes the cross product. U b and ω b are respectively the translational and angular velocity vectors of the rigid body. In the present study, the Euler forward difference method is used to solve Equations (28) and (29) for U b and ω b . The velocity vector at the surface of the rigid body can be calculated as follows:

Results and Discussions
Using the developed 2D two-phase flow model that integrates the approaches of level set (LS), immersed boundary (IB), and structural responses, selected fluid-structure interaction problems are investigated. The study cases include the following: (a) a uniform flow passing through a circular cylinder to verify the adopted IB methodology, (b) two cylinders moving against each other to test the moving boundary effect, (c) dam-break flows to examine the free-surface tracking capability, (d) dam-break wave flowing over a submerged structure, (e) wave decomposition process over a submerged trapezoid breakwater, (f) a free-falling wedge, and (g) wave packet interacting with a floating body. It should be noted the simulations carried out involve the interfacial conditions specified at the gas-liquid phase and the liquid-solid phase.

Flow Passing through a Cylinder
In order to verify the present IB method, cases with steady uniform flows passing through a circular cylinder are simulated for validation. The flow results are obtained at Re = 20, 40, 80, 100, and 200, where Re is the Reynolds number defined as Re = UD/υ. Here, υ denotes the kinematics viscosity of a fluid, U is the uniform flow velocity, and D denotes the diameter of the cylinder. In the test cases, a cylinder with a diameter of one unit is embedded in a 32 × 16 computational domain. Upstream flow velocity of one unit is assumed to be uniformly distributed in vertical direction. The accuracies of the Water 2020, 12, 543 9 of 29 predicted drag and lift coefficients, and other physical variables for the selected cases, are evaluated. For the flow conditions of Re = 20 and 40, the computed reattachment length, Lw, and drag coefficient, C D , as given in Table 1, are compared with the published experimental [42] and numerical [43][44][45][46][47] results. Additionally, to ensure the grid independence, included in Table 1 are the present results for the separately used mesh sizes of ∆h = 1/30 and ∆h = 1/40. As can be seen in Table 1 for the cases of Re = 20 and 40, the converged solutions in L w and C D are obtained, and they are in good agreements with other published results. By refining the computational grids, it is shown to be able to improve the solutions for the flow conditions at a larger Reynolds number, such as Re = 100 or 200; this confirms the global accuracy of the present method. Comparisons between the present and other numerical results [43,[46][47][48][49][50] of the drag and lift coefficients (C D and C L ), and Strouhal number (St) for the cases of Re = 100 and 200 are summarized in Table 2. The present results were computed by using a refined grid size, i.e., ∆h = 1/50. It can be seen the present results are shown to have very good agreements with others given in the literature. For the cases of Re = 100 and 200, the transient solution procedure was carried out to reach the prescribed time, i.e., t = 200. Time variations of the dimensionless drag and lift coefficients at Re = 200 are presented in Figure 3. Due to the induced vortex motions, the transition processes until the C D and C L approach to the stabilized periodic variations are clearly depicted in Figure 3. Figure 4 exhibits the instantaneous distribution patterns of pressure and vorticity at t = 200 for various Reynolds numbers of 40, 80, 100, and 200. The patterns of vortex shedding for the cases of Re = 80, 100, and 200 can be clearly observed. The predicted results of vortex shedding demonstrate the proposed method can predict accurately the transient flow patterns. As can be seen, the wake triggered at the cylinder surface is asymmetric in reference to the direction of flow. The formation of vortices on either side of the flow direction does not simultaneously occur. In fact, the shedding of vortices alternates from side to side, which had been observed in experiments. Thus, the pressure distribution around the obstacle is asymmetric about the flow direction. patterns. As can be seen, the wake triggered at the cylinder surface is asymmetric in reference to the direction of flow. The formation of vortices on either side of the flow direction does not simultaneously occur. In fact, the shedding of vortices alternates from side to side, which had been observed in experiments. Thus, the pressure distribution around the obstacle is asymmetric about the flow direction.

Two Cylinders Moving Against Each Other in Viscous Fluid
The present model is further tested to examine the hydrodynamic effect induced by two cylinders moving in opposite directions, in a domain of viscous fluid. Again, each of the two moving cylinders has a diameter of 1 unit, and the computational domain is 32 × 16. The constant velocities of the upper and lower moving cylinders are set respectively as −1 and 1 units. As shown in Figure

Two Cylinders Moving Against Each Other in Viscous Fluid
The present model is further tested to examine the hydrodynamic effect induced by two cylinders moving in opposite directions, in a domain of viscous fluid. Again, each of the two moving cylinders has a diameter of 1 unit, and the computational domain is 32 × 16. The constant velocities of the upper and lower moving cylinders are set respectively as −1 and 1 units. As shown in Figure 5, the centers of the two cylinders are initially placed at the locations of (0, − 0.75) and (16, 0.75), separately, within the domain ranges of −8 ≤ x ≤ 24 and −8 ≤ z ≤ 8. are shown in Figure 7. From Figure 7, it is noted that, for all cases, the minimum values of the drag coefficient occur roughly at t̃ 16, where t̃ = Ut/(D/2); meanwhile, for the lift coefficient, the minimum values take place at a delayed time of around t̃ = 18. In general, the results show that the drag coefficient and the positive portion of the lift coefficient decrease with an increase in the Reynolds number. The transitions of the drag and lift coefficients become particularly pronounced at the timeframe from t̃ = 14 to t̃ = 20. In addition, to demonstrate that the present model is capable of successfully simulating the flow fields, Figure 8 depicts the instantaneous vorticity contours of two moving cylinders for Re = 100 at four selected instants. The color-coded vorticity contours range from −5 to 5. In this transient simulation study, the results of hydrodynamic effect and vortex shedding phenomena under the condition of two cylinders moving against each other in viscous fluid are found to be completely different when compared to those for a single moving body at the same Reynolds number. Due to the viscous effect within the boundary layers of the solid bodies, the moving rigid bodies are noticed to play an important role in generating the vortex flow in the immersed viscous fluid domain.  To test the sensitivity of the grid size on the model solutions for this moving boundary problem, the set computational grid systems include 640 × 320, 960 × 480, and 1280 × 640. The time step was selected to satisfy the Courant-Friedrichs-Lewy (CFL) condition [51]. The CFL for n-dimensional domain can be generalized as C = ∆t is the grid size in the ith dimension, which can be a varied value in a given dimension, u x i denotes the magnitude of velocity in the ith dimension (A two-dimensional case is used here), and ∆t represents the fractional time-step. C max is typically set to be equal to 1. Figure 6 presents the time variations of the dimensionless lift coefficient for the upper cylinder at Re = 40. The results from Russel and Wang [44] are also included in Figure 6, for comparisons. Russel and Wang [44] presented an efficient method based on an underlying Cartesian grid for solving 2D incompressible viscous flows around multiple moving objects. It can be seen the present results agree well with those from [44], especially under the grid system of 1280 × 640 (i.e., ∆h = 1/40). The time varying drag coefficient (C D ) and lift coefficient (C L ) for the upper cylinder under the conditions of various Reynolds numbers are shown in Figure 7. From Figure 7, it is noted that, for all cases, the minimum values of the drag coefficient occur roughly at t = 16, where t = Ut/(D/2); meanwhile, for the lift coefficient, the minimum values take place at a delayed time of around t = 18. In general, the results show that the drag coefficient and the positive portion of the lift coefficient decrease with an increase in the Reynolds number. The transitions of the drag and lift coefficients become particularly pronounced at the timeframe from t = 14 to t = 20. In addition, to demonstrate that the present model is capable of successfully simulating the flow fields, Figure 8 depicts the instantaneous vorticity contours of two moving cylinders for Re = 100 at four selected instants. The color-coded vorticity contours range from −5 to 5. In this transient simulation study, the results of hydrodynamic effect and vortex shedding phenomena under the condition of two cylinders moving against each other in viscous fluid are found to be completely different when compared to those for a single moving body at the same Reynolds number. Due to the viscous effect within the boundary layers of the solid bodies, the moving rigid bodies are noticed to play an important role in generating the vortex flow in the immersed viscous fluid domain.

Dam-Break Problems
Simulations of selected dam-break problems are used to verify the free-surface capturing approach. The first example computes the position of the leading edge of a dam-break wave and compares the results with the experimental data of Martin and Moyce [52] and numerical results from Nomeritae et al. [53]. The other cases include a dam-break wave propagating over a submerged semicircular, or a rectangular, cylinder.
Following Martin and Moyce's [52] experimental study, the initial condition considered is a 0.057 m square water column (or H = 0.057 m) arranged behind a dam face in a 3 m long and 1 m high container. At t > 0, the dam is removed to allow the water to freely move toward downstream. Three uniform grid systems with 150 × 10 (x = z = 0.02 m), 300 × 20 (x = z = 0.01 m), and 600 × 40 grids (x = z = 0.005 m) and time step (t) of 5 × 10 −4 sec were used to simulate the dam-break flows for the determination of the nondimensionalized leading-edge position. Figure 9 shows the comparisons between the present solutions on the time varying leading-edge position and those from measurements [52] and other numerical results [53]. For the Martin and Moyce [52] dam-break flow case, the corresponding Reynolds number was about 42,792. The computed surge front positions versus the nondimensionalized time reveal a good agreement with the data from Martin and Moyce [52] and the numerical results from Nomeritae et al. [53]. Moreover, the tests of the grid size effect indicate that the converged and nicely fit solutions are obtained when the refined grid size of x = z = 0.005 m is used.

Dam-Break Problems
Simulations of selected dam-break problems are used to verify the free-surface capturing approach. The first example computes the position of the leading edge of a dam-break wave and compares the results with the experimental data of Martin and Moyce [52] and numerical results from Nomeritae et al. [53]. The other cases include a dam-break wave propagating over a submerged semicircular, or a rectangular, cylinder.
Following Martin and Moyce's [52] experimental study, the initial condition considered is a 0.057 m square water column (or H = 0.057 m) arranged behind a dam face in a 3 m long and 1 m high container. At t > 0, the dam is removed to allow the water to freely move toward downstream. Three uniform grid systems with 150 × 10 (∆x = ∆z = 0.02 m), 300 × 20 (∆x = ∆z = 0.01 m), and 600 × 40 grids (∆x = ∆z = 0.005 m) and time step (∆t) of 5 × 10 −4 sec were used to simulate the dam-break flows for the determination of the nondimensionalized leading-edge position. Figure 9 shows the comparisons between the present solutions on the time varying leading-edge position and those from measurements [52] and other numerical results [53]. For the Martin and Moyce [52] dam-break flow case, the corresponding Reynolds number was about 42,792. The computed surge front positions versus the nondimensionalized time reveal a good agreement with the data from Martin and Moyce [52] and the numerical results from Nomeritae et al. [53]. Moreover, the tests of the grid size effect indicate that the converged and nicely fit solutions are obtained when the refined grid size of ∆x = ∆z = 0.005 m is used. Water 2020, 12, x FOR PEER REVIEW 14 of 29 The present two-phase flow model is also simulated for the cases with dam-break waves over a submerged structure that is placed in an enclosed channel downstream of a dam face. Illustrated in Figure Figure 11 presents the comparisons between the present model results (velocity in x direction and pressure) and those from Flow-3D for the case of a submerged semicircular cylinder at point P1. The time variations of both the velocity and pressure, except at the time near the end of the interaction process, are in good agreements with the Flow-3D results. The results obtained from the arrangements of two different grid sizes reveal the consistent and reliable model predications. The comparisons of time varying pressures at points P2 and P3 under the cases of either a semicircular or a rectangular cylinder are presented in Figure 12. It is found the variation trends of the pressure at P2 and P3 under the influence of a submerged semicircular cylinder are similar to those under the case of a submerged rectangular cylinder. However, due to a more severe geometry change that appeared in the rectangular cylinder, the pressures at P2 and P3 under the case of a submerged rectangular cylinder are greater than those with a semicircular cylinder after the surge waves impact on the cylinders. The present two-phase flow model is also simulated for the cases with dam-break waves over a submerged structure that is placed in an enclosed channel downstream of a dam face. Illustrated in Figure 10  Two different grid sizes, ∆x = ∆z = ∆h = 0.025 m and ∆x = ∆z = ∆h = 0.033 m, and ∆t = 0.005 s are adopted for model simulations. Here, the interfacial parameter (ε) and the influential radius (r) of the Lagrangian marker points of the IB are, respectively, 2∆h and 1.5∆h. The time varying velocities and pressures at three selected locations, namely P 1 at (2 m, 0.16 m), P 2 at (4 m,0.16 m), and P 3 at (5 m, 0.26 m) (see Figure 10), are compared with the Flow-3D model results. The results obtained from the Flow-3D was based on the two-phase laminar flow model with the consideration of water and air phases. Figure 11 presents the comparisons between the present model results (velocity in x direction and pressure) and those from Flow-3D for the case of a submerged semicircular cylinder at point P 1 . The time variations of both the velocity and pressure, except at the time near the end of the interaction process, are in good agreements with the Flow-3D results. The results obtained from the arrangements of two different grid sizes reveal the consistent and reliable model predications. The comparisons of time varying pressures at points P 2 and P 3 under the cases of either a semicircular or a rectangular cylinder are presented in Figure 12. It is found the variation trends of the pressure at P 2 and P 3 under the influence of a submerged semicircular cylinder are similar to those under the case of a submerged rectangular cylinder. However, due to a more severe geometry change that appeared in the rectangular cylinder, the pressures at P 2 and P 3 under the case of a submerged rectangular cylinder are greater than those with a semicircular cylinder after the surge waves impact on the cylinders.    In terms of force computations from IB approach, the time variations of drag forces acting respectively on a semicircular or a rectangular cylinder are presented in Figure 13. As a model verification, the present results for the case of a submerged rectangular cylinder are also compared with those obtained from the Flow-3D in Figure 13b. As indicated in Figure 13a, the maximum force induced by the leading surge wave on a semicircular cylinder (600 N) is less than that on a rectangular cylinder (850 N). For the case of a dam-break wave over a rectangular cylinder, as shown in Figure 13b, the predicted maximum in-line force from IB calculations is found to be slightly greater than that from the Flow-3D model. In addition, although the variation pattern is similar, the present model predicts higher values of the forces caused by the follow-up waves due to the slightly higher water-surface level predicted (See Figure 14). Regarding the case of a dam-break wave over a submerged rectangular cylinder, the sensitivity tests for force computation, using thee different input conditions, namely ε = ∆h, r = 1.5∆h; ε = 2∆h, r = 1.5∆h; and ε = 2∆h, r = 2∆h, are presented in Figure 13b. Apparently, the difference of forces computed for cases of ε = 2∆h, r = 1.5∆h and ε = 2∆h, r= 2∆h is insignificant. In other words, the influential radius has reached its required converged value on the numerical simulation. However, the force predications have a noticeable difference when using two different interfacial thickness, i.e., ε = ∆h and ε = 2∆h, between two fluids. From results shown in Figure 13, it is evident that the values of hydrodynamic force are highly affected by the shape of the cylinder. Due to its smooth shape, the semicircular cylinder is subject to a relatively smaller hydrodynamic force than that to a rectangular cylinder. Figure 14 presents the spatial variations of the free-surface elevations at t = 0.5 s and t = 2.0 s, with the results obtained from both the present model and Flow-3D. Good agreement is noticed from the comparison plots. Due to the presence of a submerged structure, the predicted higher nonlinear growth of the free-surface flows has resulted in a more complex flow field with steep velocity gradient within the region of boundary layer. In terms of force computations from IB approach, the time variations of drag forces acting respectively on a semicircular or a rectangular cylinder are presented in Figure 13. As a model verification, the present results for the case of a submerged rectangular cylinder are also compared with those obtained from the Flow-3D in Figure 13b. As indicated in Figure 13a, the maximum force induced by the leading surge wave on a semicircular cylinder (600 N) is less than that on a rectangular cylinder (850 N). For the case of a dam-break wave over a rectangular cylinder, as shown in Figure  13b, the predicted maximum in-line force from IB calculations is found to be slightly greater than that from the Flow-3D model. In addition, although the variation pattern is similar, the present model predicts higher values of the forces caused by the follow-up waves due to the slightly higher watersurface level predicted (See Figure 14). Regarding the case of a dam-break wave over a submerged rectangular cylinder, the sensitivity tests for force computation, using thee different input conditions, namely ε = h, r = 1.5h; ε = 2h, r = 1.5h; and ε = 2h, r = 2h, are presented in Figure 13b. Apparently, the difference of forces computed for cases of ε = 2h, r = 1.5h and ε = 2h, r= 2h is insignificant. In other words, the influential radius has reached its required converged value on the numerical simulation. However, the force predications have a noticeable difference when using two different interfacial thickness, i.e., ε = h and ε = 2h, between two fluids. From results shown in Figure 13, it is evident that the values of hydrodynamic force are highly affected by the shape of the cylinder. Due to its smooth shape, the semicircular cylinder is subject to a relatively smaller hydrodynamic force than that to a rectangular cylinder. Figure 14 presents the spatial variations of the free-surface elevations at t = 0.5 s and t = 2.0 s, with the results obtained from both the present model and Flow-3D. Good agreement is noticed from the comparison plots. Due to the presence of a submerged structure, the predicted higher nonlinear growth of the free-surface flows has resulted in a more complex flow field with steep velocity gradient within the region of boundary layer.  In order to compare the flow patterns induced by either a submerged semicircular cylinder or a rectangular cylinder, the evolutions of density distributions of showing the water surface levels at various times are illustrated in Figure 15. Notably, a very similar variation trend of the free-surface elevation in cases with two tested submerged cylinders can be observed. Figure 16 presents the velocity vectors of a dam-break wave passing over the submerged structures of either a semicircular cylinder or a rectangular cylinder. The two-phase flow patterns between air and water phases are plotted particularly to show the motion of the vortices, which are generated during the process of wave-structure interaction. As seen in Figure 16, the plots also illustrate the local velocity distribution near the free-surface and around the submerged structures. The results obtained from the present two-phase flow model indicate that the shape of submerged rigid bodies has a noticeable effect on the vortex motion of the gas-phase fluid above the interface of the advancing waves. It should be noted that the results presented above for the dam-break wave problems were obtained from the simulations of the present two-phase flow model, involving the interfacial handing techniques between the air-water phase and liquid-solid phase that allow uniquely the use of a Cartesiancoordinate based mesh system. Though the single-phase flow model can be used to resolve the freesurface boundary condition of the Navier-Stokes equations, the numerical procedure adds further complexity for modeling free-surface flow problems. In order to compare the flow patterns induced by either a submerged semicircular cylinder or a rectangular cylinder, the evolutions of density distributions of showing the water surface levels at various times are illustrated in Figure 15. Notably, a very similar variation trend of the free-surface elevation in cases with two tested submerged cylinders can be observed. Figure 16 presents the velocity vectors of a dam-break wave passing over the submerged structures of either a semicircular cylinder or a rectangular cylinder. The two-phase flow patterns between air and water phases are plotted particularly to show the motion of the vortices, which are generated during the process of wave-structure interaction. As seen in Figure 16, the plots also illustrate the local velocity distribution near the free-surface and around the submerged structures. The results obtained from the present two-phase flow model indicate that the shape of submerged rigid bodies has a noticeable effect on the vortex motion of the gas-phase fluid above the interface of the advancing waves. It should be noted that the results presented above for the dam-break wave problems were obtained from the simulations of the present two-phase flow model, involving the interfacial handing techniques between the air-water phase and liquid-solid phase that allow uniquely the use of a Cartesian-coordinate based mesh system. Though the single-phase flow model can be used to resolve the free-surface boundary condition of the Navier-Stokes equations, the numerical procedure adds further complexity for modeling free-surface flow problems.

Wave Decomposition Process over a Submerged Trapezoid Breakwater
The model application is further extended to simulate wave generation, propagation, and interacting with a submerged trapezoidal bar. This case has been experimentally investigated by Beji and Battjes [54], where the measurements were performed in a 36 m long flume with a set water depth of 0.4 m. The numerical simulation was performed in a computational flume (or domain), which is 36 m long and 0.6 m high, and the water depth is 0.4 m. The spatial domain in this numerical wave tank (NWT) is divided into three zones: the wave generation zone, the wave propagation region, and the wave-absorbing layer. The length of wave generation zone is one wavelength long, and the center of the zone is located at x = 10.5 m. The numerical wave-absorbing layers are located in both the upstream and downstream domains with the ranges of 0 < x < 8 m and 28 < x < 36 m, respectively. The length of bar base is 11 m (between x = 14 m and x = 25 m), and its height is 0.3 m. The slopes of the inclined front and rear faces of the trapezoid breakwater are, respectively, 1/20 and 1/10.
The simulations start with a state of rest in water body subject to hydrostatic pressure. The free-surface elevation and velocity component for the incoming Stokes waves are specified at the inlet of the upstream boundary as the inflow conditions [55]. Separately, a total of 360,000 and 450,000 grid points together with ∆t = 0.002 s were used in the numerical calculations. The mesh sizes ∆x and ∆z for 360,000 grid points are, respectively, 0.01 and 0.005 m, whereas, for cases using 450,000 grid points, they are ∆x = 0.01 m and ∆z = 0.004 m. In the experimental study by Beji and Battjes [54], an incident wave train with a wave height, H = 2 cm, a wave period of T = 2.02 s, and a wavelength of 3.73 m was generated to propagate over a submerged trapezoidal bar. The free-surface elevations computed with the same input conditions are compared with the measurements to demonstrate the performance of the developed model.
Presented in Figure 17 Figure 17 exhibit the temporal free-surface evolution and the wave transformation processes due to the shoaling effect along the region of non-constant water depth. Because the decrease of the water depth along the upward slope of the bar, the wave elevation is shown to deform with a steepened wave amplitude and sawtoothed-like profile at x = 20.5 m and 21.5 m (See Figure 17b,c). Additionally, as shown in Figure 17d, high and sharp wave crests are formed over the bar crest at x = 22.5 m due to the wave shoaling and nonlinear interaction between waves and a submerged bar. Figure 17e reveals the decomposition of the wave with the development of higher harmonic components at regions from x = 23.5 m onward, as the wave propagates over and past the back face of the bar. It is noticed that the wave elevations obtained from the present model at the six-gauge locations match nicely with the experimental measurements [54] and the Flow-3D results. The simulated results of Flow-3D are based on a single phase with either laminar or turbulent (k − ε) flow model. The comparisons of the present and Flow-3D wave elevations indicate that both results are well predicted when compared with the measurements at the gauge locations during the period from t = 33 s to t = 39 s. However, relatively, the results from the present model are shown to have slightly better agreements with the experimental measurements than those from the Flow-3D, considering either the laminar or the turbulent (k − ε) flow module. It again demonstrates that the present two-phase flow model that combines LS and IB methods can be implemented with reasonable predictions for the study of hydrodynamic interactions between waves and structures.
gauge locations during the period from t = 33 s to t = 39 s. However, relatively, the results from the present model are shown to have slightly better agreements with the experimental measurements than those from the Flow-3D, considering either the laminar or the turbulent (k ε) flow module. It again demonstrates that the present two-phase flow model that combines LS and IB methods can be implemented with reasonable predictions for the study of hydrodynamic interactions between waves and structures.

Free Falling Wedge
In this section, with the added effect of the dynamic response of a moving solid body, a twodimensional free-falling wedge that is subject to the motions of three degrees of freedom (DOF) after its entering a fluid domain is investigated with the developed numerical model. A solid symmetric wedge with a density of b = 466.6 kg/m 3 is selected for model simulations. It is positioned with its vertex of the obtuse angle pointing downward. The length of the wedge base is 1.2 m, and the side angle is 25°. The computational results in terms of the position of wedge's bottom vertex and its velocity are compared with the experimental data collected by Yettou et al. [56]. A 2D computational domain is considered, in which the length of the domain is 20 m and the height is 4 m, including a fluid domain with water depth of 1 m. Following the experimental setup, initially, the distance from the tip of the bottom vertex to the free-surface is 1.3 m. The wedge is suddenly released into the 20 m long tank. Two grid sizes as h = 0.04 m and h = 0.02 m are used for the numerical simulations. Figure 18 shows the induced velocity fields of the fluid domains at nine selected instants under the action of a free-falling wedge and its impact on the free surface. The variations of the interface also reflect the motions of the free-surface. The induced severe interfacial oscillations and coupled air motions after the wedge entering the water body can be observed. As the present numerical approaches treats more closely the two-phase flow conditions as a process of fully nonlinear fluidstructure interaction, it is believed more accurate and physically fitted results can be generated by the present model when compared to other models considering only the single-phase flow conditions without the inclusion of the effect of air phase. The use of orthogonal mesh system is highly essential to capture the vortex generation at the regions near or around a solid structure in the cases involving the interactions of air/water/solid. The comparisons between the simulated results and measured data from Yettou et al. [56] in terms of the time varying vertical position in z coordinate and vertical velocity of the bottom vertex for a free-falling wedge described above are presented in Figure 19a,b, respectively. The results obtained by using either the grid size of h = 0.04 m or h = 0.02 m are

Free Falling Wedge
In this section, with the added effect of the dynamic response of a moving solid body, a two-dimensional free-falling wedge that is subject to the motions of three degrees of freedom (DOF) after its entering a fluid domain is investigated with the developed numerical model. A solid symmetric wedge with a density of ρ b = 466.6 kg/m 3 is selected for model simulations. It is positioned with its vertex of the obtuse angle pointing downward. The length of the wedge base is 1.2 m, and the side angle is 25 • . The computational results in terms of the position of wedge's bottom vertex and its velocity are compared with the experimental data collected by Yettou et al. [56]. A 2D computational domain is considered, in which the length of the domain is 20 m and the height is 4 m, including a fluid domain with water depth of 1 m. Following the experimental setup, initially, the distance from the tip of the bottom vertex to the free-surface is 1.3 m. The wedge is suddenly released into the 20 m long tank. Two grid sizes as ∆h = 0.04 m and ∆h = 0.02 m are used for the numerical simulations. Figure 18 shows the induced velocity fields of the fluid domains at nine selected instants under the action of a free-falling wedge and its impact on the free surface. The variations of the interface also reflect the motions of the free-surface. The induced severe interfacial oscillations and coupled air motions after the wedge entering the water body can be observed. As the present numerical approaches treats more closely the two-phase flow conditions as a process of fully nonlinear fluid-structure interaction, it is believed more accurate and physically fitted results can be generated by the present model when compared to other models considering only the single-phase flow conditions without the inclusion of the effect of air phase. The use of orthogonal mesh system is highly essential to capture the vortex generation at the regions near or around a solid structure in the cases involving the interactions of air/water/solid. The comparisons between the simulated results and measured data from Yettou et al. [56] in terms of the time varying vertical position in z coordinate and vertical velocity of the bottom vertex for a free-falling wedge described above are presented in Figure 19a,b, respectively. The results obtained by using either the grid size of ∆h = 0.04 m or ∆h = 0.02 m are concluded to have nearly identical values, confirming the converged solutions. In terms of the model performance, the present results are shown to have good agreements with the published numerical solutions [36] and experimental data [56]. The maximum vertical velocity reaches 4.8 m/s at t = 0.56 s. This computational case demonstrates the model's capability in simulating the coupled motions between the dynamic responses of a free-falling wedge and the induced free-surface waves.

Wave Packet Interacting with a Floating Body
As an extension of the modeling study given in Section 3.5, a case considering a wave-packet interacting with a floating body is numerically investigated. In this case, the numerical setup follows the experimental study performed by Hadzic et al. [57] in a small towing tank of the Berlin University of Technology. The tested floating object was a rectangular prism with the dimensions of 10 cm in length (along the x direction), 29 cm in width (along y direction), and 5 cm in high (along z direction). The density of the body was 680 kg/m 3 , and the mass of the body was 0.986 kg. As shown in Figure 20, the center of the body was situated at a distance of 2.11 m away from the wavemaker. The computational domain is set as a 13 m long and 0.8 m high channel (water depth = 0.4 m). In Hadzic et al.'s [57] experiments, a flap-type wavemaker was controlled to produce a wave packet with a focusing point at the original location of the object. The grid sizes of ∆x = 0.02 m and ∆z = 0.005 m were used for the 2D numerical computations. Temporal variations of the body motions in three DOFs are computed during the interaction process between the wave packet and the floating body. Figure 21 Figure 22. The experimentally measured data from Hadzic et al. [57] are also included in Figure 22 for comparisons. The simulated wave packet is shown to fit nicely with the measured data. The maximum values of the wave elevation occur roughly at t =

Wave Packet Interacting with a Floating Body
As an extension of the modeling study given in Section 3.5, a case considering a wave-packet interacting with a floating body is numerically investigated. In this case, the numerical setup follows the experimental study performed by Hadzic et al. [57] in a small towing tank of the Berlin University of Technology. The tested floating object was a rectangular prism with the dimensions of 10 cm in length (along the x direction), 29 cm in width (along y direction), and 5 cm in high (along z direction). The density of the body was 680 kg/m 3 , and the mass of the body was 0.986 kg. As shown in Figure  20, the center of the body was situated at a distance of 2.11 m away from the wavemaker. The computational domain is set as a 13 m long and 0.8 m high channel (water depth = 0.4 m). In Hadzic et al.'s [57] experiments, a flap-type wavemaker was controlled to produce a wave packet with a focusing point at the original location of the object. The grid sizes of x = 0.02 m and z = 0.005 m were used for the 2D numerical computations. Temporal variations of the body motions in three DOFs are computed during the interaction process between the wave packet and the floating body. Figure 21 depicts the time-based changes of the oscillating angle of the flap wavemaker that generate a wave packet. The simulated time varying wave elevations at the locations of x = 1.65 m and x = 2.66 m during the evolution of the wave packet are presented in Figure 22. The experimentally measured data from Hadzic et al. [57] are also included in Figure 22 for comparisons. The simulated wave packet is shown to fit nicely with the measured data. The maximum values of the wave elevation occur roughly at t = 5.8 s and t = 7.7 s at the locations of x = 1.65 m and x = 2.66 m, respectively. In general, the results show the increasing trend of the wave elevation until t = 5.8 s at the location of x = 1.65 but delayed to t = 7.7 s at the location of x = 2.66 m. Afterward, the wave elevation is shown to have a trend of periodic decay with time.
In terms of the three DOF motions of a floating body, it is generally defined that the rotational motion that is referenced to its longest axis (in this case, the body width in y direction) as the rolling motion. Then, in this study, the translational motion along the x direction is defined as the sway motion. The vertical motion is as usual the heave motion. The time variations of computed sway, heave, and rolling motions with values referenced to its original body center location are plotted together with the measured data from Hadzic et al. [57] in Figure 23. The comparisons suggest that the predicted body motions in 3DOFs are in reasonably good agreements with measured ones. As shown in Figure 23a, the floating body's sway motion reaches the maximum value of 0.11 m at t = 7.6 s. The observed oscillatory heave motion is indicated in Figure 23b, where the motion varies within the range of −0.04 to 0.04 m, until it damps out at roughly t = 8.5 s. Figure 23c shows the time-based variations of the rolling angles of the floating body motion in which the angle of rolling motion ranges approximately from −20° to 20°. Finally, two snapshots showing the free-surface deformation and fluid velocity vectors surrounding the floating body at t = 7.2 s and t = 7.6 s are presented in Figure  24. Overall, the good matched comparisons between the simulated fluid and body motions with the experimental measurements demonstrate, again, that the developed two-phase viscous flow model can provide reasonable predictions on waves interacting with a floating body.    In terms of the three DOF motions of a floating body, it is generally defined that the rotational motion that is referenced to its longest axis (in this case, the body width in y direction) as the rolling motion. Then, in this study, the translational motion along the x direction is defined as the sway motion. The vertical motion is as usual the heave motion. The time variations of computed sway, heave, and rolling motions with values referenced to its original body center location are plotted together with the measured data from Hadzic et al. [57] in Figure 23. The comparisons suggest that the predicted body motions in 3DOFs are in reasonably good agreements with measured ones. As shown in Figure 23a, the floating body's sway motion reaches the maximum value of 0.11 m at t = 7.6 s. The observed oscillatory heave motion is indicated in Figure 23b, where the motion varies within the range of −0.04 to 0.04 m, until it damps out at roughly t = 8.5 s. Figure 23c shows the time-based variations of the rolling angles of the floating body motion in which the angle of rolling motion ranges approximately from −20 • to 20 • . Finally, two snapshots showing the free-surface deformation and fluid velocity vectors surrounding the floating body at t = 7.2 s and t = 7.6 s are presented in Figure 24. Overall, the good matched comparisons between the simulated fluid and body motions with the experimental measurements demonstrate, again, that the developed two-phase viscous flow model can provide reasonable predictions on waves interacting with a floating body.  The promising results obtained from the present 2D two-phase flow model suggest that the numerical approaches presented in this study can be extended to the development of a 3D model to study the three-dimensional two-phase flow problems. By using similar methodologies as those proposed in the present study, the Navier-Stokes-equations-based projection method, the LS method, and the IB approach can be extended to the domain of three dimension, including using a combined Eulerian Cartesian and Lagrangian grid system. Certainly, the complexity of the problems arises from modeling the induced 3D motions of a floating body, where six degrees of freedom should be considered. The movements of a rigid body subject to the external forces and moments will be governed by the equations of motion covering three linear translational and three rotational motions. Following the similar procedure of using a combined Cartesian and Lagrangian grid system, the movements of the 3D interfacial boundaries and the dynamic responses of a moving body can be simulated. As a future research direction, the present research group plans to extend the developed 2D two-phase free-surface flow model to an applicable 3D model for solving the 3D problems.

Conclusions
A compact interfacial method by solving the Navier-Stokes equations with added source terms of external forces from immersed bodies for modeling two-phase flow problems is presented in this study. The treatments of the moving free-surface and the immersed bodies with separately using the LS and IB methods are included in the developed model. It can be noted that the present study uses a combined Eulerian Cartesian and Lagrangian grid system to avoid the re-meshing procedure for two-phase flow modeling involving coupled fluid-structure interactions. A higher-order Runge-Kutta integration in time and WENO scheme in space are included to discretize the LS equation to track the interfacial positions with severe free-surface deformations. The proposed numerical technique does not require the conformation of the grid lines onto the boundary of an immersed object, and at the same time, the effect of an immersed object can be implemented through the righthand side source term of the Navier-Stokes equations.
The study cases include both the moving objects in a viscous fluid and interfacial two-phase flows, such as dam-break flows and wave motions over a submerged structure and the motions of a floating body. Good accuracy of the computed results, as confirmed by the comparisons with published results in calculating the variations of velocity/pressure/vorticity field and the hydrodynamic forces, was achieved. It is demonstrated that the present numerical model can simulate reasonably well for the cases of a dam-break flow passing over either a submerged semicircular or a rectangular cylinder. The model was also tested successfully to simulate the periodic waves propagating over a coastal bar, where the predicted free-surface elevations reflect an improved resolution by comparing to the Flow-3D simulations, using the same grid size. The predicted results are shown to have slightly better agreements with the experimental measurements when comparing to those from Flow-3D, using either the laminar or turbulent (k ε) flow module. Moreover, the use The promising results obtained from the present 2D two-phase flow model suggest that the numerical approaches presented in this study can be extended to the development of a 3D model to study the three-dimensional two-phase flow problems. By using similar methodologies as those proposed in the present study, the Navier-Stokes-equations-based projection method, the LS method, and the IB approach can be extended to the domain of three dimension, including using a combined Eulerian Cartesian and Lagrangian grid system. Certainly, the complexity of the problems arises from modeling the induced 3D motions of a floating body, where six degrees of freedom should be considered. The movements of a rigid body subject to the external forces and moments will be governed by the equations of motion covering three linear translational and three rotational motions. Following the similar procedure of using a combined Cartesian and Lagrangian grid system, the movements of the 3D interfacial boundaries and the dynamic responses of a moving body can be simulated. As a future research direction, the present research group plans to extend the developed 2D two-phase free-surface flow model to an applicable 3D model for solving the 3D problems.

Conclusions
A compact interfacial method by solving the Navier-Stokes equations with added source terms of external forces from immersed bodies for modeling two-phase flow problems is presented in this study. The treatments of the moving free-surface and the immersed bodies with separately using the LS and IB methods are included in the developed model. It can be noted that the present study uses a combined Eulerian Cartesian and Lagrangian grid system to avoid the re-meshing procedure for two-phase flow modeling involving coupled fluid-structure interactions. A higher-order Runge-Kutta integration in time and WENO scheme in space are included to discretize the LS equation to track the interfacial positions with severe free-surface deformations. The proposed numerical technique does not require the conformation of the grid lines onto the boundary of an immersed object, and at the same time, the effect of an immersed object can be implemented through the right-hand side source term of the Navier-Stokes equations.
The study cases include both the moving objects in a viscous fluid and interfacial two-phase flows, such as dam-break flows and wave motions over a submerged structure and the motions of a floating body. Good accuracy of the computed results, as confirmed by the comparisons with published results in calculating the variations of velocity/pressure/vorticity field and the hydrodynamic forces, was achieved. It is demonstrated that the present numerical model can simulate reasonably well for the cases of a dam-break flow passing over either a submerged semicircular or a rectangular cylinder. The model was also tested successfully to simulate the periodic waves propagating over a coastal bar, where the predicted free-surface elevations reflect an improved resolution by comparing to the Flow-3D simulations, using the same grid size. The predicted results are shown to have slightly better agreements with the experimental measurements when comparing to those from Flow-3D, using either the laminar or turbulent (k − ε) flow module. Moreover, the use of an interfacial method in this study is capable of capturing directly the vortex formulations at the recirculation zone downstream of a submerged structure. With detailed calculations of the free-surface deformation, the present model can provide a better estimation of the free-surface elevation and flow patterns for strongly nonlinear interaction between waves and structures. It is concluded with certainty that the two-phase free-surface flow model developed in this study is able to successfully simulate nonlinear water wave problems, including the accurate predictions of severe wave deformation processes.
For cases with dynamic responses of moving bodies, the present model is also extended to compute with good comparison results the motions of a free-falling wedge as it is entering a water body and wave packet interacting with a 2D floating body. The results obtained from the present study suggest that a two-phase flow model with a coupled immersed body calculation can be recommended for the investigation of nonlinear water waves interacting with structures, and with the integration of equations of motion, to predict the responses of a floating body. Our findings also support the conclusions that the present model is an effective numerical tool for free-surface flow computations in a domain including gas, liquid, and solid phases, with the use of a combined Eulerian Cartesian and Lagrangian grid system.