2.1 Thin Ship Theory
Ship waves refer to the waveform of the Kelvin wake pattern generated at the stern of an advancing ship. These waves consist of divergent waves, which propagate diagonally, and transverse waves, which propagate horizontally from the stern to the rear. The wave elevation (
ζ(
x,
y)) at a certain location (
x,
y) within three-dimensional (3D) ship waves can be expressed as the real part of the overlap of unit component waves propagating at various wave angles (
θ) relative to the ship’s direction, as shown in
Eq. (1) (
Newman, 1977). In this instance,
A (
θ), a complex-valued function representing the wave elevation of unit component waves, is referred to as the complex wave amplitude function. If the wave resistance for the hull is calculated using this wave amplitude function and the energy conservation relationship, the wave resistance value (
Rw) can be obtained using
Eq. (2) (
Newman, 1977).
where
U is the forward speed of the ship, and
ρ is the density of the fluid.
k(
θ) is the wave number of the unit component wave that propagates at each advancing angle. It is expressed as
Eq. (3), using the wave number (
k02=g/U2) for the ship direction based on the dispersion relation and taking into account the advancing angle. Here,
g represents the gravitational acceleration. The wave amplitude function (
A (
θ)) can be theoretically estimated by considering the ship’s geometry, speed, and frequency. Many previous studies have proposed various methods for calculating
A (
θ), often using approximate hull shapes. In this study, Michell’s (1898) thin ship theory was applied to analyze ship waves and wave resistance. Thin ship theory assumes that the hull width is very small relative to the ship’s length. Based on this assumption, the wave amplitude function is expressed, as shown in
Eq. (4), as an integral of the product of the micro-ship waves formed on the hull surface and the x-direction slope (
ηx (
x,
z)) of the hull surface (
Michell, 1898).
The wave elevation and wave resistance around a ship moving at a constant speed can be calculated by substituting
Eqs. (1) and
(2) into
Eq. (4). Since the analytical calculation of the hull surface slope required for determining the wave amplitude function is challenging for hulls with complex geometry, the ship geometry in this study was simplified into a hexagonal column, as shown in
Fig. 1, by referring to the study by
Zhang and Maki (2023). When the dimensions of the hexagonal column in
Fig. 1 are set equal to the ship’s dimensions, the hull surface slope is expressed as different constant values depending on the location of the bow, midship, and stern, as shown in
Eq. (5). Here,
C is the hull surface slope of the bow and stern,
LP is the length where the hull surface slope is zero, and
L is the length between perpendiculars (LBP,
LWL).
The wave amplitude function can be calculated using
Eq. (6) by substituting the hull surface slope from
Eq. (5) into
Eq. (4). Substituting
Eq. (6) into
Eq. (2) yields the wave resistance calculation formula for the simplified hexagonal column hull form, as shown in
Eq. (7):
The aforementioned ship wave and wave resistance formulas can be extended to conditions in which multiple ships navigate in formations. When multiple ships are navigating, the wave field and wave resistance are altered due to the interference of ship waves. This phenomenon can be explained using wave interference effects by linearly superposition of the ship waves generated by each hull, following the thin ship theory method presented by
Tuck and Lazauskas (1998). Specifically, when the j-th ship is located at (
xj,
yj) and the wave amplitude function of the ship wave at the (
x, y) point is defined as
Aj(
θ), the total wave amplitude function for
N ships can be expressed as
Eq. (8). If each ship in the formation has identical geometry, the wave amplitude functions become equal (
Aj(
θ)=
A0 (
θ)). If the exponential term representing the interference effect due to ship positions is replaced with
F (
θ), as shown in
Eq. (9), the wave amplitude function can then be expressed as Eq. (10). Based on this, the formulas for the wave elevation in the wave field and wave resistance for multiple ships in formation are given as
Eqs. (11) and
(12), respectively.
Assuming that the separation distance between ships is identical (
yj =
w) in a parallel arrangement (
xj = 0) where a total of
N=2M+1 ships (
M ships to the left and right, respectively) are arranged transversely, the exponential term (
G (
θ)=
F (
θ)
2) based on the ship position can be expressed as
Eq. (13). By substituting this into
Eq. (12), the wave resistance for the middle ship in a parallel arrangement can be calculated using
Eq. (14) (
Zhang and Maki, 2023):
Fig. 2 shows the results of calculating the ship waves of a single ship and those of three ships moving in a parallel arrangement with a separation distance of
w/B = 3 based on the thin ship theory at a Froude number of 0.3. For the single ship, it can be observed that transverse waves are formed behind the stern and divergent waves are generated at an advancing angle of 19.47°. Meanwhile, in a parallel arrangement, the wave field geometry is altered by the interference effects between the waves generated by adjacent ships. This interference increases or decreases the transverse wave elevation, resulting in a more complex shape for the divergent waves. To verify the calculation results of the thin ship theory applied in this study,
Fig. 3 compares the wave resistance results of the middle ship in the parallel arrangement of 15 ships (seven ships to the left and right, respectively) with the calculation results of
Zhang and Maki (2023). The comparison shows that the wave resistance change characteristics with respect to separation distance and ship speed are in good agreement with the findings of
Zhang and Maki (2023)
2.2 CFD Analysis
In this study, CFD simulation was performed to numerically analyze ship waves and ship resistance, considering the accurate ship geometry and viscosity effects. STAR-CCM+ (Version 2020.2 Build 15.04.008), a commercial CFD software program, was used. A computational mesh was constructed to analyze the flow field around the ship, and solutions to the mass, momentum, and energy conservation equations were numerically calculated by applying the finite volume method.
For the governing equations of the incompressible viscous flow field, the Navier-Stokes equation (
Eq. (15)), representing the momentum conservation law, and the continuity equation (
Eq. (16)), representing mass conservation, were applied. The flow field solution was derived by imposing appropriate boundary conditions.
The Reynolds-averaged Navier-Stokes (RANS) equation was used to analyze the wave resistance of the ship moving at a constant speed, and the
k −
ϵ model was applied as the turbulence model. To examine the ship waves generated by a ship moving over a free water surface, the volume of fluid (VOF) technique, capable of simulating ideal fluids for water and air, was utilized. The ONR-Tumblehome (ONRT) hull form was selected as the target ship, as shown in
Fig. 4.
Table 1 summarizes the specifications of the full-scale and model-scale ships.
Considering the calculation efficiency of numerical analysis, CFD simulation was performed for the fluid domain that included only half of the hull, based on the
xz-plane passing through the fixed midship (half-domain). The size of the entire domain was set to 5
LWL in the flow direction, 2.5
LWL in the width direction, and 3
LWL in the vertical direction. The trim cell mesher was used for the hull surface mesh, and 10 prism layer meshers were generated from the hull surface to ensure that
y+ remained at 30 or less. To clearly observe the ship waves, a detailed fan-shaped mesh system was formed within the ship wave advancing angle of 19.47° at the stern. An additional precise mesh domain was also created in the vertical direction near the free water surface to form a dense mesh.
Fig. 5 shows the mesh system in the fluid domain. The total number of mesh cells ranged from 2.5 to 3 million, depending on the number of ships in the parallel arrangement. A time interval of △
t = 0.01 s was applied, and calculations were performed for up to 60 s.