**1. Introduction**

Stability problems associated with high-speed planing crafts have long been a notable research focus for designers, even in calm waters. It is well known that due to the existence of longitudinal or transverse instability, many kinds of dangerous accidents may occur [1]. The abrupt variation of trim causes the self-induced heave and pitch oscillations, which has been named porpoising [2]. In severe cases, the bow even suffers a violent attack and generates greater resistance in fast vessels, which threatens the safety of on-board personnel and equipment. The transverse instability causes the sudden emergence of large heeling, leads to a loss of course-keeping ability, and facilitates chine walking [3].

Considering the reasons outlined above and the associated adverse situations, it is important for engineers to control the longitudinal and transverse instability of high-speed crafts. In terms of the longitudinal instability, since the porpoising was observed during the test conducted by Clement and Blount [4], Blount and Codega [5] analyzed the dynamic instability problem, focusing on the nonlinear vibration of the planing craft through the test method, and the suggested criterion conditions for the unstable motion were given. Katayama and Yoshiho [6] conducted a series of performance tests on planing crafts, concretely involving accelerated longitudinal motion and porpoising instability, which provides crucial comparative data for many scholars. In recent years, porpoising theory has been applied in some studies based on established numerical methods [7] and Savitsky's

**Citation:** Wang, J.; Zhuang, J.; Su, Y.; Bi, X. Inhibition and Hydrodynamic Analysis of Twin Side-Hulls on the Porpoising Instability of Planing Boats. *J. Mar. Sci. Eng.* **2021**, *9*, 50. https://doi.org/10.3390/ jmse9010050

Received: 2 December 2020 Accepted: 31 December 2020 Published: 5 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

work [8,9]; those studies were primarily concerned with the inhibition of trim instability by appendages such as interceptors, trim tabs, and wedges [10–12]. The interceptor is a thin vertical plate protruding from the stern that is installed at the bottom edge. Mehran et al. [13] ascertained that the inhibition mechanism of the interceptor on porpoising in the planing boat, and Mansoori et al. [14,15] further analyzed the influence of boundary layer thickness, interceptor height, and span on the inhibition of porpoising and navigation resistance, and demonstrated that the combination of trim tab and interceptor with the same size is more beneficial to control the trim and reduce the resistance than the single interceptor. In addition, to develop a craft design utilizing an interceptor that included six heights at distinct positions of the stern bottom, Ahmet and Baris [16] tested resistance and sailing attitudes and reported that for the same-size interceptor, the effect of drag reduction and porpoising inhibition was gradually reduced when installed at the interval from keel to bilge line. Hongjie et al. [17] and Hanbing et al. [18] calculated the porpoising of a planing boat in a uniform incoming flow, indicating that moving forward of the center of gravity could reduce the resistance peak value, which is beneficial to avoid porpoising.

Regarding transverse instability research, owing to the interference of many nonlinear factors and high calculation costs, most studies merely involved model test and improving the calculation method, as discussed in Refs. [19–25]. This research is not involved in the transverse stability of the planing boat because the installing of twin side-hulls widens the transverse weight distribution of the vessel, the transverse righting moment increases, which makes the transverse stability of the craft better than that of the monomer-form.

However, well-known porpoising instability restrains the maximum speed of highspeed planing crafts, but it can be controlled by using external devices, such as trim tabs, interceptors, side-hulls, etc. The doctoral thesis of Yi [26] indicated that a high-speed trimaran planing craft was beneficial to improve longitudinal stability, but it increased the resistance at lower speeds. So a conceptual planing boat, capable of freely retracting and releasing the twin side-hulls, is proposed in this research. At lower speeds, it maintains the monomer-form state (MFS), and at higher speeds, it expands into the trimaran-form state (TFS) to inhibit the porpoising and improve the longitudinal stability. But due to the strong interference of nonlinear factors and the high test cost, there are no or few studies analyzing the inhibition of twin side-hulls on porpoising.

