1   Introduction

The importance of bow wave breaking to the hydrodynamic performance of high-speed ships has long been recognized (Wang, et al., 2020a, 2023). The main research methods of wave breaking at the bow are mainly divided into experiments and numerical simulations. Model experiments can give intuitive wave-breaking phenomena and ship resistance data. Baba (1969) confirmed that bow wave breaking has a large contribution to viscous resistance through theoretical analysis and experimental verification. Kayo and Takekuma (1981) conducted velocity field measurements and flow visualization technique on bow wave breaking and the shear flow was observed on the free surface with the use of dye particles. Dong et al. (1997) conducted experimental research on the ship model at two different speeds, with the help of particle-image-velocity measurements and free surface visualization techniques, and the formation of the bow wave sheet was intensively studied. Olivieri et al. (2007) conducted an experimental study for the scars and vortices induced by the wave breaking of model DTMB 5415, and Fr = 0.35 was intensively studied because of the large extent of plunging bow and spilling shoulder wave breaking. Wang et al. (2020b) conducted experimental investigations on KCS wave breaking at 5 different trim angles and 3 sinkages with the Froude numbers ranging from 0.26 to 0.4 in the towing tank of China Ship Scientific Research Center (CSSRC), and the influence effects of each parameter were discussed. Further, a more violent bow wave breaking was found than in the previous studies, which provided physical understanding and additional validation data like photo observation results, and the test condition selected in the present paper was based on the results.

Although experimental study on bow wave breaking has been quite in-depth, there are still shortcomings in the detailed physical description of scars and plunging breakers. With the improvement of computing power, numerical simulations become an important tool for predicting fluid patterns in ocean engineering. The application of the computational fluid dynamics (CFD) method in ship hydrodynamics is more and more extensive, and many scholars have carried out numerical studies on wave breaking of high-speed ships. The precise capture of the free surface is a difficult problem in numerical simulation due to the large deformation characteristics of the free surface caused by wave breaking at the bow of the ship. The dominant twophase flow simulation methods in grid-based CFD methods are Volume of Fluid (VOF) (Hirt and Shannon, 1968), levelset (Osher and Sethian, 1988) and CIP-based method (Hu and Kashiwagi, 2004), respectively.

Based on these free surface-capturing methods, scholars have carried out relevant numerical studies on bow wave breaking. Wilson et al. (2007) used the Reynolds averaged model (RANS) to solve the fluid governing equations, and used the level-set method to capture the free surface changes, which realized the numerical simulation of the bow rising wave of the DTMB5415 model. Based on the ship hydrodynamic solver naoe-FOAM-SJTU, Wan (2017) conducted numerical simulation with the RANS model solving the fluid equations and VOF method capturing the free surface, and the details of bow wave breaking were obtained. Yu (2020) numerically simulated the bow wave breaking of a KCS model with fixed attitude at different speeds based on the naoe-FOAM-SJTU and analyzed the effect of speed on the bow wave structure and flow field.

The phenomena of air entrainment occur during the wave plunging and breaking, and the entrained air exists in the form of bubbles, which may affect the wave-breaking patterns. Moraga et al. (2008) took the air entrainment into consideration and proposed a sub-grid model for breaking bow waves, which presented the right trend of wave breaking of the full-scale DTMB 5415. Recently, Zhang et al. (2023) proposed the bubble dynamics equation, which unifies the classical bubble equation and provides a new physical guidance for the air entrainment behavior in breaking bow waves. Despite the great success, in some aspects, there are inherent difficulties hard to deal with. Due to the characteristics of the Eulerian mesh and the limitations of the mesh size, the simulation of the wave overturning and breaking problem still faces the difficulty of the low accuracy of the free surface capture to a certain extent.

