Numerical Study on the Green-Water Loads and Structural Responses of Ship Bow Structures Caused by Freak Waves

Abstract

In recent decades, freak waves, characterized by their unusual high amplitude, sharp crest, and concentrated energy, have attracted researchers’ attention due to their potential threat to marine structures. Green-water loads caused by freak waves can be significant and may lead to local damage to the ship structures. Therefore, this paper focuses on the study of green-water loads and examines the structural responses of ship bow structures under the influence of the green-water loads caused by freak waves. Firstly, a three-dimensional numerical wave tank is established in which the superposition model is used to generate freak waves. Validations on the freak-wave generation, ship motion response and the wave loading are carried out to verify the present solvers. The simulation on the interaction between the freak wave and the ship are conducted to obtain the interaction process and green-water loads. Secondly, a finite element (FEM) model of the ship bow is built, on which the green-water loads are applied to calculate the structural responses. Finally, the displacement and stress of the deck and breakwater structures are analyzed. It is found that green water events caused by freak waves can generate enormous impact forces on the bow deck and breakwater, resulting in severe structural responses and even possible damage to the structures. The local strength of structures under freak waves needs to be considered in practical engineering applications.

1. Introduction

A ship is one main way of the transportation of goods in international trade. However, during the transportation process, ships may encounter various adverse sea conditions that endanger their safety. Freak waves, also known as rogue waves, which are rare and unpredictable waves characterized by an unusual high amplitude, sharp crest, and concentrated energy [1,2,3], are one of the most dangerous. The interaction between ships and freak waves can result in complex nonlinear phenomena including slamming and green water on deck, which may lead to severe damages on the ship structures, e.g., the ‘World Glory’ accident [4], the ‘SS El Faro’ accident [5] and many other destructions [6,7]. Therefore, investigation on the impact loads and structural response of the ship structures is crucial for ship safety.

In recent decades, many studies have been conducted to investigate the effects of freak waves on marine structures in order to understand the damage mechanisms caused by these waves. Some scholars have studied the wave force induced by freak waves. For example, Sparboom et al. [8] examined the impact of freak waves on cylinders in a large-scale physical model test and studied the influence of tilt angles. They found that the maximum impact pressure of freak waves on cylinders occurred when the incoming wave direction was opposite to the cylinder’s tilt direction. Corte and Grilli [9] utilized a wave superposition model to simulate freak waves and investigated the impact loads of freak waves on vertical rigid cylindrical piles. Kim et al. [10] investigated the impact pressure of Draupner freak waves and fifth-order Stokes waves on cylinders and found that the impact load of the freak wave was 2.8 times that of regular waves of the same scale. Li et al. [11] used an energy-focusing model in the laboratory to simulate freak waves and studied the influence of wave parameters such as amplitude, period, and spectral width on the impact load of freak waves on cylinders. Bunnik et al. [12] conducted numerical simulations of fixed offshore platforms subjected to freak-wave impacts using a numerical wave tank. They predicted the impact loads on offshore platforms caused by freak waves and designed model experiments based on simulated conditions. The experimental results were in good agreement with those of the numerical simulations. In the academic circle, focused wave groups are frequently used to generate the freak waves, and as such, some scholars have researched the interaction between focused wave groups and marine/coastal structures. Gao et al. [13] utilized OpenFOAM to examine the transient fluid resonance phenomenon occurring in a narrow gap between two adjacent boxes that were excited by incident focused waves with varying spectral peak periods and focused wave amplitudes. Liu et al. [14] investigated the transient fluid motion within a narrow gap formed by two fixed boxes under the action of an incident focused wave group, using a two-dimensional viscous flow numerical wave flume. In the study of Gao et al. [15], based on the Morlet wavelet transform and discrete Fourier transform techniques, the capability of focused transient wave groups to trigger the harbor resonance phenomenon was revealed for the first time. Some scholars have investigated the motion response and longitudinal strength of marine structures under freak waves. For example, Liu et al. [16] described how to evaluate the strength of a ship based on its nonlinear vertical bending moment (VBM), and investigated the influence of freak-wave height and speed on VBMs and deformation. Soares et al. [17] investigated the VBM in the midship induced by freak waves and explored the influence of the position in space where the freak waves were generated on the variation of the maximum bending moment. Qin et al. [18] studied the impact of freak waves and second-order Stokes waves on a fixed elastic deck using a strong coupled fluid–structure interaction method and discussed the hydroelastic responses. Luo et al. [19] conducted physical experiments to investigate the impact of a freak wave on a tension-leg platform, which suggested that a high-crest freak wave could induce violent motions of a floating platform and snap loads in tethers. Rudman et al. [20] utilized the smoothed particle hydrodynamic (SPH) method to simulate the effect of freak waves on a semi-submersible platform and compared and analyzed the motion response characteristics of the semi-submersible platform under two mooring systems: a tension leg platform (TLP) and taut spread mooring (TSM). Focusing on the local wave-structure interaction phenomena, freak waves with huge wave heights usually lead to significant green water events. Therefore, some other scholars have examined the green-water loads caused by freak waves. Hu et al. [21] found asymmetric and irregular characteristics in green water caused by freak waves when compared to regular waves. Qin et al. [22] studied the impact of green-water loads caused by Peregrine breather-type freak waves on a flat deck, and proposed an empirical formula for the quick prediction of nonlinear freak-wave forces on such deck structures. Zhang et al. [23] studied the wave height, deck pressure, and superstructure pressure of a fixed FPSO under the influence of freak waves. Liu et al. [24] investigated the mechanism and impact form of green water events under the influence of freak waves, and analyzed the ship’s motion response, water variation and pressure. Wang et al. [25] discussed the relationship between ship motion response and green water events under the influence of freak waves, and the effect of ship speed on wave-ship interaction, revealing the impact of freak-wave peaks and sequences on roll, heave, and impact pressure.

