Research on ship motion characteristics in a cross sea based on computational fluid dynamics and potential flow theory

ABSTRACT

The motion response of a ship in a cross sea is studied based on computational fluid dynamics. Firstly, according to the established numerical pool, and based on the Reynolds-averaged Navier–Stokes equations and the re-normalisation group K–ϵ turbulence model, the free surface is treated by the volume of fluid(VOF) method, and a numerical simulation method is established. The wave results obtained by the numerical simulation are compared with theoretical waves to verify the reliability of the method. Then, considering the DTMB5415 ship model, a prediction method for the ship’s motion in a cross sea is established by using an overlapping mesh and VOF technology. The influence of a cross sea on ship motion performance under different wave direction angles, different wave heights and for different periods is analyzed. In addition, a method is also established for studying the ship’s motion response based on potential flow theory, and it is found that ships sailing under superimposed waves will exhibit the phenomenon of beat vibration. Finally, by comparing the results obtained by the two methods, the consistency of the two methods is verified, which provides a strong basis for the safety assessment of ships in cross seas.

KEYWORDS:

• Computational fluid dynamics
• cross sea
• seakeeping performance
• nonlinear

1. Introduction

The security and stability of ship and ocean engineering equipment in extreme sea conditions are becoming increasingly prominent problems, especially in the North Sea, the Black Sea, South Africa, the East China Sea, the Sea of Japan, China, the South China Sea, the Taiwan Strait, and many other sea areas. The formation of the observed waves, mainly because of their instantaneous, powerful and unpredictable nature, constitute a serious threat to the safety of marine structures. When two waves in different directions generated by two weather conditions meet, a special type of wave will be generated. When the angle between the two waves is 90°, the generated wave is called a cross sea, and this cross sea is also a kind of rogue wave.

Rogue waves have strong energy and can be destructive to marine structures, especially ships. Since ancient times, there has been the risk of frequent accidents at sea. Although some of them are caused by human factors, most of them are caused by bad sea conditions (Yıldırım et al., 2019). As a special type of rogue wave, the cross sea is a perplexing condition at sea because it can prevent crew from taking the most effective action and utilizing the most effective methods to deal with it, and thus cross sea conditions are prone to cause marine accidents.

There are two most commonly used methods to study the generation of waves and the motion response of ships sailing in waves. They are simulation research and the tank ship model test. Simulation research is divided into potential flow theory and viscosity theory. Owing to the development of the computer industry, simulation research has become a hot topic in the study of waves and ships in recent years. Song et al. (2021) performed full-scale ship resistance and self-propulsion simulations using the REGAL general cargo ship with various grid numbers and wall Y+ values and compared computational and empirical results. Zhang et al. (2021) argued that computational fluid dynamics (CFD) methods are increasingly used in the early stages in the design of ships for hull shape optimization because an optimized hull shape can achieve good hydrodynamics. Lu et al. (2020) reported CFD results concerning the airwake calculated for a reduced-scale model of an isolated landing helicopter assault (LHA) model with high-quality particle image velocimetry experimental data provided by the NASA AMES Research Center. However, it is impossible to know whether the existing technology is reliable only through numerical simulation, so it is very necessary to verify the accuracy of numerical waves through error analysis of the experimental method (Ducrozet et al., 2016). However, most of the waves at sea are rogue waves, so this is particularly important for the development of wave technology in experimental pools (Xia, 2014). On the basis of existing research, Hu (2015) proposed a wave train superposition model based on probability and a second-order nonlinear solution based on a respiratory subsolution, then simulated different types of rogue wave by using a numerical wave flume, and carried out in-depth research on their interaction with structures (such as cylinders and plates). Although several theories have been proposed to explain the formation of rogue waves, the mechanism is still unclear. In the real sea, interaction between wave and current is a factor that cannot be ignored, and the existence of currents have an important influence on the parameters, propagation and shape of waves (Peng, 2014). Benetazzo et al. (2013) used a coupled ocean–atmosphere–wave–sediment transport (COAWST) modeling system to study the influence on waves of the wave–current interaction (WCI) process in the semi-enclosed Venice Bay (northern Adriatic Sea). In fact, the formation of cross seas is also caused by the interaction of waves and current.

The Reynolds-averaged Navier–Stokes (RANS) method is a common method for studying ship seakeeping (Jiao et al., 2020), and in the case of cross seas, this method is still reliable. The earliest research on intersecting waves started with Toffoli et al. (2006), who mainly analyzed the influence of s-waves on the safe navigation of ships and the causes of, and countermeasures against, safety accidents. Onorato et al. (2013) conducted a series of experiments on large ocean waves based on the nonlinear Schrödinger equation and studied in detail the abnormal phenomena generated by the collision of two groups of rogue waves. Liu (2018) studied the breaking wave phenomenon after the intersection of two wave propagation directions. By setting up a cross-two-channel device in an experimental pool, the waves from two directions of propagation can be made to propagate through their respective channels, meet and break in a specific area of the physical process. By changing the angle between the two channels, the influence of different angles on the broken wave can be studied experimentally. In the experiment, he considered the influence of 8° and 12° angles on the breaking wave. In the weak three-dimensional experiment, he systematically studied the wave surface characteristics, spectral variation and energy loss of regular and focused waves. A cross sea is the wave formed by the intersection of two groups of waves, so, this will generate breaker.

