3D Numerical Simulations of Green Water Impact on Forward-Speed Wigley Hull Using Open Source Codes

Abstract

A series of CFD RANS simulations are presented for Wigley hulls of two freeboard heights progressing with forward speed in waves. Free surface effects are captured using the Volume of Fluid (VOF) method embedded in open source software OpenFOAM. Comparisons of heave, pitch motions and added resistance of the first Wigley model against the experiments of Kashiwagi (2013) confirm the numerical validity of the hydrodynamic modelling approach. Further simulations for the lower-freeboard Wigley model reveal that the highest green water impact on decks appears in way of $\lambda /L=1.3$ and at the highest instantaneous pitch amplitude where the water propagates far downstream and across the deck. The simulations also demonstrate that the green water events are associated with air bubble entrapment.

1. Introduction

Under harsh sea conditions, vessels or offshore platforms endure large wave loads, such as bottom and bow slamming, green water on decks, etc. Severe slamming may lead to serious structural damage and may affect the safety and stability of the vessel or offshore platform. Green water may not pose a direct threat to integrity or survivability but may result in increasing expenses for repairs and renewal of damaged topside modules or outfitting equipment [1,2].

Over the years the development of non linear numerical methods for the assessment of risks associated with impact of green water on decks has been challenging. Measuring freeboard exceedance or free surface elevation relative to the deck using 3D linear diffraction theory has been extensively applied for the prediction and statistical analysis of green water risk due to its relative simplicity [3,4,5]. In reality, a green water event is highly nonlinear. This implies that linear potential flow theory and the assumption of a Rayleigh distribution for the statistics may be inaccurate [3,6,7,8,9].

Another method used is the so-called composite method. In this numerical idealisation the green water event is divided into four different stages and each of them is modeled separately using distinct methodologies. The freeboard exceedance is used as input for a dam-break model [10,11,12,13,14] or as a solution to the shallow water equations [15,16] that can simulate hydrodynamic behaviour of on-deck flows. The green water loading on deck structures is approximated either analytically by the similarity method [17] and Wagner-type analysis [18] or empirically by utilizing extensive scale model tests [3]. Computational Fluid Dynamics (CFD) is then adopted to idealise the influence of local fluid flow in way of the localized regions around the bow or to simulate local interactions between green water flows and deck structures [3,19,20]. Uncertainties associated with modeling assumptions constrain the accuracy of this method that, in any case, could be considered more suitable for preliminary design stage.

With the successful development of high performance computing, CFD simulations using the volume of fluid (VOF) method may be considered as a valuable alternative [2,21,22,23,24,25,26,27,28]. However, to capture the nonlinear evolution of wave groups and wave-structure interactions, high mesh resolution is required in such simulations. For example, for the simulation of green water events, finer mesh is required to accommodate for the influence of shallow water hydrodynamics on deck.

This paper presents a series of CFD RANS simulations for Wigley hulls of two freeboard heights progressing with forward speed in waves. Free surface effects are captured using the Volume of Fluid (VOF) method embedded in open source software OpenFOAM. Heave, pitch motions and added resistance of the first Wigley hull model are compared against the experiments of Kashiwagi (2013) [29]. Pressures on the deck bow of the lower-freeboard Wigley hull model are measured and flow observations of green water are carried out. The study provides complementary understanding of green water dynamics on deck and key physics of the green water events.

2. Mathematical Model and Numerical Method

The flow is governed by the incompressible Navier–Stokes equations

$$\begin{array}{cc}\hfill \frac{\partial u_i}{\partial x_i}& =0,\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \frac{\partial \left(\rho u_i\right)}{\partial t}+u_j\frac{\partial \left(\rho u_i\right)}{\partial x_j}& =-\frac{\partial p}{\partial x_i}+\mu \frac{{\partial }^2{u}_i}{\partial x_j\partial x_j}+g,\hfill \end{array}$$

(2)

where $u_i$ ($i=1$, 2, 3) denotes flow velocity at the x-, y- and z-axes respectively, $\rho$ represents the mixture density of the two phases separately considered, p denotes pressure, $\mu$ is dynamic viscosity of the fluid, g represents the gravitational acceleration.