Therefore, the hydrodynamic performance of twin side-hulls and their inhibiting effect on the porpoising instability in high-speed planning crafts were analyzed in this research. And a series of tests were conducted on the monomer-form models with different longitudinal locations of the gravity center (*L*cg), and the influence of the *L*cg on porpoising was determined. Then, based on the comparisons of the numerical (CFD) setup, simulations of the test model were performed, and the calculated results were verified using the whisker spray equation of Savitsky [27]. Further, the inhibition mechanism of the side-hulls on porpoising was ascertained, and a comparative analysis of the hydrodynamics of the planing craft in the MFS and TFS was performed. Finally, the influence of longitudinal and vertical side-hull locations on inhibiting porpoising was determined, and the optimal location range is provided.

#### **2. Experimental Setup**

#### *2.1. Model Design*

This study took a 1:2.5-scale test model of an actual planing hull with twin positionadjustable side-hulls as a research object. The twin side-hulls can be synchronously placed in any longitudinal, horizontal, and vertical positions of the ship broadsides within a certain scale range by the variable-structure link bridge. When on a relatively stable seas, the craft will pack up the twin side-hulls and sail forward quickly in the monomer-form state (MFS), as shown in Figure 1a. When encountering high sea conditions, it will put down the twin side-hulls and sail stably in the TFS, as shown in Figure 1b. Main geometric details of the planing hull in different navigation states are listed in Table 1.

**Figure 1.** Main geometric characteristics of the hull in different navigation states: (**a**) monomer-form state (MFS) and (**b**) trimaran-form state (TFS).

**Table 1.** Primary geometric details of the planing hull at each navigation state.


The hard chine test model (main hull), built-in wood, was square-tailed and nonstepped, had a larger knuckle line width, a plurality of spray deflectors with variant dimensions were symmetrically installed at the main hull bottom to improve its seakeeping and reduce resistance in calm water. The twin side-hulls were arranged at both broadsides of the main hull by a link bridge, were thin and sharp in cross-section, slender in whole and vertical in the bilge, provide bare buoyancy and hydrodynamic lift but improved directional stability. The molded lines of the main hull and side-hulls are presented in Figure 2. For convenient adjustment to the center of gravity (CG) position in the TFS and making the model easier to slide, the twin side-hulls were initially designed on both sides of the rear of the main hull as shown in Figure 1b, and the specific location is listed in Table 1.

**Figure 2.** Molded lines of the model: (**a**) main hull; (**b**) side-hull.

#### *2.2. Model Test*

All model tests were conducted in the towing tank (510 m × 6.5 m × 6.8 m) of the High-speed Hydrodynamic Laboratory of Special Aircraft Research Institute of China. Owing to the higher test cost, only the towing tests of the model in the MFS were completed, and considering that the twin side-hulls were far higher than the waterline when packed up in the MFS, both were not installed in the test process, the schematic diagram of the experimental setup is presented in Figure 3.

**Figure 3.** Schematic diagram of the experimental setup.

The model was attached to a carriage (0.1%), and two guide robs were also fixed in the carriage, separately inserted into the guide plates of the bow and stern to prevent yaw and roll motion; the towing points were located at the broadside and aligned with the center of gravity (CG). The cable-extension displacement sensor (0.01 mm) and the gyroscope (0.01◦) were fixed at the CG to measure the sinkage and dynamic trim, the electronic dynamometer (0.02 kg) was mounted on the carriage, pulled the tow bar for measuring the resistance, and the accuracy of the above measuring instruments is presented in the corresponding brackets. Moreover, the front and rear cameras were fixed on the forward and aft of the carriage to capture the flow phenomena of the model sailing in calm water.

In the model tests, aimed at the two displacements and three longitudinal positions of the CG, four conditions were designed, as shown in Table 2. The towed speeds of the model were 3–13 m/s (length Froude numbers *F*r = 0.63–2.72) or until porpoising occurred. The wave surface condition around the model and the stern wake in the distance were clearly captured, as presented in Figure 4, which shows the photographs from various perspectives of the model at *F*r = 1.26 under condition two, and subsequently, a Longitudinal location of CG/Main hull length (*L*cg/Lm) ratio of 0.38 occurred under condition two.