The serious wave loads during green water events may lead to severe structural responses and even local damages of the ship bow, especially to the deck and breakwater structures. However, previous research focused on the structural analysis during bottom slamming, bow flare slamming and stern slamming. The structural responses and local damage of the deck and breakwater structures were rarely examined. For example, Maki et al. [26] used computational fluid dynamics (CFD) and FEM methods to study the water entry slamming of the ship’s body and predict the structural hydroelastic response. Yang et al. [27] established a three-dimensional FEM model of a 1700 TEU container ship’s bow structure and studied the dynamic response of bow flare slamming under impact loads, using a one-fold-plate-thickness deformation criterion to determine whether the structure buckled. Ren et al. [28] studied the dynamic response of bow flare structures under impact loads using LS-DYNA software. Xie et al. [29] proposed an uncoupled CFD and FEM method to study the local dynamic response of bow flare on a 21,000 TEU container ship under extreme impact loads. Kim et al. [30] studied the slamming response of the stern of an LNG ship using numerical simulation. Qin et al. [31,32] studied the impact of green-water loads caused by freak waves on the deck and deck-house structures, which further obtained the displacement of the structures using a strong coupled fluid–structure interaction method. However, the structures they used were fixed deck models simplified as plates and beams.

From these studies, it is seen that the interaction between freak waves and marine structures has been widely studied in the aspects of wave force, motion response, longitudinal strength and green water events. Nevertheless, the green-water loads of a ship with free motions induced by freak waves and the corresponding structural analysis of the deck and breakwater structures have rarely been examined. Therefore, this paper provides a detailed study on the green-water loads induced by the freak wave and the corresponding structural responses of the deck and breakwater structures of a ship bow. To achieve this, firstly, a three-dimensional numerical wave tank is established in which the superposition model is used to generate freak waves. Validations on the freak-wave generation, ship motion response and wave loading are carried out to verify the present issues. The simulation of the interaction between the freak wave and the ship are conducted to obtain the interaction process and green-water loads. Secondly, a finite element model of the ship bow is built, on which the green-water loads are applied to calculate the structural responses. Finally, the displacement and stress of the deck and breakwater structures are analyzed, while conclusions on the safety of the ship bow is drawn. The novelty of the study lies in that a comprehensive investigation of the ship’s safety under freak-wave conditions are conducted, including the motion response, green-water impact loading and the structural response. More importantly, the structural responses and damages of the deck and breakwater structures during the green-water event induced by the realistic scale freak wave are examined for the first time, which provides a meaningful reference value for the practical engineering design of ship deck and breakwater structures.

2. Numerical Implementation

2.1. Freak Wave Based on Superposition Model

The superposition model, commonly employed as a freak-wave model, has been shown to be effective in generating freak waves experimentally and numerically. In the present study, a two-wave-train superposition model is adopted, which combines a random wave train and a convergent wave train to form the wave surface elevation. By adjusting the two trains of the waves’ energy ratio and the phase of the linear cosine component waves of the convergent wave train, a freak wave is produced at a particular time and location [33]. The wave surface elevation of the two-wave-train superposition method can be expressed as:

$$\zeta ={\displaystyle \sum _{1i=1}^{N_1}{a}_{1i}\cos \left(o_{1i}x-w_{1i}t+{\theta }_{1i}\right)}+{\displaystyle \sum _{2i=1}^{N_2}{a}_{2i}\cos \left[o_{2i}\left(x-x_c\right)-w_{2i}\left(t-t_c\right)\right]}$$

(1)

In Equation (1), the first item creates the background random waves, while the second item generates the freak-wave peaks. $N_1$ and $N_2$ represent the respective quantities of component waves present in the two wave trains. $o_{1i}$ and $o_{2i}$ represent the wave numbers of individual component waves within the two wave trains. $w_{1i}$ and $w_{2i}$ denote the circular frequency of each component wave within the two wave trains. ${\theta }_{1i}$ determines the phase of each component wave in the random wave train, while $o_{2i}{t}_c-o_{2i}{x}_c$ determines the phase of each component wave in the convergent wave train. $t$ and $x$ are the time and coordinate, and $x_c$ and $t_c$ are the position and focusing time when the freak-wave crest is formed. $a_{1i}$ and $a_{2i}$ are the amplitudes of the individual component waves within the two wave trains.

2.2. Fluid Solver

To simulate the interaction between freak waves and the ship, a numerical wave tank is required. This study utilizes a three-dimensional fluid solver to model this interaction. The incompressible N-S equations and a VOF method were used by the solver to reconstruct free surfaces. The fluid was assumed to be viscous, Newtonian, and incompressible, with governing equations consisting of continuity, momentum conservation, and volume transportation equations, written as follows:

$$\frac{\partial {\mathit{u}}^f}{\partial x}=0$$

(2)

$$\frac{\partial {\mathit{u}}^f}{\partial t}+{\mathit{u}}^f\frac{\partial {\mathit{u}}^f}{\partial x}=\frac{1}{{\rho }^f}\frac{\partial {\mathbf{\sigma }}^f}{\partial x}+{\mathit{f}}^f+\mathit{\upsilon }\frac{{\partial }^2{\mathit{u}}^f}{\partial x^2}$$

(3)

$$\frac{\partial F}{\partial t}+\frac{\partial {\mathit{u}}^{f}F}{\partial x}=0$$

(4)

where $t$ and $x$ denote the time and coordinate. ${\mathit{u}}^f$, ${\rho }^f$, ${\mathit{f}}^f$, ${\mathsf{\sigma }}^f$, $\mathit{\upsilon }$ and $F$ are the velocity, density, body force, Cauchy stress tensor of the fluid, kinetic viscosity coefficient and the transportation volume of the fluid. The standard $k-\omega$ model [34] is applied for the solution of the governing equations.

Omitting the surface tension, boundary conditions at the free surface are given as follows:

$$\frac{\partial u_n}{\partial \tau }+\frac{\partial u_{\tau }}{\partial n}=0$$

(5)

$$-p+2\mu \frac{\partial u_n}{\partial n}=-p_0$$

