A Potential Flow Theory and Boundary Layer Theory Based Hybrid Method for Waterjet Propulsion

Abstract: A hybrid method—coupled with the boundary element method (BEM) for wave-making resistance, the empirical method (EM) for viscous resistance, and the boundary layer theory (BLT) for capture of an area’s physical parameters—was proposed to predict waterjet propulsion performance. The waterjet propulsion iteration process was established from the force-balanced waterjet–hull system by applying the hybrid approach. Numerical validation of the present method was carried out using the 1/8.556 scale waterjet-propelled ITTC (International Towing Tank Conference) Athena ship model. Resistance, attitudes, wave cut profiles, waterjet thrust, and thrust deduction showed similar tendencies to the experimental curves and were in good agreement with the data. The application of the present hybrid method to the side-hull configuration research of a trimaran indicates that the side-hull arranged at the rear of the main hull contributed to energy-saving and high-efficiency propulsion. In addition, at high Froude numbers, the “fore-body trimaran” showed a local advantage in resistance and thrust deduction.


Introduction
The prediction of the waterjet thrust and thrust deduction is an important part of the hydrodynamic performance of the waterjet propulsion system.There are mainly three kinds of prediction methods: the theoretical analysis method, the model test method, and the CFD (computational fluid dynamics) method.The first method is based on experimental statistic data.The equilibrium equation composed of the Bernoulli equation and the energy loss equation determined by empirical coefficients was established to solve the equation of thrust.The model test method was proposed by the ITTC (International Towing Tank Conference).Thrust was obtained by direct measurement or by using the momentum flux method with a measured velocity field at both the inlet and the jet.The CFD method based on the momentum flux theory has been widely applied to obtain waterjet propulsion performance in recent years.
The research on waterjet propulsion carried out in the 1990s concentrates on waterjet-hull interaction [1][2][3][4].The standardization of waterjet thrust prediction and measurement, which was proposed by the Specialist Committee on Waterjets of the ITTC [5], greatly contributes to the development of the CFD method in the waterjet propulsion field.The intake duct system was researched based on the viscous CFD method [6][7][8][9][10].This research focuses on the ducting system itself without considering the influence of the ship due to the limitation of computer hardware.Later, the rapid development of the computer provided the basis for a large-scale parallel CFD numerical simulation.
Takai et al. [11,12] completed a multi-objective optimization of a ship with the waterjet propulsion system.The research of Altosole et al. [13] provides a reliable and effective method to optimize the parameters of the waterjet propulsion system for a high-speed vessel.Eslamdoost [14,15] used the potential flow module XPAN and the double model module XBOUND to study the waterjet-hull interaction problem based on the pressure jump method.
The numerical research of hydrodynamics of marine structures mainly adopts the potential flow CFD [16,17] and the viscous flow CFD method [18], which captures more details of the flow field.However, the computational efficiency of viscous flow CFD still needs to be improved, and a rapid thrust prediction method is required in the preliminary design stage for the waterjet ship.In this present study, the potential flow theory and boundary layer theory-based hybrid method is proposed to rapidly and efficiently achieve the calculation of waterjet propulsion performance.The validation process was also carried out.

The Force-Balanced Waterjet-Hull System
The waterjet-hull system consists of two parts: the waterjet system and the rest hull surface.In the waterjet duct, the fluid pressure and velocity of the fluid change is greatly due to the rotation of the pump rotor.In this study, the rotor and stator are assumed to be a virtual disk, and the pressure jump occurs when the fluid passes the virtual disk [14].The longitudinal force of the hull motion resistance is considered as follows (details of the symbols are listed at the end of the article): where N is the number of waterjet sets (waterjets units are assumed to be the same), R H is the total resistance of the hull (excluding the duct), R D,x is the longitudinal force of each duct before the pump, and R N,x is the longitudinal force of each nozzle chamber after the pump.
The pressure before and after the pump (virtual disk) is P b and P a , respectively, and the section area of duct at rotor is A i .There is a balance between the resultant force of the hull and the thrust of the waterjet system while taking the trim angle θ into account: The longitudinal force of the duct in waterjet propulsion can be calculated by using the following equation: where S n is the inner surface of the nozzle chamber.The flow passing through the pump area gets an acceleration in velocity due to the pressure jump.The pressure after the virtual disk can separate to a constant flow pressure, P c , and the pressure difference ∆p: Inside the nozzle chamber: From Equations ( 1)-( 5), we obtained: Equation ( 6) is the sum of forces except the pressure jump at the virtual disk in the whole system: Then, the expression of the pressure jump is obtained using the equation: ∆p = R RT NA n cos θ .

Iterative Solution Model
The thrust of the waterjet propulsion was solved in a virtual control volume based on the momentum flux method [19].A closed control volume consists of the capture area (1a), the duct surface (3), the virtual boundary surface (4), and the nozzle discharge (5), as shown in Figure 1.(5) From Equations ( 1)-( 5), we obtained: (6) Equation ( 6) is the sum of forces except the pressure jump at the virtual disk in the whole system: (7) Then, the expression of the pressure jump is obtained using the equation: .