Potential flow theory is also widely used in the seakeeping of ships. The results of Han et al. (2020), according to three-dimensional potential flow theory combined with simulation analysis of the frequency and time domains, show that, in grade 1∼5 sea conditions, ship wave resistance performance is good and can effectively guarantee the safety indexes, but in sea conditions at level 6, the swing angle exceeds the limit of the vessel, therefore, it is necessary to adjust the ship's attitude or add additional devices to ensure the ship's navigation safety. Combined with multi-thread parallel computing technology, the computational efficiency of potential flow theory is much better than that of CFD simulation (Chen, 2018; Ma et al., 2019). It is a practical numerical technology that can easily and quickly calculate the seakeeping performance of ships under various conditions. In addition, the CFD simulation results obtained by comparing the potential flow theory and CFD technology for each condition will be more convincing if the calculation results are consistent (Chang, 2019; Di, 2018; Li, 2014).

In summary, CFD is more widely accepted for the study of the navigation of ships at sea. This study also completes a series of studies of cross seas based on the CFD method. The difficulty to be addressed is to simulate the combined cross sea in different situations in a virtual pool. The parameters of the combined waves are relatively complex, including different wavelengths, flow rates and combinations. The wave direction angle requires a large amount of calculation. In addition, it is necessary to use the method of potential flow theory to verify the conclusions drawn by the viscous flow method, and then analyze the safety of a ship under different sea conditions according to the motion response of the ship under the combined waves.

2. Numerical simulation methodology

This section mainly provides a theoretical basis for the study of viscous flow and introduces the process of related work in detail. The numerical simulation process of computational fluid dynamics is described in detail below.

2.1 Governing equations

Fluid motion must satisfy three conservation laws: mass conservation, momentum conservation and energy conservation. In engineering practice, the viscous effect should be considered for fluid flow. The water wave problem contains two kinds of media, water and air, in which water is an incompressible fluid and air is a compressible fluid. In the water wave problem, the flow speed is much lower than the sound speed, so air can also be regarded as an incompressible fluid. For the incompressible problem, the energy conversion caused by temperature change can generally be ignored. The concrete form of the law of conservation of mass for an incompressible fluid is the continuity equation, and its form in Cartesian coordinates is (1) $\frac{\mathrm{\partial }\rho }{\mathrm{\partial }\text{t}}+\frac{\mathrm{\partial }\left(\rho u\right)}{\mathrm{\partial }x}+\frac{\mathrm{\partial }\left(\rho v\right)}{\mathrm{\partial }y}+\frac{\mathrm{\partial }\left(\rho w\right)}{\mathrm{\partial }z}=0$(1) where u, v, w are the velocity components in the Cartesian coordinate system, respectively, ρ is the density and t is the time.

The specific form of the law of conservation of momentum for an incompressible fluid is the equation of conservation of momentum, also known as the Navier–Stokes equation (the NS equation for short). The equation of conservation of momentum can be expressed in units of time as the rate of change of the momentum of a fluid cluster being equal to the sum of all external forces experienced by the fluid micelles. For incompressible fluids, the momentum conservation equation is as follows: (2) $\begin{array}{rl}\frac{\mathrm{\partial }\left(\rho u\right)}{\mathrm{\partial }t}+div\left(\rho u\overrightarrow{u}\right)& =-\frac{\mathrm{\partial }p}{\mathrm{\partial }x}+\frac{\mathrm{\partial }{\tau }_{xx}}{\mathrm{\partial }x}\\ & +\frac{\mathrm{\partial }{\tau }_{yx}}{\mathrm{\partial }y}+\frac{\mathrm{\partial }{\tau }_{zx}}{\mathrm{\partial }z}+F_x\end{array}$(2) (3) $\begin{array}{rl}\frac{\mathrm{\partial }\left(\rho v\right)}{\mathrm{\partial }t}+div\left(\rho v\overrightarrow{u}\right)& =-\frac{\mathrm{\partial }p}{\mathrm{\partial }y}+\frac{\mathrm{\partial }{\tau }_{xy}}{\mathrm{\partial }x}\\ & +\frac{\mathrm{\partial }{\tau }_{yy}}{\mathrm{\partial }y}+\frac{\mathrm{\partial }{\tau }_{zy}}{\mathrm{\partial }z}+F_y\end{array}$(3) (4) $\begin{array}{rl}\frac{\mathrm{\partial }\left(\rho w\right)}{\mathrm{\partial }t}+div\left(\rho w\overrightarrow{u}\right)& =-\frac{\mathrm{\partial }p}{\mathrm{\partial }z}+\frac{\mathrm{\partial }{\tau }_{xz}}{\mathrm{\partial }x}\\ & +\frac{\mathrm{\partial }{\tau }_{yz}}{\mathrm{\partial }y}+\frac{\mathrm{\partial }{\tau }_{zz}}{\mathrm{\partial }z}+F_z\end{array}$(4) where p is the pressure on the fluid micelles, τ_xx, τ_xy, etc. are due to the existence of viscosity and the viscous stress acting on the fluid micelles, F_x, F_y and F_z are the volume forces experienced by the fluid micelles.

The CFD method uses numerical methods to solve the flow in the computational domain. According to the distribution of dependent variables between nodes and the expression of discretization equations, the discrete methods in computational fluid dynamics can be divided into three main types, namely the finite element method, the finite difference method and the finite volume method. The discrete method used in the STAR CCM + software used in this study is the finite volume method. The turbulence model uses the realizable k–ϵ model that comes with the software.

2.2 Cross sea flow field and boundary conditions

In this study, the volume of fluid (VOF) method is used to realize the tracking and capture of a free liquid surface. The VOF method is a simple and effective method for capturing a free surface. The basic principle is to track the shape and position of the free surface by using the change of fluid volume inside the mesh cell and the volume function F of the mesh cell itself.