(6)

where $u_n$ and $u_{\tau }$ are the normal component and the tangential component of the velocity vector on the boundary, respectively, and $p$, $p_0$ and $\mu$ are the liquid pressure, air pressure and dynamic viscosity coefficient.

The governing equations are solved using the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) algorithm [35]. The solution of ship motion response is obtained by using the rigid-body six-degree-of-freedom motion equation and employing the FAVOR method to handle the dynamic boundary (Wang et al. [25], Patankar and Spalding [36]; Hirt and Sicilian [37]).

2.3. Structure Solver

According to Newton’s second law, the rate of change of structural momentum is equal to the external force load acting on the structure. Therefore, the momentum equation for the structure can be written as:

$${\rho }^s\ddot{\mathit{z}}=\nabla \cdot {\mathsf{\sigma }}^s+{\mathit{f}}^s$$

(7)

where ${\rho }^s$ is the structural density, $\mathit{z}$ is the displacement, ${\mathsf{\sigma }}^s$ is the first Piola–Kirchhoff stress tensor, and ${\mathit{f}}^s$ is the volume force acting on the structure. The expression for ${\mathsf{\sigma }}^s$ is given by the following equation:

$${\mathsf{\sigma }}^s={\mu }^s\left[\nabla z+{\left(\nabla z\right)}^T\right]+{\lambda }^s\left(\nabla \cdot z\right)\mathit{I}$$

(8)

where ${\mu }^s$ and ${\lambda }^s$ are the Lame Constants, $\mathit{I}$ is the identity tensor.

The structural responses of the deck and breakwater structures are calculated by discretizing the momentum conservation equation for the structure, using FEM for spatial discretization and the $Newmark-\beta$ method [38] for temporal discretization. The discrete form of the structural dynamic equation is:

$$\mathit{M}\cdot \ddot{\mathit{z}}+\mathit{C}\cdot \dot{\mathit{z}}+\mathit{K}\cdot \mathit{z}=\mathit{f}$$

(9)

where $\mathit{M}$ is the mass, $\mathit{C}$ is the damping and $\mathit{K}$ is stiffness matrix matrices of the structure. $\mathit{f}$ is the force vector.

Time can be discretized using the $Newmark-\beta$ method [38], which involves the calculation of displacement, velocity, and acceleration within a specific time interval.

$${\dot{\mathit{z}}}^{\mathit{n}\mathbf{+}\mathit{1}}={\dot{\mathit{z}}}^{\mathit{n}}+\left[\left(1-\beta \right){\ddot{\mathit{z}}}^{\mathit{n}\mathbf{+}\mathit{1}}+\beta \cdot {\ddot{\mathit{z}}}^{\mathit{n}\mathbf{+}\mathit{1}}\right]\Delta t$$

(10)

$${\mathit{z}}^{\mathit{n}\mathbf{+}\mathit{1}}={\mathit{z}}^{\mathit{n}}+{\dot{\mathit{z}}}^{\mathit{n}}\Delta t+\left[\left(\frac{1}{2}-\gamma \right){\ddot{\mathit{z}}}^{\mathit{n}}+\gamma \cdot {\ddot{\mathit{z}}}^{\mathit{n}\mathbf{+}\mathit{1}}\right]\Delta t^2$$

(11)

where $\Delta t$ is time step. $\beta$ and $\gamma$ are parameters that determine the accuracy and stability of the equation. Using Nastran as the structural solver in this paper, which is a widely used finite element analysis software, various dynamic loads can be defined in the time domain, and its powerful analytical capabilities have been extensively applied to structural dynamics analysis [39].

The structural boundary conditions need to satisfy the continuity of velocity and surface force at the boundary, as shown below:

$${\mathit{u}}^s={\mathit{u}}_{boundary}^s$$

(12)

$${\mathsf{\sigma }}^s{\mathit{n}}^s={\mathit{T}}_{boundary}^{}$$

(13)

here, ${\mathit{u}}_{}^s$ is the velocity of the structure, ${\mathit{u}}_{boundary}^s$ and ${\mathit{T}}_{boundary}^{}$ represent the velocity and surface force vector, respectively, given in the outer domain of the structure at its boundary, and ${\mathit{n}}^s$ is the unit normal vector at the structural boundary.

3. Numerical Validations

3.1. Validation on the Ship Motion Responses

The container ship model chosen for this validation is developed at the Institute of Ship Technology, Ocean Engineering and Transport Systems (Duisburg-Essen, Germany) as a benchmark for numerical methods by Moctar et al. [40]. The model has been widely used by scholars [41,42,43], including Mei et al. [43] who conducted physical experiments and numerical simulations to investigate the motion of this ship model in regular waves. Figure 1 presents side views of this ship;

Table 1. Main particulars of the model.

Designation	Scale
Length (m)	3.984 ± 0.001
Breadth (m)	0.572 ± 0.001
Draft (m)	0.163 ± 0.001
Displacement (kg)	245.155
Block coefficient (-)	0.661
Froude number (-)	0.139

lists its principal particulars. The regular wave is defined by a height of 0.06235 m and a period of 1.38 s.

Figure 2 illustrates the comparison between the heave and pitch motions of the container ship obtained using the present numerical method and experimental and numerical measurements conducted by Mei et al. [43]. Due to the reflection of the beach and wave maker [44], as well as the shallow-water squatting phenomenon, there are some discrepancies between the numerical simulation results and the experimental results. However, good consistency was found when comparing them with the numerical results. The simulation results obtained from the present method were in strong agreement with the numerical simulations conducted by Mei et al. [43]. This indicates that both the present solver and meshing have the ability to produce precise simulations of a ship’s heave and pitch motions in waves.

3.2. Validation of the Green Water Loading

