A PRACTICAL APPROACH TO THE ASSESSMENT OF WATERJET PROPULSION PERFORMANCE: THE CASE OF A WATERJET-PROPELLED TRIMARAN

To obtain a reasonable evaluation of the performance of waterjet propulsion at the design stage, a semi-theoretical and semi-empirical method is used to calculate the fundamental parameters of waterjet propulsion performance using an iterative approach. To calculate the ship’s resistance, a boundary element method based on three-dimensional potential flow theory is used to solve the wave-making resistance, and an empirical approach is used to evaluate the viscous resistance. Finally, the velocity and pressure of the capture area of the waterjet propulsion control volume are solved based on turbulent boundary layer theory. The iteration equation is established based on the waterjet-hull force-balance equation, and the change in the ship’s attitude and the local loss of the intake duct are considered. The performance parameters of waterjet propulsion, such as resistance, waterjet thrust, thrust deduction, and the physical quantity of the control volume, are solved by iteration. In addition, a PID-controlled free-running ship model is simulated using the RANS CFD method as a comparison. We apply the proposed approach and the RANS CFD method to a waterjetpropelled trimaran model, and the simulation process and the results are presented and discussed. Although there are some differences between the two methods in terms of the local pressure distribution and thrust deduction, the relative error in the evaluation results for the waterjet propulsion performance is generally reasonable and acceptable. This indicates that the present method can be used at the early stages of ship design without partial information about the waterjet propulsion system, and especially in the absence of a physical model of the pump.


INTRODUCTION
A waterjet-propelled vessel has characteristics such as high efficiency, strong anti-cavitation performance, low noise, and outstanding maneuverability, and waterjet propulsion performance is a significant component of the hydrodynamic performance of the waterjet system. Before the 1990s, waterjet propulsion was mainly adopted in empirical methods and model tests [1][2]. The 21st ITTC Committee on Waterjets put forward a standard procedure for performance predictions of waterjet-propelled craft, and the momentum flux method and direct thrust measurements for the evaluation of model self-propulsion (SP) tests were recommended [3]. In direct thrust measurements, the impact of the force on the intake duct is difficult to measure, so the momentum flux method is widely used, as it focuses on measurement of the mass flow in the discharge nozzle, which is more applicable. Meanwhile, based on the above momentum flux method for waterjet propulsion performance, the CFD method is widely used to simulate the detailed flow field and to compute the force components of the system.
Park performed simulations on the intake duct with and without a rotor and stator [4][5], and fully solved the flow details of the intake duct; the results for the simulated momentum flux matched the experimental data to within a relative error of 3.1%. In addition, the streamlines and pressure distribution of the duct were closely captured. The innovative development of overset mesh technology and high-performance computing enabled multi-core parallel CFD computation to be applied to the complex simulation of marine structure interactions. Takai [6][7] used an in-house viscous CFD solver called CFDSHIP-IOWA to perform barehull and self-propulsion simulations for the JHSS WJ model, generating 15 overset blocks to carefully refine the hull surface and near-field, and applied a simplified body-force model to the waterjet pump region to replace the real pump geometry and save on computational cost. The simulation results closely matched the EFD results, and demonstrated the ability to optimize the intake duct shape. Research by Altosole [8] proposed a reliable and effective method of obtaining a preliminary evaluation of the selection of waterjet units and their performance with sufficient accuracy. This method was based on the collection of extensive information about the performance of existing waterjets, and provided a nondimensional performance diagram, which was beneficial for the prediction of waterjet performance at the design stage. Based on the pressure jump method, Eslamdoost [9][10] used a potential flow calculation module called XPAN and the double model calculation module XBOUND of the SHIPFLOW software application to study the waterjet-hull interaction problem, and validated the applicability of this method in predicting waterjet thrust and thrust deduction. In recent research, Gong [11] used both experimental and viscous CFD simulations to research the waterjet-hull interaction for a four-waterjet propelled vessel. The differences in the energy/momentum velocity coefficients of the inner and outer waterjets were investigated, and the main causes were found to be the shape of the inlet duct and the presence of stabilizer fins.
The empirical method has an advantage in terms of efficiency, but provides only estimated data on waterjet propulsion performance, and further model tests or CFD based simulations are needed to obtain more details. Viscous CFD tools can simulate the flow field in more detail, but the viscous CFD method is suitable only in cases where there is sufficient physical information about the target, and the efficiency of this technique still needs to be improved. In this paper, semi-theoretical and semi-empirical iterative approaches are used to quickly and effectively predict the hydrodynamic performance of a waterjet-propelled vessel and to provide details of the wave patterns and hull pressure distribution in the partial absence of information, and especially in the absence of a waterjet pump. In the present study, simulations based on the current method and the RANS CFD method are carried out for a waterjet-propelled trimaran to evaluate the performance of waterjet propulsion, and the differences between the two approaches are analyzed and discussed.