A Reynolds averaging approach is adopted to smear out the fluctuating components of the solution variables in the instantaneous Navier–Stokes equations so as to alleviate the computational cost. Averaging the continuity and Navier–Stokes equations yields

$$\begin{array}{cc}\hfill \frac{\partial {\overline{u}}_i}{\partial x_i}& =0,\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill \frac{\partial \left(\rho {\overline{u}}_i\right)}{\partial t}+{\overline{u}}_j\frac{\partial \left(\rho {\overline{u}}_i\right)}{\partial x_j}& =-\frac{\partial \overline{p}}{\partial x_i}+\mu \frac{{\partial }^2{\overline{u}}_i}{\partial x_j\partial x_j}-\frac{\partial \left(\rho \overline{u_i^{\prime }{u}_j^{\prime }}\right)}{\partial x_j}+g,\hfill \end{array}$$

(4)

where ${\overline{u}}_i$ and $\overline{p}$ represent the ensemble-averaged velocity and pressure; the product $\rho \overline{u_i^{\prime }{u}_j^{\prime }}$ denotes the Reynolds stress component. Thereafter Reynolds-averaged Navier–Stokes (RANS) equations govern the transport of flow. The Reynolds stress in Equation (4) is modeled as a function of turbulent viscosity ${\mu }_T$ and the mean velocity gradients according to Boussinesq hypothesis

$$-\rho \overline{u_i^{\prime }{u}_j^{\prime }}={\mu }_T\left(\frac{\partial {\overline{u}}_i}{\partial x_j}+\frac{\partial {\overline{u}}_j}{\partial x_i}\right)-\frac{2}{3}\rho k{\delta }_{ij},$$

(5)

where k is is the turbulent kinetic energy and ${\delta }_{ij}$ is the Kronecker delta.

To obtain the turbulent viscosity, the shear-stress transport (SST) k-$\omega$ two-equation model proposed by Menter [30] is used as follows

$$\begin{array}{cc}\hfill \frac{\partial \left(\rho k\right)}{\partial t}+\frac{\partial \left(\rho k{\overline{u}}_i\right)}{\partial x_i}& =\frac{\partial }{\partial x_i}\left[(\mu +{\mu }_T)\frac{\partial k}{\partial x_i}\right]+\min ({\mu }_T{S}^2,10{\beta }^{*}k\omega )-\rho {\beta }^{*}k\omega ,\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \frac{\partial \left(\rho \omega \right)}{\partial t}+\frac{\partial \left(\rho \omega {\overline{u}}_i\right)}{\partial x_i}& =\frac{\partial }{\partial x_i}\left[(\mu +{\mu }_T)\frac{\partial \omega }{\partial x_i}\right]+\alpha \rho S^2-\rho {\beta }^{*}{\omega }^2+2(1-F_1)\rho {\sigma }_{\omega ,2}\frac{\partial k}{\partial x_i}\frac{\partial \omega }{\partial x_i},\hfill \end{array}$$

(7)

where $\omega$ represents the specific dissipation rate,

$${\mu }_T=\frac{\rho k}{\omega }/\max \left[\frac{1}{{\alpha }^{*}},\frac{S{F}_2}{{\alpha }_1\omega }\right],$$

(8)

and $F_1$ and $F_2$ are the blending functions [30].

The volume of fluid (VOF) method with artificial bounded compression technique proposed by Hirt and Nichols (1981) [31] was used to capture the free surface. The free surface was considered as a mixture of water and air and accordingly the VOF transport equation was expressed as

$$\frac{\partial \alpha }{\partial t}+\nabla ·\left({\mathit{U}}_r(1-\alpha )\alpha \right)=0,$$

(9)

where ${\mathit{U}}_r={\mathit{U}}_{water}-{\mathit{U}}_{air}$ is the relative velocity between the water and air and termed as “compression velocity”; $\alpha$ represents the volume fraction defined as the relative volume proportion of water in a cell, namely