Rosetti et al. [45] conducted research on the green-water loads of various regular waves using a fixed FPSO model in a wave tank at the University of Sao Paulo. The experimental investigation was conducted in the wave flume of the Department of Naval Architecture and Ocean Engineering, University of São Paulo. As shown in Figure 3, the FPSO model employed in the experiment had a slightly inclined deck to facilitate water drainage after green-water events. A wave probe (WM01) was located on the upstream surface of the deck; the load cell used for measuring the impact force was located on a 20 mm × 20 mm plate. The regular wave height was 0.088 m, the wave period was 0.97 s, and the wavelength was 1.46 m.

Figure 4 illustrates the comparison between the water elevations detected by WM01 and impact forces among the results obtained from the present solver, and the experimental data by Rosetti et al. [45]. They noted that CFD simulations tend to slightly overestimate the maximum value, which could result from monitoring errors during testing, but such results remain satisfactory. These findings indicate that the present simulation methodology is highly effective in accurately simulating green-water events and loads.

4. Wave-Ship Interaction and Green-Water Loads

4.1. Numerical Wave Tank

A three-dimensional numerical wave tank was constructed to simulate the interaction between container ships and freak waves. As depicted in Figure 5, the wave tank features a length of 1000 m, width of 300 m, and water depth of 180 m. A container ship model, shown in Figure 6, was situated 150 m from the left end of the tank. It should be mentioned that an on-deck breakwater structure was considered in this ship model, which is usually used for on-deck water blocking. A Piston wave-maker was located at the left boundary to generate freak waves, while a wave absorber zone was positioned near the right boundary to eliminate wave reflections based on the sponge layer relaxation method proposed by Mayer et al. [46]. As the present study focused solely on the head wave condition, the ship was limited to a rigid container ship model with heave and pitch motions, which is similar to the condition of Zhao and Hu [47].

4.2. Generation of the Freak Wave

The target freak wave was similar to the one recorded at Yura fishery in the Japanese sea area by Mori et al. [48]. The measured freak wave had a height of 11.2 m, crest height of 7.4 m, and significant wave height of 4.6 m, occurring in water with a depth of 43 m. Considering that the container ship model used in present study possesses dimensions of 398 m in length, 57 m in width, and 34 m in height, the freak wave was scaled with a scaling factor of 1:3.3. A wave surface elevation monitoring point was set 150 m away from the wave-making boundary to record the wave surface elevation at that location within the numerical wave tank. Subsequently, the free surface elevation of the freak wave was transformed into the displacement of the Piston. As mentioned by Cui et al. [49], the displacement $S\left(t\right)$ of the Piston wave-maker for freak waves can be expressed as:

$$S\left(t\right)=\frac{{\displaystyle \sum _{1i=1}^{N_1}{a}_{1i}\cos \left(o_{1i}x-w_{1i}t+{\theta }_{1i}\right)}}{W_{1i}}+\frac{{\displaystyle \sum _{2i=1}^{N_2}{a}_{2i}\cos \left[o_{2i}\left(x-x_c\right)-w_{2i}\left(t-t_c\right)\right]}}{W_{2i}}$$

(14)

$$W_i=\frac{4{\sinh }^2\left(o_{i}d\right)}{{\sinh }^2\left(o_{i}d\right)+2o_{i}d}$$

(15)

where $W_i$ is the transfer function associated with the ith component of the propagating wave, $d$ is the water depth, while the other variables have the same definitions as in Equation (1). Additionally, the theoretical value of the surface elevation at this location was calculated using Equation (1) as a reference for comparison with the numerical simulation results.

Figure 7a displays the free surface elevations of freak waves generated using three different grid groups: coarse, medium, and fine. The coarse grid had a size of 1 m × 1 m × 1 m, the medium grid had a size of 0.5 m × 0.5 m × 0.5 m, and the fine grid had a size of 0.25 m × 0.25 m × 0.25 m. Figure 7b gives a comparison between the theoretical and medium-grid results. As can be observed, the numerical simulation results are in excellent agreement with the theoretical values. Despite some small differences in the random waves used as background waves, the heights of the freak-wave crest and peak are highly consistent.

Klinting and Sand [50] put forward $\alpha$, ${\beta }_1$, ${\beta }_2$ and $\mu$ to describe the features of freak waves. These characteristic parameters are defined as follows: (1) $\alpha =H_f/H_a$; (2) ${\beta }_1=H_f/H_{f-1}$; (3) ${\beta }_2=H_f/H_{f+1}$; (4) $\mu ={\eta }_f/H_f$. Here, $H_f$ is the freak wave’s height, $H_{f+1}$ and $H_{f-1}$ are the wave heights of adjacent waves before and after the freak wave, $H_a$ is the significant wave height of the wave train, and ${\eta }_f$ is the peak height corresponding to $H_f$. As shown in

Table 2. Comparison of parameters between the theoretical result, numerical result and the measured freak wave in Yura fishery.

Parameters	α	β1	β2	μ
Measured by Mori et al. [48]	2.42	4.03	2.02	0.67
Theoretical result	2.46	3.65	1.96	0.70
Numerical result	2.48	3.70	1.84	0.71

, it can be observed that the characteristic parameters of the generated freak wave largely agree with those of the measured freak wave, indicating the effectiveness of the numerical wave tank and freak-wave model in generating such waves.

4.3. Measurement of the Impact Pressures

In order to obtain the green-water loads induced by the target freak wave, pressure monitoring points were installed on both the deck and breakwater structures of the ship bow. As illustrated in Figure 8 and Figure 9, D100 to D811 were set to obtain the pressures on the deck, while B100 to B805 were set to obtain the pressures on the breakwater. Upon obtaining the pressure time histories during the green water event, the green-water loads could be approximated and loaded to the FEM model of the ship bow structures.

4.4. Ship Motion Responses and Green-Water Loads

4.4.1. Ship Motion Responses Induced by the Freak Wave