DEFINITIONS USED IN WATERJET PROPULSION
The global coordinate system satisfies the right-hand rule, and its origin O is located at the transom on the free surface. The X-axis points to the bow, along the longitudinal section of the main hull, while the Y-axis points to the port side and the Z-axis is perpendicular to the XOY plane, pointing upwards. The inflow speed is U (Ux, Uy), where Ux points to the negative direction of the X-axis, and Uy points to the positive direction of the Y-axis.
As shown in Fig. 1, the waterjet-hull system includes two parts: the waterjet system, which contains the virtual control volume and the waterjet pump; and the hull system, which contains the remaining hull surface. According to the widely used definitions of van Terwisga [1], the gross thrust (T g ) of the waterjet system is calculated based on the flux momentum method in the control volume, which is equal to the horizontal component of the moment flux change. Moreover, the net thrust (T net ) is defined as the force vector acting on the physical material boundaries of the waterjet system, and is directly passed to the hull.
where ρ is the fluid density; u x is the fluid velocity vector in the x component; u k and n k are the fluid velocity vector and unit normal vector, respectively; k=1, 2, 3 are the surface stresses in the x direction; and F px is the unit pump force in the x direction. Q, u 0 , NVR, and α are the mass flow rate, ship speed, nozzle velocity ratio, and momentum velocity coefficient, respectively. A in , A out , A duct , A virtual , and V pump correspond to the surface of the capture area, nozzle discharge, duct, virtual surface, and volume of the pump region, respectively. During the numerical simulation of self-propulsion to predict the self-propulsion (SP) point of the real ship, the difference in the frictional resistance between the model and the real ship needs to be corrected. The correction force F D is the so-called towing force or rope force, and is expressed as: where S w is the area of the wetted hull surface, λ is the scale ratio, and C fm , Cf s , and ΔC f are the frictional coefficients of the ship model and the real ship, and the roughness correction coefficient, respectively. When the barehull (BH) model test and SP test or simulations have been carried out, and the ship resistance and gross thrust are obtained, the thrust deduction fraction of the waterjet propulsion system can be obtained using the following formula:

ITERATIVE SOLUTION MODEL
As shown in Fig.1, the longitudinal force on the hull (with trim and sinkage) is as follows: where N is the number of waterjet units, R H is the resistance of the hull except for the intake duct, and R D,x and R N,x are the longitudinal forces on the intake duct and the nozzle chamber.
In the present research, the region of the waterjet pump is assumed to be a virtual disk, and the pump force is regarded as the body-force that accelerates the fluid. The circumferential rotation of the waterjet rotor is not accounted for in the simplified body-force method, as the complicated rotational process needed to generate waterjet thrust is not the primary research target here. We therefore simplified the process and took the final generated axial pump force into account. The pressure difference in the virtual disk is the reason for the variation in the pressure and the velocity of the fluid passing through the pump area. This pressure difference can be written as [9][10]: where θ is the trim angle, p c is the constant flow pressure without the waterjet pump, τ is the shear stress, and S n is the surface of the nozzle chamber. The Bernoulli equation is applied to the center of the capture area and the nozzle, as follows [12] The ratio between the nozzle velocity and the inlet velocity can be obtained as: where in u and out u , in p and out p , in h and out h represent the average velocities, average pressures and average heights above the baseline of the capture area and above the nozzle, respectively. K is the local loss due to frictional retardation, changes in the flow direction, and contraction of the duct.

Fig. 2. Flow chart for the proposed method and the RANS CFD method for comparison
Eq. (9) gives the iterative solution model for the waterjet propulsion once the parameters can be calculated. As shown in the upper part of Fig. 2, the present method for calculating the waterjet propulsion performance contains three components that need to be solved: wave-making, viscous resistance, and the physical quantities of the capture area. The velocity ratio k is the initial input value for the calculation, and after solving the three parts of the equation, a new velocity ratio k can be obtained from the formula, a process which is iterated until the calculation converges. The waterjet thrust and thrust deduction can then be identified.