**Table 2.** Model test conditions in calm water.

**Figure 4.** Photographs of the monomer-form model.

#### *2.3. Experimental Results*

The measured parameters mainly include total resistance *R*T, sinkage *Z*cg, and trim angle *τ* when sailing stably; the *R*<sup>T</sup> was treated into the dimensionless form of *R*T/Δ. The test results on the *R*T, *Z*cg, and τ of the monomer-form model under different conditions are shown in Figure 5.

**Figure 5.** Test results of the monomer-form model: (**a**) total resistance; (**b**) sinkage and trim angle.

Combing test phenomena and results, we found that at low speeds (*F*r < 1.05), the bow was gradually lifted upwards as the speed increased. At the semi-planing state (*F*r = 1.05), the sinkage and trim angle increased significantly. When crossing the resistance peak, entering the planing regime (*F*r > 1.05), the *R*<sup>T</sup> appeared to be notably reduced; however, after descending to a certain threshold, the *R*<sup>T</sup> rose and exceeded the previous peak as the speed further increased. When *F*r > 2.31, the *Z*cg increased slowly, but the *τ* always reduced as the speed increased during the planing stage.

For the model with equal displacement during the semi-planing state, as the *L*cg moved backward, the resistance peak increased. When *F*r > 1.05, as speed further increased, the more the *L*cg moves backward within a certain scale, the more the *R*<sup>T</sup> decreased, but the *Z*cg and *τ* increased. In addition, when the *L*cg remained at the same and speed increased, the *R*T, *Z*cg, and *τ* of the small-displacement boat were all less than those of the largedisplacement boat, especially for the small-displacement boat, *R*<sup>T</sup> is lesser when crossing the resistance peak.

The porpoising of the model under conditions two and three occurred at *F*r = 2.73 (*v* = 13 m/s) and 2.52 (*v* = 12 m/s), respectively. From the videos recorded in the experiment, we observed that the coupled heave and pitch oscillation amplitudes under condition three were extremely larger compared with condition two, which shows that when *L*cg becomes

increasingly backward, porpoising occurs more easily, and the oscillation amplitudes are even larger it limits the maximum speed.

#### **3. CFD Setup**

#### *3.1. The Numerical Method*

To simulate the viscous flow field around the sailing vessel, the governing equations of viscous incompressible fluid were introduced, and that were solved based on the Finite Volume Method (FVM) described in Ref. [28]; the main CFD solver is the platform of Star-CCM+, the Reynolds-averaged Navier–Stokes (RANS) and continuity equations jointly constitute the governing equations, and as follows

$$\frac{\partial(\rho u\_i)}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_i}(\rho u\_i u\_j) = -\frac{\partial p}{\partial \mathbf{x}\_i} + \frac{\partial}{\partial \mathbf{x}\_j} \left( \mu \frac{\partial u\_i}{\partial \mathbf{x}\_j} - \rho \overline{u\_i' u\_j'} \right) + \mathcal{S}\_i \tag{1}$$

$$\frac{\partial \rho}{\partial t} + \frac{\partial (\rho u\_i)}{\partial x\_i} = 0 \tag{2}$$

where *u*<sup>i</sup> and *u*<sup>j</sup> are the time mean of the velocity component, (*i*, *j* = 1, 2, 3), *p* is the pressure mean, ρ is the fluid density, μ is the coefficient of dynamic viscosity, *ρu i u <sup>j</sup>* is the Reynolds stress term, *S*<sup>i</sup> represents the generalized source term.

To close the governing equations, the Shear Stress Transport turbulence (SST) k–ω model [29] was adopted to calculate the Reynolds stress term in this research. Despite the existence of the strong adverse pressure gradient, the RANS model has been shown to be inaccurate and has some limitations as discussed in Refs. [30–32], but the SST k–ω turbulence model is still widely used to deal with high Reynolds number flow problems and has higher precision on solving the flow field around the high-speed planing craft [33].