Iterative Solution Model
The thrust of the waterjet propulsion was solved in a virtual control volume based on the momentum flux method [19].A closed control volume consists of the capture area (1a), the duct surface (3), the virtual boundary surface (4), and the nozzle discharge ( 5), as shown in Figure 1.Considering the local losses caused by frictional retardation, flow direction changes in the duct, duct expansion, and contraction, the pressure head balance applied in the center of the capture area and the nozzle discharge, there is: The nozzle velocity ratio was obtained using the following equation: ( -)cos Considering the local losses caused by frictional retardation, flow direction changes in the duct, duct expansion, and contraction, the pressure head balance applied in the center of the capture area and the nozzle discharge, there is: The nozzle velocity ratio was obtained using the following equation: where u 0 is the advance velocity of the ship, IVR = u 1a /u 0 is the inlet velocity ratio, u 1a and u 5 , p 1a and p 5 , and h 1a and h 5 represent the average velocity, average pressure, and average height above the baseline of the capture area and nozzle respectively.K is the local loss coefficient of the duct, including friction loss, flow direction loss, and duct contraction loss.Xia provides more details in his book [20] about the local loss coefficient.NVR was set as the initial input value in the calculation and a new velocity ratio NVR can be obtained by using Equation (9), once the force acting on ship and the pressure jump ∆p are solved.Therefore, the size correlation between the ratio |∆NVR|/NVR and the residual ε we set was the criterion used to determine whether the stop criterion was reached.This concept is illustrated in Figure 2. The waterjet thrust could be computed once the stop criterion was reached.
NVR was set as the initial input value in the calculation and a new velocity ratio can be obtained by using Equation ( 9), once the force acting on ship and the pressure jump are solved.
Therefore, the size correlation between the ratio and the residual we set was the criterion used to determine whether the stop criterion was reached.This concept is illustrated in Figure 2. The waterjet thrust could be computed once the stop criterion was reached.

Wave-Making Resistance
Based on the basic assumption of the potential flow theory, the fluid is in-viscid, in-compressible, and ideal.The total velocity potential without considering the unsteady disturbances is written as: (10) where represents the velocity potential of uniform flow, and is the perturbation velocity potential caused by the uniform velocity motion of the ship, there is: (11) The free surface boundary condition is: The hull boundary conditions should be specially treated.An additional boundary condition for the nozzle discharge was added to the waterjet propulsion calculation.The hull and duct

Wave-Making Resistance
Based on the basic assumption of the potential flow theory, the fluid is in-viscid, in-compressible, and ideal.The total velocity potential without considering the unsteady disturbances is written as: where −u 0 x represents the velocity potential of uniform flow, and ϕ is the perturbation velocity potential caused by the uniform velocity motion of the ship, there is: The free surface boundary condition is: The hull boundary conditions should be specially treated.An additional boundary condition for the nozzle discharge was added to the waterjet propulsion calculation.The hull and duct surface, except the nozzle discharge, were similar to the wall without penetration.The nozzle discharge was set as the velocity boundary with a fixed velocity of NVR • u 0 .The two boundary conditions were defined as follows: The infinite boundary condition and transom boundary condition [21] were written as: where n x is the outward normal component of the ship hull surface in the x direction and NVR is the nozzle velocity ratio.In the transom condition, h T is the height of the transom stern edge.x T and y T are the longitudinal and transverse coordinates of the transom stern, respectively.∆x is a the finite distance on the free surface after the transom stern, which can be evaluated as half the longitudinal scale of the panel.The Rankine source panel method [22][23][24] was used to obtain the perturbation potential in the flow field using the discretization of the above equations.The Rankine source Green function can be expressed as G(p, q) = 1 r(p, q) , where p(x, y, z) is the field point and q(x , y , z ) is the source point.
Then, the potential velocity at p can be denoted as: where r(p, q) is the distance between p and q.The total number of panels (N t ) can be divided into the panels of the hull (N s ), panels of the nozzle discharge (N z ), and panels of the free surface (N f ).This can be represented as N t = N s + N z + N f .S i indicating the number i panels.Applying Equations ( 13)-( 15), we can obtain the discrete hull surface boundary equation: Free surface and transom boundary condition are converted to the new discrete form: N i=1, j i To fulfill the infinity radiation condition and avoid the singularity in the integration of the free surface panels, the source points of the free surface panels should manually rise up for 0.5-2.5 ∆x and the field points should move forward about ∆x.By combining the discrete Equations ( 15)- (18), the source strength could be solved to get the potential velocity.
Then the hull surface pressure distribution was obtained according to the Bernoulli equation: The fluid force in the x directional was obtained by integration of the pressure distribution P, which is the wave-making resistance of the hull: The trim and sinkage of the ship were resolved though the force and moment balance equations as follows [25,26]: x 2 y(x, 0)dx (22) where C p = P/(0.5ρu 2 0 ) is the pressure coefficient of the control point on the hull panel, s represents the sinkage, and y(x, 0) is the width of the waterline at x under floating in an even keel state.The initial value of s and θ was set as 0. The continuous iteration was used to calculate a new s and θ value along with the mesh being updated until the balance equations converged.