NUMERICAL CALCULATION OF THE RESISTANCE AND CAPTURE AREA HULL RESISTANCE
The basic assumptions of potential flow theory are that the fluid is inviscid, incompressible, and ideal. The total velocity potential is expressed as , without considering the unsteady disturbances, x u 0 represents the velocity potential of uniform flow, and ö is the constant disturbance potential caused by the uniform velocity motion of the ship. We take Zhou's settings [13] as references to set the free surface boundary condition, the hull and duct boundary conditions, the infinite boundary condition, and the transom boundary condition. These boundary conditions are defined in Eqs. (10)- (15), where n x is the outward normal component of the hull surface in the x-direction, and k is the velocity ratio. In the transom condition, h is the height of the transom stern edge, and x T and y are the longitudinal and transverse coordinates of the transom stern, respectively. ∆x is a finite distance on the free surface after the transom stern, which can be evaluated as the half longitudinal scale of the panel.
The boundary element method is used to obtain the perturbation potentialϕ in the flow field by discretization of the above equations. The wave-making resistance is then obtained in the x-direction of the hull surface pressure distribution, according to the Bernoulli equation: Using the balance equations of the hull force and moments, the trim and sinkage are iteratively converged: where C p , s, and y(x,0) represent the pressure coefficient, the sinkage, and the width of the waterline at position x.
The hull friction resistance is calculated using the 1957ITTC formula for the friction resistance coefficient, and the final expression is: where the Reynolds number is , L wl is the waterline length, and ν is the kinematic viscosity of water.
The viscous pressure resistance coefficient is based on the approximate formula: 20) and the hull viscous pressure resistance is: where A m is the midship cross-sectional area, and L r is the run length, which should fulfill the condition 4.08

CAPTURE AREA
Based on the turbulent boundary layer theory, the velocity distribution at the wall is as follows: where y is the vertical distance away from the hull surface at the section of the capture area, u is the velocity of the point, and the value of the index n is 7 when the model scale is calculated, and 9 in the calculation of the real ship [14]. δ is the boundary layer thickness, which is calculated by the Weighardt formula at high Reynolds number as we could also use the plate flow boundary layer formula In general, a semi-elliptical shape is used to represent the capture area of the control volume, i.e. 2 2 0 where W is the width of the capture area, and Y 0 is the effective inflow thickness. The relationship between the hull boundary layer thickness and Y 0 determines the boundary layer influence coefficient, the flow velocity, and the pressure distribution. The specific formula used to calculate Y 0 is based on Liu's method [14][15]. The momentum velocity coefficientα is defined as: The average velocity in u of the capture area can be calculated after determining the relationship between δ and Y 0 , and at the scale of the model, the average velocity is determined as follows:

APPLICATION OF THE PROPOSED METHOD TO A TRIMARAN MODEL
The proposed method for numerical calculation of the waterjet propulsion performance is applied to a model of a waterjet-propelled trimaran [12,16] with a transom stern, as shown in Fig. 3. Three sets of waterjet units are installed at the stern. The main dimensions are listed in Table 1.

HULL MODIFICATION AND MESH DISTRIBUTION
Since the physical surface area of the waterjets exceeds the hull transom area in the waterjet-propelled vessel, if the free surface panel and hull panels are directly generated for the wave-making resistance calculation, numerical errors will be obtained in the influence coefficient matrix in the potential flow calculation, due to the penetration between the duct panel and the free surface panel. To avoid this situation, while maintaining the duct model unmodified, the main hull of the trimaran should be lengthened appropriately to ensure that the calculation is correct. As shown in Fig. 4, the stern is extended by two amounts, 1.25%L pp , and 2.5%L pp . Since the 1.25%L pp lengthened model is still partly outside of the transom, the 2.5%L pp lengthened trimaran model is more suitable and is therefore applied in the remaining calculations. The increase in the total resistance due to the lengthened part is corrected, and this is expressed as where S add is the wetted surface area of the lengthened part and C t is the calculated total resistance coefficient. Due to symmetry characteristics, a half mesh is generated for the free surface and the trimaran. The wetted modified hull is divided into 1520 panels, and a local mesh refinement approach is applied. The free surface is divided into 1430 panels (Fig. 5).