In the VOF method, a function should be defined to reflect the ratio of the volume of liquid in any mesh to the volume of the mesh, and this function should be defined as the fluid volume function F. As shown in Figure 1, the value of F is one for a fully liquid-filled mesh cell and zero for a fully air-filled mesh cell. In the simulation process, the F value of each unit is updated by constantly solving the governing equation of the F value, so as to realize the real-time capture of the free liquid surface and the simulation of the wave’s liquid surface.

Figure 1. Schematic diagram of the VOF method.

Based on the propagation and flow characteristics of a cross sea, a numerical wave pool is created. A schematic diagram and the size of the wave pool are shown in Figure 2.

Figure 2. Schematic diagram of the numerical flume model.

The length of the numerical wave pool is 24 m, the width is 24 m, the length and width of the working area are 18 m; the height is 11.5 m and the depth is 8 m. The coordinate system is the right-handed Cartesian coordinate system, as shown in Figure 3. The X- and Y-axes point to the propagation direction of the wave, and the Z-axis is vertically upward. Here, the wave flume is divided into four regions, namely the wave-making region, the front-end wave elimination region, the working region and the tail-end wave elimination region.

Figure 3. Numerical flume coordinate system.

The boundary condition settings are shown in Figure 4.

Figure 4. Schematic diagram of the boundary conditions of the numerical pool.

The mesh used in the numerical wave pool model is called a ‘trimmed mesh’. This mesh is suitable for computational fluids and is of better quality than normal tetrahedral meshes. Local encryption takes place near the free surface, as shown in Figures 5 and 6.

Figure 5. Schematic diagram of global meshing.

Figure 6. Mesh division of free liquid surface.

In order to ensure the accuracy of numerical simulation, the numerical pool is first verified. Stokes superimposed waves with a wave height h of 0.08 m and wavelengths λ of 4.576 and 6.864 m were selected for simulation. The propagation direction vectors of superimposed waves were (1,0,0) and (0,1,0), respectively. The pool velocity was u = 0, v = 0 m/s, and u = 0.5, v = 0 m/s. Table 1 shows the numerical simulation conditions of a deep water regular Stokes wave.

Table 1. Numerical simulation conditions of Stokes wave in deep water rules.

No.	h (m)	λ (m)	v
1	0.08	4.576	u = 0 m/s, v = 0 m/s
2	0.08	6.864	u = 0 m/s, v = 0 m/s
3	0.08	4.576	u = 0.5 m/s, v = 0 m/s
4	0.08	6.864	u = 0.5 m/s, v = 0 m/s

After the simulation is stabilized, the free liquid surface is shown in Figure 7. Figure 8 shows the diagonal wave height monitoring of the cross sea free liquid surface.

Figure 7. Cross sea free surface diagram.

Figure 8. Wave height monitoring of cross sea free liquid surface on the diagonal.

As can be seen from Figure 7, the superimposed wave surfaces of high and low frequencies have significant characteristics and are highly differentiated, and the superimposed wave peaks and troughs are evenly arranged, showing obvious cross sea properties. At the same time, it can be observed from Figure 8 that the numerical calculation results can simulate the shape of the Stokes waveform well, and the details of the free liquid surface at the peaks and valleys are captured in the mesh refinement area. The mesh resolution within the range of wave height and wavelength is fine, indicating that the mesh division strategy of the numerical pool is excellent.

The origin of coordinates (as shown in Figure 9) was selected to monitor its wave height, and the change of wave height with time under different flow rates u = 0, v = 0 m/s, and u = 0.5, v = 0 m/s, was compared. See Figure 9.

Figure 9. Curve of wave height over time.

According to the results in Figure 9, the maximum amplitude of the monitoring point is basically the same under both conditions of flow velocity and no flow velocity. The changing period of wave height in the condition of current velocity is different from that in the condition of no current velocity, and the period is shorter, and there are two to three small amplitude fluctuations in each period, and the same regular characteristics are presented in the condition of two wave parameters.

In addition, this paper also tests the accuracy of cross sea surface elevation results based on numerical simulation, and carries out wave height temporal monitoring at the origin of coordinates (as shown in Figure 2). Here, it is given that u = 0, v = 0 m/s, and u = 0.5, v = 0 m/s, under two different flow rate conditions at infinite water depth. A comparison between the wave height monitoring results based on numerical simulation and the wave height target value curve directly calculated based on a theoretical formula is shown in Figure 10.

Figure 10. Curve of cross sea height with time in infinite water depth.

As shown in Figure 10, the wave surface elevation of a cross sea obtained based on numerical simulation is basically consistent with the wave height target value directly calculated based on a theoretical formula, indicating that the numerical simulation method adopted in this study to complete the simulation of a cross sea has high accuracy and authenticity.

2.3 Prediction method of ship motion response

In this paper, a DTMB5415 model ship with a scale of 1:24.832 is used for calculation, and the length of the model ship is 5.72 m. The geometric model is shown in Figure 11. The horizontal line in Figure 12 shows the position of the free liquid surface. The specific parameters of the ship model are shown in Table 2.

Figure 11. DTMB5415 ship geometric model.

Figure 12. Schematic diagram of computing domain.

Table 2. Parameters of ship model.

Parameter	Values
Scale	Scale ratio = 24.832
Length between perpendiculars	L = 5.72 m
Beam	B = 0.80 m
Draft	T = 0.2248 m
Wetted surface area	SDWL = 4.861 m^2
Block coefficient	CB = 0.506
Inertia moment	(Ixx,Iyy,Izz) = (37.8741, 1911.993, 1911.993) kg-m^2
Gravity center	G = (2.86, 0.0, 0.155) m

