Fluid Structure Interaction of Buoyant Bodies with Free Surface Flows: Computational Modelling and Experimental Validation

In this paper we present a computational model for the fluid structure interaction of a buoyant rigid body immersed in a free surface flow. The presence of a free surface and its interaction with buoyant bodies make the problem very challenging. In fact, with light (compared to the fluid) or very flexible structures, fluid forces generate large displacements or accelerations of the solid and this enhances the artificial added mass effect. Such a problem is relevant in particular in naval and ocean engineering and for wave energy harvesting, where a correct prediction of the hydrodynamic loading exerted by the fluid on buoyant structures is crucial. To this aim, we develop and validate a tightly coupled algorithm that is able to deal with large structural displacement and impulsive acceleration typical, for instance, of water entry problems. The free surface flow is modeled through the volume of fluid model, the finite volume method is utilized is to discretize the flow and solid motion is described by the Newton-Euler equations. Fluid structure interaction is modeled through a Dirichlet-Newmann partitioned approach and tight coupling is achieved by utilizing a fixed-point iterative procedure. As most experimental data available in literature are limited to the first instants after the water impact, for larger hydrodynamic forces, we specifically designed a set of dedicated experiments on the water impact of a buoyant cylinder, to validate the proposed methodology in a more general framework. Finally, to demonstrate that the proposed numerical model could be used for a wide range of engineering problems related to FSI in multiphase flows, we tested the proposed numerical model for the simulation of a floating body.