In recent years, with the continuous progress and development of GPU (Graphics Processing Unit) technology, meshless methods have been more widely used and promoted, the most representative of which is smoothed particle hydrodynamics (SPH) (Gingold and Monaghan, 1977; Lucy, 1977; Liu and Zhang, 2019; Sun et al., 2021; Lyu et al., 2023). By discretizing a continuous fluid into a set of particles and simulating the fluid motion through inter-particle interactions, SPH has a unique advantage in dealing with complex flow fields and including physical effects, such as free-surface flow (Liu and Zhang, 2019; Shao et al., 2012; Sun et al., 2021), multiphase flow (Sun et al., 2019a and b), etc., and can simulate the details and motions of the fluid in a more fine-grained way (Liu and Zhang, 2019; Xu et al., 2023). Landrini et al. (2012) applied a hybrid BEM-SPH method to study the bow wave splashes. Colagrossi (2001) proposed a 2D+t SPH model that gave a good qualitative description of the breaking phenomenon. Marrone et al. (2011, 2012) proposed two different models, 2D+t and 3D models, both based on SPH, to resolve bow wave breaking, and robustness and accuracy were validated by comparing with the experimental data.

When the bow wave breaking phenomenon occurs, the air is entrapped and mixed into the water and the mixture flow goes down to the ship wake flow. However, the size of the air entrainment is very small. Taking the air phase into simulation will lead to a significant increase in computing costs. Therefore, the effects of air are not considered in the present study, and air entrainment will be considered in subsequent research.

In the present work, the utilized SPH method is described and deals with the prediction and analysis of bow wave breaking induced by benchmark ship model KCS. Subsequently, the results of SPH are compared with those of the grid-based commercial software STAR-CCM+, which is based on a RANS solver coupled with a VOF method. The main framework of this paper goes as follows. The first part is the numerical methods, where SPH schemes and numerical techniques are presented. The second part is the geometry model and test conditions. Then comes the simulation part, where bow wave breaking simulations are present at different Froude numbers. In this part, extensive comparisons are performed against the simulations obtained by STAR-CCM+ and experiments. Finally, a conclusion of this paper is drawn.

2   Numerical method

2.1   Governing equations and their discrete forms

Based on the weak compressibility assumption, the fluid governing equations (Liu and Liu, 2003; Sun, 2018) are as follows:

where the first line is the mass conservation equation and the second line is the momentum conservation equation. D/Dt is the material derivative concerning time. ρ denotes the density of the fluid, u represents the velocity vector of the fluid, p signifies the isotropic pressure of each fluid microcluster, μ is the coefficient of dynamic viscosity of the fluid, and g denotes the vector of gravitational acceleration. c₀ denotes the artificial speed of sound and ρ₀ denotes the reference density of the fluid. In order to satisfy the hypothesis of weak compressibility, the c₀ is usually obtained through (Sun, 2018)

where U_max and P_max, respectively, represent the maximum expected velocity and fluid pressure during the simulation.

The governing equations are treated in the SPH method employing the kernel approximation and the particle approximation (Liu and Liu, 2003). In the kernel approximation, the continuity function f (r) and its derivatives ▽f (r) can be obtained by interpolating the kernel approximation of the function in the tightly branched domain around it. The selection of kernel function should be very careful, and the conditions, such as Dirac function properties, symmetry conditions, normalization conditions, and tight branching conditions (Liu and Liu, 2003; Liu and Liu, 2010; Sun et al., 2019a), should be satisfied. In the present SPH framework, the Wendland C2 kernel function is applied. In the particle approximation, the flow field is discretized into a series of Lagrangian particles carrying information about the position r, velocity u, pressure p, and density ρ, etc. The solution of the governing equations is achieved through the interaction between the Lagrangian particles. In summary, the discrete form of the SPH fluid governing equations under the weak compressibility assumption can be obtained as follows (Liu and Liu, 2003; Sun et al., 2018):

where subscript i denotes the variable associated with the i-th particle while j represents the neighboring particles located within the kernel support of i. The W_ij denotes the kernel function about the relative position ‖r_i − r_j‖ and smoothing length h. ▽i represents the gradient taking at the position vector r_i. The V_j denotes the volume of the particle j. The $\mathcal{T}_i^{(D)} $ is known as the density diffusion term to avoid the high-frequency noise of pressure, which is the most critical superiority of the δ-SPH model. The $\mathcal{T}_i^{(D)} $ is given as (Sun et al., 2019a)

where the coefficient term δ is usually adopted as 0.1. $ \langle\nabla\rangle_i^L$ is the renormalized spatial gradient (Sun et al., 2019a).