A square computing domain and overlapping mesh technology are used in the computing domain. In the overlapping mesh method, the parent mesh covering the whole computing domain and the sub-mesh of the moving object are superimposed to describe the relative motion of the objects. The computational mesh in each region is generated independently, and the flow field information is matched and coupled at the boundary of the overlapping region by an interpolation function. This mesh technology has strong adaptability to complex shapes, and can ensure that the spatial topology of the mesh in the computing domain remains unchanged during the process of object motion. It does not need to regenerate the mesh when dealing with flow around the object surface in relative motion. The mesh generation quality is good, and the calculation accuracy is good.

The calculation domain size and boundary conditions are shown in Figure 12. The origin of the computed domain falls on the free surface and passes through the forward perpendicular, with the X-axis pointing aft, the Y-axis pointing starboard and the Z-axis pointing straight up. The calculation domain extends about twice the length of the ship in front of the ship, twice the length at the back of the ship and twice the length along the side of the ship. The bottom of the tank is about 1.5 times the length of the free surface, and the top of the pool calculation domain is about 0.4 times the length of the free surface.

Figure 13 shows the computational mesh used in this study, with local encryption near the hull surface, the free liquid surface and the Kelvin wave system. The mesh at the free surface is encrypted inn order to ensure accurate simulation of the wave flow and ensure that the mesh scale ΔX in the X-direction is not greater than 1/100 times the wavelength and ΔZ in the Z-direction is not greater than 1/20 times the wave height. Since this project needs to simulate the superimposed wave with an intersection angle of 90°, the mesh cell scale ΔY in the Y-direction is equal to the mesh scale ΔX in the X-direction. The scale ratio of adjacent meshes is set to two. A boundary layer mesh is set on the hull surface, y+ is 40 and a wall function is used.

Figure 13. Hull calculation mesh and its refinement and boundary layer mesh.

In the calculation condition, the velocity is set as v = 0.5 m/s, the velocity direction vector is (1,0,0), the included angle of the superimposed wave is always 90°, and the combined wave direction angle of the superimposed wave is (0°; 90°), (45°; 135°), (90°; 180°), (135°; 225°), the wavelength length ratios are taken as 0.8, 0.9, 1.0, 1.1 and 1.2, respectively, and the wave height h is 0.08 m. Table 3 lists the calculation conditions. The definitions of flow velocity and wave direction are shown in Figure 14.

Figure 14. Schematic diagram of wave direction angle and velocity direction.

Table 3. Calculation condition table.

Case number	Wave direction angle	Wavelength–ship length ratio
1	(0°; 90°)	0.8
2		0.9
3		1.0
4		1.1
5		1.2
6	(45°; 135°)	0.8
7		0.9
8		1.0
9		1.1
10		1.2
11	(90°; 180°)	0.8
12		0.9
13		1.0
14		1.1
15		1.2
16	(135°; 225°)	0.8
17		0.9
18		1.0
19		1.1
20		1.2

In the following, cross sea (0°: 90°), cross sea (45°: 135°), cross sea (90°: 180°) and cross sea (135°: 225°) are referred to as cross sea 1, cross sea 2, cross sea 3 and cross sea 4, respectively, for the sake of simplicity.

2.4 Mesh convergence verification

In this section, the mesh convergence is verified in the conditions of Case 3, i.e. the combined angle of the cross sea is (0°, 90°), the wavelength-to-length ratio is one and the speed is v = 0.5 m/s. A total of four sets of meshes are selected, in which the number of meshes are 2.3, 3.24, 4.59 and 5.54 million. The physical quantity calculated and compared is the maximum amplitude of the ship’s heave, and the specific results are shown in Figure 15.

Figure 15. Comparison results of the maximum amplitude of ship heave under different meshes.

It can be seen from Figure 15 that, as the number of meshes increases, the maximum amplitude of ship heave gradually decreases. It also shows that the more meshes there are, the more ideal are the calculation results. Compared with mesh 4, the error of mesh 3 is only 1.4%, and the change of heave amplitude is very small. Therefore, in order to balance the calculation efficiency and the accuracy of the results, this paper selects mesh 3 as the calculation mesh for this study.

2.5 Time step convergence verification

This paper selects five different time steps from long to short to verify the convergence of time steps, i.e. 0.05, 0.025, 0.01, 0.005 and 0.0025 s, respectively. In this section, they are defined as Step1, Step2, Step3, Step4 and Step5, respectively, in which the wave direction angle of the combined wave is (0°, 90°), the wavelength length ratio is 0.7, the wave height is 0.04 m, the velocity direction of the flow field is (1.0.0) and the velocity is 0.5 m/s.

It can be seen from Figure 16 that the heave amplitude of the ship increases with the shortening of the time step, but the results of Step 4 and Step 5 are relatively close. The results are compared in pairs, with the latter as the benchmark, the difference of the results for each time step is 75.3%, 76.2%, 40.5% and 6.8%, respectively. Therefore, considering the accuracy of the calculation results and the limitation of the calculation resources, the time step selected in this paper is Step 4, i.e. 0.005 s.

Figure 16. Comparison of maximum amplitude of ship heave at different time steps.

2.6 Numerical verification

In this section, the numerical verification of the DTMB5415 ship model is carried out according to Zong’s experiment (Zong, 2021). In the verification, the reference physical quantity is the ship roll motion response, where the wave period is 1.42 s, and the basic parameters of the ship model are consistent with the references. After the numerical calculation is completed, the stable phase of the ship’s roll motion response is compared with the experimental data, as shown in Figure 17.

Figure 17. Comparison of experimental roll motion curve and simulation roll curve.