$$\left\{\begin{array}{cc}\alpha =0\hfill & \mathrm{air}\hfill \\ 0<\alpha <1\hfill & \mathrm{interface}\hfill \\ \alpha =1\hfill & \mathrm{water}\hfill \end{array}\right.$$

(10)

The spatial variation of fluid density and dynamic viscosity can be then expressed as

$$\begin{array}{c}\hfill \rho =\alpha {\rho }_{water}+(1-\alpha ){\rho }_{air},\end{array}$$

(11)

$$\begin{array}{c}\hfill \mu =\alpha {\mu }_{water}+(1-\alpha ){\mu }_{air},\end{array}$$

(12)

where the subscripts “water” and “air” represent the corresponding fluid property of water and air respectively.

Rigid body motion equations, in which all other motions are constrained apart from heave and pitch, were used to model the dynamic response of the ship as follows

$$\begin{array}{c}\hfill m\ddot{z}=F_z,\end{array}$$

(13)

$$\begin{array}{c}\hfill I_{yy}\ddot{\theta }=N,\end{array}$$

(14)

where m denotes the ship mass; z is displacement in heave; $F_z$ represents the total force in the vertical direction; $I_{yy}$ denotes moment of inertia in pitch; $\theta$ represents angle in pitch, N is moment of force in pitch. The Newmark-$\beta$ method [32] with $\gamma =0.25$ and $\beta =0.5$ was used to solve Equations (13) and (14).

Mesh deformation was achieved by moving the node points of the mesh according to the ship motion and without changing the mesh topology. The displacements of the node points were then calculated by using the equation

$$\nabla ·(\frac{1}{r}{D}_i)=0$$

(15)

where $D_i$ represents the displacements of the node points; r denotes the distance of the cell centre to the nearest moving wall boundary.

3. Problem Setup

At first instance Equation (9) was solved for a fraction of volume and then integrated across the entire flow field. Then, the Reynolds-averaged Navier–Stokes equations using the SST $k-\omega$ model were solved by the OpenFOAM interFoam subroutine as follows

$$\frac{2y}{B}=\left[1-{\left(\frac{z}{d}\right)}^2\right]\left[1-{\left(\frac{2x}{L}\right)}^2\right]\left[1+0.6{\left(\frac{2x}{L}\right)}^2+{\left(\frac{2x}{L}\right)}^4\right]+{\left(\frac{z}{d}\right)}^2\left[1-{\left(\frac{z}{d}\right)}^8\right]{\left[1-{\left(\frac{2x}{L}\right)}^2\right]}^4,$$

(16)

where B denotes the waterline breadth; L denotes the length between perpendiculars; d represents the draft. To observe the green water occurrence, two blunt Wigley hull models with different freeboard heights were considered (see Figure 1 and

Table 1. Principle parameters of the ship model.

Length between perpendiculars, L (m)	2.5
Waterline breadth, $B$ (m)	0.5
Draft, $d$ (m)	0.175
Displacement volume, V (m${}^3$)	0.139
Height of gravity center, $KG$ (m)	0.145
Gyration radius in pitch, $k_{yy}/L$	0.236
Freeboard height, f (m)	0.175 (model 1), 0.0875 (model 2)

).

A top view of the layout of the computational domain according to Chen et al. (2019) [2] is shown in Figure 2. The first-order Stokes regular head wave was generated in the wave tank for Froude number ($Fn$) of 0.2. Assuming that the wave propagates along the negative x-axis direction, the wave elevation $\zeta$, the encounter frequency ${\omega }_e$ and the incident wave potential $\phi$ were written as

$$\begin{array}{cc}\hfill \zeta & =Acos(kx+{\omega }_{e}t),\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill {\omega }_e& =\omega +kU,\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill \phi & =\frac{gA}{2\omega }\frac{cosh\left[k\right(z+H\left)\right]}{sinh\left(kH\right)}(kx+{\omega }_{e}t),\hfill \end{array}$$

(19)

the velocity components of the water were then defined as follows

$$\begin{array}{cc}\hfill u_x(x,y,z,t)& =U+A{\omega }_e{e}^{kz}cos(kx-{\omega }_{e}t+\theta ),\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill u_y(x,y,z,t)& =0,\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill u_z(x,y,z,t)& =A{\omega }_e{e}^{kz}sin(kx-{\omega }_{e}t+\theta ),\hfill \end{array}$$

(22)

where A represents the wave amplitude; k is the wave number; $\omega$ denotes the natural frequency; U represents the ship speed; H is the water depth.

For all cases, the free surface was set at $z=0$ with the draft $d=0.175$. Equations (20)–(22) were then used to define the boundary condition at the inlet. At the outlet, a damping zone of wave length $\lambda$ was set to damp the vertical motions of an interface by applying a force ($d\left(m{u}_z\right)/dt=-\kappa m{u}_z$) to the momentum equation proportional to the momentum of the flow in the direction of gravity. Fluid damping effects lead to exponential decay of the vertical velocity ($u_z=u_{z0}{e}^{-\kappa t}$), where $\kappa$ is set based on the desired level of damping and the residence time of a perturbation through the damping zone. At the mid plane of the hull, a symmetry boundary condition was used.

Unstructured mesh including hexahedral and tetrahedral cells was generated using snappyHexMesh tool of OpenFOAM. Top view and side view of the mesh in the computational domain are shown in Figure 3. Refinements with 5 different levels were used to refine the mesh around the hull and free surface. The finest mesh size was set to satisfy $\lambda /\Delta x>250$ in the longitudinal direction and $d/\Delta z=20$ in the vertical direction. In addition, five layers of thin cells were added for the boundary layer around the ship hull. The first layer of cells around the ship hull was set to satisfy the nondimensional thickness of the cell ${\eta }^{+}\approx 1000$, where ${\eta }^{+}=U\eta /\nu$ (U is the ship speed, $\eta$ denotes the thickness of the first layer of cell,$\nu$ represents kinematic viscosity of the fluid). Meanwhile, standard wall functions were applied to computing ${\nu }_T$ at the first off-wall nodes (${\nu }_T={\mu }_T/\rho$ represents turbulence kinematic viscosity).

4. Results and Discussions

4.1. Code Validation

To carry out a validation for the wave model, the same wave parameters to Kashiwagi (2013) [29] were used (see

Table 2. Wave parameters for the simulations.

Froude Number,	Wave Amplitude,	Wave Length,	Wave Steepness,	Encounter Frequency,
$\mathit{Fr}$	$\mathit{A}/\mathit{L}$	$\mathit{\lambda }/\mathit{L}$	$\mathit{H}/\mathit{\lambda }$	${\mathit{\omega }}_{\mathit{e}}$
		0.5	0.048	11.998
		0.7	0.0343	9.49
		0.9	0.02667	7.996
0.2	0.012	1.1	0.0218	6.994
		1.3	0.01846	6.267
		1.6	0.015	5.48
		2.0	0.012	4.75

). Seven cases with different wave lengths were then implemented for model 1. The nondimensional total resistance of the hull was defined as

$$C_t=\frac{{\overline{R}}_{wave}}{\rho g{A}^2{B}^2/L},$$

(23)

and was computed using 1.51, 2.33 and 3.02 million cells of meshes and three time steps ($\Delta t$) to test the grid and time step convergence (see Figure 4). In Equation (23) ${\overline{R}}_{wave}$ denotes the mean resistance of the hull in waves. The figure shows the convergence and a mesh with about 2.33 million cells and $\Delta t=0.00025$ was adopted for all the further simulations.

In the simulation of each case, 15 encounter periods ($T_e=2\pi {\omega }_e$) of computation, which takes about 50 CPU hours of high performance computing with eight Intel Xeon Gold processors, were undertaken and 10 periods of data were used to calculate the results of the wave-induced motion and resistance after the system reached a steadily periodic state. Time histories of heave and pitch were expanded into Fourier series

$$\zeta \left(t\right)=\frac{a_0}{2}+\sum _{n=1}^N{a}_{n}cos(n{\omega }_{e}t+{\theta }_n)$$

(24)

where $a_n$ is the nth harmonic amplitude and ${\theta }_n$ is the corresponding phase angle. The nth dimensionless harmonic amplitudes of the heave (${\zeta }_3^{\left(n\right)}=a_3^{\left(n\right)}/A$) and pitch (${\zeta }_5^{\left(n\right)}=a_5^{\left(n\right)}/\left(kA\right)$) were then obtained using Fourier transformation. Profiles of the first harmonic amplitude of heave ${\zeta }_3^{\left(1\right)}$ over a range of dimensionless wave length ($\lambda /L$) are shown in Figure 5. Experimental data, computed results by Enhanced Unified Theory (EUT), New Strip Method (NSM) and Rankine Panel Method (RPM) originally introduced by Kashiwagi (1995 [33], 2013 [29]), Iwashita and Ito (1998) [34] are also included in this comparison. Both EUT and RPM are frequency-domain linear calculation methods, in which the disturbed flow is divided into the steady (${\varphi }_s$) and unsteady (${\varphi }_u$) flows. In EUT, the steady disturbance potential ${\varphi }_s$ is ignored. In RPM it is computed numerically as double-body flow and its effects on the body and free-surface boundary conditions are considered. In NSM hydrodynamic actions on and motions of the ship are determined by integrating the results from two-dimensional hydrodynamic coefficients over the ship length. As shown in Figure 5, EUT and NSM are not in good agreement with the experimental results. In general comparisons against RPM are are good with the exception of $\lambda /L=1.3$.

Variations of the first harmonic amplitude of pitch ${\zeta }_5^{\left(1\right)}$ with the wave length are shown in Figure 6. Both EUT and NSM results do not coincide well with the experiments and visible underpredictions appear at large $\lambda /L$. Results obtained by the RPM model are less accurate than those obtained by the novel CFD method presented in this paper for $\lambda /L<1.3$.

Added resistance is defined as a difference between the mean (time-averaged) resistance in waves and in calm water at the same speed U. The nondimensional added resistance was defined as

$$C_{aw}=\frac{{\overline{R}}_{wave}-{\overline{R}}_{cw}}{\rho g{A}^2{B}^2/L},$$

(25)

where ${\overline{R}}_{wave}$ and ${\overline{R}}_{cw}$ represent the mean resistance of the hull in waves and in calm water respectively; $\rho$ is the density of fluid; g is the gravitational acceleration; A is the wave amplitude; B and L are waterline breadth and length between perpendiculars of the ship respectively. Comparisons of the distributions of the nondimensional added resistance obtained by experiments, EUT and CFD are shown in Figure 7. It is shown that EUT underpredicts added resistance. CFD results are generally in better agreement with experiment. However, differences at low $\lambda /L$ become more obvious.

4.2. Green Water on Decks

To observe the green water events, simulations associated with seven wave lengths were implemented in Model 2 (see Figure 1). The impact of green water on decks was examined by measuring the pressure at five monitoring points with coordinates (1.25, 0), (1.1, 0), (1.0, 0), (0.9, 0) and (0.8, 0) on the deck along the centreline of the bow (see Figure 8).

In green water events free surface exceeds the ship freeboard and hydrodynamic flow propagates over the deck with high-velocity. In way of impact the fluid flow pressure sharply increases. Figure 9 shows the variation of time histories of pressure coefficients over nondimensional time ($t/T_e$) in way of points 1 and 2 (see Figure 8). For each case the pressure coefficient was defined as

$$C_p=\frac{p}{0.5\rho U^{2}LB}.$$

(26)

It was found that the green water impact at point 1 increased from $\lambda /L=0.5$ to 0.7 and decreased from $\lambda /L=0.7$ to 0.9. On the other hand, the pressure at point 2 is larger than that in way of point 1 at $\lambda /L=0.9$. This is because at $\lambda /L=0.9$ green water becomes stronger and the water flow over the freeboard rapidly propagates downstream while an air layer is formed underneath. Eventually, an air bubble is entrapped inside the on-deck water. The pressure at point 2 increases from $\lambda /L=0.9$ to 1.3 and reaches its highest value at $\lambda /L=1.3$. This implies that the highest influence of green water on decks occurs at $\lambda /L=1.3$. The green water event then starts to decline until it disappears at $\lambda /L=2.0$. Occurrence of the strongest green water impact is associated with the highest 1st harmonic amplitudes of the pitch.

Figure 10 shows the pressure coefficient distributions over the time at five points for $\lambda /L=0.9$ and 1.3. The profiles show that for $\lambda /L=0.9$ pressures at all the monitoring points are smaller than those of $\lambda /L=1.3$. This confirms that the green water at $\lambda /L=0.9$ is relatively weaker than at $\lambda /L=1.3$. This demonstrates that while a large amount of water exceeds the freeboard and keeps propagating along the deck, pressures remain high at the five monitored points. The time lag in pressure peaks experienced between point 1 and points 2 to 5 reflects that water first impacts point 1, then impacts point 2 as it climbs up the freeboard. The water rapidly spreads on deck (points 2–5) without major pressure fluctuations.

4.3. Flow Visualization

Flow visualization was achieved by plotting velocity contours on the free surface and pressure contours in way of the x-z plane at $y=0$ at seven instants in one motion period corresponding to $\lambda /L=1.1$. Free surface was defined at $\alpha =0.5$. As shown in Figure 11 and Figure 12 on the top left side of each flow velocity contour, the pressure profile at point 1 is attached to highlight the pressure magnitude. At first instance, the hull is almost at a horizontal state and the pressure at point 1 is low (see Figure 11a,b). At the second instant, a trim by the bow starts and wave-structure interaction causes climb up of the water above the bow height (see Figure 11c,d). Consequently, water impact with high-velocity flow causes sharp increase of pressure at point 1. At the third time instant, it can be seen that the pressure at point 1 decreases (see Figure 11e,f). This is because the water flow has very high velocity. As the vessel trims by the bow, a low pressure air layer grows between the water and the deck (see Figure 11f). At the fourth time instant, a trim by the aft starts and the wave-structure interaction becomes stronger (see Figure 11g,h). This causes the pressure at point 1 to reach its highest magnitude. When the trim by aft starts, the water over the deck rolls down and entraps the air underneath to form an air bubble in a form of a low-pressure zone on the deck bow (see Figure 11h). At the fifth time instant, the air bubble propagates downstream along the deck and air bubbles are formed downstream (see Figure 12a,b). At the sixth time instant, the hull trims by the aft, the bow is no longer submerged in waves and the on-deck water is starting to decline, thus leading to decrease of the pressure at point 1 (see Figure 12c,d). At the seventh instant, the volume of water at the deck bow continues to decline and some of on-deck water propagates downstream on the deck. This leads to lower pressure at the other four monitoring points (see Figure 12e,f). The free surface ahead of the bow is falling and the pressure at point 1 approaches its lowest magnitude.

4.4. Motion Quantities and Added Resistance

During the green water event, the wetted surface also decrease at the moment of occurrence of a trim by aft. This may lead to changes in the hydrodynamic forces on the hull. The green water impact would also have an effect on the hydrodynamic forces on the hull, variations in motions and added resistance. Figure 13 illustrates variations of heave, pitch and added resistance for both Wigley models (see Figure 1). It may be concluded that because of the higher free-board height of Model 2 in comparison to Model 1 the impact of green water on decks is low and does not influence that much the first harmonic heave/pitch amplitudes at $\lambda /L=0.5$, 0.7 and 1.6. For both models from $\lambda /L=0.9$ to 1.3, the green water event remains strong, thus both the decrease of the wetted surface and the green water impact cause decreases of the first harmonic amplitudes of the heave and pitch as well as added resistance.

5. Conclusions

This paper presented 3D RANS CFD simulations for low and high freeboard Wigley hull models. Free surface effects were solved by the VOF method in the openFOAM numerical wave tank. Comparisons between computations and experimental results helped to validate fluid structure interaction modeling procedures, realise fluid flow phenomena and pressure peaks of relevance to the phenomenon of green water on decks. Results revealed that CFD simulations provide improved physics representation. The highest green water impact appears at $\lambda /L=1.3$ at which the highest pitch is also present. For $\lambda /L=1.1$ to 1.3, the water propagates far downstream and an air bubble is entrapped in a form of a small low-pressure zone.