To ensure result accuracy, grid independence testing must be conducted beforehand. Given the complexity of wave-ship interactions, it is preferable to select larger computational domains while maintaining a fine simulation of the flow field, thereby effectively reducing the impact of wave reflection caused by the boundary. As such, the present study utilized nested mesh methods, dividing the mesh region into two parts: background and free surface regions, with the free surface regions being denser. Through these methods, computational speed was maximized while maintaining result accuracy. Three types of meshes were tested, including coarse mesh (background region mesh size of 2 m × 2 m × 2 m, free surface region mesh size of 1 m × 1 m × 1 m), medium mesh (background region mesh size of 1 m × 1 m × 1 m, free surface region mesh size of 0.5 m × 0.5 m × 0.5 m), and fine mesh (background region mesh size of 0.5 m × 0.5 m × 0.5 m, free surface region mesh size of 0.25 m × 0.25 m × 0.25 m). Simulation time steps were determined according to the Courant number criterion (Anderson and Wendt [37]), with values of 0.1 s, 0.05 s, and 0.025 s used for the coarse, medium, and fine meshes, respectively. As shown in Figure 10, the results of the grid-independence tests indicate that the medium mesh size provides a balance between simulation accuracy and computational speed.

Since present study focuses on the interaction between freak waves and the ship, the following content examines only the period from when a freak wave appears until the interaction with the ship concludes. Specifically, this refers to the time interval between 530 s and 565 s in Figure 7 (hereinafter, 0 s corresponds to 530 s in Figure 7).

As shown in Figure 11, at 13.4 s, the freak wave comes in touch with the ship bow, causing the wave to roll over due to the large bow flare. Subsequently, a substantial amount of water surges over the deck at an extremely high speed. During this process, the ship bow is lifted by the freak wave. At 15.7 s, the on-deck water slams the breakwater and rapidly moves upwards along its surface. At 16.6 s, the water climbs to the highest point, with some of it overtopping the breakwater. Following this, the water drains off the deck, marking the end of the first green water impact. At this moment, the ship bow is greatly lifted, with increased heave motion reaching the maximum value, and the pitch motion also increasing to the first peak value. Additionally, then, the freak wave moves towards the middle of the ship, and the amplitude of the heave and pitch motions decreases as the ship gradually returns to its normal floating state. As the freak wave approaches the stern of the ship, the stern is gradually lifted, causing the ship bow to sink into the water. At 27.0 s, the bow is submerged below the free surface. At 29.20 s, a small amount of water surges over the deck, leading to the second green water impact but far less than the first. At this moment, the absolute value of the pitch motion reaches its maximum. As the freak wave propagates across the ship, the interaction between the freak wave and the ship ends at about 34.0 s.

4.4.2. Green-Water Loads Induced by the Freak Wave

As the amount of water on the deck during the second green-water impact was significantly less than the first impact, the impact pressure generated from the second impact on the deck was also much smaller than that of the first. Therefore, the following sections only study the impact pressure generated by the first green water impact.

Figure 12 and Figure 13 show the selected time histories of the impact pressure on the deck and breakwater caused by the freak wave. For the deck, the impact pressure starts at 13.8 s and gradually returns to zero after 21.5 s. The maximum impact pressure on the deck is located near the breakwater. The impact pressure at the front of the deck increases slowly, while the area of the deck near the breakwater reaches its peak value at a faster rate. These phenomena may be due to the gradual reduction in water thickness as the water flows over the deck, resulting in a decrease in pressure. However, the velocity of the water flow also increases, and when the water reaches the area of the deck near the breakwater, the flow velocity is the highest. At the same time, a large amount of water accumulates near the breakwater, leading to a large impact pressure in a short period of time. At around 19.6 s, a small impact pressure appears again due to the falling water. Compared with other literature, such as Zhang et al. [23], the impact pressure generated from the falling water was relatively smaller here. This is because other models used in previous studies applied higher superstructures rather than the lower breakwater in the present study. Due to the lower height of the breakwater, most water overtops over the breakwater and only a small portion of the water falls and impacts the deck, resulting in a smaller impact pressure from the falling water.

The breakwater starts to experience impact pressure at 15.5 s and gradually returns to zero at 20.5 s. The maximum impact pressure on the breakwater occurs near the bottom close to the deck. For the breakwater, the impact pressure it experiences is greater than that of the deck, generating a large pressure peak during the impact moment before decreasing at a slightly slower rate. This phenomenon may be due to the fact that when the water flow impacts the surface of the breakwater, the flow velocity reaches its maximum, carrying a huge amount of energy to impact the surface of the breakwater. At this moment, the breakwater mainly experiences dynamic pressure, reaching the peak pressure in an instant. After the peak pressure, two small peaks appear on the surface of the breakwater due to the water overturning under the influence of gravity, resulting in another impact on the breakwater. It is worth noting that the breakwater was impacted by green-water loads before the deck area near it, as has been elaborated by Gomez-Gesteira et al. [51] and Hu et al. [52] in their articles, and will not be reiterated here.

5. Structural Responses of the Deck and Breakwater

5.1. Modelling of the Ship Bow Structures

To investigate the structural response of the deck and breakwater structures under freak waves, it is essential to establish a finite element model of the ship. The dimensions of the container ship were identical to those detailed in Section 4.1. To reduce computational costs, the bow section was selected for the model. The slamming areas subjected to green-water wave loads lack a clear definition among classification societies across various countries. Zhang et al. [23] and Liu et al. [24] conducted research on the slamming load characteristics from the bow to the superstructure area. According to the guidelines of DNV [53], BV [54], GL [55], and KR [56], the impact zone of the bow flare and bottom is situated at a range of 0.1 L from the FP (forward perpendicular) [29,57]. Based on this information, the truncated region extending 0.1 L from the FP is selected as the slamming region to construct an FEM model of the ship’s truncated bow, which is illustrated in Figure 14.