The data in Figure 17 show that the ship roll angle of the numerical simulation is basically consistent with the experimental results, 17(a) is the experimental result and 17(b) is the numerical simulation result. It can be seen from the figure that the ship’s motion response curves of the two are basically consistent. According to the eight cycles in the figure, the average roll amplitude of the ship model simulation is 12.77°, and the value of the experimental results is 13.669°. Compared with the experimental results, the error of the numerical simulation is 6.6%, which is acceptable. According to the above numerical verification results, it can be considered that the numerical simulation method is feasible.

3. Potential flow theory

At present, the numerical calculation method based on three-dimensional potential flow theory is also one of the methods for calculating the motion response of a floating body in waves in the field of ship and ocean engineering. In this section, potential flow theory is used to study the motion response of ships in cross seas, and the accuracy of CFD numerical simulation is verified by this method.

3.1 Potential flow calculation method for the motion and load of a floating body in the frequency domain

As a branch of potential flow hydrodynamics calculation, the frequency domain analysis method can be used to study the steady motion of a floating body and its frequency response characteristics. The main research object of this section is a floating structure at sea. In addition to the basic hypothesis of space potential flow, the following hypotheses are added in the research process: the water depth is infinite, and the floating body is rigid. At zero speed, the motion of the floating body is in a stable state.

The velocity potential of the flow field in the whole fluid domain satisfies the Laplace equation, namely (5) ${\mathrm{\nabla }}^2\mathrm{\Phi }\left(X,Y,Z,t\right)=0$(5) As shown in Figure 18, the spatial fixed coordinate system O-XYZ and the body coordinate system o-xyz coincide at the origin when at rest, and both conform to the right-hand rule.

Figure 18. Schematic diagram of reference coordinate system.

The time and space variables of the velocity potential are separated: (6) $\mathrm{\Phi }\left(X,Y,Z,t\right)=\text{Re}\left\{\mathrm{\Phi }\left(X,Y,Z\right){\text{e}}^{-\text{i}\omega t}\right\}$(6) where the space velocity potential can be expressed as the superposition of the incident potential, the radiation potential and the diffraction potential: (7) $\begin{array}{rl}\mathrm{\Phi }\left(X,Y,Z\right)& ={\mathrm{\Phi }}_{\text{I}}\left(X,Y,Z\right)+{\mathrm{\Phi }}_{\text{R}}\left(X,Y,Z\right)\\ & +{\mathrm{\Phi }}_{\text{D}}\left(X,Y,Z\right)\end{array}$(7)

Because of the large inertia of the floating bodies in ships and ocean engineering, it will produce a large inertia force when moving in waves, so the inertia characteristics of the floating bodies are very important. In hydrodynamic analysis, the mass model of a floating body should be established first. The basic format of the mass model is the distribution of mass in space, i.e. the mass value and the location of the centroid are given as follows: (8) $m=m\left(x,y,z\right)$(8) In general, the mass distribution model is discrete, i.e. the block mass is given. For the selected coordinate system, the floating mass matrix can be expressed as (9) $M=\left[\begin{array}{cccccc}M& 0& 0& 0& S_z& -S_y\\ 0& M& 0& -S_z& 0& S_x\\ 0& 0& M& S_y& -S_x& 0\\ 0& -S_z& S_y& I_x& -I_{xy}& -I_{xz}\\ S_z& 0& -S_x& -I_{yx}& I_y& -I_{yz}\\ -S_y& S_x& 0& -I_{zx}& -I_{zy}& I_z\end{array}\right]$(9) where $M$ is the mass of the floating body, S_x, S_y and S_z are the static moments of the floating body mass with respect to the coordinate plane, and I_x and I_xy are the moments or products of inertia of the floating body mass with respect to the coordinate axes.

When calculating the inertia of a part of a floating body, it is necessary to synthesize the mass distribution of this part, and the inertia matrix is a function of that section of the floating body. Another characteristic of inertial matrices is that they are insensitive to small changes in mass and its distribution, mainly because the mass of floating bodies is generally large, which is why sometimes better results can be obtained by using cruder mass distribution models.

There are three free-swing modes for free-floating structures with six degrees of freedom, namely heave, roll and pitch. Taking a floating structure with two symmetric planes as an example, its static water restoring force coefficient is non-zero only corresponding to the above three motion modes, and the static water restoring force stiffness matrix is given as follows: (10) $C=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& C_{33}& 0& 0& 0\\ 0& 0& 0& C_{44}& 0& 0\\ 0& 0& 0& 0& C_{55}& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]$(10) where $C_{33}=\rho gA$, $C_{44}=\rho gV{H}_x$, $C_{55}=\rho gV{H}_y$, $\rho$ is the fluid density, $g$ is the acceleration due to gravity, A is the waterline surface area, V is the drainage volume, H_x is the transverse metacentric height and H_y is the longitudinal metacentric height. It is obvious that the static water restoring force does not have a coupling relationship in this case, which is different from the general ship-type floating body.

3.2 Wave superposition principle and beat vibration motion

In the case of a cross sea, it actually consists of two simple harmonics propagating and meeting in the same medium. In the case of potential flow, in the region where several waves meet, the displacement of a particle is equal to the vector sum of the displacements caused by each of the waves propagating separately, which is the superposition principle of waves. While the wave height of a cross sea is at the level of a small or medium amplitude, the wave intensity is small, and the two waves propagating vertically meet the superposition principle, and the motion response also meets the superposition principle. Therefore, after obtaining the downward motion response of each wave in the frequency domain, the motion response can be superimposed to obtain the motion response result under the combination of multiple wave directions.