THE ITERATIVE CONVERGENCE PROCESS
In the iterative calculation used in the present approach, there are two iterative processes. One is the inner iteration of attitude, based on Eq. (17) and Eq. (18), when calculating the wave resistance, and the other is the outer iteration of the nozzle velocity ratio, based on Eq. (9). Fig. 6 shows the inner iterative convergence process for the trim and sinkage for Fr = 0.30, 0.45, and 0.60 (the trim at the stern is positive). At speeds of Fr = 0.30 and 0.45, the hull motion converges within four to five iterations, and the convergence criterion is that the difference between the last two iterations is less than 1.5%. When Fr=0.6 at high speed, the convergence criterion is satisfied after seven iterations.
It can be seen that as the speed increases, the numerical fluctuations in the calculation process also become larger, and the number of iterations needed gradually increases. In the outer iteration, substantial convergence for the nozzle velocity ratio can be achieved after four or five iterations. When using a 3.6 GHz processor for single-core computing, each iteration takes about 1.5 minutes; it therefore takes about half an hour to calculate a low-speed point, and several hours to calculate a series of speed points. If multi-threaded parallel computing is used, the time cost will be greatly reduced.

PID CONTROLLER AND BODY-FORCE METHOD
For the SP problem, the traditional method in viscous CFD is to fix the ship model in the longitudinal position to encounter the water flow at a specified speed; then, the speed of the propulsor is gradually increased with a fixed step size, and the calculated resistance and thrust curves are finally interpolated to obtain the SP point. In the present RANS CFD comparison simulation, a free-running model is used that can travel freely in real space, rather than using the conventional CFD self-propelled method to encounter water flow in a fixed position. To maintain the required hull speed in self-propulsion, the PID controller [17] is used to automatically adjust the rotating speed of the waterjet rotor. The PID-controlled rotor revolution is as follows: is the error term, 0 u , and M u , and i u are the target hull speed, the hull speed at simulated time t, and the hull speed at the iteration i, respectively. M is the total simulated iteration number, Δt is the time step, and a is the hull acceleration. P, I, and D are the coefficients of the proportion, integration, and differentiation terms, respectively. A widely used propulsion approach, the body-force method uses the concept of an actuator disk rather than a real propeller, and has been shown to be a reasonable approach to the simulation of maneuvering and self-propulsion [18][19]. In the current viscous CFD model, the waterjet rotor is used to simulate the open water performance and obtain convergent KT and KQ curves. These performance curves are applied to the body-force region. As shown in Fig. 8, the original rotor is replaced by a virtual disk of the same radius.

GRID CONVERGENCE ANALYSIS
The computational domain is cuboid, with a size of 4.5L pp ×3.0L pp ×1.5L pp . The prism layer grid is 2.0 cm thick and contains six layers, the growth factor is 1.414, and the y+ value of the hull surface mesh is less than 100. The volume mesh is generated using prism and trimming volume techniques. To carry out the unstructured mesh convergence analysis, Baek [20] and Gong [11] used the total number of grids to define the refinement ratio, and the mesh independence of the RANS CFD was validated using the following equations: where N is the total number of the grids, r represents the refinement ratio, S is the solution to the corresponding mesh density, and F s is the safety factor, which is equal to three.
To carry out a grid convergence analysis for the RANS CFD contrast model, three different mesh densities are applied to simulate the BH model and the SP performance for Fr=0.30. The total resistance for the BH model and the total resistance and RPM for the SP model are validated as shown in Tables 2  and 3, where Extr. is the linear extrapolated predicted value of the finer mesh adopting the same refinement ratio. The final uncertainty is normalized by Extr. The results show that with an increase in the mesh density, the uncertainty in the total resistance drops from 16.72% to 4.56% and from 6.53% to 1.05% for the BH and SP models, respectively, and the uncertainty in the RPM for the SP model drops from 8.93% to 4.65%. The uncertainty results for the fine mesh are reasonable and reliable, and thus the fine mesh density is applied to the comparison models. The final mesh and the computational region used are shown in Fig. 9.

Fig. 9. Computational domain and mesh for the RANS CFD model
Three different speed inputs (Fr=0.30, 0.45 and 0.60) are used in the comparison models. The towing force at different speeds is also added to the hull. Fig. 10 shows the convergence process for the rotor and hull speeds at different Froude numbers. We can see that the increased resistance rapidly reduces the hull speed at all three values of Fr after starting the simulations. Through the PID control, the rotor RPM is automatically adjusted to the maximum speed within a simulation time of 0.2-0.3 s; the RPM is then gradually reduced, and the hull speed is gradually increased. The rotor RPM, hull speed, thrust, resistance, and attitudes all start to fluctuate slightly, and converge steadily within 2-4 s of simulation time. The traveling wave of the ship converges slowly compared to the above physical quantities, and the simulation time needed for stabilization of the wave pattern is at least 10 s.
The PID controller technique used in the rapid response of the ship's speed apparently reduces the fluctuations in the physical quantities, which is beneficial for convergence of the simulation. The body-force method used for the propulsion also contributes to the time saving. Although the above techniques are used in the RANS CFD simulations, for the self-propulsion simulation at the specified speed, the time required for a convincing simulation result is still several tens of hours using parallel computing.  Fig. 11 shows the boundary layer velocity distribution for the hull, the capture area, and the control volume streamlines for the three waterjet propulsion ducts in the RANS CFD simulation at Fr=0.30. The thickness of the boundary layer reaches a maximum in the middle bottom of the hull, while it reaches a minimum in the bilge. The capture area varies with the hull surface curvature, and the shape of the middle duct capture area is a regular semi-ellipse, as it is located at the flat hull bottom, while the shapes of the other two ducts are inclined semi-ellipses as they are located close to the bilge area where the curvature is considerably changed.