Moreover, the volume of fluid (VOF) method, proposed by Nichols and Hirt [34,35] in 1981, was also used to track the change of the free surface in this research. In the VOF method, the most critical aspect is obtaining the volume ratio function (*F*) of the specified fluid occupancy in a grid cell. When the calculation of the *F* values in each grid cell is completed, the motion interface of the liquid-gas, two-phase flow can be tracked.

Then, to acquire the hull position, based on velocity and pressure in the flow field, the centroid motion theorem and centroid moment of motion theorem were used as follows

$$\begin{array}{c} \stackrel{d\overrightarrow{B}}{\frac{d\overrightarrow{B}}{dt}} + \stackrel{\rightarrow}{\Omega} \times \stackrel{\rightarrow}{\dot{B}} = \stackrel{\rightarrow}{F} \\ \stackrel{d\overrightarrow{K}}{\frac{d\overrightarrow{K}}{dt}} + \stackrel{\rightarrow}{\Omega} \times \stackrel{\rightarrow}{K} + \stackrel{\rightarrow}{V} \times \stackrel{\rightarrow}{\dot{B}} = \stackrel{\rightarrow}{\dot{M}} \end{array} \tag{3}$$

where <sup>→</sup> *<sup>B</sup>* is the momentum of the model, <sup>→</sup> <sup>Ω</sup> (*p*, *<sup>q</sup>*, *<sup>r</sup>*) is the angular velocity, <sup>→</sup> *F* (*X*, *Y, Z*) is the combined force, <sup>→</sup> *<sup>K</sup>* is the momentum moment relative to the CG, <sup>→</sup> *V* (*u, v, w*) is the speed and <sup>→</sup> *M* (*L*, *M*, *N*) is the combined moment.

→ *<sup>F</sup>* and <sup>→</sup> *M* can be solved as below

$$\stackrel{\rightarrow}{\vec{F}} = \int\_{S} ([\pi] - P[I]) \cdot \stackrel{\rightarrow}{n} ds - \stackrel{\rightarrow}{G} \tag{4}$$

$$\stackrel{\rightarrow}{M} = \int\_{S} \left( \stackrel{\rightarrow}{r} - \stackrel{\rightarrow}{r}\_{G} \right) \times \left( [\pi] - P[I] \right) \cdot \stackrel{\rightarrow}{n} ds \tag{5}$$

where [*τ*], *<sup>P</sup>*[*I*] and <sup>→</sup> *G* are the shear stress, pressure, and gravity, respectively. *S* is the hull surfaces. <sup>→</sup> *r* is the displacement of mesh nodes, and <sup>→</sup> *r <sup>G</sup>* is the displacement of CG.

The overset-grid method, as described in Ref. [36], was adopted due to the complex hull motion. The solver procedure of the numerical method is shown in Figure 6.

When variations of the forces and moments were less than the tolerance (*ε*) or the total iteration time (*T*) reaches, the calculation was terminated.

**Figure 6.** Solver procedure of the numerical method.

In the post-processing, the pressure, shear force, unit node coordinates at each unit are known; thus, the longitudinal moment of the side-hulls to CG can be calculated by integrating the element moments to the CG, and the dimensionless form as follows

$$\mathbf{C}\_{\text{M}\_{\text{d}}} = \frac{\iint\_{\text{S}\_{\text{d}}} \left[ P \cdot \stackrel{\rightarrow}{I} \cdot \stackrel{\rightarrow}{n} \right] ds + \iint\_{\text{S}\_{\text{d}}} \left[ \tau \cdot \stackrel{\rightarrow}{I} \cdot \stackrel{\rightarrow}{n} \right] ds}{\text{Ag}\text{L}\_{\text{m}}} \tag{6}$$