Viscous Resistance
The 1957ITTC friction resistance coefficient formula was used to calculated the friction resistance of the ship hull: The Reynolds number was defined using the equation Re = u 0 L wl /ν, where L wl is the waterline length, ν is the kinematic viscosity of water, and the roughness allowance is ∆C f = 0.4 × 10 −3 in the process of converting the model-scale results to the full-scale ship.
According to the approximate formula of the viscous-pressure resistance coefficient [27]: The ship's viscous pressure resistance is: where A m is the cross section area at midship, and L r is the run length, which should meet the condition: L r ≥ 4.08

Physical Quantities of the Capture Area
The velocity and pressure distribution at the capture area were calculated using the boundary layer theory.The velocity distribution formula in the turbulent boundary layer theory is as follows: where y is the vertical distance from the hull, and u is the local velocity.The index n was 7 in the model-scale, and 9 in the full-scale [19].δ is the thickness of the boundary layer, which is calculated using the Weighardt formula, δ = 0.27x(u 0 x/υ) −1/6 , at high Re or the plate flow formula, δ = 0.37x(u 0 x/υ) −1/5 [28].
The experiments used to measure the capture area with the series IVR (inlet velocity ratio) ranged from 0.25 to 0.71 in front of the inlet point of the tangency carried out by Roberts [7].It showed that, for a low IVR of <0.39, the capture area occurs as a concave phenomenon.For higher values (i.e., >0.58), the capture area resembles a semi-elliptical.This agrees with the elliptic capture area phenomenon noticed by Alexander [2].
The semi-elliptical shape was used to represent the transverse section of the capture area: where W is the width of capture area recommended to be 50% larger than the geometrical intake width [19].Y 0 is the effective inflow thickness of the capture area.
There is a basic definition as follows: where, in the equation i = 1, 2, 3. Q 1 , Q 2 , and Q 3 represent the inflow volume, inflow momentum, and inflow kinetic energy, respectively.The correlation between the hull boundary layer thickness δ and Y 0 determined the above parameters and the inlet momentum coefficient.These parameters were deduced in Liu's paper [28].While Y 0 > δ, the capture area was divided into two parts, as shown in Figure 3. Liu took the upper part as a rectangle to do the integration.Under this assumption, as the difference between δ and Y 0 becomes smaller, the numerical error increases.In the present work, the accurate integral formulas were expressed to improve the calculation while Y 0 > δ: When Y 0 ≤ δ: where The average velocity of the capture area can be calculated as .The average pressure of the capture area was determined using the formula . is the average dynamic pressure and is the average static pressure.
Y0 is obtained by solving Equations ( 29) or (30) when NVR or Q1 is given.Then, IVR and the inlet momentum coefficient can be calculated at capture area station 1a, as shown in Figure 1: The relationship between effective thickness Y 0 and the boundary layer.
The average velocity of the capture area can be calculated as u 1a = Q 1 /A 1a .The average pressure of the capture area was determined using the formula p 1a = p d + p s .p d is the average dynamic pressure and p s is the average static pressure.
Y 0 is obtained by solving Equations ( 29) or (30) when NVR or Q 1 is given.Then, IVR and the inlet momentum coefficient c m1a can be calculated at capture area station 1a, as shown in Figure 1:

Athena Model and Barehull Validation
As an early developed standard ship model, the Athena ship with a transom stern and V-shape stem, can be used as a platform to test ship resistance [29], motions [30], wave field [31], and wake flow [32,33] over a large Froude number range both at a model-scale and full-scale.The ITTC waterjet committee recommends the Athena ship as the baseline hull for waterjet propulsion research [19].The principal dimensions are listed in Table 1.The ITTC provided the barehull experimental data, including total resistance, trim angle, and sinkage from seven different respondents for the 1/8.556scale Athena ship model with the inlet covered.The EFD (experimental fluid dynamics) data originally named Set A to Set I by the ITTC.
In the present method, a half hull and free surface were taken into account to decrease the number of panels used in calculating wave-making resistance.The half hull surface was formed from 740 panels (Figure 4).The free surface was a 3.5L pp × 1.0L pp ranged area with 1776 panels distributed which meets the panelization requirements mentioned by Janson [34] and Eslamdoost [14].