In particular, the artificial viscosity term $\boldsymbol{\mathcal{T}}_i^{(V)} $ is applied instead of the physical viscosity term for numerical stability, which has very little effect on the results. The artificial viscosity term is given as (Sun et al., 2019a)

where α is set to 0.02.

2.2   Solid boundary conditions

The breaking of bow waves is related to the separation point between the fluid and the hull surface, which is directly related to the implementation of the boundary condition.

Due to the fact that the fixed ghost particle technique is more suitable for complex boundaries (Adami et al., 2012), as it does not require the evaluation of the normal direction on the solid surface, in the present work, the fixed ghost particle technique is implemented. The pressure of the ghost particle is interpolated from the surrounding fluid particles using the Shepard formula (Shepard, 1968). In the particle approximation, the pressure gradient term is described in the form of $ \frac{1}{\rho_i} \sum\limits_j\left(p_i+p_j\right) \nabla_i W_{i j} V_j$. In the numerical simulation of fluid-solid interaction using the SPH method, the pressure summation of a fluid particle and a solid particle may become below zero, which manifests that the structure has an attractive effect on fluids, resulting in a disturbing separation process of fluids and solids. A numerical treatment is implemented when the i and j do not belong to the same particle set (fluid particles or solid particles) and the summation of the pressures of the i-th and j-th particles is negative (i.e., p_i + p_j < 0). The formula is given as follows:

where, Ω_Fs and Ω_S denote the particles belonging to the free surface and the ghost particles, respectively. A numerical study demonstrating the effectiveness of the proposed numerical treatment will be presented in Section 4.

2.3   Inflow and outflow boundaries

To reduce computational complexity, the concept of relative flow is adopted to simulate the motion of ships in still water. By introducing inflow and outflow boundaries (Sun et al., 2019a) in the SPH algorithm, a uniform flow numerical tank is established. As shown in Figure 2, the simulation domain is divided into three parts, which are the inflow region, the SPH computational region and the outflow region. The inflow region is located at the left end of the numerical tank, in which the velocity of the particles is assigned at the given speed. The outflow region is located at the right end of the numerical tank, maintaining computational stability. Special operations on some physical variables are required in the inflow and outflow regions, which are different from the SPH computational region. In the inflow region, the particle pressure is calculated using the hydrostatic pressure formula. Meanwhile, when particles are in the inflow or outflow region, the density derivative and acceleration are frozen, i.e., Dρ/Dt = 0 and Du/Dt = 0.

Figure  1  The phenomenon of ship bow wave breaking

Download: Full-Size Img

Figure  2  Schematic diagram of inflow and outflow boundary conditions

Download: Full-Size Img

During the simulation, the positions of these particles, in the inflow or outflow regions, are updated based on the given velocity without being updated by the SPH calculation. When a particle moves from the inflow region to the SPH computational region, it is immediately converted into an SPH particle whose physical quantities start to be updated by the SPH calculation. At the same time, when the leftmost particle at the inflow region (see Figure 2) has a movement of more than a particle spacing after it is generated, a new particle from the candidate set will be added to the inflow, being placed at a particle spacing to the left of the leftmost particle, which is used to allow for a constant stream of particles from the inflow region into the SPH calculation region and ensure the integrity of the support domain of the fluid particles. Conversely, when a particle moves out of the outflow region, it will be deleted from the simulation and supplied into the candidate set.

3   Geometry and test conditions

3.1   Geometry model

The geometry of interest is the benchmark ship model KCS, the standard model for the 2023 International Hydrodynamics Conference and the 8th National Conference of Marine Computational Fluid Dynamics, shown in Figure 3. The scale ratio of the model is 1∶37.89, the length between the vertical lines of the ship model is 6.0702 m, and the width of the ship model is 0.8498 m. The principle parameters of the ship model KCS are shown in Table 1.

Figure  3  Three-view drawings of KCS hull model

Download: Full-Size Img

3.2   Test conditions