where <sup>→</sup> *I* and <sup>→</sup> *n* are the local normal vector of the grid element and the displacement relative to the CG, respectively, *S*<sup>d</sup> represents the surface area of twin side-hulls.

Likewise, adopting the same integral strategy, the dimensionless forms of the resistance and lift for the twin side-hulls are acquired as follows:

$$C\_{R\_d} = \frac{\iint [\pi] ds}{\Delta \mathbf{g}} \tag{7}$$

$$C\_{N\_{\rm d}} = \frac{\iint [P] ds}{\Delta \mathbf{g}} \tag{8}$$

#### *3.2. Computational Domains and Boundary Conditions*

To prevent the reflection of waves in the computational domain and obtain a better precision, the domain should be no less than five times the hull length (Lm) [37]; thus, it was designed to a cuboid region with dimensions of 7.5 Lm × 2 Lm × 2.6 Lm, and its specific dimensions and the boundary conditions were depicted in Figure 7. Considering

the symmetry of the model and flow field, only half a domain of the model was established to reduce the simulation duration.

**Figure 7.** Dimensions and boundary conditions in the domain.

To obtain a better simulation on the sailing attitudes of model, the overset region was embedded in the background region; the boundary conditions were set as follows: The background region, incoming-flow inlet, top, bottom, and side were all set to the velocity inlet, the outlet was the pressure outlet, and the middle longitudinal section was a symmetry plane; for the overset region, the mid-ship section was still a symmetry plane, the ambient planes were set to the overset mesh boundaries, and the hull body was defined as the non-slip wall.

#### *3.3. Mesh Generation*

The grids in the background and overset region were automatically generated, as shown in Figure 8, and the grids of the two regions defined as follows: The cutting hexahedron grids were mainly used to discretize the background region, both the cutting hexahedron grids and the prismatic layer grids were selected to fill the overset region due to the complex geometrical details of the hull body. To obtain the more accurate flow field around the hull, the circumambience of the hull was refined by the smaller hexahedron grids in the volume control region; for the grids around the free surfaces in the two regions, that were set to more than twenty layers in the vertical direction to clearly capture the change of free surface; and the boundary layer grids were imposed on the surface of hull body.

#### *3.4. Wall Non-Dimensional Distance of the First Layer Grid (y+) and Time Step Set-Up*

To obtain the accurate stress and pressure of the flow field around the hull, wall functions, and boundary layer grids were required, referred to the Ref. [37]. The *y*+ value was ascertained as below

$$y^{+} = \frac{y}{\upsilon} \sqrt{\frac{1}{2}} \mathcal{U}^2 \frac{0.074}{\text{Re}\_L \frac{1}{2}} \tag{9}$$

where *y*+ and *y* separately represent the non-dimensional distance and the height of the first layer grid, *U* is the speed of vessel, *L* is the waterline length of the vessel, *υ* is the fluid viscosity coefficient and Re represents the Reynolds number.

In this research, five boundary layers were adopted with a grid growth rate of 1.3. The generated boundary layer grids on the hull body are presented in Figure 8d. The details of when the *y*+ = 250, the *y* values, and total boundary layer grid reached thicknesses Δ*h* at different speeds are shown in Figure 9.

In addition, considering the time step of iteration (Δ*t*) depends on flow characteristics in the implicit unsteady simulation, as discussed in Ref. [38], the range of Δ*t* is initially defined as follows:

$$
\Delta t = (0.005 \sim 0.01) \cdot \text{Lm/} \Omega \tag{10}
$$

The maximum number of internal iterations and total calculation time were set to five times and 15 s, respectively.

**Figure 8.** Detailed grid partition of the entire computational domain: (**a**) axis view; (**b**) front view; (**c**) vertical view; (**d**) partial view.

**Figure 9.** The height of the first layer grid (*y*) values and total thicknesses of the boundary layer grid at different speeds.

#### **4. Numerical Simulation Verification**

#### *4.1. Grid Parameters*