Athena Model and Barehull Validation
As an early developed standard ship model, the Athena ship with a transom stern and V-shape stem, can be used as a platform to test ship resistance [29], motions [30], wave field [31], and wake flow [32,33] over a large Froude number range both at a model-scale and full-scale.The ITTC waterjet committee recommends the Athena ship as the baseline hull for waterjet propulsion research [19].The principal dimensions are listed in Table 1.In the present method, a half hull and free surface were taken into account to decrease the number of panels used in calculating wave-making resistance.The half hull surface was formed from 740 panels (Figure 4).The free surface was a 3.5Lpp × 1.0Lpp ranged area with 1776 panels distributed which meets the panelization requirements mentioned by Janson [34] and Eslamdoost [14].For barehull resistance, this method showed a good agreement with the EFD data at Fr = 0.3-0.5 in Figure 5.While Fr > 0.5, the computed resistance scatters were distributed near the Set C, D, H, and I curves.The mean error of the other three sets of EFD data was about 7% at Fr = 0.5-0.7.
For Fr = 0.3-0.6, the calculated trim angle (Figure 6) had almost the same value as Set H, but was larger than the other data sets.The trim distribution trend line was fit with the mean EFD data.The computed sinkage reached the lowest point of −15 mm at Fr = 0.45 and matched the trend line of Set E, H, and I, as shown in Figure 7.
The noticeable distinction between the EFD data from the seven institutes mainly resulted from the non-uniform test conditions and differences in data measurements.The model differences of weight and temperature ranged from −2.6% to +3.9% and −20.7% to +14.7% compared with the mean values, respectively.Temperature determines the water density and viscosity, while weight determines the ship shape immersed in water.Furthermore, towing tanks of different sizes have been used by each institute.Different blocking effects significantly contributed to the observed differences.The joint effect of different factors led to a scatter in the experimental data.Despite the difference among the EFD data sets, their good agreement with the majority of the EFD data validates the present method for calculating ship resistance with free trim and sinkage.
differences.The joint effect of different factors led to a scatter in the experimental data.Despite the difference among the EFD data sets, their good agreement with the majority of the EFD data validates the present method for calculating ship resistance with free trim and sinkage.differences.The joint effect of different factors led to a scatter in the experimental data.Despite the difference among the EFD data sets, their good agreement with the majority of the EFD data validates the present method for calculating ship resistance with free trim and sinkage.differences.The joint effect of different factors led to a scatter in the experimental data.Despite the difference among the EFD data sets, their good agreement with the majority of the EFD data validates the present method for calculating ship resistance with free trim and sinkage.

Validation for Capture Area Parameters
The validation of the capture area parameters was carried out for two real ships mentioned in Jin's [35], Zhang's [36], and Liu's papers [28].The data for this calculation is listed in Table 2.The results in Table 3 show that the momentum coefficient c m1a agrees well with the other methods for both two ships.For Ship 1, Y 0 and δ had a value closer to Liu's results.For Ship 2, Y 0 was very similar to Liu's results, and δ was closer to Zhang's.The differences in the values of δ that occurred between the three methods was mainly due to the variety of dynamic viscosity and δ formulas.The validation results indicate that the present method gave a reasonable and acceptable calculation of parameters for the capture area.

Validation of Self-Propulsion for the Athena Ship
The ITTC waterjet self-propulsion experimental test for the waterjet Athena model was carried out at Fr = 0.6.The measured data included waterjet thrust, rope force, trim, sinkage, IVR, NVR, thrust deduction, etc.Among the data, only the reported scatter points were displayed in the following comparison figures.The numerical results from Eslamdoost [14] were taken as contrast for the present hybrid method in the self-propulsion validation procedure.
For the waterjet Athena model in the experiment, water flowed into the duct, passed through the transom, and then sprayed out from the nozzle.This process in the wave-making resistance calculation could not be achieved by using the boundary element method, which is based on the potential theory.Panel penetration between the free surface and nozzle at the stern results in a numerical error in the computation.To avoid this situation and proceed with the calculation, the Athena waterjet model was modified by cutting the outside part of the nozzle, as Figure 8 shows.Figure 9 displays the mesh of the free surface and the Athena hull with an intake duct.

Validation for Capture Area Parameters
The validation of the capture area parameters was carried out for two real ships mentioned in Jin's [35], Zhang's [36], and Liu's papers [28].The data for this calculation is listed in Table 2.The results in Table 3 show that the momentum coefficient agrees well with the other methods for both two ships.For Ship 1, Y0 and had a value closer to Liu's results.For Ship 2, Y0 was very similar to Liu's results, and was closer to Zhang's.The differences in the values of that occurred between the three methods was mainly due to the variety of dynamic viscosity and formulas.The validation results indicate that the present method gave a reasonable and acceptable calculation of parameters for the capture area.