The finite element model employs three types of elements: plate elements are primarily utilized to simulate the shell structures of various components, including the deck, side shells, longitudinal bulkheads, top transverse bulkheads, deck beams, transverse bulkheads, panels, transverse webs, hatch coaming and its panel, and elbow plates. The majority of plate elements are quadrilateral elements, with a minimal number of triangular elements used for component connections and arc transitions. Beam elements are mainly used to simulate large and continuous longitudinal members, such as stringers, stiffeners, and supporting materials, while accounting for beam cross-sections and eccentricities based on actual conditions. Rod elements are mainly used to simulate small-sized reinforcements, such as openings in panels and discontinuous stiffeners. With a total count of 44,476, the deck and breakwater were set to a grid size of 0.7 m, which is fine enough for finite element analysis [29,57,58,59]. In addition, mesh convergence testing were conducted. Mesh sizes ranging from 0.3 m to 1.0 m were chosen, and Figure 15 shows the variation in the peak von Mises stress at Deck D111 and Breakwater B105 (the positions of D111 and B105 are shown in Figure 8 and Figure 9) under different mesh sizes. It can be seen that the results converge when the mesh size is less than 0.7 m, and as the size is further reduced, the peak stress difference is less than 2%. Excessively small mesh sizes would significantly increase computational time, which means that a mesh size of 0.7 m can be considered appropriate for simulating the impact problem in this paper. The time step used for finite element analysis was 0.0045 s. The computer used for the paper is equipped with two AMD EPYC 7402 2.8 GHz CPUs and 128 GB of memory. The time required for mesh independence verification was 624 min, and the average convergence time required per time step was 30 s when performing finite element analysis with a mesh size of 0.7 m.

Reasonable boundary conditions are critical for calculating the dynamic response of structures subjected to slamming loads. Xie et al. [29] used a six-degree-of-freedom constraint to determine the truncation position. Yang et al. [58] indicated that only the aft-side of the model was fixed in six degrees of freedom to ensure that all forces could be balanced with minimal boundary effects on the six brackets. The approach described above was fixed in this article to constrain the ship bow, as illustrated in Figure 16.

Parameters of the finite element model are as follows: elastic modulus $E$ of 2.06 × 10⁵ MPa, Poisson’s ratio $v$ of 0.3, and material density $\rho$ of 7.8 × 10³ kg/m³. When subjected to slamming pressure, the constitutive behavior of the material differs from that observed under static loading conditions. To account for this, the dynamic constitutive model Cowper–Symonds model [60] commonly used in structural impact analysis was introduced.

$${\sigma }_{dy}={\sigma }_y\left[1+{\left(\frac{\dot{\epsilon }}{D}\right)}^{\frac{1}{e}}\right]$$

(16)

where $\dot{\epsilon }$ represents the plastic strain rate, while $D$ and $e$ are coefficients that must be determined based on test data for various materials. Paik et al. [61] determined a set of Cowper–Symonds coefficients for high-tensile steel, determining $D$ = 3200 s⁻¹ and $e$ = 5.

5.2. Loading on the Ship Bow Structures

In this study, the slamming load was first obtained in Section 4.4.2. Thereafter, the impact loads were applied to the FEM model and the resulting structural response was calculated. In problems where structural deformations are significant and impact the fluid field, a coupled fluid–structure interaction method is necessary. However, given that the deformations of the structure can be considered to have a negligible effect on the fluid field, an uncoupled CFD-FEA method was employed in the present study. This approach involved independent utilization of the wave loads simulation and the finite element analysis. For the sake of result accuracy, impact pressure values obtained from the wave load simulation were collected over the hull surface grids every 0.05 s. Subsequently, the time series data were directly applied to the finite element nodes in the same location. In the following sections, the impact pressure obtained from 13.75 to 18.25 s in Section 4.4.2 were applied to the finite element model for structural response calculation. The loading method at different times is shown in Figure 17. The structural solver Nastran [39] was applied for the structural response calculation.

5.3. Structural Response Analysis

5.3.1. Structure Responses of the Deck and Breakwater

This paper focuses on the structural strength assessment under slamming loads, and the evaluation process mainly involves verifying the yield strength of the structure. For plate elements, the Mises stress ${\sigma }_{vm}$ is:

$${\sigma }_{vm}=\sqrt{{\sigma }_x^2-{\sigma }_x{\sigma }_y+{\sigma }_y^2+3{\tau }_{xy}^2}$$

(17)

where ${\sigma }_x$ and ${\sigma }_y$ are the normal stresses in the x and y directions of the element, and ${\tau }_{xy}$ is the shear stress of the element.

The yield factor ${\lambda }_y$ is used to determine whether there is a yielding phenomenon in the ship’s structure, and the specific judgment formula is as follows [29]:

$$\mathrm{Deck}: {\lambda }_y=\frac{k{\sigma }_{vm}}{225}$$

(18)

$$\mathrm{Breakwater}: {\lambda }_y=\frac{k{\sigma }_{axial}}{235}$$

(19)

where $k$ is the material coefficient, which for the high tensile steel used in this model is 0.92 [62].

Using the dynamic analysis method, structural responses of the deck and breakwater structures were obtained. Figure 18 shows the contour plots of total displacement at different moments of the ship bow. During 14.0 to 15.5 s, the green water gradually flows from the deck front to the breakwater, causing the increase of the displacement. The area of the deck affected by green-water loads becomes larger over time. At 15.5 s, most of the deck is subjected to the impact of the green-water loads, but the water has not yet struck the breakwater. At 15.7 s, the impact load is applied to the breakwater, which results in a significant displacement of the breakwater. The maximum displacement occurs above the centerline near B105, and the value of the maximum displacement reaches as high as 1 m. This indicates that the area is likely to have been damaged or even fractured due to the impact. Over time, the green water gradually flows out of the bow, and the deck displacement caused by the impact of green-water loads gradually decreases to zero. It can be preliminarily inferred that the deck remains undamaged during the green water event. However, residual displacement of different degrees still exists in the breakwater area, indicating that the breakwater is likely to be damaged under the green-water loads.