Suppose that the steady-state responses of a certain degree of freedom movement in the frequency domain under two separate waves are (11) $\begin{array}{rl}& A_1\sin \left({\omega }_{1}t+{\phi }_1\right)\end{array}$(11) (12) $\begin{array}{rl}& A_2\sin \left({\omega }_{2}t+{\phi }_2\right)\end{array}$(12) where ${\omega }_1$ and ${\omega }_2$ are, respectively, the encountering frequencies of the two waves downwards, and it is assumed that ${\omega }_1$ > ${\omega }_2$

Then the sum of the two is (13) $\begin{array}{rl}& A_1\sin \left({\omega }_{1}t+{\phi }_1\right)+A_2\sin \left({\omega }_{2}t+{\phi }_2\right)\\ & =\left(A_1+A_2\right)\sin \frac{{\omega }_1+{\omega }_2}{2}t\\ & +\frac{{\phi }_1+{\phi }_2}{2}\cos \frac{{\omega }_1-{\omega }_2}{2}t+\frac{{\phi }_1-{\phi }_2}{2}\\ & +\left(A_1-A_2\right)\cos \frac{{\omega }_1+{\omega }_2}{2}t\\ & +\frac{{\phi }_1+{\phi }_2}{2}\sin \frac{{\omega }_1-{\omega }_2}{2}t+\frac{{\phi }_1-{\phi }_2}{2}\end{array}$(13) and then to $\sin \text{(}{\omega }_1+{\omega }_2/2)t+({\phi }_1+{\phi }_2/2)$ and $\cos \text{(}{\omega }_1+{\omega }_2/2\text{)}t+({\phi }_1+{\phi }_2/2)$ as a benchmark, others are considered to be coefficients and the combined type can be written as (14) $A\sin \frac{{\omega }_1+{\omega }_2}{2}t+\frac{{\phi }_1+{\phi }_2}{2}+\varphi$(14) The above is the synthetic expression of the response of a certain degree of freedom body under two waves, where (15) $\begin{array}{rl}& A=\sqrt{A_1^2+A_2^2+2A_1{A}_2\cos \left[\left({\omega }_1-{\omega }_2\right)t+\left({\phi }_1-{\phi }_2\right)\right]}\end{array}$(15) (16) $\begin{array}{rl}& \tan \varphi =\frac{A_1-A_2}{A_1+A_2}\tan \left(\frac{{\omega }_1-{\omega }_2}{2}t+\frac{{\phi }_1-{\phi }_2}{2}\right)\end{array}$(16) It can be seen from the above equations that the resultant response function contains the vibration frequency ${\omega }_1+{\omega }_2/2$ and is modulated by the low frequency ${\omega }_1-{\omega }_2/2$. The maximum value of the synthesized response is (17) $A_{max}=\sqrt{A_1^2+A_2^2+2A_1{A}_2}=A_1+A_2$(17) When different frequencies are encountered, the resultant motion response will appear as the phenomenon of ‘beat vibration’, as shown in Figure 19.

Figure 19. Schematic diagram of beat vibration.

The amplitude of the resultant waveform will change slowly from time to time periodically.

3.3 Ship motion response under unidirectional waves based on potential flow theory

Based on DTMB5415, which was used in the numerical simulation of ship motion response in a cross sea in the theoretical calculation of potential flow, this computational ship was continued to be used to build a hydrodynamic mesh of the whole ship, and the cubic spline difference method was used to determine the mesh nodes along the longitudinal and vertical directions, respectively. The final hydrodynamic mesh model is shown in Figure 20.

Figure 20. Hydrodynamic mesh model of the DTMB5415 ship.

In order to verify the accuracy of viscous flow numerical simulation, the motion response of ships under unidirectional waves should be predicted according to the principles of potential flow theory. In the following, ship motion response curves under unidirectional waves obtained by potential flow theory can be combined to obtain the ship’s motion response curves under a cross sea, which can be compared with the results obtained by viscous flow numerical methods.

Consistent with the calculation conditions of numerical simulation, the single wave directions based on potential flow theory are 0°, 45°, 90°, 135° and 180°, and the waves’ natural frequencies are determined according to the wavelength–length ratios (0.8, 0.9, 1.0, 1.1 and 1.2) to be 0.74, 0.69, 0.66, 0.63 and 0.60 rad/s, respectively.

Figures 21–23 show the spectrum analysis results of the motion response under the calculation conditions of each single wave direction. The motion responses of heave, pitch and roll, the three degrees of freedom, are mainly investigated here. According to the above calculation conditions, five wave direction angles, 0°, 451°, 90°, 135° and 180°, were selected and their wavelength–length ratios were one.

Figure 21. Heave RAO curve.

Figure 22. Pitch RAO curve.

Figure 23. Rolling RAO curve.

4. Results

In a cross sea, there must be a certain difference between the hull disturbance and the absence of hull disturbance, and there must be a certain difference in the wave height change; under the same cases, the calculation results of CFD numerical simulation and potential flow theory are also contrasted for significance of the difference.

4.1 Prediction results of wave height for a ship sailing in a cross sea

Considering the disturbance of ship while sailing, the wave height of a cross sea is bound to be affected, but the disturbance situation will be different in different places. Therefore, it is necessary to compare and verify the wave height values of the two. Here, wave height changes at five points around the ship are monitored, as shown in Figure 24, in which the wave direction angles of the cross sea are (0°; 90°) for cross sea 1, the wavelength–ship length ratio is 1.0, the wave height H is 0.08 m and the flow rate V is 0.5 m/s.

Figure 24. Position of wave height monitoring point of ship sailing in a cross sea.

Based on the wave height monitoring points in Figure 24, the curve of wave height change with time at each monitoring point was obtained by numerical simulation, as shown in Figure 25. The wave height target value directly calculated based on the theoretical formula is also given in the figure, which represents the wave height change value without ship rigid body disturbance.