Validation of Self-Propulsion for the Athena Ship
The ITTC waterjet self-propulsion experimental test for the waterjet Athena model was carried out at Fr = 0.6.The measured data included waterjet thrust, rope force, trim, sinkage, IVR, NVR, thrust deduction, etc.Among the data, only the reported scatter points were displayed in the following comparison figures.The numerical results from Eslamdoost [14] were taken as contrast for the present hybrid method in the self-propulsion validation procedure.
For the waterjet Athena model in the experiment, water flowed into the duct, passed through the transom, and then sprayed out from the nozzle.This process in the wave-making resistance calculation could not be achieved by using the boundary element method, which is based on the potential theory.Panel penetration between the free surface and nozzle at the stern results in a numerical error in the computation.To avoid this situation and proceed with the calculation, the Athena waterjet model was modified by cutting the outside part of the nozzle, as Figure 8 shows.
Figure 9 displays the mesh of the free surface and the Athena hull with an intake duct.Data at Fr = 0.6 show that the computed trim was closer to the Set H data point, with a 0.26% error.
The sinkage was closer to Set A and was 1.7 mm shallower.For the self-propulsion hull, the hydrodynamic force acting on the duct generated a forward moment which reduced the trim.
Meanwhile, the low pressure distribution on the hull surface around the intake duct resulted in a deeper immersion.The waterjet self-propulsion was computed over the Fr range of 0.3-0.7 based on the present hybrid method.The trim and sinkage of the barehull and self-propulsion were compared in Figure 10.The trim under two conditions stayed almost the same at Fr < 0.4.However, when Fr > 0.4, the trim difference increased as the Fr increased, and the maximum gap reached 0.135 • at Fr = 0.7.The two sinkage curves deviated from Fr = 0.4.The observed maximum difference was at Fr = 0.625.Data at Fr = 0.6 show that the computed trim was closer to the Set H data point, with a 0.26% error.The sinkage was closer to Set A and was 1.7 mm shallower.For the self-propulsion hull, the hydrodynamic force acting on the duct generated a forward moment which reduced the trim.Meanwhile, the low pressure distribution on the hull surface around the intake duct resulted in a deeper immersion.Data at Fr = 0.6 show that the computed trim was closer to the Set H data point, with a 0.26% error.
The sinkage was closer to Set A and was 1.7 mm shallower.For the self-propulsion hull, the hydrodynamic force acting on the duct generated a forward moment which reduced the trim.
Meanwhile, the low pressure distribution on the hull surface around the intake duct resulted in a deeper immersion.The wave cut profiles of the barehull and waterjet hull were generated in y/L pp = 0.106, 0.246 at Fr values of 0.424 and 0.608.The experimental data was referenced from Fu's model test report for the 1/8.25 model-scale Athena ship [29].The four wave cut profiles showed the same trend with the measured results (The computed data has been converted using the coordinates of the model test), as illustrated in Figure 11.The peak and trough points matched well with the EFD curves.Even the wave fluctuation was observed after the stern, as shown in Figure 11b.As the nozzle was cut off outside the hull to ensure wave-making resistance computation went smoothly, the wave patterns at the stern showed a similar trend between the barehull and waterjet hull.
IVR and NVR are the important parameters of waterjet thrust in the momentum flux method.A sharp contrast of IVR and NVR values between the present method and the comparison results is illustrated in Figure 12.The reason for this difference is the consideration of the local loss coefficient in the present method, which resulted in the decrease of NVR in Equation ( 8).The present results had a lower value at Fr values of 0.35-0.7.Both the two NVR curves increased steadily and reached their peak at Fr = 0.5.The IVR value based on the boundary layer theory appears to be underestimated compared with the results of Eslamdoost [14], based on the viscous CFD method.The relationship between the δ and Y 0 curves illustrates that the capture area was located inside the boundary layer, which resulted in an underdeveloped velocity distribution in the capture area and the relatively lower position of IVR, as shown in Figure 12.
A good agreement between Q 1 and Y 0 was observed between the prediction and Set H data at Fr = 0.6 in Figure 13.This was despite the fact that the predicted Q 1 curve had an average 4.67% lower than the comparison curve.The relationship between the and Y0 curves illustrates that the capture area was located inside the boundary layer, which resulted in an underdeveloped velocity distribution in the capture area and the relatively lower position of IVR, as shown in Figure 12.
A good agreement between Q1 and Y0 was observed between the prediction and Set H data at Fr = 0.6 in Figure 13.This was despite the fact that the predicted Q1 curve had an average 4.67% lower than the comparison curve.Considering the influence of the boundary layer, in theory, the waterjet thrust calculated by the momentum flux method in the virtual control volume was: where and are the average momentum velocity at the capture area and nozzle discharge, respectively.
Hu and Zangeneh [37] researched the waterjet pump flow using the RANS CFD method.They calculated the momentum velocity of the pump nozzle discharge.was found to be 0.4% higher than the average volumetric velocity, .Regardless of the tiny difference between the two average velocities, was used instead of in Equation (33).The waterjet thrust formula was: While the assumption was equally applied to the capture area, was taken as .
Equation (34) becomes:   Considering the influence of the boundary layer, in theory, the waterjet thrust calculated by the momentum flux method in the virtual control volume was: where and are the average momentum velocity at the capture area and nozzle discharge, respectively.
Hu and Zangeneh [37] researched the waterjet pump flow using the RANS CFD method.They calculated the momentum velocity of the pump nozzle discharge.was found to be 0.4% higher than the average volumetric velocity, .Regardless of the tiny difference between the two average velocities, was used instead of in Equation (33).The waterjet thrust formula was: While the assumption was equally applied to the capture area, was taken as .
Equation (34) becomes: Considering the influence of the boundary layer, in theory, the waterjet thrust calculated by the momentum flux method in the virtual control volume was: where u m1a and u m5 are the average momentum velocity at the capture area and nozzle discharge, respectively.Hu and Zangeneh [37] researched the waterjet pump flow using the RANS CFD method.They calculated the momentum velocity of the pump nozzle discharge.u m5 was found to be 0.4% higher than the average volumetric velocity, u 5 .Regardless of the tiny difference between the two average velocities, u 5 was used instead of u m in Equation (33).The waterjet thrust formula was: While the assumption was equally applied to the capture area, u m1a was taken as u 1a .Equation (34) becomes: Equations ( 34) and ( 35) are the two widely used waterjet thrust formulas.In the present work, c m1a was 1.49% larger than IVR on average, as shown in Figure 12.The difference between the two formulas led to a 1.64% average discrepancy in waterjet thrust, as illustrated in Figure 14.The two waterjet thrust curves were more consistent with the Set B and C data at Fr = 0.6.As Eslamdoost uses a correction coefficient to correct the wave-making resistance during computation, the correction greatly increased the total resistance at high Froude numbers.Meanwhile, no correction was made in the present work as the calculated resistance was within the range of the experimental scatter, which resulted in the gap between the two resistance curves.A highly similar result for rope force was reported within the two curves and average EFD data.
Equations (34) and (35) are the two widely used waterjet thrust formulas.In the present work, was 1.49% larger than IVR on average, as shown in Figure 12.The difference between the two formulas led to a 1.64% average discrepancy in waterjet thrust, as illustrated in Figure 14.The two waterjet thrust curves were more consistent with the Set B and C data at Fr = 0.6.As Eslamdoost uses a correction coefficient to correct the wave-making resistance during computation, the correction greatly increased the total resistance at high Froude numbers.Meanwhile, no correction was made in the present work as the calculated resistance was within the range of the experimental scatter, which resulted in the gap between the two resistance curves.A highly similar result for rope force was reported within the two curves and average EFD data.The thrust deduction fraction prediction formula for the full-scale ship is: (36) Two thrust deduction fraction curves were formed from the two waterjet thrust curves in the present method, which had a 1.65% difference in their averages, as shown in Figure 15.Negative thrust deduction occurred at Fr = 0.3-0.45,as the thrust deduction fraction increased steadily.The peak point was reached at Fr = 0.55.The thrust deduction fraction was underestimated in the present method compared with the average EFD data at Fr = 0.6, but this was within the scope of scatter.
Although some differences in the forces and deduction fraction were reported between the present hybrid method and the contrast data, in general the computation results generally met the experimental data well.The algorithm was validated to be applicable in a wide range of Froude number values, which could be applied to other ship types.The thrust deduction fraction prediction formula for the full-scale ship is: Two thrust deduction fraction curves were formed from the two waterjet thrust curves in the present method, which had a 1.65% difference in their averages, as shown in Figure 15.Negative thrust deduction occurred at Fr = 0.3-0.45,as the thrust deduction fraction increased steadily.The peak point was reached at Fr = 0.55.The thrust deduction fraction was underestimated in the present method compared with the average EFD data at Fr = 0.6, but this was within the scope of scatter.
Although some differences in the forces and deduction fraction were reported between the present hybrid method and the contrast data, in general the computation results generally met the experimental data well.The algorithm was validated to be applicable in a wide range of Froude number values, which could be applied to other ship types.