PHYSICAL QUANTITIES OF THE CAPTURE AREA
In this paper, the results of the capture area calculation obtained from the turbulent boundary layer method are compared with the results from the RANS CFD method, as shown in Table 4. The width used in the present method is fixed at W=1.2D, where D=0.06 m is the maximum diameter of the intake duct. The RANS CFD results show that the width W and height Y 0 decrease with an increase in Fr. The results of the proposed method give a similar trend for the variation in Y 0 , while the value of Y 0 is 29% larger than the viscous results on average, and the maximum difference in the width W is 8.9%. The velocity distribution of the capture area determines the momentum velocity coefficient á . Different computed values of á are observed at the three capture areas of the left, center, and right intake duct. Since the outer waterjets are affected by the inlet shape of the bilge, the momentum velocity coefficients are both reduced compared with that of the center waterjet, but this difference is reduced with an increase in Fr. The maximum difference between the present approach and the viscous CFD method is 3.9%. In the range Fr=0.30-0.60, the range of á is 0.782-0.795, and this increases gradually with Fr.   Fig. 12 shows a comparison of the wave patterns in the present approach and the RANS CFD method, where the y/ L pp =0.06 wave-cut profile lies between the main body and the side hull of the trimaran, and the y/L pp =0.24 wave-cut profile is located outside of the side hull. In general, the ship traveling wave has an obvious speed effect. The wave peak of the bow increases by nearly a factor of two between Fr=0. 3 and Fr=0.6. The wave crest at the stern increases and moves backwards, and the wave trough formed between the main hull and the side hull also moves backwards and becomes deeper. A comparison of the two methods shows that except for the first wave trough at the rear of the ship with Fr=0.6, the value calculated using the present method is smaller, and the result for y/L pp =0.24 far from the hull shows good agreement with the RANS CFD method, including the small fluctuations in the wave patterns and the amplitude of the peak and troughs.

WAVE PATTERNS
The main reason for the differences in the wave patterns between the two methods is that since the potential flow theory is based on the boundary element method, it cannot simulate the physical phenomenon of the actual water spray from the nozzle discharge at the stern, so the interaction between the flow field and water flow from the waterjets cannot be modeled. This effect is not apparent while the nozzle velocity is small at low speeds (Fr=0.30), but when the nozzle velocity is increased sufficiently at high speeds (Fr=0.60), a more direct effect occurs, and the flow field around the stern is changed. This is also the reason why the two wave-cut profiles show a relatively large difference at the stern and rear of the ship when Fr=0.60, which gradually decreases with the spread of the wave patterns (from y/L pp =0.06 to y/L pp= 0.24).