Figure 25. Time-history variation curves of wave height at different monitoring points in a cross sea: (a) monitoring point 1 (0°); (b) monitoring point 2 (45°); (c) monitoring point 3 (90°); (d) monitoring point 4 (135°); (e) monitoring point 5 (180°).

The distribution of the five wave height monitoring points in Figure 24 is based on the intersections of 0°, 45°, 90°, 135° and 180°, which is the wave height observation point, with the center of the ship as the circle and 3.86 m as the radius. As shown in Figure 25, the numerical results at wave height monitoring point 1 are different from those directly calculated based on the theoretical formula, because the ship sails in a cross sea. The bow wave will sure be different from the direct calculation result obtained by the theoretical formula in the pure wave case. The difference is more obvious at wave height monitoring point 5 because the stern wave will disturb the flow of the cross sea. Observation monitoring point 2, 3 and 4 of the numerical results of wave height, owing to its monitoring distance, so the disturbance is small, and the results of wave height are basically consistent, in line with the real physical phenomena. It is further proved that the reliability and accuracy of the numerical simulation method in simulating the cross sea. It is further proved that the reliability and accuracy of the numerical simulation method in simulating the cross sea have a certain guarantee.

4.2 Consistency between numerical simulation and potential flow theory

According to the calculation settings in Table 3, the initial wave surface of cross sea ship motion and the wave surface of a fully developed flow in the various cases are shown in Figure 26.

Figure 26. Contrast diagram of initializing wave surface and fully developed wave surface of a ship sailing in a cross sea for various cases. Case n: (a) the initial wave surface; (b) the fully developed wave surface.

As shown in Figure 26, at the initial time of calculation, the intersection positions of the first wave front of the superimposed waves were all set at the origin of coordinates (the position of the free liquid surface of the ship’s bow). With the increase of calculation time, waves gradually developed fully in the calculation domain. Finally, wave surface characteristics of different wave direction angles under the same wavelength calculation conditions were basically the same without significant difference. For different wavelengths, the wave characteristics are naturally different, and will greatly affect the motion of the ship.

It can be observed from Figure 27(a) that the amplitude of heave is positively correlated with wavelength, and the law is basically consistent under different wave direction conditions. In Figure 27(b), the amplitude of ship pitching response under cross seas 1 and 3 has a linear relationship with wavelength, and the overall amplitude is smaller than that under cross sea 2 and cross sea 4. It can be seen from Figure 27(c) that the amplitude of the roll response does not vary significantly with wavelength, and only presents a certain linear relationship under cross sea 2. Additionally, the maximum amplitude of the roll response under all wave directions is around 5°, except under cross sea 4.

Figure 27. Amplitude of motion response: (a) heave motion response; (b) pitch motion response; (c) roll motion response.

The maximum motion amplitude can only reflect part of the characteristics of ship motion, and the addition of a time variable can more intuitively represent a ship’s motion response. Figure 28 shows the graph curve of ship motion changing with time under each cross sea parameter.

Figure 28. Curve of ship’s heave motion over time in various cross seas: (a) cross sea 1; (b) cross sea 2; (c) cross sea 3; (d) cross sea 4.

It can be seen from Figure 28 that the curve of heave response value changing with time shows different variation rules under the conditions of various waves: under the conditions of cross sea 1 and cross sea 3, the response value of ship heave motion changes relatively little with time; under the conditions of cross sea 2, the response value changes most dramatically with time. Under the conditions of cross sea 4, the heave response period is consistent with the wave period, and the amplitude changes periodically with time.

Figure 29 shows the time-varying curve of ship pitch motion response in a cross sea with different combined wave angles.

Figure 29. Curve of a ship’s pitch motion over time in various cross seas: (a) cross sea 1; (b) cross sea 2; (c) cross sea 3; (d) cross sea 4.

It can be seen from Figure 29 that the variation law of a ship’s pitch response under the conditions of various waves is basically consistent with the heave response result, and the pitch amplitude decreases with time to a certain extent. Under the conditions of cross sea 4, the ship’s pitch response period increases.

Figure 30 shows the time-varying curve of ship rolling motion response in a cross sea with different combined wave angles.

Figure 30. Curve of ship’s roll motion over time in various cross seas: (a) cross sea 1; (b) cross sea 2; (c) cross sea 3; (d) cross sea 4.

It can be seen from Figure 30 that, under the conditions of cross seas 1 and 3, the ship’s rolling amplitude changes with time in a similar way, and the response period is consistent with the wave period, while under the conditions of cross sea 4, the roll amplitude is basically zero. Combined with the potential flow theory results, it can be seen that the combined wave is symmetrical with the sailing direction of the ship at this time, and the effect of mutual cancelation is generated in the direction (0, 1, 0), so the roll amplitude is almost zero.

In Section 3.3, three response amplitude operator (RAO) curves of the ship under the conditions of potential flow theory have been obtained. The RAO of the wave parameter motion response with a wavelength–length ratio of one is selected for superposition, and the heave, pitch and roll motion response time curves of the ship under various cross sea conditions can be obtained. In order to compare with the ship motion characteristics obtained from the potential flow theory, the data consistent with the wave parameters were selected from the CFD simulation results for data collation, and the corresponding motion curves obtained.

Figures 31 and 32 show the comparison of results obtained by the two methods under the conditions of cross sea 1 and cross sea 3, wherein the abscissas of the left-hand sides of the figures is the number of iteration steps (NOI), and the abscissas of the right-hand sides of the figure is time, the corresponding relationship is that 100 iteration steps are equal to ones.