Application in the Side-Hull Configuration Research for a Trimaran
The target ship was a 120 m length waterjet-propelled trimaran with a service speed of 40 kn (Fr = 0.6).A 1/30 scale ship model was constructed, as shown in Figure 16.The previous research on the side-hull configuration was focused on added mass research [23,24] and resistance based optimization [38,39].As the hybrid method was validated on the Athena mono ship, the waterjet performance based side-hull configuration research for the trimaran model was numerically performed at Fr = 0.3-0.6 using the discussed method.The side-hull layout parameters and main data of the trimaran are provided in Figure 16 and Table 4.

Application in the Side-Hull Configuration Research for a Trimaran
The target ship was a 120 m length waterjet-propelled trimaran with a service speed of 40 kn (Fr = 0.6).A 1/30 scale ship model was constructed, as shown in Figure 16.The previous research on the side-hull configuration was focused on added mass research [23,24] and resistance based optimization [38,39].As the hybrid method was validated on the Athena mono ship, the waterjet performance based side-hull configuration research for the trimaran model was numerically performed at Fr = 0.3-0.6 using the discussed method.The side-hull layout parameters and main data of the trimaran are provided in Figure 16 and Table 4.
Depending on the design, four positions were used in both the longitudinal (a/L pp ) and transverse (b/L pp ) directions, resulting in a total of 16 different side-hull layouts.Both the barehull and waterjet hull were used in the same side-hull configuration to generate a mesh and compute the waterjet performance.