PRESSURE DISTRIBUTION AT THE INTAKE DUCT
Since a lengthened hull is used in the calculation of wave-making resistance to avoid numerical error, and the viscous CFD contrast model uses the actual model without modification, the pressure distribution on the intake duct surface between the above two methods will show some differences, as illustrated in Fig. 13. The pressure coefficient near the nozzle is negative (C p =−0.5 to −1.0) in the present method, while the pressure coefficient in viscous CFD is close to zero. The above phenomenon is mainly due to the adoption of the velocity exit boundary condition in the calculation of potential flow, which implies a negative pressure distribution at the nozzle and means that the peripheral surface element pressure distribution is also negative. Another area with an apparent pressure difference is the stator region in viscous CFD, where the stator model is retained due to the need to eliminate the rotational energy of the accelerated fluid in the actual waterjet propulsion system. When the fluid passes through the stator after rotational acceleration, the circumferential rotating energy of the fluid is eliminated by the rectification of the stator, leaving only the axial velocity to eject. Hence, in the stator region, due to the impact of the rotating fluid, a sizeable local pressure coefficient (up to a maximum of 2.5) will occur. At the same time, a simplified body-force method is adopted in the pump area (compared with the body-force method used in RANS CFD, the simplified method neglects the circumferential rotation of the rotor and only considers the axial force). The stator is not set, so the pressure coefficient distribution in the stator area is lower, but the discrepancy between them is relatively small in the pump area.  Fig. 14 shows the pressure coefficient curves on the longitudinal section, in the center plane of the intake duct and the nearby hull. The difference in the pressure distribution between the two methods can be intuitively observed. For the waterjet propulsion system, the value of the pressure coefficient in the viscous CFD method is close to zero, and the value for the present method is near −0.9, while for the stator region, the value in the present method is between −0.9 and 0.3, and that in the viscous CFD method ranges from −0.7 to 2.5; for the oblique ascending part of the duct, the values of pressure coefficients in both methods are relatively small, and are between −0.3 and 0.4. For the hull system area, the pressure distribution on the hull surface in front of the intake duct is kept the same (Fore_hull area). The part of the hull that is below the lower surface of the duct has a longer distribution curve, due to the 2.5% lengthening of the hull in the proposed method, and value from the viscous CFD method is slightly higher in the Aft_hull region, where the pressure coefficients tend to be close to each other. In general, the pressure distribution on the intake duct of the waterjet propulsion system in the two methods is different due to the different physical models used, while the pressure distribution on the hull surface of the hull system tends to be the same.

WATERJET THRUST AND THRUST DEDUCTION
The distribution of the gross waterjet thrust and the hull resistance of the trimaran are shown in Fig. 15. The difference between the barehull resistance in the present method and the RANS CFD method is less than 3%, and the difference percentage in the gross thrust is less than 5%. The difference in the maximum hull resistance between the present approach and the viscous CFD method is 4.2% at Fr=0.3 for the waterjet hull. The hull resistance and waterjet thrust obtained via the proposed approach are in good agreement with the RANS CFD results.
The variation curve for the nozzle velocity ratio with Fr, as calculated by the iterative method, is shown in Fig. 16. For Fr=0.30-0.60, the trend in k is downward, while the calculated values are distributed around the value of 1.5. Compared with the results of the RANS CFD method, the difference in k between the two methods is less than 3%.
The total thrust deduction fraction (shown in Fig. 16) can be calculated using Eq. (4). Both methods show a positive thrust deduction within a range of Froude numbers of 0.30-0.60, and the general trend in t is downward with an increase in the Froude number, except for the peak value (marked in red in Fig. 16) at Fr=0.375 for the proposed method. We can find the reason for this inflection point in Eq. (4), as the fraction of total thrust deduction is determined by the difference between T g and the R BH −F D term. If the difference is negative, a negative thrust deduction is obtained, and contrary positive thrust deduction is obtained for a positive difference. For the present calculated positive value of t, the thrust deduction curve follows the same trend as the difference term variation curve. The difference term reaches a maximum at Fr=0.375, as shown in Fig. 15, and thus the maximum thrust deduction fraction is obtained as shown in Fig. 16.

CONCLUSIONS
The applicability of the present approach to the prediction and evaluation of the performance of waterjet propulsion was confirmed by comparing simulation results with the widely used viscous CFD method. For both of these methods, the mathematical and physical models used, the convergence process, and the results were described and discussed in detail, and the following conclusions can be drawn: 1) Although there were some differences in the local pressure distribution and the results of thrust deduction between the proposed method and the CFD method, based on a comparison between the capture area quantities, wave patterns, hull surface pressure distribution, resistance, thrust and thrust deduction, it can be seen that the calculated results for waterjet propulsion are in good agreement with the results of the converged RANS CFD method. 2) The time-saving advantage of the present approach to the assessment of waterjet propulsion performance compared with the viscous CFD method is very large, and is generally a factor of a few tens in terms of efficiency. 3) In the proposed numerical approach, the boundary element method used in the calculation of wave-making resistance can provide more details of the wave patterns and hull pressure distribution, and offers significant time savings. The simplified body-force model is used in the area of the waterjet pump to extend the applicability of the present method, especially in the absence of details about the waterjet pump; an empirical approach is used to calculate the viscous resistance; and the boundary layer theory method is used to calculate the physical quantities of the capture area. This saves on resource costs while giving reasonable accuracy. A comparison indicates that the present research provides a practical approach to the assessment of waterjet propulsion performance.