Figure 31. Curve of ship’s motion over time in cross sea 1. (a) Heave: (1) CFD, (2) potential flow; (b) pitch: (1) CFD, (2) potential flow; (c) roll: (1) CFD, (2) potential flow.

Figure 32. Curve of ship’s motion over time in cross sea 3. (a) Heave: (1) CFD, (2) potential flow; (b) pitch: (1) CFD, (2) potential flow; (c) roll: (1) CFD, (2) potential flow.

Figures 31 and 32 show that the results of the CFD numerical simulation are in good agreement with the results obtained by the potential flow theory. The movement trend of the ship has a strong consistency, and there is an obvious phenomenon of beat vibration, i.e. the ship shows a regular range of motion in multiple cycles. The superposition of waves will lead to the beat vibration of the waves themselves, and the same phenomenon will be produced by the ship through energy transfer. It also shows that validity verification has been achieved. The specific comparison results are shown in Table 4.

Table 4. Validity verification of ship motion response calculation results in a cross sea based on potential flow theory and numerical simulation.

	Heave(m)			Rolling (°)			Pitch (°)
Cross sea	Potential flow	CFD	Discrepancy	Potential flow	CFD	Discrepancy	Potential flow	CFD	Discrepancy
(0°; 90°)	0.065	0.058	10.8%	4.20	4.80	14.3%	1.45	1.25	13.8%
(45°; 135°)	0.035	0.030	14.3%	4.40	4.20	4.5%	2.80	2.40	14.3%
(90°; 180°)	0.045	0.040	11.1%	4.50	4.80	6.7%	1.00	0.90	10.0%
(135°; 225°)	0.042	0.040	4.8%	0.045	0.042	6.7%	2.40	2.10	12.5%

It can be seen from Table 4 that the difference of the calculated results in different cross seas is in the range 4.5%–14.3%, which is acceptable, indicating that the numerical simulation results of ship motion response are in good agreement with those based on potential flow theory, and the numerical simulation results obtained are acceptable. At the same time, the above is the result of ship motion response under unit wave amplitude. Multiply the amplitude of the actual cross sea as a correction to get the result of ship motion response under real-time sea state wave height.

5. Conclusion and discussion

In this paper, the motion prediction of ship model DTMB5415 is successfully achieved by CFD simulation and potential flow theory. According to different crossing parameters, the motion response curves of heave, pitch and roll are obtained, and the results obtained by the two methods are compared. The results show that the CFD simulation results are consistent with the potential flow theory results. In the above comparison data, the error values of roll, pitch and heave results are all about 10%. In particular, the roll error values for cross sea 2, cross sea 3 and cross sea 4 are about 5%. In general, the difference between the results of the two methods is small, so both methods can be used as reliable methods for ship seakeeping prediction. To sum up, the reliability and accuracy of the ship motion response in the CFD numerical simulation of cross seas have also been effectively verified by the proof of the theoretical potential flow results.

The common phenomenon obtained by both methods is that the heave of a ship is always proportional to the wavelength. That is, no matter what the wave direction angle of the cross sea is, the larger the wavelength is within a certain range, the larger the ship’s heave amplitude will be. Next, it can be seen in Figure 29 that, in the case of cross sea 2, the heave frequency of the ship is faster than in the other cases. Therefore, from the perspective of ship navigation safety, in any cross sea field, the heave stability of a ship needs to be paid attention. Especially when the synthetic wave wavelength is too large, the heave motion of the ship needs to be observed when encountering the situation of cross sea 2, and it is necessary to pay attention to the speed of the ship because reducing the speed can make the ship heave at a relatively low frequency. The ship’s pitch law is similar to that of heave, but in the case of cross sea 2 and cross sea 4, the ship’s pitch amplitude is too large, and the former’s motion response period is short, so the ship’s pitch stability needs to be focused on at this time. Therefore, when encountering the situations of cross sea 2 and cross sea 4, if the pitch frequency and amplitude are too large at this time, in addition to reducing the speed, adjusting the sailing direction of the ship is also a very important maneuvering method, and the ship should try to sail obliquely to the cross sea. In the case of cross sea 4, owing to the symmetrical distribution of the combined waves and the ship’s sailing direction, the ship’s roll amplitude is minimum at this time. Except for this case, the ship’s roll period is relatively short, the ship’s motion is unstable, and the relationship with wavelength is not obvious. Therefore, in a cross sea field that is not symmetrical with the ship’s sailing direction, the roll stability of the ship needs to be focused on. Therefore, in such cases, if the rolling frequency and amplitude are too large, the ship should run in the same direction as the cross sea, so that the wave and the ship’s running direction are symmetrical.

The conclusions drawn in this study have some limitations. For the combined wave direction angle, this paper only considers several typical cases, and the combined wave conditions directly ahead are ignored here, i.e. the (−45°, 45°) combined wave. In addition, the flow velocity of the combined wave itself is also considered more conventionally; the cross wave in reality should be more complex, so this study is not comprehensive enough in this regard. In the authors’ forthcoming work, they will consider the relative flow velocity and direction of arrival of more combined waves, and will also consider the case of three or more combined waves.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

Additional information

Funding

This work was funded by the National Natural Science Foundation of China General Program Project [Grant No. 52171293]; the National Natural Science Foundation of China Youth Project [Grant No. 51809029]; the Ph.D. Scientific Research Fund of the Natural Science Foundation of Liaoning Province [Grant No. 2019-BS-025]; the high-level talent innovation support program of Dalian [Grant No. 2020RQ009]; Key discipline project of Dalian Science and Technology Innovation Fund [Grant No. 2020JJ25CY016]; the Fundamental Research Funds for the Central Universities [Grant No. 3132022112].