Application in the Side-Hull Configuration Research for a Trimaran
The target ship was a 120 m length waterjet-propelled trimaran with a service speed of 40 kn (Fr = 0.6).A 1/30 scale ship model was constructed, as shown in Figure 16.The previous research on the side-hull configuration was focused on added mass research [23,24] and resistance based optimization [38,39].As the hybrid method was validated on the Athena mono ship, the waterjet performance based side-hull configuration research for the trimaran model was numerically performed at Fr = 0.3-0.6 using the discussed method.The side-hull layout parameters and main data of the trimaran are provided in Figure 16 and Table 4.
Depending on the design, four positions were used in both the longitudinal (a/Lpp) and transverse (b/Lpp) directions, resulting in a total of 16 different side-hull layouts.Both the barehull and waterjet hull were used in the same side-hull configuration to generate a mesh and compute the waterjet performance.Meanwhile, the smaller the lateral position of the side-hull, the closer the side-hull was to the main hull.The closer the side-hull was to the main hull, the better the drag reduction.This low-resistance layout distribution is well understood.When a/Lpp < 0, a part of the side body beyond the main hull  The barehull resistance was computed based on the barehull model.Figure 17 shows the resistance distribution of different side-hull layouts.Evidently, in the range of Fr = 0.3-0.6, the side-hull located at a/L pp = −0.25 had an obvious advantage of low drag in the longitudinal position.Meanwhile, the smaller the lateral position of the side-hull, the closer the side-hull was to the main hull.The closer the side-hull was to the main hull, the better the drag reduction.This low-resistance layout distribution is well understood.When a/L pp < 0, a part of the side body beyond the main hull transom is equivalent to an increase in the "effective ship length", which is beneficial to the reduction of the wave-making resistance.When the body is laterally close to the main hull, the "effective ship width" of the trimaran is reduced and the trimaran becomes more slender, which is also contributes to the wave-making resistance reduction.It is worth noting that, at a high speed of Fr = 0.6, the trimaran with a side-hull located at a/L pp = 0.5 near the bow (i.e., a "fore-body trimaran" [40]), shows the characteristics of reduced resistance.This is shown in Figure 17, which is consistent with the conclusions obtained in the model test in Jia's paper [40].
Waterjet thrust and other parameters were calculated using the waterjet hull models.Then, the thrust deduction was obtained using Equation (36), and the thrust deduction with the Fr distribution curve is shown in Figure 18.
Similar to what was shown regarding resistance distribution, Figure 18a shows that the trimaran with a side-hull located at a/L pp = −0.25,b/L pp = −0.1 had an advantage in thrust deduction, which means this side-hull layout has an efficient propulsion performance.As displayed in Figure 18c,d, when b/L pp > 0.2, the "fore-body trimaran" had the lowest thrust deduction at Fr > 0.35.distribution curve is shown in Figure 18.
Similar to what was shown regarding resistance distribution, Figure 18a

Conclusions
The hybrid method for determining wave-making resistance, viscous resistance, and capture area physical quantities was proposed to calculate waterjet propulsion performance using iterative computation.The validation process for the barehull, capture area, and self-propulsion was applied.The differences of the IVR, NVR, resistance, waterjet thrust, and thrust deduction fraction values between the present method and contrast results were discussed.In the results of the present study, the lower NVR was obtained by considering the local loss in the iteration equation.This shows a 0.2% error compared to the EFD Set B data.As the location of the computed capture area was inside the boundary layer, the lower IVR was obtained, which had a 0.4% difference compared to the EFD Set B data.The error of waterjet thrust was 2.5% compared with the EFD Set B data.Although the differences exist, the present work showed similar tendencies in the comparison curves and reached a good agreement with the experimental data.A comparison between two widely used waterjet thrust formulas was investigated while using the volumetric inlet velocity instead of the momentum inlet velocity.On average, their difference and that of the waterjet thrust and thrust deduction was 1.49%, 1.64%, and 1.65%, respectively.
In addition, the application of the hybrid method in the side-hull configuration research further verified the applicability of the present method.While the side-hull position of the trimaran was located on the rear side of the main hull, its layout was conducive to energy-saving and efficient waterjet propulsion performance.Furthermore, if a large deck area was required and the side-hull was relatively far from the main hull, the layout of the "fore-body trimaran" could be selected to reduce the drag and thrust deduction at a high Froude number.
The hybrid method presented in this paper was preliminarily verified as an effective way to research and predict the characteristics of waterjet propulsion without considering the impeller rotation and cavitation.Research on more waterjet hull types and applications based on the hybrid approach will be undertaken in further works.

Figure 1 .
Figure 1.Definitions of the reference stations of the waterjet vessel.

Figure 1 .
Figure 1.Definitions of the reference stations of the waterjet vessel.

