1. Introduction
Moonpool-type floaters were initially proposed to be applied either as artificial islands or as structures enabling activities in their inner area, as they used their solid boundaries as protecting walls to the wave action. Nevertheless, in recent years, they have been associated with their use as wave energy converter (WEC) devices suitable for wave energy absorption, as well as floating pontoons for semi-submersible designs.
Several types of bottomless cylindrical hulls (known as moonpools) have been examined in the literature. The most common is a vertical cylinder with a vertical opening from the deck to keel through the hull [
1], whereas square floaters with square moonpools have been also studied [
2,
3]. Additional types concern toroidal configurations consisting of a circular (core) section centered upon a large circle (ring) with an infinitesimally small [
4] and finite [
5] slenderness core-to-ring-radius ratio, as well as two coaxial surface piercing truncated circular cylinders [
6]. The latter geometry was encompassed by an exterior partially immersed toroidal structure of finite volume and an interior coaxial free surface-piercing truncated cylinder. Thus, an internal annular free surface was developed, which was totally enclosed between the solid boundaries of the cylinders and open to the exterior fluid domain beneath the concentric bodies.
Hydrodynamics of the isolated truncated hollow cylindrical bodies have been investigated in the past. Garrett [
7] presented an analytical solution of the wave scattering problem by a vertical bottomless partially immersed circular cylinder with infinitesimal wall thickness using matched axisymmetric eigenfunction expansions. Analytical computations of the exciting wave forces and moments on the hollow cylindrical shell structure, together with sway and pitch hydrodynamic coefficients, were derived by Miloh [
8]. He also derived analytic expressions for the amplification ratios of heave (pumping) and pitch (seiche) modes of water motions inside the pond. The corresponding hydrodynamic characteristics of an open-ended vertically floating circular cylinder, as well as the induced water motion in the interior basin, accounting for the effect of the finite wall thickness, were presented by [
1,
9]. Furthermore, Liu et al. [
10] investigated the second-order sum frequency wave loads on a hollow cylindrical shell structure in long crested irregular seas. Additionally, in [
11,
12] the solutions of the diffraction and radiation problems of a moonpool cylindrical body of infinite thickness and of a partially bottom-opened moonpool were presented. A moonpool with a restricted entrance was also studied. In [
13], it was shown that the motion of the water inside the body was amplified as the size of the moonpool entrance was reduced. In addition, Liu et al. [
14] measured the viscous effect of a moonpool cylindrical structure with a restricted entrance by introducing a quadratic dissipation assuming an additional dissipative disk at the moonpool entrance.
Regarding the two coaxial surface piercing truncated circular cylinder configuration, the solutions of the linearized diffraction and radiation problems around two independently moving concentric truncated cylinders were presented in [
6,
15]. Furthermore, the mean- and the time-dependent second-order wave loads on such a type of floating structures were evaluated in [
16,
17]. Recently, Chao and Kim [
18] investigated the hydrodynamic performance of a two-body WEC by applying analytical solutions, whereas Kong et al. [
19] examined a moonpool platform-wave energy buoy composed of an inner cylindrical buoy and an outer toroidal-cylinder for wave power extraction by the relative heave motions between the inner and outer buoys.
Apart from wave energy applications, concentric surface-piecing circular cylindrical bodies have also been investigated, combined with the trapped wave modes that they can develop, i.e., the localized oscillations of the enclosed formed free surface area, whereas the motion in the exterior fluid region decayed to zero. Abramson [
20] presented a comprehensive review of liquid sloshing problems in moving containers, while Shipway and Evans [
21] examined the wave trapping phenomena by two concentric circular cylindrical shells with a vertical axis and a zero thickness. As an extension of [
21]’s work, McIver and Newman [
22] studied the case of vertically axisymmetric trapping structures, of finite volumes, composed by two concentric interior free surfaces.
All the aforementioned research is based on potential flow methodologies concerning the hydrodynamic solutions of several types of bottomless cylindrical hulls. Nevertheless, over the last decade, there has been a significant interest in computational fluid dynamics (CFD) modelling due to its detailed results. CFD modelling solves the fluid flow problem using Navier–Stokes equations, where the hydrodynamics of random-shape floating structures can be calculated in detail, providing representative results of the flow physics. Several authors have studied the importance of including viscosity in moonpool problems using CFD modelling. Lu et al. [
23] investigated the possibility of finding an artificial, empirically-based damping coefficient to the free surface condition inside a moonpool device applied to potential flow models. The latter exhibit overestimated magnitudes of wave forces as the fluid resonance takes place; hence, by introducing an artificial damping term, the potential flow model may work as well as the viscous fluid model. Lu et al. [
24] extended the results from [
23] on the wave elevation inside a moonpool body; whereas, in [
25,
26], the piston-mode in a 2D moonpool under forced heave oscillations was evaluated using a linear potential solver coupled to a viscous Navier–Stokes modelling. In addition, Fredriksen et al. [
27] investigated, numerically and experimentally, how a low forward/incoming current speed influences the resonant piston-mode resonance in a moonpool. They concluded that the moonpool behavior was marginally reduced with a low forward velocity, whereas this reduction was dependent on the heave forcing amplitude. Recently, in [
28] a viscous damping model for energy losses of a 2D fluid response in a moonpool due to wall friction and flow separation was proposed. Furthermore, the limitations of boundary element method solvers when viscous effects near floater’s sharp edges were taken into account and were exploited in [
29] in the case of a moonpool type floater.
In the present article, two different geometric configurations of a moonpool-type floater were investigated; a floater encompassed by two coaxial, free-surface piercing toroidal cylinders with vertical symmetry axes (Configuration 1, see
Figure 1a) and a floater consisting of an external toroidal cylinder and an internal coaxial surface piercing truncated circular cylinder (Configuration 2, see
Figure 1b). In the annular and the cylindrical fluid areas formed between the cylinders’ vertical walls, oscillations of the enclosed water columns were developed. A semi-analytical model was used to solve the wave diffraction problem in the context of linear potential theory [
30]. The method of matched axisymmetric eigenfunction expansions was applied to solve the relevant hydrodynamic problem, according to which the flow field around and inside the examined body was subdivided in coaxial ring-shaped fluid regions. Hence, appropriate series representations of the velocity potential can be established in each region, which were matched at the boundaries of adjacent fluid domains by enforcing the continuity conditions of the hydrodynamic pressure and the radial velocity in order to determine the unknown coefficients in the expressions of the diffraction potential. The exciting wave forces/moments and the wave elevation on both case studies were calculated and compared focusing on the resonant wave frequencies. Furthermore, the accuracy of the semi-analytical formulation was assessed in comparison to high fidelity CFD simulations.
The present work is structured as follows: in
Section 2 the solution of the diffraction problem is formulated, along with the evaluation of the hydrodynamic forces and wave elevation around/inside the examined case studies.
Section 3 is dedicated to the presentation of the used CFD solver, whereas, in
Section 4, the outcomes of the two methods (i.e., semi-analytical and CFD solvers) are compared for the two moonpool arrangements. Finally, conclusions are drawn in
Section 5.
3. Computation Fluid Dynamics—MaPFlow Solver
In this section, a brief description of the CFD modelling applied to the present bodies’ configurations is presented. For a more detailed review of the solver’s methodology, the reader may refer to [
29,
33]. In the previous mentioned publications, a series of test cases have been examined, including a moonpool configuration excited by incident waves that prove the liability of the solver in this type of simulations.
MaPFlow treats free surface flows as flows of two immiscible and incompressible fluids. The presence of the two fluids is described using the Volume of Fluid (VoF) approach [
34], while the requirement of the incompressibility is satisfied using the artificial compressibility method (ACM) [
35]. The governing equations are presented in Equation (46). The system of equations, in 3 dimensions, consists of 5 scalar equations. The equations are augmented by the pseudo-time derivatives of the variables. The aim of the numerical procedure is to drive these derivatives to zero; thus the original unsteady system of equations will be retrieved. The coupling of the equations is performed during the pseudo-time, where a relation between the density and the pressure field is assumed. The coupling is controlled through the relation
, where
β is a free parameter that in typical free surface flows ranges from 5 to 15 [
33,
36].
In the above integral equation,
is a control volume with boundary
,
is the vector of the unknown variables (pressure
, velocity
and volume fraction
), vector
contains the various source terms of the equations such as the gravitational field, while
and
denote the real and fictitious time, respectively. The original ACM when applied to two-phase flows scales locally with the density field and this results to a stiff numerical behavior. In order to remedy this dependency, the preconditioner
of Kunz [
37] is used; it is presented in Equation (47). Furthermore, to recover the conservative form of the equations, the unsteady terms of Equation (46) are supplemented with the matrix
, which is also presented in Equation (47).
Finally,
is the vector of the convective fluxes and
the vector of the viscous fluxes. The two vectors are given by Equation (48).
In Equation (48),
is the velocity difference between the contravariant velocity
and the grid face velocity due to the mesh motion
. In the present study, because the structures were considered fully restrained, the mesh was stationary; thus,
. Furthermore, the viscous fluxes are computed through the viscous stresses
, which in turn were expressed using the Boussinesq approximation, by Equation (49).
In the above equation, is the viscosity of the mixture, is the turbulent viscosity, the turbulent kinetic energy and the Kronecker’s symbol. Moreover, since here only small amplitude linear waves were investigated, the turbulence modelling was switched off ().
The generation and damping of the free surface waves were performed by adding source terms to the momentum equations in user-specified regions near the boundaries of the computational domain. The form of this source term is described by Equation (50). The source terms scale according to the coefficient
, which is controlled though the factor
and the function
defined inside the wave generation or damping zone.
has typically an exponential or polynomial form [
38];
is a free parameter that depends on the flow itself (e.g., the waves entering the damping zone) and the discretization process. In [
33], a parametric analysis of the effects of the coefficient
is presented. The source term
drives the solution to the desired velocity field
. The target velocity in the case of generation zone was the velocity field provided by an analytical solution of free surface waves (linear, Stokes
nth order waves, etc.); in the case of damping zone, the usual strategy was to set the normal velocity to the boundary equals to zero.
The above system of equations was discretized using the finite volume method. In each cell of the computational mesh, a control volume
was defined, where the unknown variables had a constant value and the control volume was represented by the geometric center of the computational mesh (cell-centered approximation). The surface integrals were decomposed into a sum of constant surface integral over the faces that make up the cell. The semi-discretized form of the equations is presented in Equation (51). The source terms of the equations are considered constant in each cell.
is the spatial residual of the discretization process (see Equation (52)).
The convective fluxes are computed using the approximate Riemann solver of Roe [
39,
40], whereas a second order interpolation scheme supplemented by a directional derivative to account for the skewness of the mesh has been applied for the computation of the viscous fluxes.
For the time discretization, a second order backwards differentiation scheme was used [
41] for the real time derivatives, while, for the pseudo-time derivatives, a first order explicit Euler scheme was utilized. In order to facilitate convergence, the pseudo-timestep scales locally in each cell based on a given Courant-Friedrichs-Lewy (CFL) number and the characteristics of the flow. For an analytic description of the time discretization process the reader can refer to [
33].
Finally, regarding the boundary conditions, at the wall boundaries, the no-slip boundary condition was applied; the symmetry plane was treated as a slip boundary, while the far field boundaries were approximated using ghost cells. At every far field boundary face, an approximate Riemann problem was solved between the state inside the domain and the far field conditions imposed by the ghost cell. It should be noted that, in the case of wave generation, the ghost cell’s values were given by the current wave solution.
The discretization of the equations led to a non-linear system of equations. The linearization process took places in pseudo-time. The resulting linear system was solved using a Gauss–Seidel iterative method, while it was permutated using the Reverse Cuthill-Mckee reordering scheme to accelerate convergence.
The above methodology was implemented as part of the CFD solver MaPFlow [
42,
43]; MaPFlow was developed in NTUA, solved the URANS equations in unstructured meshes and is capable of running up to 1000 parallel processes utilizing the MPI protocol.
4. Results
In accordance with the procedure presented in the previous sections, it is obvious that the accuracy of the method and the CPU time required was affected by the evaluation procedure of the Fourier coefficients in each fluid domain. In the present calculations, 80 terms were used for the series expansions of the velocity potential in the fluid domains of the I, III and V types, whilst 150 terms are retained for the velocity representation in the II, IV types.
Two different geometric configurations of a moonpool-type floater were investigated; a floater encompassed by two coaxial, free-surface piercing moonpool cylinders with vertical symmetry axes (Configuration 1, see
Figure 1a) and a floater consisting of an external toroidal cylinder and an internal coaxial surface piercing truncated circular cylinder (Configuration 2, see
Figure 1b) [
30]. Specifically, for the Configuration 1 (i.e., coaxial hollow cylinders) it held:
α = 13 m,
b = 12 m,
c = 9 m,
l = 6.083 m, and
h1 = 14 m,
h2 = 5.5 m (see
Figure 2a), whereas, for the Configuration 2 (i.e., coaxial truncated cylinder), it was assumed:
α = 13 m,
b = 12 m,
c = 9 m,
h1 = 14 m and
h2 = 5.5 m (see
Figure 2b). Both arrangements were considered restrained to the wave impact at a water depth of
d = 70 m.
Regarding the numerical results that were produced by the in-house CFD solver, the two different configurations were excited by linear waves with wave frequencies ranging from 0.6 rad/s to 1.2 rad/s, with a constant wave height of 0.6 m. Based on the linear wave theory of Airy, the velocity field of the wave was imposed on the one end of the domain, while, at the other end, wave damping was used. The waves were propagated along the x-axis, while the free surface was placed perpendicular to the z-axis. In order to save computational resources, only the half domain was simulated by setting symmetry boundary conditions to the xz-plane. Due to the large deviation of the wave lengths, two different computational domains were used, for each case. For wave frequencies between 0.6 and 0.8 rad/s (which corresponds to wave lengths between 170 to 95 m, respectively) the domain had a total length of 600 m, while, for the wave frequencies between 0.85 and 1.20 rad/s (where the wave lengths range between 85 and 45 m), a domain of total length of 300 m is used.
The domain discretization is depicted in
Figure 3. The discretization of the domain upwind of the configuration was based on the wave characteristics. Near the free surface, the domain was uniform along the z-axis, with a length that corresponded to 18 cells per wave height. The mesh along the wave propagation was also uniform, with a discretization that corresponded to at least 100 cells per wavelength in each domain for the studied wave frequencies. Furthermore, special care was taken near the sharp edges of the structure, as illustrated in
Figure 3a, to account for the flow separation that was expected in these regions. Finally, the mesh coarsened as it reached the lateral and the outflow boundary, to make sure that no wave reflections occurred. In order to capture the viscous effects close to the wall, the first point was placed at 0.5 mm normal to the wall.
Figure 4 presents the grid independency study. The M2 grid has the aforementioned mesh characteristics and a total size of approximately 5 million cells. Additionally, three more grids were examined. The M1 grid consisted of 1.6 million cells. The main difference compared to the M2 grid is that it had 8 cells along the wave height and 50 cells at the smallest wavelength. Furthermore, the characteristic size of the mesh close to the wall is two times larger, compared to the previous mesh. The M2-w mesh had the same characteristics as the M2 grid; however, the first node from the wall was placed at 1 mm normal to the wall. Finally, the M3 had a total size of 8.5 million cells. The M3 grid differed from M2 in the characteristics of the mesh near the structure, as well as upwind of the structure in the direction of the wave propagation. The cells per wave height were the same as M2; the total cells per smallest wavelength were 200; the mesh near the structure had two times smaller characteristic size. From the independency study, the M2, M2-w and M3 produced similar results. The characteristics of the M2 grid were chosen for all simulations. The time discretization was chosen as 1000 timesteps for the corresponding wave period [
33].
Figure 5 and
Figure 6 depict the horizontal and vertical exciting forces on the moonpool floater of the 1st examined Configuration, derived by the presented potential flow formulation. The results are also compared with the corresponding exciting forces on the outer moonpool floater (i.e., Body 2) (see
Figure 2a), without the presence of the inner moonpool floater (i.e., Body 1), as well as on Body 1, without the presence of Body 2 (i.e., isolated bodies). The results of the exciting forces are normalized by the term
πρgα2(
H/2).
The peculiar behavior of the horizontal exciting forces (see
Figure 5) on the Body 1 and Body 2, in the neighborhood of ω ~ 1.74 rad/s and ω ~ 1.23 rad/s, respectively, was due to the resonance pitch oscillations of the interior water basin of each Body 1 and 2. The same behavior has been reported in [
1] regarding moonpool floaters. It should be further noted that the peculiar behavior of the horizontal exciting force on the total Configuration 1 (i.e., floater of two coaxial, free-surface piercing toroidal cylinders) was moved at a lower wave frequency, i.e., ω ~ 1.13 rad/s. This transpose could be attributed to the annular water area that formed between the two toroidal cylinders, which resonated in a pitch at a lower wave frequency than the corresponding one of the independent moonpool floaters (i.e., Body 1 and Body 2). In addition, a secondary resonance of the horizontal exciting force on the total Configuration 1 was observed in the neighborhood of ω ~ 1.74 rad/s. In the latter wave frequency, the resonance pitch oscillations of the interior water basin inside the Body 1 occurred.
With regard to the vertical forces (see
Figure 6) on Body 1 and Body 2, resonances appeared in the vicinity of ω ~ 1.05 rad/s and ω ~ 0.7 rad/s, respectively. These peaks were due to the pumping resonance of the fluid motion in the interior water area of the moonpools Body 1 and Body 2. Concerning the vertical exciting force on the total Configuration 1, it can be seen (
Figure 6) that its pumping resonant wave frequency was located in between the corresponding pumping wave frequencies of Body 1 and Body 2 (i.e., ω ~ 0.8 rad/s).
Next, the horizontal and vertical wave exciting forces on the examined Configuration 2 (see
Figure 2b), derived from the theoretical formulation, are presented. Specifically, in
Figure 7 and
Figure 8 the values of F
x/(
πρgα2(
H/2)) and F
z/(
πρgα2(
H/2) on the Configuration 2 are depicted, respectively. The results are also compared with the corresponding exciting forces on the isolated Body 1 and Body 2, when the latter were considered alone in the wave field (i.e., Body 1 was assumed alone in the wave field without the presence of Body 2, and vice versa for Body 2).
The resonance of the horizontal exciting force on the Configuration 2 (see
Figure 7), which occurred in the neighborhood of ω ~ 1.125 rad/s due to the pitch fluid motion in the annular fluid area between the interior cylinder and the exterior toroidal body, was moved to ω ~ 1.24 rad/s for the Body 2 case, when the latter was assumed isolated in the wave field. As far as the forces on the isolated Body 1 are concerned, it can be seen that they followed a smooth variation pattern at the examined wave frequency range, without any sharp maximizations.
Regarding the vertical exciting forces on the Configuration 2 (see
Figure 8), the pumping resonance of the fluid motion in the annular area was depicted at ω ~ 0.85 rad/s [
20], whereas, for the isolated Body 2, this pumping resonance lied near ω ~ 0.7 rad/s (see discussion of
Figure 5). It should also be noted that the smooth variation pattern of the vertical exciting force on the isolated cylinder follows (i.e., Body 1).
Comparing
Figure 5 and
Figure 7, it can be concluded that both configurations 1 and 2 attained, in the non-resonance regime, comparable results concerning the horizontal exciting wave forces. Furthermore, their resonance pitch oscillations occurred at the vicinity of ω ~ 1.13 rad/s. Hence, it could be derived that the wave motion inside the pond of Body 1 in the first configuration did not seem to affect the floater’s horizontal forces. On the other hand, the secondary resonance of the horizontal exciting wave force on the Configuration 1 (i.e., ω~1.74 rad/s) was not observed in the case of Configuration 2. Regarding the vertical exciting forces F
z/(
πρgα2(
H/2)), see
Figure 6 and
Figure 8; in the case of Configuration 2, due to the truncated inner cylindrical body, they attained higher values compared to the vertical loads on Configuration 1. Hence, assuming both Configurations as heaving point absorbers for wave energy exploitation, Configuration 2 seems to be more efficient compared to Configuration 1.
In the sequel, the results from the theoretical formulation were supplemented by numerical ones, which took into consideration the viscous effects near the sharp edges of the examined moonpool-type floaters. The numerical results presented follow the configuration terminology described in the previous section. A grid consisting of 5 million cells was employed with symmetry conditions (half of the moonpool is modelled). Since the moonpool was considered fixed and the waves were essentially in the linear region, the flow is assumed laminar [
29,
44].
Regarding Configuration 1 the comparison of the non-dimensional horizontal and vertical exciting forces could be seen in
Figure 9 and
Figure 10. The illustration of the horizontal forces (
Figure 9) suggested that the comparison was relatively good for both methods. CFD predictions follow the analytical method’s trends with minor discrepancies. In the proximity of the peak excitation frequency (between 1.0–1.2 rad/s), the differences between the predictions became larger. More specifically, in the near peak frequency region, CFD results predicted earlier the peak of the Body 1 than the analytical method, leading to a rise in total force. At 1.13 rad/s, both methods predicted the maximum horizontal excitation force. However, it was evident that the amplitude of the CFD predictions was significantly smaller. This could be justified by both the phase difference of the peak excitation loads between the two bodies in the CFD predictions (discussed in the following) induced by the viscous effects that were taken into account in the CFD computations.
Regarding the vertical excitation force in
Figure 10, the comparison between the methodologies was fair, with the CFD simulations predicting smaller peak loads and shifted towards the lower frequencies. This has also been noted in [
29]. Nevertheless, the qualitative comparison for both the total force and the individual body forces was considered to be in good agreement.
The respective comparison for Configuration 2 (external toroidal cylinder with internal coaxial truncated cylinder) can be found in
Figure 11 and
Figure 12. The agreement between the numerical results is good in both the horizontal and vertical forces. As
Figure 11 suggests, CFD predictions followed the analytical method trends; nevertheless, the peak excitation frequency was predicted earlier. As opposed to Configuration 1, the individual body forces peak excitation frequencies were predicted in the same frequency, suggesting that the forces on the two bodies were in phase. This was also the case for the analytical method’s predictions. Finally, regarding the vertical force in
Figure 12, again the CFD predicted the peak excitation frequency earlier; however, as in Configuration 1, the predicted trends were the same.
In order to gain a better insight of this phase difference between each body’s horizontal force in Configuration 1, instantaneous free surface elevation contours at four snapshots of the wave period are presented in
Figure 13 at ω = 1.1 rad/s.
As
Figure 13 suggests, the sloshing mode was excited between the outer (Body 1) and the inner (Body 2) cylinder, which explains the phase difference between the horizontal force of each individual body.
Figure 14 presents the free surface elevations inside the moonpool (between Body 1 and 2) for both configurations with respect to the incident wave frequency. The space averaged free surface elevation inside the moonpool is presented, as well as in three circumferential positions at 0, 90 and 180 degrees (as shown in
Figure 13a). In both configurations up to ω = 1 rad/s, it was evident that the free surface elevation did not greatly vary between the various spatial positions, suggesting that the piston mode dominated the free surface elevation inside the moonpool. As the frequency increased beyond 1 rad/s the gauge at 90 degrees followed the spaced averaged elevation; however, the other two locations at 0 and 180 degrees clearly deviated (they increase in amplitude), thus indicating the appearance of a sloshing mode inside the moonpool. As the radial frequency was further increased, this mode dampens out. Furthermore, it is worth noting that peak frequencies of the piston mode resulted in an increase of the vertical forces, while the sloshing mode resulted in the resonance of the horizontal excitation forces.
This is further illustrated in
Figure 15, where the phase angle of the signal in the three wave gauges is presented. In smaller frequencies, the three signals had no phase difference; however, in larger frequencies, the phases started to deviate. In the case of ω = 1.1 rad/s, the signal in the upwind gauge (0 deg) and in the downwind gauge (180 deg) had a 90 degrees phase difference compared to the middle gauge (90 deg) and 180 degrees between them. In higher frequencies, there were still some phase differences; however, the amplitude of free surface was relatively small (see
Figure 13).
Lastly,
Figure 16 presents the y-vorticity field in the symmetry plane in the case of the Configuration 1 for different wave frequencies. Close to the resonance of the piston mode, excessive flow separation was noted near the sharp edges of the structure (
Figure 16b). Due to the loss of energy in this region of recirculating flow, the viscous solver predicted smaller amplitude of the piston mode (and thus of the vertical forces) compared to the one of the analytical solution. Contrary to
Figure 16b, the vorticity field in larger frequencies (where the amplitude of the piston mode is reduced significantly), the vortex shedding was not significant; thus, minor discrepancies were noted between the two approaches. It was also noted that the performance of inviscid potential models is influenced by the radiation effect of the inner water area, leading to discrepancies between the viscous solver and the analytical solution.