Figure 19 shows the contour plots of the Von Mises stress at different moments of the ship bow. To facilitate the visualization of the stress distribution beneath the deck, half of the deck plate was hidden in the Mises stress contour plot. Similar to the displacement contour plots, during 14.0 to 15.5 s, the green water event happens as time increases, resulting in larger impact forces and affected deck areas over time. At 15.7 s, the green water impacts the breakwater with a sudden and enormous force. Further investigation reveals that the point of maximum stress response occurs at the stiffeners of the breakwater B101. The stress response value exceeds the material’s yield limit, suggesting that this area may have been damaged. Subsequently, the impact force gradually decreases, but residual stresses still exist in the breakwater area. It is noteworthy that the stress response of the deck near the breakwater is greater than that of the forward deck, and the stress response near the longitudinal centerline of the deck and breakwater is greater than that on both sides. Therefore, particular attention should be paid to the above regions when evaluating the safety of the ship’s bow.

Based on the previous description, the stress response increases as the deck approaches the position of the breakwater. Figure 20 presents the time histories of the stress response and total displacement at locations D108, D110, and D111 on the deck (the positions of D108, D110, and D111 are shown in Figure 8). These three locations experience a greater impact load compared to other areas of the deck, resulting in a correspondingly greater stress response and displacement. In terms of numerical patterns, compared with the results obtained in Section 4.4.2 for the historical impact load, it can be seen that the stress response results of the deck exhibit synchronous variations with the impact load applied to the deck, and the shape of the curves and the time of occurrence of the maximum values are basically the same, with a relatively short duration of the maximum values. However, after the maximum peak disappears, there is still a slight fluctuation in stress, which may be due to the influence of other positions after the first impact ends, and a slight fluctuation in stress response occurs subsequently. The total displacement at the deck also exhibits synchronous changes with the impact load applied to the deck.

The stress response in the region near the longitudinal centerline of the breakwater was also significant and requires particular attention. The contour plots in Figure 18 and Figure 19 show that the breakwater suffered severe damage due to the impact of the green water. The stress response peak at the breakwater near the longitudinal centerline exceeds its yield limit and may have undergone fracture deformation. Other areas of the breakwater also experience varying degrees of damage, with some residual stress and displacement remaining even after the reduction in the impact loads. When further analyzing the failure mechanism of the breakwater, it was found that the maximum stress response occurs at the bottom of the stiffener behind the breakwater near the longitudinal centerline (i.e., the stiffener located behind breakwater B101), reaching a peak of 315.5 MPa. The areas with the largest displacement were located in the upper part of the breakwater, with the maximum displacement occurring near the longitudinal centerline of the upper part of the breakwater (at breakwater B105). Figure 21 shows the stress response and displacement time histories of the areas with the maximum stress response and displacement. The stress response of the stiffener behind the breakwater was greater than that of the breakwater itself, indicating that the stiffener helped to absorb a considerable amount of the impact force, improving the overall strength and bearing capacity of the structure. However, due to the excessive impact load, the stress response of the stiffener behind breakwater B101 exceeds its yield limit, and the stiffener may fracture under the impact of the green water load. The breakwater near the longitudinal centerline experiences significant displacement as it may have lost support. The upper part of the breakwater is less fixed compared to the lower part, resulting in the largest displacement at that location. The contour plot results at 18.25 s show that in addition to the area near the centerline of the breakwater, where residual stress and displacement exceeds the limit, other parts of the breakwater also exhibit residual stress and displacement. This indicates that other areas of the breakwater also suffered damage under the impact of green water.

Through comparing the stress response with the yield limit, it was found that the stiffeners near the longitudinal centerline of the breakwater may be damaged during green-water events. However, residual stresses and displacements are still present in other areas of the breakwater after the water drains off the deck, indicating potential damage to these areas as well. Therefore, the safety evaluation of the ship’s bow structure cannot solely rely on comparing peak stress responses with the yield limit. In

Table 3. Safety evaluation results of the ship’s bow.

Part	Von Mises Stress (MPa)	Allowable Stress (MPa)	λy	Safety Assessment Results
Deck	95.7	244.5	39.1%	qualified
Deck stiffener	155.4	244.5	63.6%	qualified
Breakwater	264.3	255.4	103.5%	unqualified
Breakwater stiffener	293.0	255.4	114.7%	unqualified

, the peak stress responses of four key areas on the ship’s bow, namely the deck, deck stiffener, breakwater, and breakwater stiffener (excluding the areas that may be damaged due to exceeding the yield limit), are compared against allowable stress. The standard for allowable stress is taken from the Rules for Classification of Sea-going Steel Ships [62] published by the China Classification Society (CCS). The maximum stress responses of the deck and deck stiffener are within the yield limit, indicating that the dynamic response evaluation of the deck and deck stiffener is qualified. However, the dynamic response evaluation of the breakwater and breakwater stiffener is unqualified. For these areas, increasing the thickness of the breakwater and the number of breakwater stiffeners to increase their strength and avoid structural damage is recommended.

5.3.2. Comparison of Dynamic and Static Analysis Methods

Based on the previous evaluation, the green-water loads caused by the freak wave are likely to cause damage to the ship bow structures. However, the results calculated in Section 5.3.1 are based on the dynamic analysis. Generally, equivalent static-analysis methods are computationally efficient and widely applicable in engineering, but their results are generally conservative, which can easily cause material waste (Yang Bin et al. [57]). The concept of equivalent static analysis is to consolidate the peak values of each component at different times into a single moment, which results in more conservative outcomes, and makes it suitable for the preliminary analysis of ship structures. For ships with high structural requirements, dynamic response analysis is needed (Xie et al. [29]). In the present study, structural responses using dynamic and static analyses at deck location D111 and breakwater stiffener location B101 are compared in

Table 4. Comparison of structural response between equivalent static method and dynamic analysis method.

Part	Von Mises Stress (MPa)
	Equivalent Static Analysis Method	Dynamic Analysis Method
D111 of the deck	232.9	95.7
B101 of the stiffener	735.0	315.5

.