To verify the convergence of the grid, four grids were designed, whose sizes on the hull surface and in the background region were the same, but the grid sizes in the overset region, volume control region, and around the free surface increased with the refinement

ratio (*<sup>r</sup>* <sup>=</sup> <sup>√</sup>2, Lg = 1 m), the grid parameters of four grids are shown in Table 3, and the grids around the hull and the free surface are presented in Figure 10. Then, the sailing of the model in calm water was simulated, and the results of resistance were compared with the test. The numerical setting is consistent with the Section 3, and the towing speed *F*r = 1.26 (*v* = 6 m/s) was selected for verification.

**Table 3.** Parameters of the four grids.


**Figure 10.** The four designed grids: (**a**) grid 4, (**b**) grid 3, (**c**) grid 2, (**d**) grid 1.

Table 4 and Figure 11 indicate that the resistance curves of the four grids have good convergence; comparing the numerical results of the four grids, the coarse grid 4 had a larger deviation compared with test results, and as the grid was refined, the deviation gradually decreased. However, when the grid was fine enough (grid 2), the further refinement of grid 1 did not greatly improve the simulation accuracy; inversely, the calculation efficiency dropped significantly due to too many grids.

**Table 4.** Resistance comparisons of the test (EFD) and numerical (CFD) results for different grids.


Summing up the above, considering the high accuracy and the less computing time at the same time, the parameters of grid 2 were selected to be the more suitable grid input setting, and subsequent grid settings were all based on grid 2.

**Figure 11.** Resistance curves of the different grids.

#### *4.2. Y+ Values*

In general, the value of *y*+ has a profound influence on the calculation accuracy. In 2005, Wang et al. [39] analyzed the influence of y+ on the turbulence problem, indicating that the *y*+ value near the hull should be restrained between 30 and 300. Thus, to ascertain the influence of *y*+ value on the calculation accuracy, six different *y*+ values of 50, 100, 150, 200, 250, and 300 were selected, and as in the previous CFD set up, simulations of different *y*+ values were conducted, and the wall *y*+ and resistance after calculation were obtained, which are shown in Figure 12 and Table 5, respectively.

**Figure 12.** Wall non-dimensional distance of the first layer grid (*y*+) after calculation: (**a**) *y*+ = 50; (**b**) *y*+ = 100; (**c**) *y*+ = 150; (**d**) *y*+ = 200; (**e**) *y*+ = 250; (**f**) *y*+ = 300.

Figure 12 indicates that the setting of boundary layer under different *y*+ values was reasonable, and the amount waterline of the vessel arose changed due to the increase of the trim angle; the reduction of the waterline length made the *y*+ value less than the theoretical value.

Table 5 and Figure 13 indicate that the resistance curves of different *y*+ values had the same convergence. However, a smaller *y*+ value made the calculation deviation larger due to the sharp reduction of the waterline length. When *y*+ = 250, the deviation was smallest, so subsequently, the *y*+ value of 250 was adopted in this research.


**Table 5.** Resistance comparisons of the test and numerical results at different *y+* values.

**Figure 13.** Resistance curves of the different *y*+.

#### *4.3. Time Step of Iteration*

The determination of time step (Δ*ts*) is usually based on the incoming flow characteristics, and it must satisfy not less than the distance of the designed minimum grid at the same time. Thus, with the aim of high accuracy and solution efficiency, the influence of the Δ*ts* on the deviation was analyzed in this section.

To ascertain the appropriate value of Δ*ts*, five values of 0.01 s, 0.008 s, 0.006 s, 0.004 s, 0.002 s were used to compute the resistance, and the numerical results of different Δ*ts* values were compared with the test, and are presented in Table 6 and Figure 14.

**Table 6.** The numerical results of resistance at different determination of time steps (Δ*ts*).