Figure 2 .
Figure 2. Iteration illustration of the waterjet thrust and thrust deduction.

Figure 2 .
Figure 2. Iteration illustration of the waterjet thrust and thrust deduction.

J 22 Figure 3 .
Figure 3.The relationship between effective thickness Y0 and the boundary layer.

Figure 4 .Figure 4 .
Figure 4. Mesh of the free surface and the Athena barehull model.

Figure 5 .
Figure 5. Barehull resistance of the present method and the International Towing Tank Conference (ITTC) experimental fluid dynamics (EFD) data.

Figure 6 .
Figure 6.Trim angle of the present method and the ITTC EFD data.

Figure 5 .
Figure 5. Barehull resistance of the present method and the International Towing Tank Conference (ITTC) experimental fluid dynamics (EFD) data.

Figure 6 .
Figure 6.Trim angle of the present method and the ITTC EFD data.

Figure 5 .
Figure 5. Barehull resistance of the present method and the International Towing Tank Conference (ITTC) experimental fluid dynamics (EFD) data.

Figure 6 .
Figure 6.Trim angle of the present method and the ITTC EFD data.

Figure 7 .
Figure 7. Sinkage of the present method and the ITTC EFD data.

Figure 9 .
Figure 9. Mesh of the free surface and the Athena waterjet model.

Figure 9 .
Figure 9. Mesh of the free surface and the Athena waterjet model.

J
. Mar. Sci.Eng.2019, 7, x FOR PEER REVIEW 13 of 22 had a lower value at Fr values of 0.35-0.7.Both the two NVR curves increased steadily and reached their peak at Fr = 0.5.The IVR value based on the boundary layer theory appears to be underestimated compared with the results of Eslamdoost[14], based on the viscous CFD method.

Figure 14 .
Figure 14.Comparison of the resistance, thrust, and rope force.

Figure 14 .
Figure 14.Comparison of the resistance, thrust, and rope force.

Figure 16 .
Figure 16.The trimaran model and side-hull configuration parameters.

Table 4 .
The main data of the trimaran and side-hull layout.was computed based on the barehull model.Figure17shows the resistance distribution of different side-hull layouts.Evidently, in the range of Fr = 0.3-0.6, the side-hull located at a/Lpp = −0.25 had an obvious advantage of low drag in the longitudinal position.

Figure 16 .
Figure 16.The trimaran model and side-hull configuration parameters.
shows that the trimaran with a side-hull located at a/Lpp = −0.25,b/Lpp = −0.1 had an advantage in thrust deduction, which means this side-hull layout has an efficient propulsion performance.As displayed in Figure 18c,d, when b/Lpp > 0.2, the "fore-body trimaran" had the lowest thrust deduction at Fr > 0.35.

Author Contributions:
Conceptualization, L.Z., J.-N.Z. and Y.-C.S.; Methodology, L.Z., J.-N.Z. and Y.-C.S.; Computation-Validation and Analysis, L.Z., and Y.-C.S.; Original Draft Writing, L.Z.; Writing-Review and Editing, L.Z., J.-N.Z. and Y.-C.S.; Resourses, L.Z. and J.-N.Z.; Supervision, J.-N.Z.; Project Administration, J.-N.Z.; Funding Acquisition, J.-N.Z.Funding: This research was funded by the High-tech Ship Research Project of Ministry of Industry and Information Technology of China; National and International Scientific and Technological Cooperation Special Projects of China, grant number 2013DFA80760, and the Fundamental Research Funds for the Central University of China, grant number 3132016339.Conflicts of Interest:The authors declare no conflict of interest.NotationA 1a Area of the capture surface (m 2 ) a/L pp Logitudinal distance of side hull A i Duct section area at the impeller (m 2 ) b/L pp Transverse distance of the side hull A m Cross section area at midship (m 2 ) c m1a Inlet momentum velocity coefficient A n Area of the nozzle discharge (m 2 ) g Gravity acceleration (m/s 2 ) C p Pressure coefficient of hull panel h T Height of the transom edge (m) C pv Pressure viscous coefficient h 1a Centroid height of the capture area (m) Run length of the ship (m) p 1a Average pressure of Capture area (Pa) N Number of the waterjets units p 5 Average pressure of nozzle discharge (Pa) N f Panel number of free surface p d Average dynamic pressure (Pa)

Table 1 .
Data of the full-scale ship and 1/8.556model-scale ship.

Table 1 .
Data of the full-scale ship and 1/8.556model-scale ship.
The ITTC provided the barehull experimental data, including total resistance, trim angle, and sinkage from seven different respondents for the 1/8.556scale Athena ship model with the inlet covered.The EFD (experimental fluid dynamics) data originally named Set A to Set I by the ITTC.

Table 2 .
The main data of the real ship for validation.

Table 3 .
Comparison of results between different methods.

Table 2 .
The main data of the real ship for validation.

Table 3 .
Comparison of results between different methods.

Table 4 .
The main data of the trimaran and side-hull layout.