According to existing literature (Wang et al., 2020b), the Froude number has a significant impact on the wavebreaking phenomenon at the bow of the KCS benchmark model ship. When the Froude number is greater than 0.35, significant wave-breaking phenomena were observed in the experiment. Therefore, a series of approaching speeds corresponding to Froude numbers of 0.35, 0.425, 0.5 and 0.6 are taken into account to investigate the wave-breaking phenomenon, which is much higher than the commonly used Fr = 0.26, corresponding to the hypothetical fullscale ship design speed, and not possible in a real KCS voyage. In addition to the Froude numbers, the trim of the KCS model ship also has a strong effect on the bow wave breaking, especially when the trim occurs in the bow, which will cause violent wave breaking. Based on the above research (Wang et al., 2020b), the KCS ship model is fixed with a bow inclination of 1 degree during the simulation, which means a greater draft at the bow.

In SPH simulation, to reduce computational costs, mirror processing is performed along the longitudinal profile of the KCS ship model, and a half ship is selected for numerical calculation. In the calculation, only the influence of water on wave breaking is considered, without considering the effect of air.

A uniform flow numerical tank is established by setting inflow and outflow boundaries. The KCS ship model is located on one side of the numerical tank, with a bow distance of 0.3 L_pp from the inflow boundary, a static water surface distance of 4B from the bottom of the tank (noted in literature (Marrone et al., 2011), and a numerical tank width of 1.0 L_pp.

According to the results in reference (Yu and Wan, 2020), it can be seen that the amplitude variation of KCS hull surface waves at Fr = 0.35 ranges from −0.06 m to 0.12 m. To capture this flipping process, the spacing between SPH particles is set to 0.01 m, and the total number of particles is about 13 million. In the commercial computational software STAR-CCM+, the same calculation example is set as a comparison, and the grid size around the calculation area and the ship body is consistent with the particle spacing in the SPH method. Since the air area is also meshed, the final calculation requires a grid number of around 20 million. In the experiment conducted by Wang et al. (2020b), the violent wave-breaking area is recognized, which is located at the bow and shoulder of the ship. Therefore, 15 typical positions in the area are selected to set up wave gauge measurement points to detect wave height changes in the process of bow wave breaking, as shown in Figure 4. The position parameters of the wave gauges are shown in Table 2.

Figure  4  Schematic diagram of the measurement points of the wave height

Download: Full-Size Img

4   Simulation results and analysis

4.1   Validation of the proposed surface separation treatment

In this section, the effectiveness of the numerical treatments mentioned in Section 2.2 is verified. Fr = 0.425 is selected for the comparative study due to the large extent of bow wave breaking. The numerical simulation is conducted within the SPH scheme described in Section 2, and the difference between the two cases is whether the proposed surface separation treatment is applied or not. In the former case, no treatment is done with the pressure gradient term, while in the latter case, the pressure summation of a fluid particle and a solid particle at the free surface is forced to zero when it becomes negative. Figure 5 shows the results of bow wave overturning at different time instants at x = 0.5 m obtained from two cases. The negative pressure has a significant impact on the bow wave breaking phenomena. The ghost particles within the kernel radius of any free surface particle have an attractive effect on the free surface particles, which will lead to the accumulation of fluid particles on the surface of the ship.

Figure  5  Comparison of results of solid boundary negative pressure treatment. The left plot shows the SPH result allowing negative pressure on the ship surface close to the free surface; the Right plot shows the result using the negative pressure prevention on ghost particles adjacent to the free surface

Download: Full-Size Img

With the treatment, the attractive effect of the solid particles, coming from the negative pressure term in momentum equations, is suppressed, and the phenomena of wave overturning at the bow is smoother, which is consistent with the conclusions in the references (Wang et al., 2020b; Sun et al., 2018). Therefore, in the following, the treatment is implemented for all the simulations.

4.2   Validation of the wave breaking profiles

The wave profile on the KCS hull surface at Fr = 0.35 is depicted in Figure 6 and compared against the results obtained by STAR-CCM+. It can be seen from the figure that near the bow of the ship, the wave generation shows a peak valley trend, and this trend evolves towards the stern. The first peak appears around x/L_pp = 0.2, and the first valley appears around x/L_pp = 0.3. From the trend of the curve and the position of peak and valley values, a fairly good agreement is obtained with STAR-CCM+.

Figure  6  Comparison of wave profile on KCS hull surface at Fr = 0.35 obtained by SPH and STAR-CCM+

Download: Full-Size Img

Snapshots from the experiment (Wang and Wan, 2019) and SPH simulation are compared in Figure 8 for different Froude numbers 0.35 and 0.425, respectively, and similar wave overturning phenomena can be observed. It should be noted that, since the experimental data under the same test conditions used in the paper have not yet been published, in the selected experiment, the trim angle is zero, which means the breaking phenomena will be less violent than that in the numerical simulation. Therefore, the experimental photographs are only used for the description of wave-breaking phenomena. It can be seen that the wave overturning phenomenon is observed in both the experiment and SPH numerical simulation at Fr = 0.35 and 0.425. Comparing the results of Fr = 0.35 and Fr = 0.425, the phenomenon of wave breaking at the bow of the ship becomes increasingly severe as the Froude number increases.

At Fr = 0.425, the scars induced by the bow wave breaking are observed and are consistent with the experimental phenomenon, which indicates that the utilized SPH method can be robust and reliable in predicting the bow wave breaking of high-speed ships.

The results of the SPH method are analyzed qualitatively in the above. The wave-breaking phenomenon simulated by the SPH method obtains consistent results with the experiment. To quantitatively analyze the phenomenon of wave breaking at the bow of the ship, wave height obtained from 3 wave gauges near the bow position is presented and compared with the results of STAR-CCM+.

Before the comparison, the wave height detection methods implemented in SPH and software STAR-CCM+ are introduced, shown in Figure 7. In the STAR-CCM+, the wave height is obtained by the top outline of the free surface at measurement points. While, in the SPH framework, wave height detection requires 2 steps. Firstly, a series of wave gauge particles are generated according to the position of the measuring point, and the spacing between these particles is smaller. Then, the wave height is obtained by the interpolation of the surrounding SPH fluid particles to the wave gauge particles with the formula, written as

Figure  7  Schematic diagram of the difference between wave height detection in SPH and STAR-CCM+

Download: Full-Size Img

Figure  8  Comparison of free surface in bow region obtained by SPH simulation and experiment (Wang and Wan, 2019) at Fr = 0.35 and 0.425

Download: Full-Size Img

where the Ω_F denotes the particles belonging to the fluid. The difference between the two detection methods can explain the discontinuity in the curves obtained by STARCCM+ and the discrepancy in the curve of the wave heights noted below.

The time history results of the wave height at typical positions obtained by the SPH method and STAR-CCM+ are compared in Figures 9 and 10 for Fr = 0.35 and Fr = 0.425, respectively.

Figure  9  Comparison of the wave height time history at typical positions at Fr = 0.35 obtained by the SPH method and STAR-CCM+

Download: Full-Size Img

Figure  10  Comparison of the wave height time history at typical positions at Fr = 0.425 obtained by the SPH method and STAR-CCM+

Download: Full-Size Img

The vertical axis in the figure represents the wave height nondimensionalized by L_pp. The evolutions of the time history show a fair agreement between the SPH method and STAR-CCM+. Although there may be a slight discrepancy in the wave height on some time instants, this can be explained. The wave height error calculated by the SPH method and STAR-CCM+ software in the figures is within 0.01 m, which is an order of magnitude with particle spacing in SPH simulation. The slight discrepancy in the wave height can be reduced by setting a smaller particle spacing, indicating that this error is acceptable. As for the discontinuity in the wave height curve obtained by STAR-CCM+, it is due to the wave height detection method, which is applied by capturing the top outline of the free surface at measurement points.

From time history results of the wave height of 3 wave height gauges at Fr = 0.35, it can be observed that there is a stable uplift of the wave surface at the bow breaking point (see Figure 9(a)). Due to the attend of the KCS bow, the wave has a propagation perpendicular to the direction of flow, with peaks and valleys propagating outward (see Figure 9(b)). After leaving the bow area, the height of the free surface decreases (see Figure 9(c)). As the Froude number rises from 0.35 to 0.425, the wave height is much larger which indicates that the bow wave breaking phenomena is more intense at Fr = 0.425. Figure 11 shows the comparison of the mean wave height obtained by SPH and the experimental (Wang et al., 2020b; Liu et al., 2022) at Fr = 0.35. The upper part of the figure is the mean wave height contour obtained by SPH simulation, and the lower part is the results obtained by the experiment. As can be seen from Figure 11, the results of SPH are consistent with the experiment in terms of the distribution of the mean wave height, which demonstrates that the SPH method can predict the change in wave surface height when the bow wave breaking occurs.

Figure  11  Comparison of the mean wave height obtained by SPH and the experiment (Wang et al., 2020b; Liu et al., 2022)

Download: Full-Size Img

4.3   Analysis of velocity field and wave elevation

To further investigate the wave characteristics of the bow wave breaking at the bow of the KCS ship, velocity fields in 3 directions at Fr = 0.35 and Fr = 0.425 are depicted in Figure 12. It can be obviously seen that due to the attend of the bow of the ship, the x-direction velocity of the fluid is significantly weakened, and the y-direction velocity of the fluid reaches its maximum at the highest point of the wave crest, which is consistent with the characteristics of wave overturning. Once the wave overturning is formed, it develops in the direction of flow. Comparing the wavebreaking results at the bow at Fr = 0.35 and Fr = 0.425, it is not difficult to see that the wave breaking at Fr = 0.425 is more severe, with higher wave height and larger wavebreaking region, and there are multiple obvious wave scars induced by wave overturning at bow region.

Figure  12  Comparison of velocity contour at Fr = 0.35 and Fr = 0.425

Download: Full-Size Img

Figure 13 shows free surface evolutions at Fr = 0.425 and detailed information on the process of the bow wave breaking at high-speed KCS ship is provided. The wave overturning at the KCS bow and the formation of the scars are observed. In the initial stage of wave breaking, the fluid around the bow of the ship climbs along the surface of the ship. After reaching the highest point, the climbing water begins to move downward due to gravity (see Figure 13(a)). The wave overturning develops constantly until the wave crest collapses and hits the water surface like the jet flow (see Figure 13(b)). At the same time, since the fluid has a propagation perpendicular to the direction of flow, a scar is formed at the wave contact point with evidence of the wave jet (see Figure 13(c)). Subsequently, the area of the bow wave breaking extends to a larger range in the vertical direction and more scars and plunging waves are generated with the continuous evolution of the bow wave breaking (see Figure 13(d)).

Figure  13  The process of bow wave breaking on KCS model ship at Fr = 0.425 s

Download: Full-Size Img

Further, to show the effect of the Froude number, the wave breaking results obtained by SPH at four different Froude numbers, namely 0.35, 0.425, 0.5, and 0.6 are given in Figure 14. With the increase of the Froude number, the range of the area of bow wave breaking becomes larger, and the wave breaking phenomenon becomes more violent.

Figure  14  Comparison of the bow wave breaking at different Froude numbers

Download: Full-Size Img

5   Conclusions and prospects

In the present paper, the phenomena of the bow wave breaking are simulated using the SPH method. The inflow and outflow boundary is implemented to improve the numerical efficiency. Numerical treatment for free surface separation is conducted to allow a smooth separation as observed in experiments.

The wave overturning at the KCS bow and the formation of the scars are observed in the SPH simulation, which is consistent with the experimental photograph, demonstrating that the SPH method can be robust and reliable in predicting the bow wave breaking of high-speed ships.

The time history results of the wave height are compared. The wave breaking at Fr = 0.425 is more severe than the Fr = 0.35, with higher wave height and larger wavebreaking region, and there are multiple obvious wave scars induced by wave overturning at the bow region at Fr=0.425. Furthermore, different types of wave breaking such as spilling breaker, plunging breaker and scars are captured in SPH simulation and the evolution of the bow wave breaking is provided.

Next, adaptive particle refinement technique can be introduced into the GPU parallel framework to achieve higher resolution accurate simulation of bow wave breaking on ships. Further, the air phase could be taken into consideration to study the effect of air entrainment on the bow wave breaking.