Table 6 and Figure 14 shows that the five resistance curves had the same convergence trend, the larger Δ*ts* (0.01 s) had a weaker calculation accuracy, and the deviation became smaller with the shortening of the time step, which indicates a smaller time interval is conducive to improve the solution precision. However, when the Δ*ts* was less than 0.004 s, further reducing the Δ*ts* value greatly increased the convergence duration, which caused a waste of calculation efficiency, so the time step Δ*ts* = 0.004 s was adopted because of its higher accuracy and lesser calculation duration.

Summing up the above, the influences of the grid parameters, *y*+ values, and time step on the calculation accuracy and efficiency were verified. The results show that the numerical method had good convergence, high accuracy, and appropriate efficiency in simulating the navigation of the model in calm water.

**Figure 14.** Resistance curves of the different time steps.

#### *4.4. Comparison of Numerical and Experimental Results*

As in the aforementioned CFD setup, simulations were performed for the model at *F*r = 0.63–2.52 or until the porpoising appears when *L*cg/Lm = 0.38 (condition two) and 0.35. The comparison between experimental data and simulated values for total resistance (*R*T), sinkage Zcg, and trim angle *τ* are shown in Figure 15, and the deviations are listed in Table 7.

**Figure 15.** Comparison between EFD and CFD results for (**a**) total resistance *R*<sup>T</sup> and the *R*<sup>T</sup> added to the whisker spray resistance, (**b**) sinkage and trim angle.

Whether the ratio *L*cg/Lm was 0.38 or 0.35, the changing trend of the CFD results for the *R*T, *Z*cg, and *τ* with speed was consistent with the test, as shown in Figure 15. During the planing stage (*F*r ≥ 1.05), for the *Z*cg and τ when *L*cg/Lm = 0.38, the maximum deviations were 10.06% and −6.6%, respectively; when *L*cg/Lm = 0.35, the maximum deviations were 10.42% and −4.87%, respectively, most deviations of the remaining were below 10%, as shown in Table 6, that indicates the CFD setup could accurately forecast the sailing attitudes (*Z*cg and *τ*) of the model in the planing regime. But for the *Z*cg and *τ* in

the draining state (*F*r < 1.05), the deviations of *Z*cg were relatively larger, the deviations of τ were smaller, controlled at approximately 12% or lower.

In addition, for the computation of *R*<sup>T</sup> in the draining and semi-planing state (*F*r ≤ 1.05), the deviations were controlled below 5%, which implies that the CFD setup had better forecast precision for the *R*<sup>T</sup> at that stage. However, when entering the planing state, as speed increased, the deviations of *R*<sup>T</sup> under the two conditions got larger, especially at *F*r = 2.52 when *L*cg/Lm = 0.38, the maximum deviation attained was 18.15%, which indicates the numerical method could weakly forecast the *R*<sup>T</sup> of the model during the high-speed planing stage (*F*r > 2.1). In this stage, the spray resistance accounted for a large proportion of the total resistance, but the size of the refined mesh was unable to finely capture the tiny splash and jet. Concurrently, the nonlinear factors such as splash and jet were not effectively solved by the turbulence model, which caused the forecast precision to descend.

To obtain more accurate results of the *R*<sup>T</sup> during the high-speed planing stage, based on the trim angle obtained by the CFD method and combining the whisker spray equation reported by Savitsky [27], the spray resistance (*R*S) was computed and further added to the *R*<sup>T</sup> solved by the CFD method, the amendatory *R*<sup>T</sup> are summarized in the brackets of the third column in Table 7 and Figure 15a. Comparing the amendatory *R*<sup>T</sup> with test results, the maximum deviation was −6.53% at *F*r = 1.89 when *L*cg/Lm = 0.35, which proves the calculated results of the *R*<sup>T</sup> could be remarkably improved after the correction of the whisker spray equation of Savitsky [27].

In addition, under conditions two and three, porpoising respectively occurred at *v* = 13 m/s and 12 m/s in the simulation process, which accorded with the observed phenomena during the test. In brief, the combination of the CFD method and the whisker spray equation of Savitsky [27] achieved good forecasting of the total resistance and sailing attitude of the vessel during the high-speed planing stage.