Introduction
Fluid structure interaction (FSI) studies the physical phenomena related to the motion of a solid body, possibly compliant, within an encompassing fluid [1]. The development of reliable models for FSI is relevant for several engineering fields [2,3], including naval or marine engineering [4][5][6][7][8][9][10][11][12], damping, a proper simulation of fluid-structure interaction is crucial and in some cases involves solid bodies interacting with a free surface (i.e., wave energy harvesting). In such a case, the interaction between the solid and the water may be responsible for large impulsive loadings, with the humbling presence of the free-surface, which represents a jump of density of approximately three orders of magnitude and whose dynamics is very complex to be predicted. Therefore, the development of reliable numerical methods and experimental techniques to characterize such impulsive loading is crucial for the design of several kind of engineering devices, such as marine vessels and hydraulic structures. Several works are available in literature dealing with the simulation of fluid-structure interaction in the presence of a free surface [6,9,[35][36][37]. Different numerical models, alternative to Navier-Stokes, may be also employed to simulate the fluid flow with a free surface, such as Smoothed Particles Hydrodynamics (SPH) [54,55], Lattice Boltzmann Method (LBM) [8,56] and constrained interpolation profile (CIP) method [57,58]. However, only few reliable numerical methods with strong coupling modelling techniques have been proposed in the scientific literature so far and only part of them refer to buoyant rigid bodies (e.g., References [59,60]), a problem relevant to the design of marine vessels [7,[61][62][63] and aerospace structures [64] that are often designed to resist to large displacements and impulsive accelerations.
In this paper, we propose a computational model for the FSI of buoyant rigid bodies in the presence of free surface flows and we validate it with literature and custom designed experimental data for different fluid structure interaction problems. We adopt tightly coupled algorithm through a Dirichlet-Neumann partitioned approach with fixed-point iterations. The incompressible free surface flow is modelled through the volume of fluid approach (VOF), the partial differential flow equations are discretised thorugh the finite volume method and solid displacement is described by the Newton-Euler equations. The proposed algorithm is implemented within the foam-extend-3.2 open source framework [65], allowing a substantial code reutilization. Specifically, the flow solution is based on the interDyMFoam solver that utilizes the Arbitrary Lagrangian Eulerian method (ALE) and the PIMPLE algorithm [66] to solve the multiphase flow equations in a deformable domain, while we implement the solution algorithm for the Newton-Euler equations and the coupling procedure.
The rest of the paper is organized as follows. In Section 2 we describe the proposed methodology. In particular, in Sections 2.1 and 2.2 we discuss on the fluid and the structure modeling, respectively, while in Section 2.3 we describe the proposed coupling algorithm. The methodology is experimentally validated in Section 3, by comparing numerical results to literature ones and ad-hoc experiments. Conclusions are drawn in Section 4.

Methodology
The problem in study is schematically represented in Figure 1. Solving an FSI problem requires the physical modeling of two non-overlapping domains, the fluid, F and the solid, S, ones, separated by an interface Σ.
Herein we consider two immiscible fluid phases, namely phase 1 and phase 2 and we assume that the flow is incompressible. The solid domain is rigid and free to move subject to fluid, gravity or other external forces.
We consider an inertial reference frame Oxyz and a body fixed frame Gxȳz, which defines the orientation of S with respect to Oxyz (see Figure 1). G is the center of mass of S and the superimposed bar (·) identifies the vectors expressed in body fixed coordinates.
We solve the FSI problem through a partitioned Dirichlet-Neumann approach with fixed-point iteration [2,45]. Thereafter separate models are utilized for F (see Section 2.1) and S, (see Section 2.2).

Fluid Flow
The incompressible flow of two immiscible fluids is described by the following continuity and momentum conservations: where u is the fluid velocity, p is the pressure, µ is the dynamic viscosity, ρ is the density, g is the gravitational acceleration and D(·)/Dt denotes the material derivative. Despite the flow is incompressible, ρ and µ vary in space and time according to the fluid phase. The force exerted by the fluid on the solid body is, where I is the identity matrix, Σ is the body surface and n is the unit normal to the solid body. The moment of F with respect to G is being d = r P − r G the distance of each surface element from the center of gravity of the body. The flow is numerically approximated utilizing the interDyMFoam algorithm included within the foam-extend-3.2 software package [65]. It implements the volume of fluid (VOF) [67][68][69][70] and models the incompressible flow of two immiscible fluids through a single set of the Navier-Stokes and continuity equations (Equation (1)). The interface between the two phases is tracked by solving Equation (4) [69]: where α is the water volume fraction and u r is the relative velocity of the two phases. The second term of Equation (4) is introduced to contrast the numerical diffusion of the fluid interface [69].
Equations (1) and (4) are discretized through a generalized unstructured finite volume approach with a collocated variable arrangement [71]. The pressure velocity coupling and the non linearity of Equation (1) are solved using the PIMPLE algorithm [72]. The discrete version of Equation (4) is solved through the Multidimensional Universal Limiter with Explicit Solution (MULES) procedure [66].
Turbulent flows can be modeled with different levels of approximations including Reynolds averaged Navier Stokes (RANS), large eddy simulation (LES) and direct numerical simulation (DNS) [66]. The most common turbulence models for RANS and LES are already included in the foam-extend-3.2 environment.
As the motion of the solid body distorts the fluid domain, the computational nodes are re-allocated after each time-step to avoid depletion of the mesh quality. Several procedures for the mesh deformation are implemented within the foam-extend library [45,[73][74][75]. The numerical solution of Equations (1) and (4) accounts for the mesh deformation through the arbitrary Lagrangian Eulerian (ALE) procedure [45,74].

Rigid Body Dynamics in Space
The position (r p ) of the generic point p in Oxyz can be expressed as where r G is the position vector of G and A is the rotation matrix that defines the orientation of S with respect to Oxy. The motion in space of an unconstrained body follows the Newton-Euler equations [76], being m the mass, F ext the external force,J the inertia matrix,M ext the external moment, andω the angular velocity of S. The system of Equation (6) can be regarded as a system of differential equations in r andω. Sinceω cannot be directly integrated, we introduce the four component vector containing the Euler parameters to obtain the orientation of the body as follows [76,77]: where e T = [e 1 , e 2 , e 3 ] andẽ is the skew matrix associated to e. The following equation provides the relation betweenω and p:ṗ where G = [−e,ẽ + e 0 I].
The system of ordinary differential equations from (5) to (9) is numerically integrated through the Euler explicit method.

Coupling Algorithm
The FSI problem is solved through a partitioned (or segregated) approach, which means that the mathematical models for the fluid flow (i.e., Equations (1) and (4)) and for the solid motion (i.e., Equations (5)-(9)) are integrated separately [2]. To accomplish the interaction between the fluid and the structure, F communicates F and M to the solid and S communicates the displacement to F . A strong or tight, coupling is required to ensure convergence when the density of the solid is comparable or lower, than the fluid density [1,2]. Being the algorithm segregated, such a coupling is accomplished through the iterative procedure reported in Figure 2   At the beginning of each time step, τ, the time step duration t(τ) is calculated according to the Courant condition, assuming that the fluid flow determines the most stringent stability constraint. Then, the interface fields are initialized. In particular, the residuals of the FSI sub-cycle and the interface displacement are set to zero (i.e., R 0 = 0, ∆ 0 = 0) and the relaxation factor φ 1 is initialized 0.25, with the superscripts referring to the FSI subcycle.
After initialization, the FSI sub-cycle, identified in the followings by the superscript n, begins with the calculation of the interface forces. Fluid force and moment are calculated through Equations (2) and (3) and then utilized to solve the system of Equations (5)- (9). The displacement of the points on the fluid-solid interface is then estimated as: The algorithm proceeds with the calculation of the interface displacement. The displacement of the generic node i on Σ with respect to the previous time-step is estimated as: Under-relaxation is necessary to guarantee the convergence of the segregated algorithm, as widely recognized in the literature [1,2]. Having determined the displacement of Σ, the computational mesh for the fluid domain can be adjusted to meet the new boundary configuration. For each FSI sub-cycle, the flow equations have to be integrated within the same time-step. The ALE procedure calculates the convective fluxes accounting for the displacement of the faces of each computational cell within a time-step. A correct estimation of such fluxes requires that each computational node is displaced from its position at time-step (τ − 1) to the final position at sub-cycle n of time-step τ. Within the foam-extend-3.2 framework, this is accomplished by first moving back the mesh to its configuration at time-step τ and then smoothing the nodes to their final position. The procedure described in Section 2.1 is then utilized to solve the flow equations. At any step, mass is updated guaranteeing continuity. The residuals are defined as: The convergence of the FSI sub-cycle is assessed by comparing the magnitude of the residuals to a user-defined threshold R max ∑ where N Σ is the number of computational nodes in the fluid solid interface. If the inequality reported in Equation (13) is not verified, the relaxation factor is estimated before cycling the FSI sub-cycle starting from the flow solution (see Figure 2), according to the Aitken's procedure [1,2]: where φ max is a limiting value that can be varied for each application.

Validation
The methodology described in Section 2 is validated by comparison against experimental data. Both independent literature data [40,78,79] and custom designed experiments are considered in this validation.
Despite the proposed methodology is able to simulate 3D hydrodynamics and 6-DOF motions, for validation we concentrate on 2D hydrodynamics and 1-DOF structure motion to facilitate the comparison between numerical results and experimental data.

Water Entry of a Non-Buoyant Wedge
Zhao et al. [78] experimentally studied the water impact of a rigid wedge characterized by a total mass of 241 kg, a deadrise angle of 30 • , a width of 0.5 m and a breadth of 1 m. The specimen is bounded to move in the vertical direction and impacts the water at a velocity of 6.15 m/s. Following the indications in Reference [80], we assume that the flow is two dimensional (2D), being the breadth of the specimen larger than its width. Thereafter the fluid domain is discretized through a 2D structured mesh with about 5200 computational cells. Moreover, according to References [10,32], we assumed that turbulence modeling is not required for impulsive impact phenomena.
Temporal discretization, is performed through the first order accurate Euler implicit scheme [66]. The Gauss upwind [66] scheme is utilized to discretize the convective terms in Equation (1) and the Gauss interface compression scheme is employed for the convective term in Equation (4) [69]. Other spatial discretizations are performed through the Gauss linear procedure.
In Figure 3 we compare the numerical and experimental results for the vertical velocity of the wedge. We note that numerical results compare very well with the experiments. Specifically, the average relative error is being N sample the number of the considered experimental samples, v CFD the numerical approximation vertical component of the velocity of the specimen center and v EXP its experimental evaluation.

Water Entry of a Buoyant Wedge
The test performed in the previous section is not particularly challenging, as the specimen is relatively heavy and not buoyant, being the maximum mass of the displaced water 144 kg much lower than the mass of the specimen (241 kg). In this section, the proposed numerical model is verified against the water impact of buoyant bodies performed by Panciroli et al. [79]. Therein, the water impact of several buoyant wedges with different average deadrise angles, curvatures and impact velocities are reported. All the wedges have a square base of width 0.2 m. Here, we focus on the specimens with a flat wetted surface, a constant deadrise angle of 35 • and a mass of 0.412 kg. The impact velocity varies between 1.87 m/s and 3.82 m/s as a function of the drop height h and is reported in Figure 4 together with its theoretical value, for all the considered experiments.  Analogously to the numerical test described in Section 3.1, we consider a 2D flow [80] and we utilize a laminar modeling approach [10,32]. A hexa-dominant unstructured mesh, with about 15,000 computational cells is used to discretize the fluid domain. We utilized the same numerical setting described in Section 3.1 except for the discretization of the convective term of Equation (1b) that is now performed through the linear-upwind procedure.
Numerical results are compared to experimental measurement from Reference [79] in terms of vertical displacement of the center of mass. We denote with δ EXP and δ CFD the experimental and numerical displacement, respectively. Figure 5 and Table 1 show the good agreement between the proposed numerical model and the experimental data for all the considered cases. The accuracy of CFD is quantitatively assessed by the average relative error, evaluated as and reported in Table 1. Notably, ε δ is always below 4%. Figure 5 shows that CFD slightly underestimates the wedge penetration. This behavior can be explained by the 2D approximation that overestimates the fluid forces on the specimen. In fact, border effects, which are not considered in the 2D simulation, are expected to decrease the pressure close to the wedge borders, thus decreasing the fluid reaction [80].

Water Entry and Exit of a Buoyant Cylinder
The test cases proposed in Sections 3.1 and 3.2 focus on the water entry problem. Experimental data are, in fact, limited to the first instants after the water impact, when hydrodynamic forces are larger. To validate the proposed methodology in a more general framework, we specifically designed a set of dedicated experiments on the water entry and exit of a buoyant cylinder [81].

Experimental Setup
The experimental setup consists of a water tank made of stainless steel and Lexan R , with an internal volume of 1.5 m × 1.85 m × 0.7 m. The tank is filled with water up to a level of 0.4 m. Two transparent sides grant optical access to the tank. The specimen, a hollow PVC cylinder filled with polystyrene with a diameter d = 0.16 m and a breadth l = 0.295 m, is rigidly connected to a wood sledge that is allowed to move along two vertical aluminum rails. We consider four drop heights of 0.25 m, 0.50 m, 0.75 m and 1.00 m and two specimen masses, m = m 0 and m = m 0 + 1 kg, being m 0 = 2.214 kg.
The specimen position over time is measured through two Spectrasymbol thinpot linear potentiometers embedded in one of the rails. A wiper with tip diameter of 2.7 mm is connected to the sledge and actuates the position sensor.
The optical acquisition of the experiments is performed via a Phantom R Miro 110 monochromatic high speed camera operating at a frame-rate of 2500 Hz for m = m 0 and of 2700 Hz for m = m 0 + 1 kg.
A capacitive accelerometer Adafruit ADXL335-5V is rigidly connected to the sledge and measures the specimen acceleration before water impact. Acceleration data are acquired only before the water impact because the measurement range of the accelerometer is ±3 g, while acceleration ensuing from the water impact is expected to be much larger.
The analogue signals from the accelerometer and the displacement sensor are synchronously acquired by a National Instrument NI USB-6009 acquisition board at a sampling rate of 48 kHz at a resolution of 14 bits. The acquisition system is controlled through an ad-hoc realized LabVIEW [82] dashboard.
Each drop test is repeated three times to assess the data fidelity against random and systematic errors. We average the data from the three repetitions and apply a smoothing average filtering to obtain the effective position and acceleration as functions of time. We set t = 0 when the specimen impacts the water free surface.
The impact velocity is of paramount importance in water entry problems, as discussed in Reference [83] and we obtain a robust estimation of it averaging the values yielded by three procedures: (i) Numerical differentiation of the displacement measurement through central difference approximation; (ii) Numerical integration of the measured acceleration through the trapezoid rule; (iii) Least square fitting of the specimen displacement to a second order function and algebraic differentiation.
The obtained impact velocity is reported in Figure 7, as a function of the specimen mass and drop height. For more details on the experimental setup, the reader can refer to References [81,84,85].

Numerical Setup
We consider here a 2D flow [80] and we utilize a laminar modeling approach [10,32]. We adopt the same numerical setup described in Section 3.2. In this case the hexa-dominant computational mesh is characterized by about 19,000 computational cells.

Results and Dicussion
Numerical and experimental evaluations of the cylinder displacement are reported in Figure 8 for all the combinations of h and m. We observe a general good agreement between CFD and experiments. A visual estimation of the cylinder position through the high speed camera is also reported in Figure 8 as a further confirmation of the accuracy of the proposed methodology. The reliability of the present simulations is quantitatively assessed by the percent relative errors estimated through Equation (16) and reported in Table 2. Remarkably, ε δ < 11% for all the considered tests.   We note that the differences between numerical results and experimental data may be generally related to an overestimation of the fluid forces by the numerical model. In fact, numerical simulations underpredict the maximum penetration of the specimen and larger fluid forces also cause a higher rebound after motion inversion (see Figure 8). As already highlighted in Sections 3.1 and 3.2, the 2D approximation may explain such an underprediction. Moreover, we discarded any friction dissipation on the experimental setup after water impact. Such a dissipation also concur in reducing the specimen momentum, in particular after the motion inversion, thus smearing the amplitude of the specimen oscillations. Figure 9 represents the free surface dynamics ensuing from the water impact, by comparing the high speed camera images to the CFD evaluation of the water volume fraction for different time instants. The proposed numerical model correctly captures the most relevant features of the hydrodynamics related to the water impact of the cylinder, such as the pile-up and water jet formation (Figure 9b), the air cavity formation (Figure 9c) and closure ( (Figure 9d,e) in the proximity of the maximum penetration. Slight differences between CFD and experiments are detected for the free surface pattern in Figure 9b,c, before the closure of the air cavity. Images from the high speed camera show a larger cavity compared to CFD. Moreover, the cavity closure is slightly anticipated in the CFD with respect to the experiments, as evidenced in Figure 9c, where the experimental air cavity borders points away from the cylinder, while the proposed numerical model predicts that the air cavity borders are directed inward. In this respect, we comment that such a discrepancy may be related to a poor prediction of the air flow. As a matter of facts, we focused on the hydrodynamics ensuing from water impact, that largely influence the cylinder motion, rather than on the modeling of the air flow. In particular, initial conditions are not realistic for the air flow, as the fluid domain is initialized with a constant null velocity and the cylinder is impulsively started at the impact onset. On the contrary, in the experiments the air flow at the impact onset is already completely developed, as a result of the specimen free fall path in air. We also note that the lack of turbulence modeling may also hinder the correct evaluation of the air flow behind the cylinder. Figure 10 shows the convergence history for h = 0.5 m and m = m 0 . For sake of clarity we reported only 25 equally spaced time-steps from t = 0 to t = 2 s. Convergence is relatively fast being reached with a number of FSI subcycles comprised between 4 and 7. We also comment that φ is equal to 0.25 for the first 3 iterations (see Equation (14)) and is then increased up to 0.6 thanks to the Aitken's procedure. Similar convergence histories are obtained also for the other cases.

Roll Motion of a Rectangular Structure
The previous tests are all focused on hull slamming problems. However, the proposed numerical model could be used for many other engineering problems related to FSI in multiphase flows and, among them, for the simulation of floating structures and their response to wave actions, which have attracted much interest in the hydraulic and ocean scientific community in the last decades [86][87][88][89].
In this section, we test the proposed numerical model for the simulation of a floating body. Jung et al. [40] experimentally studied the roll motion of a rectangular acrylic structure 0.3 m large, 0.1 m high and 0.9 m wide (normal to the motion plane). The structure is hinged on a water tank 35 m long, 1.2 m deep and 0.9 m wide and is free to roll about its center of gravity. Translations, as well as yaw and pitch, are not permitted. The water depth is 0.9 m. At equilibrium the width of the rectangle is horizontal and the draft of the body is 0.05 m. The moment of inertia with respect to an axis normal to the motion plane through G is J = 0.236 kg m 2 .
The free roll decay test was performed by releasing the structure from an initial roll angle of 15 degrees in calm water and measuring the roll angle variation as a function of time with a rotary sensor [40].
We utilize the numerical setup described in Section 3.2 and the same mesh structure. In this case the domain discretization results in about 7000 computational cells.
As evidenced in Reference [41], contrarily to the hull slamming problems described in the previous sections, the turbulent eddy dissipation is expected to significantly impact the amplitude of the roll motion. Thereafter, a proper turbulence treatment is necessary to capture the hydrodynamics induced by the solid body motion and, in turn, the rotation amplitude. To this aim, we adopt a RANS modeling approach with a k − ε standard closure [90,91]. Standard wall functions [91] are utilized for the solution of the turbulent kinetic energy and turbulent kinetic energy dissipation velocity equations at solid boundaries.
The reliability of the proposed methodology is assessed in Figure 11 that reports the rotation angle (β) as a function of time. The agreement between CFD and experimental results from Reference [40] is remarkable. Specifically, the average relative error is To demonstrate the potentiality of the proposed methodology in terms of flow analysis we represent the flow field generated by the solid body motion in Figure 12. Therein, the velocity vectors are represented only within water to allow the visualization of the vortex structures ensuing from the rotation of the rectangle. We comment that, is this case, the air motion is expected to have a minor impact on the structure dynamics, with respect to water. As expected, the fluid velocity is maximized at the onset of the motion, when the angular velocity of the solid structure is larger. Two counter-rotating vortexes are generated at the corners of the rotating box. Such vortexes subtract mechanical energy from the structure contributing to the reduction of the oscillation amplitude. This mechanical energy is dissipated through a turbulent dissipation mechanism [41].

Conclusions
In this paper we develop a computational methodology to model the FSI of rigid bodies in free surface flows. We adopt a tightly coupled partitioned approach to cope with structures whose density is comparable or lower to the fluid one.
The methodology is implemented within the foam-extend-3.2 open-source framework that allows a significant model re-utilization. The solution of the flow equation relies on the interDyMFoam solver that implements the VOF model and approximates Equations (1) and (4) through a generalized finite volume method. The 6-DOF solid body motion is described through the Newton-Euler equations that are numerically integrated utilizing the Euler explicit method. Quaternions or Euler parameters, are employed to parametrize the solid body orientation with respect to an inertial reference frame. Simplified solid body motions are obtained by nullifying the proper acceleration components in the Newton-Euler equations. Within the Dirichlet-Newmann approach, the tight coupling is obtained through fixed point iterations, whose convergence is boosted utilizing the Aitken's dynamic under-relaxation.
We validate the proposed methodology by comparing numerical results to experiments in four different simplified 2D geometries (i.e., 2D flow field and 1-DOF solid body motion), yet computationally demanding test cases. Specifically, we select three different water impact problems and the free decay of the roll motion of a buoyant rectangle.
Both literature and custom designed experiments are considered. We first studied the water impact of elongated wedges varying the deadrise angle, the mass and the impact velocity. Comparison with experiments from the literature is very satisfactory with an average relative error below 4%. Nevertheless, literature data were limited to a relatively short time after the water impact. Thereafter, we performed a dedicated set of experiments to assess the capabilities of the proposed methodology also for the water exit. By contrasting CFD with those experiments, we note that the average relative error is around or below 10%. Moreover, the method is capable of correctly describing the initial water impact, the motion inversion and the successive oscillations about its equilibrium position of a buoyant cylinder. Finally, since rotations were not included in the previous tests, we simulated the natural decay of the roll motion of a buoyant rectangle to assess the capability of the proposed methodology to simulate floating bodies. Also in this case the agreement between numerical predictions and experimental data is very good, being the average relative error below 4%. Summarizing, the proposed numerical model performs remarkably well for all the tests cases. Even in the worst case, the discrepancies fall within the confidence interval expected for the adopted numerical models. The simplifying hypotheses in the computation, such as 2D flow, laminar modeling or the eddy viscosity model for the turbulent closure may explain the differences between CFD and experiments.
We comment that, despite we considered only 1-DOF cases for validation, due to the availability of experimental results, the combinations of the four proposed test cases includes both pure translations and rotations. Thus, the proposed methodology can be considered reliable also in the general 6-DOF case.
Besides demonstrating the reliability of the proposed methodology, the test cases described throughout the paper also evidence its relevance for engineering applications. In fact, such problems are important in several engineering fields such as naval, civil and hydraulic engineering. Moreover, FSI combined with multiphase flow modeling can be useful for other applications such as bio-engineering or energy harvesting.
Future work will be devoted to implementing kinematic constraints that limit the solid body motion and to modeling compliant bodies.