Abstract
Potential (inviscid-irrotational) and Navier-Stokes equation-based (viscous) flows were numerically simulated around a 2D floating rectangular body with a moonpool; particular emphasis was placed on the piston mode through the use of finite volume method (FVM) and boundary element method (BEM) solvers. The resultant wave height and phase shift inside and outside the moonpool were compared with experimental results by Faltinsen et al. (2007) for various heaving frequencies. Hydrodynamic coefficients were compared for the viscous and potential solvers and sway and heave forces were discussed. The effects of the viscosity and vortex shedding were investigated by changing the gap size, corner shape, and viscosity. The viscous flow fields were thoroughly discussed to better understand the relevant physics and shed light on the detail flow structure at resonant frequency. Vortex shedding was found to account for the most of the damping. The viscous flow simulations agreed well with the experimental results, showing the actual role and contribution of viscosity compared to potential flow simulations.
1. Introduction
A moonpool is a vertical well directly connecting the free surface to the deck of offshore vessels for drilling, diving, installation, or other offshore support works. The discontinuous openings on the bottom of ships not only increase resistance but also cause harmful wave motions that degrade the safety and working ability of crews. The behavior of the free surface inside the moonpool has two types of modes: piston and sloshing. The piston mode involves the vertical movement of water, and in this mode, a violent motion at resonant frequency is generated. Horizontal motions are observed in the sloshing mode, which results in increased resistance.
Most previous research can be divided into two groups: resistance increase and violent water motions. Studies on resistance increase mainly considered drillships. English [1] measured the resistance and wave height for different moonpool shapes of a drillship. Yoo and Choi [2] examined detailed moonpool shapes with appendages. Veer and Tholen [3] surveyed the principal dimensions of moonpools for several drillships and discussed the mode of the free surface related to the increase in resistance.
Fukuta [4] and Aalbers [5] studied the water motion and its resonant frequency for drillships. They both employed a spring-mass-damper system. Fukuta [4] proposed an empirical formula to estimate the resonant frequency. Aalbers [5] measured the damping coefficient from free-decay tests for the spring-mass-damper model. Molin [6] derived analytic formulas of the resonant frequency for 2D and 3D moonpools based on the potential theory.
Two-dimensional moonpools—for example, catamarans and side-by-side moored vessels—have the same resonant behavior. Faltinsen [7] presented a summary on the piston-mode resonance for twin hulls. As new concepts for floating structures are developed, such as liquefied natural gas floating production storage offloading (LNG-FPSO) and liquefied natural gas floating storage and regasification unit (LNG-FSRU) with a side-by-side mooring arrangement, piston mode resonance has drawn more attention. The free-surface elevation inside the gap and the motion of the two floating bodies are important design factors. In general, linear potential theory in the frequency domain is widely used in the design process. However, this approach inherently overestimates the wave elevation near resonance. Koo and Kim [8] and Hong et al. [9] studied side-by-side moored systems and addressed the discrepancy in the piston mode. The overestimation of the free-surface elevation can be suppressed by employing a damping lid or artificial damping [10–12]. However, the exact value of the artificial damping is not available prior to model tests, and improper use of the lid or artificial damping coefficient may decrease reliability in the numerical estimation.
Faltinsen et al. [13] investigated the piston mode occurring in a square-shaped catamaran hull form under forced heave motion by using experiments and linear numerical analysis. Kristiansen and Faltinsen [14] compared their simulation results by using the vortex-tracking method with experimental results. They recently used a finite volume method-based domain decomposition method, in which the viscous domain for the body region communicates with the outside potential region [15].
The aim of this study is to determine the effects of eddy damping and viscosity on the piston mode by varying the gap size, corner shape, and numerical solver. Fully nonlinear potential flow solver and Navier-Stokes equation solver for viscous flows were employed for computations using different gap sizes and corner shapes to determine the contributions of viscosity and vortices to the gap flows. In this paper, we present the governing equations and numerical methods of each flow solver. The computational and boundary conditions along with the grid convergence test results are described. Numerical results are presented and compared with the experimental results of Faltinsen et al. [13]. Vortical flow fields according to computational fluid dynamics (CFD) are discussed in detail to better understand the relevant physics. The force response amplitude operator (RAO) and hydrodynamic coefficients are also compared for different flow solvers. Finally, conclusions are drawn based on the computations along with remarks on the eddy damping.
2. Numerical Modeling
A fully nonlinear potential code [16] and a commercial CFD software (FLUENT ver. 6.3) were used. The governing equations and numerical schemes are summarized as follows.
2.1. Potential Flows
Assuming the computational domain is filled with homogenous, inviscid, incompressible, and irrotational fluid, the Laplace equation is applied as a governing equation using the velocity potential :
On the free surface, the fully nonlinear dynamic and kinematic boundary conditions are, respectively, applied as follows: where is the time, is the free surface elevation, is the gravitational acceleration, is the water density, is the pressure on the free surface assumed to be zero, and is the vertical coordinate.
No normal-flux condition is applied to rigid boundaries:
To solve the Laplace equation with all boundary conditions, the Green function is applied, and the governing equation can be transformed to the boundary integral equation: where is the basic Green function satisfying the Laplace equation and is a solid angle (α = 0.5, when singularities are on the boundary).
By discretizing all boundary surfaces and placing a node on each panel—which is called the constant panel method—and integrating the discretized boundary integral equation, the velocity potential and its normal derivative on each segment can be obtained at each time step. For 2D problems, the Green function () as a simple source is given by where R1 is the distance between source and field points.
The velocity potential and free-surface elevation at each time step are obtained by using the Runge-Kutta fourth-order time integration scheme and the mixed Eulerian-Lagrangian method to solve (2). The total derivative can be used to modify the fully nonlinear free-surface conditions in the Lagrangian frame.
The nonlinear body pressure can be calculated from Bernoulli’s equation:
Finally, the nonlinear wave forces on a body can be obtained by integrating the nonlinear pressure over the instantaneous wetted body surface at each time step: where is the unit normal vector of the instantaneous body surface.
An artificial damping zone is used in front of the end wall to dissipate the outgoing wave energy as much as possible; that is, the wave energy gradually dissipates in the direction of wave propagation. This is necessary to avoid wave reflection from the rigid end wall. Both and -type damping terms are added to the fully nonlinear dynamic and kinematic free-surface conditions.
A Chebyshev five-point smoothing scheme is applied along the free surface during time marching to avoid the so-called saw-tooth numerical problem, which is not a physical phenomenon but a numerical instability. Detailed mathematical formulations for the numerical schemes used here are given by Koo and Kim [16, 17].
2.2. Viscous Flows
The governing equations for viscous, incompressible, and rotational fluids are the continuity and Navier-Stokes equations: where is the stress tensor, which is defined as follows: where is the molecular viscosity and is the unit tensor.
The governing equations can be written in the following integral form for a general scalar on an arbitrary control volume whose boundary is moving at the mesh velocity of the moving mesh : where is the control surface, is the area vector normal to , and is the source term.
The integral equation is discretized by the FVM into a linear algebraic system. The second-order upwind and central difference schemes are used for the convection and diffusion terms, respectively. The pressure implicit with splitting of operators (PISO) algorithm is used to couple the pressure and velocity for the continuity equation.
The free surface is represented by an explicit volume of fluid (VOF) method with a modified high-resolution interface capturing scheme. The range of the volume fraction of the air and water varies from 0 to 1, where the free surface is defined at 0.5.
3. Boundary and Computational Conditions
The computational domain was set up based on experiments by Faltinsen et al. [13], who performed tests in a 2D wave flume tank and measured the wave height inside and outside the moonpool. Figure 1 shows the schematic definition of the experimental setup.
In this research, we carried out numerical simulations for two moonpool geometries: narrow and wide gaps. The narrow moonpool had square and rounded corner shapes. Principal particulars are summarized in Table 1. Details on the boundary conditions are explained in the subsection for each flow solver.
3.1. Potential Flow Computations
The water depth was set to 1.03 m, and the width of the entire computational domain was 13.5 m; these are very similar to the wave tank conditions in the experiments. The bottom and the two bodies were treated as solid walls, that is, no flux through the boundary. A numerical damping zone of 2 m was placed on both sides of the free-surface area to prevent reflection from the side boundary walls.
In total, there were 265 elements and about 20 panels distributed over one wavelength . The number of panels per wavelength was chosen from a set of panel convergence tests; this is discussed in the next section.
The time step was T/128, where T is the period of prescribed body motion. Because the simulation was fully nonlinear, the instantaneous wetted body surface was updated and regridded at each time step. The fully nonlinear free-surface condition was also satisfied at the instantaneous free-surface positions.
3.2. Viscous Flow Computations
The CFD computations were carried out with half of the numerical wave tank; a symmetry boundary condition was applied to the center line of the moonpool. The width of the computational domain was 15.0 m, consisting of a 7.0 m wave zone and 8.0 m artificial damping zone. The horizontal size of the mesh in the damping zone was exponentially increased for artificial damping to prevent reflected waves [18]. On the bottom boundary, 1.03 m deep, no-slip wall was applied for the viscous computations. Hydrostatic pressure was imposed, and a gradient condition for the velocity and volume fraction was applied at the outlet.
The entire computational domain with the grid system periodically moved up and down during the computation according to the prescribed body motion, as shown in (11); the size of the cells remained unchanged except for the top and bottom boundary cells. Only the closest rows of cells to the top and bottom boundaries changed their height to absorb the imposed heaving motion at every time step: where , of 0.0025 m, and represent the vertical velocity, heave motion amplitude, and angular frequency, respectively.
The viscous flow was assumed to be laminar because the prescribed oscillating velocity was low with a Reynolds number at the approximate order of 103. A total of 176,000 rectangular cells, which is a highly dense grid for 2D computation, were used after a grid dependency test. Such a high-resolution grid may be capable of capturing flow structures including shedding vortices.
The minimum height and width of the cell were /20 and 0.0042 B, respectively. Equally spaced regular rectangular cells were allocated inside the moonpool and on the free surface to capture the precise vortex structure and wave elevation.
The time step was set to 0.005 s, which is around T/250. A linear ramp duration was employed during the initial 10 s to prevent sudden transient fluctuations.
4. Numerical Results
4.1. Grid Convergence
The number of panels per wavelength is an important variable for potential flow computations. Tests were performed on the narrow moonpool for 10, 20, 30, and 40 panels per wavelength. Figure 2(a) shows the wave height RAOs inside and outside the moonpool. The case of 10 panels per wavelength showed a discrepancy from the others; at least 20 panels per wavelength were required.
(a) Panel convergence
(b) Grid convergence
For the CFD simulations, three grid sets were compared: coarse, medium, and fine grids. Table 2 lists the details of each grid. The medium grid was only denser than the coarse one in the horizontal direction near the free surface. The fine one was also twice as dense vertically over the height of the body. Figure 2(b) shows the time history of the wave elevation for the narrow moonpool at the peak frequency. The three grid sets showed fully converged agreement with each other; the fine grid set was chosen to capture vortex shedding more precisely.
4.2. Wave Elevation
The free-surface elevations reached a periodically steady state within 60 s for all runs. Computations were performed which took longer than 60 s; data sampling and analysis were carried out on elevations obtained after 60 s. Figure 3 shows the free-surface elevation of viscous flows around the peak frequency. The black-dashed dot line represents the forcing amplitude for reference, and (red line) and (blue dotted line) are the wave elevation at the symmetry line inside the moonpool and 700 mm outside the body, respectively. The tendency of the time series of wave elevations observed in the CFD computations was very similar to that of the experiment results of Faltinsen et al. [13].
| (a) , narrow |
| (b) , wide |
The wave height inside the moonpool—a major concern in moonpool operations and side-by-side moorings—was compared with the experimental results shown in Figure 4. In the narrow moonpool, the fully nonlinear viscous results predicted the resonance frequency and its amplitude well, whereas the fully nonlinear potential flow solver predicted higher wave elevations than the experiments, which is the same as the linear results of Faltinsen et al. [13]. In the present example, the small forcing amplitude did not bring about appreciable nonlinearity in the potential flow solvers. Therefore, the discrepancy between the potential results and experiments can be attributed to the viscous effects. The viscous flow computations actually produced excellent correlation against the measurement and did not show any overestimation at the peak frequency.
(a) Narrow
(b) Wide
Interestingly, such overestimation by the potential flow solvers almost disappeared in the wide-gap case. The potential and viscous results are very similar, which means that viscous effects in this case are not significant. This will further be discussed in the next section. Compared to the narrow gap case, the free-surface elevation became lower, and the resonance occurred at a higher frequency. The difference in peak frequency relative to the experiment results was because of the reflected waves of the experiments, as pointed out by Kristiansen and Faltinsen for their experiments [14].
4.3. Viscous Flow Fields
CFD is an expensive computational tool compared to panel methods. Nevertheless, it can provide valuable clues to understanding the relevant physics. The velocity, vorticity, and pressure fields were investigated to better understand the flow fields.
Figure 5 shows the velocity vectors and vorticity contours at two frequencies off the resonant frequency: low and high. Weak vortices were generated at both sides. Vortices were attached around the corners and the flow fields were not disturbed at the low frequency. In contrast, the vortices at the high frequency were shed into a circular trail and slightly disturbed the vertical water motion like a water column inside the moonpool. The difference in inner and outer velocity magnitudes was not significant; thus, a small amount of sway force was induced.
| (a) |
| (b) |
Stronger vortices were generated at the resonant frequency, as shown in Figure 6. Because the water velocity inside the moonpool was faster than that on the outside, more vigorous vortex shedding was observed in the vicinity of the inner corner. The vortex changed the direction of rotation at the corner as the body changed the direction of its vertical motion. The shed vortices trailed into a circular motion, for which the diameter reached approximately half the size of the moonpool. The vortices interrupted the vertical motion and disturbed most of the water column, leaving only the upper part undisturbed. Shear-stress-driven vortices were observed on the side wall of the moonpool; they become another source of viscous damping.
The resonant vertical movement of the water column affected the pressure field, as shown in Figure 7. The inner pressure was periodically higher or lower than the outer pressure corresponding to the free surface elevation. The difference in pressure between inside and outside the moonpool induced a sway force. When the free surface reached its lowest level, the inner pressure became minimum; thus, a sway force inward the moonpool and heave force downward were generated and vice versa.
Figure 8 shows the velocity vectors and vorticity contours in the wide gap at the peak frequency. A similar vortex behavior as that in the narrow gap was observed. Vortices generated at the inner corner showed lower strength than those at the narrow moonpool owing to its weaker resonance. The shed vortices moved along a big circular trace; however, they could not block the movement of the water column as much as in the narrow gap. This is one reason for the small discrepancy between the potential and viscous flows in the wave height RAO shown in Figure 4.
4.4. Corner Shapes
In the previous section, the shed vortices from the square corner were shown to disturb the vertical water motion. We changed the corner shape to be rounded to investigate its effect on vortex shedding. A bilge radius of 5% of the draft was applied. Although the vortices did not completely vanish, the strength was weaker than for the square corners, as shown in Figure 9. The vortices in the boundary layer along the wall by shear stress still remained. The circular track of the shed vortices from the square corner was no longer observed. Instead, the water column moved up and down very freely like a solid column owing to the weaker vortex strength and the disappearance of the circular trace of the shed vortices.
The effect of weakened vortex shedding due to rounded corner is shown in Figure 10(a) with the wave height RAOs. In the fully nonlinear potential flows, the change in corner shapes little affected the maximum wave amplitude, as expected, and they showed similar RAOs. However, applying the bilge radius significantly increased the wave height in the viscous flows. This fact may be very important in the design of oscillating water column (OWC) wave energy converter (WEC). At low and high frequencies, the wave height was not affected by the corner shapes. As the frequency approached the resonant frequency, the wave height increased by up to 32%. This significant increase in the wave height explains the important role of vortex shedding associated with the interaction with free-surface motion. Vortex shedding from the hull or appendages such as damping plates may contribute to the suppression of harmful free-surface behavior. It is also notable that the damping effects are dominant around the peak frequency; thus, artificial treatments to suppress the wave elevation might be taken at the resonance.
(a) RAOs
(b) Phase shift
The phase shift was not related to the corner shapes and vortex shedding as shown in Figure 10(b). The phase shift of the two corner shapes in potential and viscous flows coincided with each other and showed good agreement with the experimental results.
5. Hydrodynamic Forces
5.1. Force RAOs
As discussed in the previous section, the movement of water inside the moonpool generates a periodic pressure change that causes hydrodynamic forces. Figure 11 shows the sway and heave force RAOs, respectively, calculated using viscous CFD solvers. The wave height RAO is included for reference. It is definitely observed that the wave height RAO is in the great agreement with the force RAOs, showing a similarity. The phase of the two forces also became close to that of the wave elevation and approached .
(a) Sway force RAOs
(b) Phase shift of sway force
(c) Heave force RAOs
(d) Phase shift of heave force
Table 3 shows the effects of the corner shape which resulted in the vortex damping near the peak frequency. In the narrow gap case, the corner shape can reduce the wave elevation and sway and heave forces by about 30%. This remarkable amount of damping should be accounted for while designing offshore structures. Taking the similarity between the wave elevation and the forces into account, accurate estimation of the wave height is highly required to estimate hydrodynamic forces acting on the bodies or mooring lines if they are moored.
5.2. Hydrodynamic Coefficients
The linear equation for heave motion is given in (12). With Newton’s second law in (13), we can use the least square method to calculate hydrodynamic coefficients [16]: where A33, M, B33, C33, f(t), F(t), and are the added mass, body mass, linear damping coefficient, restoring coefficient, wave exciting force, total force on the body, and heave displacement, respectively. The quadratic damping term may be added to the equation of motion.
By subtracting the restoring force from F(t), the added mass and linear damping can be derived as
If we consider the quadratic damping term, another term with the quadratic damping coefficient emerges. The linear and quadratic damping terms are coupled to each other. In the present study, we considered only linear damping.
Figure 12 shows the hydrodynamic coefficients for the narrow moonpool. The present results showed typical characteristics of the piston mode in the moonpool well, such as the negative added mass and linear damping approaching zero around the peak frequency.
The hydrodynamic coefficients from the potential-flow simulations were appreciably higher than those from the viscous flows, especially for sharp corners. This result is related to the overestimation of the free-surface elevation by the potential-flow solvers. Since forced-motion amplitude is small in the present example, the fully nonlinear potential-flow simulation is not significantly different from that of Faltinsen’s linear calculation. The effects of the corner shape are small, as intuitively expected, in the potential flows; in contrast, it has a great effect on the viscous flow, which can take the rotational fluid motions into account.
The hydrodynamic coefficients were found to be very sensitive to the hull forms; only the viscous flow solver provided reliable results. It is difficult to derive a general formula that consists of the main particulars of the moonpools. Further study on the eddy damping and its quantification is needed for practical purposes. We carefully propose to quantify the difference in hydrodynamic coefficients between the viscous and potential results. This might be implemented as nonlinear viscous damping in time domain motion analysis to account for the eddy and viscous damping.
6. Conclusion
A two-dimensional moonpool by a rectangular section under forced heave motion, including piston mode, was successfully simulated by fully nonlinear solvers for potential and viscous flows.
All numerical methods of the present study predicted the resonant frequency for different gap sizes well. The viscous flow CFD solver estimated the most reliable and acceptable resonance amplitudes compared to the experimental results. The potential flow solver tends to overestimate the free-surface elevation near the resonance frequency. Outside the resonance region, potential solutions are also reliable. In the wide-gap case, the moonpool resonance is relatively weak and their differences are not appreciable. The viscous solver produces significant differences in moonpool elevation between the narrow-gap cases of sharp and rounded corners, while the phenomenon cannot be reasonably predicted by the potential simulations.
The free-surface elevation inside the moonpool was found to be related to the hydrodynamic forces. To accurately estimate the hydrodynamic forces and hydrodynamic coefficients, the free-surface elevation should be correctly predicted first.
Vortex shedding and surface viscosity are the main sources of the damping that differentiated the viscous and potential flow solvers. The vortical flow fields at the resonant frequency clearly showed the mechanism and contribution of vortex shedding, and it should be taken into account for moonpool problems. The amount of eddy damping was found to depend on the sharpness of corners; therefore, free-surface elevation may be controlled by changing the corner shape or appended devices.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.