Both for the undamaged deck and the breakwater stiffeners that may be damaged, the results obtained from the equivalent static analysis method are significantly higher than those from the dynamic analysis method. As analyzed in Section 5.3.1, the impact load varies in space and time, and the peak time at different locations is also different. The equivalent static analysis method concentrates the peak loads of all components into one moment, resulting in a conservative estimate. The dynamic analysis method can reflect the temporal and spatial differences of the impact load on the ship, and can consider the structural dynamic response and nonlinear behavior at different times, thus providing more realistic and accurate results.

5.3.3. Comparison of Linear and Nonlinear Analysis Methods

The finite element analysis method used in this paper was a nonlinear dynamic analysis; however, there is another commonly used finite element analysis method, namely, linear dynamic analysis (Wang et al. [63]). Linear analysis has a broader range of applications and can solve various dynamic-response problems. Additionally, it has a simpler calculation process and faster computational speed. However, linear analysis involves a series of assumptions, such as assuming that materials remain in an elastic state, stress–strain relationships are always proportional, loads are proportional to structural responses, and structures undergo only small displacements and deformations. Nonlinear dynamic analysis, on the other hand, can consider factors such as large deformations, large displacements, and material nonlinearity under load, resulting in the capability to simulate more complex structural responses, more accurately reflect real-world conditions, and provide more reliable design results. To compare the differences between linear and nonlinear dynamic analysis, this paper compared the structural response results obtained by the two methods. Figure 22 presents a comparison of the stress response and total displacement results obtained at deck location D111 and breakwater stiffener location B101 using nonlinear dynamic response analysis and linear dynamic analysis methods.

Under the same mesh and time step, for places where the stress response does not exceed the material yield limit, such as the deck, the results of nonlinear dynamic analysis method and linear dynamic analysis method were basically consistent. Figure 22a,b show the analytical results at deck D111. However, for places where the stress response exceeds the material yield limit, such as the breakwater stiffener, there was a significant difference between the results of the nonlinear dynamic analysis method and the linear dynamic analysis method. As shown in Figure 22c,d, in the nonlinear dynamic analysis method, when the impact load causes the failure of the breakwater stiffener and enters the plastic stage, the structure undergoes significant displacement, and the rate of increase in the stress response of the breakwater stiffener begins to slow down. After the impact load decreases, the stress response of the breakwater stiffener does not immediately decrease but decreases slowly. However, the results of the linear dynamic analysis method show a significant stress response and small displacement. This indicates that the nonlinear dynamic analysis method can accurately simulate plastic behavior or instantaneous damage that cannot be described by the linear dynamic analysis method.

Therefore, in the assessment of the slamming strength of the ship hull, if the slamming load is small and does not cause the structure of the ship to exceed its yield limit, the linear dynamic analysis method can be used to simplify the calculation process. However, if the stress response obtained from the linear analysis method is significantly higher than the material’s yield limit, then the nonlinear dynamic analysis method should be used to obtain more accurate structural stress-response and displacement results. The nonlinear dynamics analysis method is more complex compared to the linear dynamics analysis method, but the results are generally more accurate. It is more suitable for analyzing ships with higher structural requirements.

6. Conclusions

The purpose of this paper was to study the green-water loads and structural responses of ship bow structures induced by freak waves. To achieve this, a three-dimensional numerical wave tank was established to generate freak waves and to obtain the wave-ship interaction and green-water loads. A finite element model of the ship bow was built, on which the green-water loads were applied to obtain the displacement and stress of the deck and breakwater structures. The main conclusions are as follows:

(1)

It can be concluded from numerical simulation results that freak waves have a crucial impact on ship motion response. It should be noted that wave crests and troughs can cause significant heave and pitch motion. Therefore, it is necessary to consider the influence of extreme sea conditions on ship stability in the design and verification of ships.

(2)

Freak waves have significant impacts on green-water loads. The green-water event caused by freak waves generates a huge impact force on the deck and breakwater of the ship bow, with a greater slamming load on the breakwater than on the deck. The tremendous impact loads may pose a threat to the safety of personnel on the deck and cause damage to the equipment on deck. Therefore, in engineering applications, it is necessary to pay attention to checking the local strength of the bow deck and the breakwater to ensure the safety of the ship and personnel.

(3)

Freak waves are able to cause severe structural responses in the deck and breakwater structures of the ship bow. Under the impact of green water caused by the freak wave, the breakwater and breakwater stiffener suffered relatively severe damage. Therefore, in the design of ships, it is necessary to strengthen the structure of the breakwater to prevent damage to the breakwater caused by green-water events resulting from freak waves. Increasing the thickness of the breakwater and the number of breakwater stiffeners is recommended to guarantee local strength and avoid structural damage.

(4)

Equivalent static analysis methods are computationally efficient and widely applicable in engineering, but tend to produce conservative results, which may result in material waste. As such, they are best suited for the preliminary design stage of ship hulls. However, this study focuses on the high-speed water impact caused by freak waves, which is a transient dynamic problem. The equivalent static analysis methods cannot reflect the temporal and spatial differences in the green water impact phenomenon. Therefore, the dynamic analysis methods are more suitable for the analysis of such problems.

(5)

The linear dynamic analysis method can capture the spatial and temporal characteristics of slamming loads, but fails to account for geometric and material nonlinearity. It is best suited for scenarios where the slamming load is small and does not result in material failure. However, these methods cannot accurately reflect the state of structural damage or accurately describe the areas of the structure that have been damaged. The damaged areas determined by these methods are often overestimated. Therefore, if the stress response of the linear analysis methods exceeds the yield limit of the material, nonlinear dynamic analysis methods should be used. Nonlinear dynamic analysis methods can not only more accurately display the damaged areas and states of the structure, but also obtain more precise stress responses and displacement results of the structure.

This study may contribute to the better understanding of ship safety design under freak-wave conditions. The green-water loads and corresponding structural analysis may be of reference value from an engineering point of view. Limitations of this study include the lack of physical experiment verification and the assumption of the uncoupled fluid–structure interaction calculation, which could be further investigated.