CFD Body Force Propeller Model with Blade Rotational Effect

: The purpose of this study is to consider propeller geometry and blade rotation in the propeller model in a CFD code. To predict propeller performance, a body force propeller model was developed based on blade element theory and coupled with URANS (unsteady Reynolds-averaged Navier–Stokes) solver CFDSHIP-IOWA V4.5 both implicitly and interactively. The model was executed inside the ﬂow solver every inner iteration. The grid points inside each 2D blade geometry were identiﬁed by a numerical search algorithm. To calculate the lift coefﬁcient, the total ﬂow velocities at 25% foil chord length were obtained using the inverse distance weighting interpolation from the RANS solution. The body forces were distributed linearly along the chord length with the maximal value located at the leading edge and zero at the trailing edge. The main achievements are: (1) for a KP505 propeller in an open water condition, the error of the thrust coefﬁcient generally is around or less than 3%, which is a better prediction than the previous model. (2) For a behind-hull condition, the error is about 1%. (3) For an E1619 propeller in an open water condition, the error is around 6%. (4) The blade-to-blade effect and unsteady ﬂow ﬁeld between blades are sufﬁciently resolved by the model.


Introduction
Over the decades, CFD (computational fluid dynamic) has been widely applied to solve ship hydrodynamics problems. In this paper, CFD is also referred as the RANS (Reynolds-averaged Navier-Stokes) solver or URANS (unsteady RANS). CFD not only predicts the global values but also provides flow field details. The capability of the stateof-the-art CFD to consider discretized real propeller geometry behind a ship hull was demonstrated in the latest CFD workshop on ship hydrodynamics [1] (abbreviated as T2015). Three types of ship hull forms were suggested, with their own propellers and rudders: JBC (Japan Bulk Carrier), KCS, and ONRT (ONR tumblehome).
In some CFD studies, the dynamic overset grid method provides the domain connectivity information among overlapping grid blocks of ship hull and propeller. In T2015, it was implemented in OpenFOAM by Shen and Korpus [2] and the naoe-FOAM-SJTU group [3][4][5][6][7] for unstructured grids. For structured grids, Mofidi et al. [8] performed selfpropulsion and course keeping for ONRT with twin propellers using overset grid in REX.
Another method is to exchange flow field data through an interface between the rotating and fixed part of grid, the so-called sliding mesh approach. Many CFD codes equip this function including commercial, academic, and open-source code, such as the following works presented at T2015: Liu et al. [9] and Qiu et al. [10] for Fluent; Bugalski et al. [11], Park and Jun [12], Kim and Jun [13], and Tenzer et al. [14] for Star-CCM+; Wu et al. [15] and d'Aure et al. [16] for FineMarine; Schuiling et al. [17] for ReFRESCO; Deng et al. [18] for ISIS-CFD; Arai et al. [19] and Abbas et al. [20] for OpenFOAM.
Regardless of whether the overset grid or sliding mesh approach is used, CFD simulations considering realistic and rotating propellers are still expensive and time-consuming.
Meanwhile, several body force propeller models have been developed to save computational time and grid generation effort. Generally, propeller programs are based on the small disturbance method and potential flow theory. Thus, a numerical procedure or coupling method is required to subtract the propeller-induced velocity or obtain the effective wake from a viscous flow solution. Force distribution is predicted by propeller program and then transferred to RANS as body force, which is the source term of momentum equations. Normally, a propeller/actuator disk or a specific propeller plane is utilized to exchange those velocities and body forces between RANS and propeller code.
A body force propeller model can be a prescribed model as well, e.g., a Hough and Ordway circulation distribution [21]. Based on known propeller open water data, a given thrust and torque can be transformed to a body force distribution in CFD grids. In T2015, the model was implemented by Broglia et al. [22] in CNR-INSEAN code, Sadat-Hosseini et al. [23] for CFDSHIP-IOWA, and Lidtke et al. [24] in OpenFOAM and Star-CCM+.
In fact, body force propeller models coupled with RANS code have been under development, have a long history of use, and have been described in many publications. Until now, it has been very useful, effective, and efficient for ship applications. Stern et al. [25] coupled PUF-2 propeller code with a RANS code based on the finite analytical method. Simonsen and Stern [26] introduced a linear-induced velocity subtraction between the Yamazaki propeller model [27] and CFDSHIIP-IOWA. The propeller model was based on the lifting line method. PUF-14 code was coupled with CFDSHIP-IOWA later, as well in Chase et al. [28]. Lee and Chen [29] implemented MPUF3A propeller code inside a Chimera RANS solver to study propeller cavitation in propeller-hull interaction. Viscous nominal wake was obtained first, and then, the wake was updated by solving body forces in RANS. In the end, an effective wake for MPUF3A would be obtained. In a similar manner, Chou et al. [30] used MIT-PSF-2 propeller code in a RANS code UVW. The propeller inflow for MIT-PSF-2 is circumferentially averaged from the RANS solution. Using the circumferentially averaged inflow and updating the nominal wake to an effective one, Chen et al. [31] coupled an unsteady propeller panel code developed by Hsin [32] with Star-CCM+. Widen et al. [33] applied the blade element momentum theory in OpenFOAM to model a propeller behind a ship. Basically, those models were based on the quasi-steady assumption that the propeller rotates much faster than the ship's advance speed. Those models only could predict time-averaged propeller flow field.
More advanced propeller models have also been proposed to distribute body force on 3D propeller geometry to consider the blade-to-blade effect or blade rotation. Kerwin et al. [34] combined a RANS solver and a vortex lattice method. The blade geometry was described by the B-spline surface for easy manipulation in design purpose. The RANS cells inside the 3D blade shape were first identified. The force per volume solved by the vortex lattice method was imposed at the cell center. The propeller panel code panMARE was utilized by Greve et al. [35] in RANS code FreSco+. A plane was located at an upstream propeller to transfer the flow velocities along the propeller's angular and radial directions. A cell search algorithm was proposed to find the RANS cell center closest to the panel for panMARE. Later, the free surface effect was considered [36]. The singularities were distributed on the free surface geometry extracted by Fourier transform from RANS simulations with the propeller model. In the presented work, the cell was only searched in the 2D propeller plane. The side/lateral geometry of the propeller blade was neglected.
In the above-mentioned propeller models, generally, the strength of the free and bound vortex was solved inside the propeller program. The bound vortex was distributed on the propeller surface or along the foil sections. The free vortex was shedding from the trailing edge of the propeller or foil to the far downstream. Therefore, not only the propellerinduced velocity subtraction was required, but the direction of the free vortex was also assumed. However, in CFD simulations and real situations, especially when the propeller operates behind the ship hull, the free vortex is not necessarily guaranteed to be the same as assumed in the propeller program. It is much more complicated and involves the propeller tip vortex and the influence from hull-propeller-rudder interaction. To resolve this issue, a propeller model based on a quasi-steady assumption and blade element theory (BET) had been proposed [37] and has had some applications [38][39][40][41][42]. In the following discussion, it is called the "average model", i.e., the body forces are distributed in average on the circular plane or actuator disk. The total velocity from the CFD flow solution is directly used as propeller inflow, so that the propeller-induced velocity subtraction is no longer needed. As a result, only the free vortex is solved in the propeller model. The bound vortex imbedded inside the CFD flow solution is used in the propeller computation, i.e., the propeller model and the CFD share the same bound vortex.
The purpose of the presented work is to develop the average model [37] further to consider blade rotation and more realistic geometry. It is still based on BET using total velocity and coupled with the RANS solver CFDSHIP-IOWA V4.5 implicitly and interactively. The innovation of the proposed propeller model is summarized as follows. The blade geometry in the 2D propeller plane is constructed by chord length, pitch, and skew distribution. The grid points inside the blade's leading and trailing edges are identified by a numerical search algorithm. The propeller's local thrust and torque are calculated using the lift and drag coefficients. The total velocities at 25% of chord length along foil section are provided via the inverse distance weighting interpolation from a RANS solution. The body force is distributed linearly along the chord length with a peak value at the leading edge and a value of zero at the trailing edge. Several actual applications of the proposed model are conducted successfully. A five-bladed KP505 propeller is simulated in open water for advance coefficient J = 0.3 to 1.0. A grid independence test is conducted for J = 0.5, 0.7, and 0.8. The predicted thrust and torque coefficients show good agreement with experimental values. The unsteady flow phenomena of the blade-to-blade effect and vortex shedding are simulated successfully. The E1619 propeller in open water and KP505 propeller behind KCS (KRISO container ship) hull are investigated for one J.

CFDSHIP-IOWA
The simulations were performed by the URANS (unsteady Reynolds-averaged Navier-Stokes) solver CFDSHIP-IOWA V4.5 [43] to model the viscous flow field around ship and propeller complex geometries, ship motions and unsteady propeller rotational flow. A multi-block structured grid with overset grid capability was used. The solver featured semi-captive, full-6DOF (degree of freedom), parallel, and high-performance computing (HPC) capabilities. The non-dimensional governing equations in the solver were [41]: where t is time, u is the vector form of velocity field, u G is the grid velocity for considering ship motions, and P S is the pressure term. The gravity term is included in P S . S b is the source term, which will represent the body forces produced by propeller model (Sections 2.2 and 2.3) . Re e f f is the effective Reynolds number defined by ship speed U, ship length L, kinematic viscosity of fluid ν, and turbulent viscosity ν t solved by the SST (shear stress transport) k-ω turbulence model without wall function [44]: The SST k-ω model is widely used in ship hydrodynamics application, and its robustness and global reliability had been concluded at T2015 [45] and the previous workshop [46]. Appl. Sci. 2022, 12, 11273 4 of 27 To model the free surface, the governing equation of the single-phase level set method [47] is: where φ is the level set function or distance function. The free surface is defined along the iso-surface of φ = 0, and the undisturbed free surface is located at z/L = 0. φ > 0 indicates that the grid is full of water. In the region of φ < 0, i.e., air, URANS is not solved. The numerical methods to discretize and solve the above-mentioned governing equation are based on the second-order finite-difference method. The Euler backward method is used for the time derivatives, the upwind method is used for the convection term, and the central difference is used for the diffusion term. The velocity and pressure field were coupled by the projection method. The details have been addressed in [41].

Average (Body-Force Propeller) Model
In the average model [37], the flow total velocities (u/U, v/U, and w/U) provided by the URANS solver were directly input into the inflow to calculate the lift and drag on each blade element as shown in Figure 1. The lift coefficient C L was obtained by where k 1 is an empirical correction for finite blade width and α is the AOA (angle of attack) against the zero-lift line. Drag coefficient is assumed to be constant: The local thrust dT and torque dQ, at each grid point inside the propeller disk (cylindrical block) are calculated from the C L and C D on the corresponding blade element. Propeller pitch is considered to calculate α. Chord length is needed when computing dT and dQ from C L and C D . Finally, the axial and tangential body forces (fb x , fb θ ), i.e., force per volume, are obtained as: As implied in Equation (7) of [37] and Equation (8), the model is based on the average effect of dT and dQ among N blades, propeller rotation perimeter 2πr at certain radius r, and thickness of the propeller actuator disk ∆x. Appl. Sci. 2022, 12, x FOR PEER REVIEW 4 of 28 where is the level set function or distance function. The free surface is defined along the iso-surface of = 0, and the undisturbed free surface is located at z/L = 0. 0 indicates that the grid is full of water. In the region of 0, i.e., air, URANS is not solved.
The numerical methods to discretize and solve the above-mentioned governing equation are based on the second-order finite-difference method. The Euler backward method is used for the time derivatives, the upwind method is used for the convection term, and the central difference is used for the diffusion term. The velocity and pressure field were coupled by the projection method. The details have been addressed in [41].

Average (Body-Force Propeller) Model
In the average model [37], the flow total velocities (u/U, v/U, and w/U) provided by the URANS solver were directly input into the inflow to calculate the lift and drag on each blade element as shown in Figure 1. The lift coefficient CL was obtained by where k1 is an empirical correction for finite blade width and α is the AOA (angle of attack) against the zero-lift line. Drag coefficient is assumed to be constant: The local thrust dT and torque dQ, at each grid point inside the propeller disk (cylindrical block) are calculated from the CL and CD on the corresponding blade element. Propeller pitch is considered to calculate α. Chord length is needed when computing dT and dQ from CL and CD. Finally, the axial and tangential body forces (fbx, fbθ), i.e., force per volume, are obtained as: As implied in Equations (7) and (8), the model is based on the average effect of dT and dQ among N blades, propeller rotation perimeter 2πr at certain radius r, and thickness of the propeller actuator disk Δx. The average effect, i.e., the quasi-steady assumption, is useful for ship applications because propellers rotate at high speeds in many scenarios such as self-propulsion in calm water or waves and seakeeping and manoeuvring problems. The momentum equations The average effect, i.e., the quasi-steady assumption, is useful for ship applications because propellers rotate at high speeds in many scenarios such as self-propulsion in calm water or waves and seakeeping and manoeuvring problems. The momentum equations in the URANS solver are solved along with those body forces as the source terms. The six components of total propeller forces and components are considered in the equations of motions in a 6DOF solver. The flow field is solved in an earth-fixed coordinate, 6DOF motions are solved in ship-fixed coordinate, and the propeller model computation is in the propeller-shaft coordinate. Thus, several coordinate transformations are required in the code.

Rotating Blade (Body Force Propeller) Model
The numerical procedure to consider 2D blade geometry and rotation is explained as follows: According to the non-dimensional chord length C/ D P , pitch P/ D P , and skew θ s distribution along the propeller radius, which are generally provided by the propeller offset table, the propeller geometry on an axial plane (x: shaft axis) or a y-z plane can be constructed. For an N-bladed propeller, the mid-chord line would be: herein, n b = 1, 2 . . . , N. Here, k = 1 is to produce mid-chord line. The leading and trailing edge can then be generated from the mid-chord line: At certain time t, the coordinate transformation due to blade rotation is: in which n is propeller rotation rate (RPS, revolution per second). If a 3D model of propeller geometry is available, the leading and trailing edge points (x L , y L , z L ) and (x T , y T , z T ) relative to the propeller center (0,0,0) can be extracted first. Next, the propeller skew (θ s ), pitch P, and chord length C can be calculated by Equations (19)-(21), respectively. Finally, the propeller blade outline can be constructed by following the previous Equations (9)- (15).

Cell Search Algorithm
Searching for and identifying the RANS grids inside the propeller geometry outline is required. In this study, the search algorithm is based on dot product and angle summation. If the propeller geometry outline from Section 2.3.1 is split into several convex cells, the vector between the RANS grid point and each cell edge point can be defined as shown in Figure 2. By calculating the dot product, the angle between a pair of two adjacent vectors is obtained. There are four pairs of vectors corresponding to four angles. If Equation (22) is satisfied, the RANS grid point is inside the cell (propeller outline).
Instead, if Equation (23) is satisfied, the RANS grid point is outside the cell (propeller outline). Appl. Sci. 2022, 12, x FOR PEER REVIEW 6 of 28

Cell Search Algorithm
Searching for and identifying the RANS grids inside the propeller geometry outline is required. In this study, the search algorithm is based on dot product and angle summation. If the propeller geometry outline from Section 2.3.1 is split into several convex cells, the vector between the RANS grid point and each cell edge point can be defined as shown in Figure 2. By calculating the dot product, the angle between a pair of two adjacent vectors is obtained. There are four pairs of vectors corresponding to four angles. If Equation (22) is satisfied, the RANS grid point is inside the cell (propeller outline).
Instead, if Equation (23) is satisfied, the RANS grid point is outside the cell (propeller outline).

Blade Element Theory at ¼ Chord Position
Once the RANS grid point is identified inside one blade, the propeller computation is conducted on the corresponding foil section/blade element at that propeller radius. The total velocities are obtained for the 25% chord position ( , , ), which is the usual location of the aerodynamic center for a foil section.
( , , ) can be generated by Equations (3) and (4) using k' = 0.5. The cell search algorithm in Section 2.3.1 is performed again to locate the RANS grid points surrounding the point ( , , ), such as for points A, B, C, and D as shown in Figure 2. Accordingly, the total velocities ( , , ) at ( , , ) are obtained by the inverse distance weighting interpolation from the velocities on point A, B, C, and D:

Blade Element Theory at 1 /4 Chord Position
Once the RANS grid point is identified inside one blade, the propeller computation is conducted on the corresponding foil section/blade element at that propeller radius. The total velocities are obtained for the 25% chord position x q , y q , z q , which is the usual location of the aerodynamic center for a foil section.
x q , y q , z q can be generated by Equations (3) and (4) using k' = 0.5. The cell search algorithm in Section 2.3.1 is performed again to locate the RANS grid points surrounding the point x q , y q , z q , such as for points A, B, C, and D as shown in Figure 2. Accordingly, the total velocities u q , v q , w q at x q , y q , z q are obtained by the inverse distance weighting interpolation from the velocities on point A, B, C, and D: where r i is the distance between x q , y q , z q and grid point. The grid point A, B, C, and D are assigned to i = 1, 2, 3, 4, respectively. The power parameter p is a positive real number, and is set as 1.51294159474 which is a common value used in CFD-SHIPIOWA V4.5. The main part of propeller computation is the same as in the average model in Section 2.2. The only difference in this case is the body force distribution. In the rotating blade model, the linear body force distribution is specified along the chord length C with the maximal value at leading edge, and a value of zero at the trailing edge. If x q , y q , z q is at x'/C portion of the chord length C, the body forces are distributed axially and tangentially as below: Figure 3 illustrates the linear distribution of the axial and tangential body force.
where is the distance between ( , , ) and grid point. The grid point A, B, C, and D are assigned to i = 1, 2, 3, 4, respectively. The power parameter p is a positive real number, and is set as 1.51294159474 which is a common value used in CFD-SHIPIOWA V4.5.
The main part of propeller computation is the same as in the average model in Section 2.2. The only difference in this case is the body force distribution. In the rotating blade model, the linear body force distribution is specified along the chord length C with the maximal value at leading edge, and a value of zero at the trailing edge. If ( , , ) is at x'/C portion of the chord length C, the body forces are distributed axially and tangentially as below: Figure 3 illustrates the linear distribution of the axial and tangential body force.

Computational Domain and Boundary Conditions
Two simulation conditions are considered in the presented study. The first condition is the simulation of a single propeller in open water conditions. Two different propellers, KP505 and E1619, are simulated in open water using the same CFD setup and grid system. The second simulation condition is the behind hull condition, i.e., the simulation of a KCS model ship hull appending an operating KP505 propeller simplified by the propeller models (Sections 2.2 and 2.3). The computational domain and boundary conditions for the two conditions are explained separately in the following Sections 2.4.1 and 2.4.2.

Open Water Propeller Test
The overset grid system and computational domain are the same as in the previous work [39] for the wave condition with the deepest propeller immersion depth. Here, it is only for the calm water condition, but the free surface is considered. Although only the propeller is simulated here, for the further application of the behind hull condition, the CFD computational domain is non-dimensionalised by the model length of KCS: L = Lpp (length between perpendiculars) = 3.2 m. The domain size is −0.24 < x/L < 1.75 from upstream to downstream, −1 < y/L < 1 from side to side, and −0.79 < z/L < 0.22 from top to bottom. The propeller plane is at x/L = 0.0. The domain consisted of two overset grid

Computational Domain and Boundary Conditions
Two simulation conditions are considered in the presented study. The first condition is the simulation of a single propeller in open water conditions. Two different propellers, KP505 and E1619, are simulated in open water using the same CFD setup and grid system. The second simulation condition is the behind hull condition, i.e., the simulation of a KCS model ship hull appending an operating KP505 propeller simplified by the propeller models (Sections 2.2 and 2.3). The computational domain and boundary conditions for the two conditions are explained separately in the following Sections 2.4.1 and 2.4.2.

Open Water Propeller Test
The overset grid system and computational domain are the same as in the previous work [39] for the wave condition with the deepest propeller immersion depth. Here, it is only for the calm water condition, but the free surface is considered. Although only the propeller is simulated here, for the further application of the behind hull condition, the CFD computational domain is non-dimensionalised by the model length of KCS: L = Lpp (length between perpendiculars) = 3.2 m. The domain size is −0.24 < x/L < 1.75 from upstream to downstream, −1 < y/L < 1 from side to side, and −0.79 < z/L < 0.22 from top to bottom. The propeller plane is at x/L = 0.0. The domain consisted of two overset grid blocks: the propeller cube and far-field background, as shown in Figure 4. Inside the propeller cube, the grid number is 50 × 50 on the x/L = 0 plane to cover the propeller diameter. Table 1 shows the grid number detail. The total grid number is 3.5 M (million). The uniform inflow boundary condition is specified on the inlet, the exit condition for the outlet, the zero-gradient condition on both sides, and the far-field condition on the top and bottom face of the domain. The non-dimensional inflow velocity U is set as 1, and the open water test for different J is conducted by adjusting propeller rotation rate n. The mathematical formulae for the boundary conditions are listed in Table 2 and specified on the boundaries indicated in Figure 4a. There is no solid surface in the open water test here because no physical propeller exits inside the domain. The body force propeller models simplify the propeller. The grid system and CFD setup are the same for the rotating blade and average propeller model, and for KP505 and E1619 propeller.
blocks: the propeller cube and far-field background, as shown in Figure 4. Inside the propeller cube, the grid number is 50 × 50 on the x/L = 0 plane to cover the propeller diameter. Table 1 shows the grid number detail. The total grid number is 3.5 M (million). The uniform inflow boundary condition is specified on the inlet, the exit condition for the outlet, the zero-gradient condition on both sides, and the far-field condition on the top and bottom face of the domain. The non-dimensional inflow velocity U is set as 1, and the open water test for different J is conducted by adjusting propeller rotation rate n. The mathematical formulae for the boundary conditions are listed in Table 2 and specified on the boundaries indicated in Figure 4(a). There is no solid surface in the open water test here because no physical propeller exits inside the domain. The body force propeller models simplify the propeller. The grid system and CFD setup are the same for the rotating blade and average propeller model, and for KP505 and E1619 propeller.

Boundary
Condition

Boundary
Condition

. Propeller behind Hull Condition
The propeller KP505 behind the KCS hull with a rudder in calm water is also tested in the presented work. A propeller disk (cylinder grid block) is required for applying both the average and rotating blade body force propeller models. The same CFD grid system including the propeller disk, and the setup is shared between the two propeller models. Figure 5 explains the computational domain and boundary conditions detailed in Table 2 accordingly. Figure 6 shows the propeller disk is attached on the ship stern bulb with the propeller boss and in front of the rudder. The disk has a central tunnel, which is also the propeller boss surface. The rotating hub velocity, i.e., propeller rotational velocity, is specified on the inner surface of the disk central tunnel. In this section, the ship vertical motions including heave and pitch are considered to predict the sinkage and trim in CFD simulation using dynamic overset grid performed by Suggar. An identical grid system was used in the previous study [41] as well.
propeller boss and in front of the rudder. The disk has a central tunnel, which is also the propeller boss surface. The rotating hub velocity, i.e., propeller rotational velocity, is specified on the inner surface of the disk central tunnel. In this section, the ship vertical motions including heave and pitch are considered to predict the sinkage and trim in CFD simulation using dynamic overset grid performed by Suggar. An identical grid system was used in the previous study [41] as well.
The entire grid system is composed of nine blocks to describe the complicated geometry: starboard and portside for hull, cover (hull top), stern (including shaft tube and dummy boss), and rudder (attached under transom) and one background. The total grid number is around 13 M (million). To capture turbulence, the minimum size of the grid approximately is 1 × 10 −6 normal to the solid surface corresponding to y+ < 1. Fine grids are generated in the x-direction near the ship and the z-direction near the free surface to capture and maintain the propagating waves.

KP505 Propeller in Open Water Test
KP505 is a five-bladed propeller with a skew angle designed for a KCS container ship. The main propeller dimensions are available on the T2015 workshop website [1]. The geometric details are shown in Figure 7 for the technical drawing and Table 3 for the offset table. Because the offset table was available, the 2D geometry outline was extracted using Equations (3)-(9). The entire grid system is composed of nine blocks to describe the complicated geometry: starboard and portside for hull, cover (hull top), stern (including shaft tube and dummy boss), and rudder (attached under transom) and one background. The total grid number is around 13 M (million). To capture turbulence, the minimum size of the grid approximately is 1 × 10 −6 normal to the solid surface corresponding to y+ < 1. Fine grids are generated in the x-direction near the ship and the z-direction near the free surface to capture and maintain the propagating waves.

KP505 Propeller in Open Water Test
KP505 is a five-bladed propeller with a skew angle designed for a KCS container ship. The main propeller dimensions are available on the T2015 workshop website [1]. The geometric details are shown in Figure 7 for the technical drawing and Table 3 for the offset table. Because the offset table was available, the 2D geometry outline was extracted using Equations (3)-(9).

KP505 Propeller in Open Water Test
KP505 is a five-bladed propeller with a skew angle designed for a KCS container ship. The main propeller dimensions are available on the T2015 workshop website [1]. The geometric details are shown in Figure 7 for the technical drawing and Table 3 for the offset table. Because the offset table was available, the 2D geometry outline was extracted using Equations (3)-(9).  (3) and (4). Red and green lines are the leading and trailing edges, respectively, generated by Equations (5) and (6).   (3) and (4). Red and green lines are the leading and trailing edges, respectively, generated by Equations (5) and (6). As shown in Figure 8, the propeller thrust and torque coefficient (KT, KQ) predicted by the proposed rotating blade model were compared with the average model result and experimental data provided by NMRI (National Maritime Research Institute, Tokyo, Japan) for the T2015 workshop [1]. The result of both propeller models shows good agreement with the experiment, especially for the high J (advance coefficient) region. The KT values were from J = 0.3 to 1.0 in the experiment, and the rotating blade and average models are listed in Table 4. In the table, the error E%D for the two models are included as well, and evaluated as follows: wherein D represents the experimental data and S is the simulation value. For the rotating blade model, the KT in low J (J ≤ 0.5) is under-predicted and the KT is over-predicted at high J (J > 0.5). The absolute error |E%D| is less than or around 3%, except for at the two extreme values of J, i.e., 0.3 and 1.0. The smallest error is less than 1%, at J = 0.5. Zero error is expected between J = 0.5 and 0.6. From low to high J excluding J = 1.0, the average model provides under-predicted errors increasing from around 3% to more than 10%. The |E%D| of the average model is generally larger than that of the rotating blade model (except for at J = 0.7 and 1.0). As shown in Figure 8, the propeller thrust and torque coefficient (KT, KQ) predicted by the proposed rotating blade model were compared with the average model result and experimental data provided by NMRI (National Maritime Research Institute, Tokyo, Japan) for the T2015 workshop [1]. The result of both propeller models shows good agreement with the experiment, especially for the high J (advance coefficient) region. The KT values were from J = 0.3 to 1.0 in the experiment, and the rotating blade and average models are listed in Table 4. In the table, the error E%D for the two models are included as well, and evaluated as follows: wherein D represents the experimental data and S is the simulation value. For the rotating blade model, the KT in low J (J 0.5) is under-predicted and the KT is over-predicted at high J (J > 0.5). The absolute error |E%D| is less than or around 3%, except for at the two extreme values of J, i.e., 0.3 and 1.0. The smallest error is less than 1%, at J = 0.5. Zero error is expected between J = 0.5 and 0.6. From low to high J excluding J = 1.0, the average model provides under-predicted errors increasing from around 3% to more than 10%. The |E%D| of the average model is generally larger than that of the rotating blade model (except for at J = 0.7 and 1.0).    For J = 0.1, KT is close to zero, leading to prediction difficulty. Furthermore, it is easy to obtain larger error value because the very small D value is the denominator of Equation (30). In the following discussion for KT, J = 0.1 will not be included.
As J decreases, the KT and KQ of both models predict lower values and deviate more clearly from the experimental data. Especially at low J value, the performance of both models is limited by the blade element theory which cannot consider foil section geometry. Equation (5) may not be suitable for high-angle attack conditions. Because C D is assumed as a small constant, C L is the dominated component which is computed in the current propeller computation. With the consideration of 2D blade geometry, more realistic body force distribution and blade-to-blade interaction, the rotating blade model can provide a larger C L , producing better KT prediction. In contrast, the average model is a highly simplified model, which is merely an actuator disk. In other words, due to lack of the foil geometry consideration, the C L of the average model is generally too low. The KT of the rotating body force is around 7% larger than the average model one, on average.
The KQ values of J = 0.3 to 1.0 of the experiment, rotating blade, and average model are listed in Table 5. In the table, the error E%D for the two models are included as well and evaluated as in Equation (3). Both models show a similar trend and error level. As J decreases, |E%D| becomes worse. The error raises close to 10% when J ≥ 0.5. This issue is caused by the assumption of constant C D = 0.01, i.e., Equation (6). With low values of J corresponding to a high-angle attack condition, the real C D is larger than 0.01, so both models predict lower a lower value KQ. With high values of J, e.g., J > 0.6, C D = 0.01 is very close to a real situation with a small angle of attack. Consequently, both models could provide considerably small |E%D| with under-or over-prediction. However, the rotating blade model does not perform better than the average model for KQ. For only J = 0.7, 0.8, and 0.9, the rotating blade model provides better error, and the |E%D| of J = 0.9 is the lowest amongst all data. In contrast to KT, The KQ of the rotating body force is around 2% smaller (on average) than that of the average model. For propeller torque computation, the foil section drag is more relevant, but C D only can be assumed as a constant in both propeller models currently. Because C D = 0.01 is only applied inside the blade region, the KQ is all smaller for the rotating blade model. Instead, C D = 0.01 was specified on the whole propeller disk area in the average model, so that the KQ is all larger.

Grid Independence Test
Based on the Section 3.1.1 result, the condition J = 0.5, 0.7, and 0.8 are selected to check grid independency of the propeller model and CFD method of the presented work. J = 0.5 has the minimal KT error, which is under-predicted less and around 1% in Table 4. J = 0.7 has the minimal KQ error, which is over-predicted nearly only 0.2% in Table 5. The KT and KQ predictions are performed in balance at J = 0.8. with error about 1.7-1.8%. The averaged absolute error for J = 0.8 is 1.75% which is second lowest among all J values. It is only slightly higher than the error of J = 0.7 (1.57%), because the KQ error of J = 0.7 is really low.
The analysis of grid independence and uncertainty here is based on the verification and validation theory suggested by ITTC 7.5-03-01-01 guideline [48]. A different grid size and density are required. The grid in the previous section and Table 2 is served as medium grid, and then the grid number is increased by (N i , N j , N k ) × √ 2 to become a fine grid. Reducing the grid number, the coarse grid is built by (N i , N j , N k )/ √ 2. Finally, the total grid number is 1.2 M (1,217,280) and 9.8 M (9,777,246) for the coarse and fine grid, respectively. The simulation values of fine, medium, and coarse grid are assigned in this section as S 1 , S 2 , and S 3 , respectively. The grid dependence test results for KT and KQ are presented in Tables 6 and 7, respectively. All S 2 values are the same as the S values in Section 3.1.1. The D value and E%D calculation follows Section 3.1.1 as well. The grid independence is judged by the R G value, which is calculated as below: Grid independence is accomplished for |R G | < 1, meaning that the value difference between the medium and fine grid is smaller than the difference between medium and coarse grid. As shown in Tables 5 and 6 which is presented as percentage of D in Tables 5 and 6. The result achieving a U G > |E%D| of S 1 would be considered as valid. As a consequence, for a KT of J = 0.5 and 0.8 and a KQ of J = 0.7 and 0.8, validation is achieved. Once grid convergence is confirmed, it strongly depends on if S 1 could produce small enough error, i.e., how accurate the simulation can reach as the grid grows finer. For J = 0.8, KT and KQ both are validated. This indicates that the current CFD method and setup is the most suitable for the J = 0.8 condition.
In actuality, for KT of J = 0.7, the U G is only slightly smaller than |E%D| of S 1 with a less than 1% difference. Thus, the comparisons between the CFD results themselves are only slightly less confident than comparisons with the experiment. This is also because the KT is slightly deviated from D somehow from S 2 to S 1 , i.e., when refining the grid from medium to fine grid. The |E%D| increase is very small, also less than 1%. In fact, if comparing U G with |E%D| of S 2 , it can be considered to be validated, i.e., U G > |E%D| of S 2 .
For KQ of J = 0.5, a very small U G value is found, which reflects that S 1 and S 2 are extremely close to each other. Thus, the confidence level for the comparison of the CFD results themselves is much higher than that of comparison with the experiment. As Table 6 reveals, E%D rises to close to 9% as the grid number increases from S 3 to S 2 , and then to S 1 . This is related to the fact that constant C D = 0.01 is unsuitable for low J condition, previously discussed. A larger drag coefficient would increase propeller torque in low J condition because the high angle of attack. A severe increase in error can be observed in Figure 8 and Table 5 as well. Future development of the propeller model, e.g., a more sophisticated C D correction, is required. reach as the grid grows finer. For J = 0.8, KT and KQ both are validated. This indicates that the current CFD method and setup is the most suitable for the J = 0.8 condition.

Body Force Distribution
In actuality, for KT of J = 0.7, the UG is only slightly smaller than |E%D| of S1 with a less than 1% difference. Thus, the comparisons between the CFD results themselves are only slightly less confident than comparisons with the experiment. This is also because the KT is slightly deviated from D somehow from S2 to S1, i.e., when refining the grid from medium to fine grid. The |E%D| increase is very small, also less than 1%. In fact, if comparing UG with |E%D| of S2, it can be considered to be validated, i.e., UG > |E%D| of S2.
For KQ of J = 0.5, a very small UG value is found, which reflects that S1 and S2 are extremely close to each other. Thus, the confidence level for the comparison of the CFD results themselves is much higher than that of comparison with the experiment. As Table  6 reveals, E%D rises to close to 9% as the grid number increases from S3 to S2, and then to S1. This is related to the fact that constant CD = 0.01 is unsuitable for low J condition, previously discussed. A larger drag coefficient would increase propeller torque in low J condition because the high angle of attack. A severe increase in error can be observed in Figure 8 and Table 5 as well. Future development of the propeller model, e.g., a more sophisticated CD correction, is required.  Figure 10 shows the velocity profile behind the propeller for J = 0.8. The flooded contour is the axial velocity, and the vector field includes both lateral and vertical velocities. For the rotating blade model, the blade-to-blade effect to the flow field is clearly observed at x/L = 0.005 (Figure 10a). In the radial direction, the high-velocity region is located around r/R = 0.4 to 0.8. It is not axisymmetric along the circumferential direction. The  Figure 10 shows the velocity profile behind the propeller for J = 0.8. The flooded contour is the axial velocity, and the vector field includes both lateral and vertical velocities. For the rotating blade model, the blade-to-blade effect to the flow field is clearly observed at x/L = 0.005 (Figure 10a). In the radial direction, the high-velocity region is located around r/R = 0.4 to 0.8. It is not axisymmetric along the circumferential direction. The velocity drops slightly between two adjacent blades, i.e., the trailing edge of the current blade and the leading edge of the next blade. The area shrinks as well. The flood contour near the blade leading edge has a clear concave shape. However, at just slightly downstream (e.g., x/L = 0.01 in Figure 10b) the flow field of rotating blade model becomes almost axisymmetric. The contour pattern is almost identical to the average model one in Figure 10c, which is steady, i.e., it would not change by time.

Propeller Wake Field
the same timepoint. It decays more obviously and faster in the far downstream region as well. This phenomenon is more obvious for J = 0.6, and the waves become more laterally scattered downstream in the rotating blade model. This is because the body forces are only distributed on the blades, producing discontinuous influence in the lateral direction on the free surface. Comparing the different J values, the lower J shows a higher axial velocity on the propeller plane, a larger wave amplitude, and a shorter effective pitch indicated by the streamlines.  In Figure 11 (J = 0.8) and Figure 12 (J = 0.6), the free surface waves caused by the propeller excitement have a very small amplitude, because the propeller immersion depth is in deep water. In the comparison of both propeller models, the wave amplitude of the rotating blade model at a specific timepoint is smaller than that of the average model at the same timepoint. It decays more obviously and faster in the far downstream region as well. This phenomenon is more obvious for J = 0.6, and the waves become more laterally scattered downstream in the rotating blade model. This is because the body forces are only distributed on the blades, producing discontinuous influence in the lateral direction on the free surface. Comparing the different J values, the lower J shows a higher axial velocity on the propeller plane, a larger wave amplitude, and a shorter effective pitch indicated by the streamlines. Appl

E1619 Propeller in Open Water Test
In contrast to KP505 and Section 3.1, this section will investigate a propeller with a short chord length, i.e., a small EAR (expanded blade area ratio) and available for a 3D geometry model instead of an offset table. The seven-bladed and highly skewed submarine propeller E1619 was chosen. The main particulars are shown in Table 8. Its 3D model was provided by INSEAN (Italian Ship Model Basin), and the 2D geometry outline was extracted using Equations (10)- (15). The CFD setup and grid system are exactly the same as the one used in Section 3.1. This is also proof of the high flexibility of the presented propeller model; it is convenient and easy to switch between not only two body force propeller models but also different propeller geometries. For J = 0.8, the CFD using the presented blade rotating blade propeller model predicts KT = 0.205. Against the experimental KT = 0.218, the error E%D is around 6% under-predicted. Figure 14 illustrates the axial velocity distribution directly following the propeller plane. To compare the result in Figure 14, please refer the Figure 4.13 on p. 36 in [49]. The area of high velocity behind each blade body is described well for real propeller simulations (CFDSHIP-IOWA using overset grid) [49], LDV experiments [49], and the proposed rotating blade propeller model. The velocity drops between the two blades are captured as well.

E1619 Propeller in Open Water Test
In contrast to KP505 and Section 3.1, this section will investigate a propeller with a short chord length, i.e., a small EAR (expanded blade area ratio) and available for a 3D geometry model instead of an offset table. The seven-bladed and highly skewed submarine propeller E1619 was chosen. The main particulars are shown in Table 8. Its 3D model was provided by INSEAN (Italian Ship Model Basin), and the 2D geometry outline was extracted using Equations (10)- (15). The CFD setup and grid system are exactly the same as the one used in Section 3.1. This is also proof of the high flexibility of the presented propeller model; it is convenient and easy to switch between not only two body force propeller models but also different propeller geometries. For J = 0.8, the CFD using the presented blade rotating blade propeller model predicts KT = 0.205. Against the experimental KT = 0.218, the error E%D is around 6% underpredicted. Figure 14 illustrates the axial velocity distribution directly following the propeller plane. To compare the result in Figure 14, please refer the Figure 4.13 on p. 36 in [49]. The area of high velocity behind each blade body is described well for real propeller simulations (CFDSHIP-IOWA using overset grid) [49], LDV experiments [49], and the proposed rotating blade propeller model. The velocity drops between the two blades are captured as well.
However, when using the rotating blade propeller model the lower velocity area outside the propeller radius next to each blade tip is not clear for the E1619 propeller ( Figure 14). In contrast, it is clearer for the KP505 propeller in Figure 15. It can also be observed clearly in the real propeller simulation [49] and the LDV experiment [49]. Because of the small chord length and expanded blade area ratio of E1619, the grid numbers across its narrow tip are less than can be covered by that of the KP505. The grid number is insufficient to adequately resolve the flow field phenomena. This issue will be investigated by a finer grid in the future. However, when using the rotating blade propeller model the lower velocity area outside the propeller radius next to each blade tip is not clear for the E1619 propeller ( Figure  14). In contrast, it is clearer for the KP505 propeller in Figure 15. It can also be observed clearly in the real propeller simulation [49] and the LDV experiment [49]. Because of the small chord length and expanded blade area ratio of E1619, the grid numbers across its narrow tip are less than can be covered by that of the KP505. The grid number is insufficient to adequately resolve the flow field phenomena. This issue will be investigated by a finer grid in the future.  However, when using the rotating blade propeller model the lower velocity area outside the propeller radius next to each blade tip is not clear for the E1619 propeller ( Figure  14). In contrast, it is clearer for the KP505 propeller in Figure 15. It can also be observed clearly in the real propeller simulation [49] and the LDV experiment [49]. Because of the small chord length and expanded blade area ratio of E1619, the grid numbers across its narrow tip are less than can be covered by that of the KP505. The grid number is insufficient to adequately resolve the flow field phenomena. This issue will be investigated by a finer grid in the future.  The propeller downstream vector field was presented in Figure 16 for the rotating blade model. Due to lack of an experimental vector field, the PIV measurement behind a submarine with an X-rudder [50,51], e.g., Figure 13 in [51], could be used a reference for a brief comparison. A submarine is a slender body and the rudder is outside the propeller radius. Both results show that the vector directions are twisted in front of the blade's leading edge in the propeller's rotational direction. The vectors pointing to the propeller center change to the propeller's rotational direction around 0.7 R. The axial velocity increases behind the blade body.
The propeller downstream vector field was presented in Figure 16 for the rotating blade model. Due to lack of an experimental vector field, the PIV measurement behind a submarine with an X-rudder [50,51], e.g., Figure 13 in [51], could be used a reference for a brief comparison. A submarine is a slender body and the rudder is outside the propeller radius. Both results show that the vector directions are twisted in front of the blade's leading edge in the propeller's rotational direction. The vectors pointing to the propeller center change to the propeller's rotational direction around 0.7 R. The axial velocity increases behind the blade body.   Table 9 shows the results of the propeller thrust coefficient KT for the experiments and CFD using the average and the rotating blade body force model for the behind-hull condition. The propeller rotational speed is 16.3 RPS and was obtained using a selfpropulsion test in a towing tank at Osaka University for a 3.2 m long ship model at Froude number Fr = 0.26 [41]. Both models under-predicted the KT by around 1%. The rotating blade result was obtained by restarting from the average model result to save computational time. The time step is required to be very small for the rotating blade model to obtain a converged and stable result. In the presented work, the time step for the average model was around 0.0055 s, while it was 0.00022 s for the blade rotating model (around 1.3 degrees of propeller rotation).  Figure 17 presents the flow field plotted by the axial velocity contour (u/U) and vector field for cross flow (v/U, w/U) after the rudder, i.e., x/L = 1.025. U is the ship model speed.

KP505 Propeller behind KCS Hull
x/L = 0.0 is at the ship FP (front perpendicular). Figure 17a is the PIV (particle image velocimetry) measurement. Figure 17b is the result of the CFD using the average model. Figure 17c is the result of the CFD using the rotating blade model. The three results show a similar flow pattern. Because of the right-rotating propeller, high axial velocity is located mainly in the starboard. The highest u/U values in (a) and (c) are in a more similar position near z/L = −0.04, where is the lower part of the starboard side flow field. The highest u/U in Figure 17b is a longer band region in the middle of the starboard side flow field. The clockwise-rotating flow is interfered by the rudder in the middle, so the portside cross flow was twisted to form an upward jet and it is pointing downward in the starboard. The hub vortex is separated, becoming two cores as well: the portside core is in the lower position and starboard core is in the higher position. A vortex was induced by outer upward flow and starboard downward cross flow around (y/L, z/L) = (0.01, −0.02). This is the corner of the high u/U contour (in yellow color, i.e., u/U = 1.2-1.3). In Figure 17a,c, the corner shape is sharper than the one in Figure 17b. Basically, the flow field away from the propeller is nearly steady even though the propeller rotational flow is unsteady. This can also be seen in Figure 18 for the flow field between the propeller and rudder, i.e., x/L = 0.99. The propeller plane is at x/L = 0.9825 (the ship AP is at x/L = 1.0 and the ship FP is at x/L = 0). In Figure 18, both models show a similar flow field pattern, but the average model has a larger high-axial-velocity area in the starboard side.  When the wake field is near the propeller, the rotating blade model reveals an unsteady phenomena. The x/L = 0.985 section is right behind the propeller plane (x/L = 0.9825). Figure 19 indicates the average model result. Because of the right-rotating propeller, the region of high axial velocity is observed in the starboard side. This is a time-averaged and steady flow field. Instead, the results of the rotating blade model plotted in Figure 20 for different blade positions show an unsteady flow field changing by time (the force computation is still under a quasi-steady assumption). In Figure 20, the flow field around every 18 degrees of the index blade rotation is plotted. Four time instants are divided, from the index blade near approximate 0 degrees (tip pointing up vertically) to 54 degrees. The index blade at 72 degrees is the next blade at 0 degrees again because it is a five-bladed propeller (360/5 = 72 degrees). At a different time instant, the high-axial-velocity area still occurs in the starboard side but at a lower position, i.e., the fourth quadrant ( Figure 19). As the blade rotates into the fourth quadrant one after another, high axial velocity is generated on the blade area inside the quadrant and moves along with the blade until it leaves the quadrant.  When the wake field is near the propeller, the rotating blade model reveals an unsteady phenomena. The x/L = 0.985 section is right behind the propeller plane (x/L = 0.9825). Figure 19 indicates the average model result. Because of the right-rotating propeller, the region of high axial velocity is observed in the starboard side. This is a time-averaged and steady flow field. Instead, the results of the rotating blade model plotted in Figure 20 for different blade positions show an unsteady flow field changing by time (the force computation is still under a quasi-steady assumption). In Figure 20, the flow field around every 18 degrees of the index blade rotation is plotted. Four time instants are divided, from the index blade near approximate 0 degrees (tip pointing up vertically) to 54 degrees. The index blade at 72 degrees is the next blade at 0 degrees again because it is a five-bladed propeller (360/5 = 72 degrees). At a different time instant, the high-axial-velocity area still occurs in the starboard side but at a lower position, i.e., the fourth quadrant ( Figure 19). As the blade rotates into the fourth quadrant one after another, high axial velocity is generated on the blade area inside the quadrant and moves along with the blade until it leaves the quadrant. When the wake field is near the propeller, the rotating blade model reveals an unsteady phenomena. The x/L = 0.985 section is right behind the propeller plane (x/L = 0.9825). Figure 19 indicates the average model result. Because of the right-rotating propeller, the region of high axial velocity is observed in the starboard side. This is a time-averaged and steady flow field. Instead, the results of the rotating blade model plotted in Figure 20 for different blade positions show an unsteady flow field changing by time (the force computation is still under a quasi-steady assumption). In Figure 20, the flow field around every 18 degrees of the index blade rotation is plotted. Four time instants are divided, from the index blade near approximate 0 degrees (tip pointing up vertically) to 54 degrees. The index blade at 72 degrees is the next blade at 0 degrees again because it is a five-bladed propeller (360/5 = 72 degrees). At a different time instant, the high-axial-velocity area still occurs in the starboard side but at a lower position, i.e., the fourth quadrant ( Figure 19). As the blade rotates into the fourth quadrant one after another, high axial velocity is generated on the blade area inside the quadrant and moves along with the blade until it leaves the quadrant. Figure 19. Wake field after propeller (x/L = 0.985) for the average model. Figure 19. Wake field after propeller (x/L = 0.985) for the average model. Appl

Conclusions
The rotating blade propeller model has been developed successfully based on the blade element theory and using total velocity from RANS solution. A cell search algorithm was implemented to identify the RANS grid points inside the 2D propeller geometry. The propeller computation is based on the 25% chord length position along a foil section. The body forces were distributed linearly across the chord line with a maximal value at the trailing edge and a value of zero at the linear edge. Simulations for the KP505 and E1619 propellers in open water, and the KP505 propeller behind the KCS hull were successfully conducted.
For KP505 propeller, the results showed good agreement with experimental data and had a small error mainly in high-load conditions. The propeller thrust coefficient was predicted well by the rotating blade model. For the KP505 propeller in open water conditions,

Conclusions
The rotating blade propeller model has been developed successfully based on the blade element theory and using total velocity from RANS solution. A cell search algorithm was implemented to identify the RANS grid points inside the 2D propeller geometry. The propeller computation is based on the 25% chord length position along a foil section. The body forces were distributed linearly across the chord line with a maximal value at the trailing edge and a value of zero at the linear edge. Simulations for the KP505 and E1619 propellers in open water, and the KP505 propeller behind the KCS hull were successfully conducted.
For KP505 propeller, the results showed good agreement with experimental data and had a small error mainly in high-load conditions. The propeller thrust coefficient was predicted well by the rotating blade model. For the KP505 propeller in open water conditions, the absolute error of the thrust coefficient was generally around or less than 3%. The lowest error was less than 1% under-prediction for the advance coefficient J = 0.5. Prediction of torque coefficient was not as successful because of the constant drag coefficient assumption. Grid independence was achieved for J = 0.5, 0.7, and 0.8. The CFD result was also validated for J = 0.8, i.e., the comparison between the converged CFD result and the experiment was more confident than the comparison between the CFD results themselves.
In open water and behind-hull conditions, the unsteady phenomena and blade-toblade effect occur only when the propeller wake field is very close to the propeller. This also supports the conclusion that the body force propeller model can be applied to many ship applications.
For KP505 appended and operating behind the KCS hull in calm water, the thrust coefficient of the rotating blade model error was about 1%. For the downstream wake field, the rotating blade model provided a steady result closer to the experimental result than the average model's result. The high-velocity distribution was more scattered and had a lower value because the body forces were only distributed on the blades and the thrust was under-predicted. For the wake field right after the propeller, the unsteady phenomena were obvious. As the blades rotated into the fourth quadrant in sequence, the starboard area of high axial velocity appeared on the incoming blade and moved together with the blade until the blade left the quadrant. For the average model, the area of high axial velocity simply remained on the starboard side constantly and did not change over time.
For E1619 in open water condition, the error was under-predicted by around 6%. The high-velocity areas behind each blade and velocity drops between the two blades are sufficiently captured by the rotating blade model. However, the area of lower velocity outside the propeller tip was not clear, because the grid number was not sufficient around the narrow blade tip.
The contribution of the presented study is the implementation of the total velocity concept and blade rotation with 2D blade geometry in a CFD body force propeller model. Discretized real propeller simulation is still time consuming. Grid generation is difficult and complicated. The rotating blade body force model will be a great alternative for considering unsteady flow field phenomena and blade-to-blade effect. Only an additional cylinder disk is required. Grid generation is easy with a lower grid number. The same grid system can be applied for different propellers. In this study, it was also proven that the propeller model works in open water and behind-hull conditions.
In future work, the propeller computation will be developed further to solve the vortex lattice on the grid cell found on propeller plane. Additionally, the linear body force distribution will be replaced with data from literature or a more realistic curve to represent the pressure difference distribution on a foil. Furthermore, side/lateral geometry will be included to consider propeller chamber and thickness. To do so, a 3D cell search algorithm will be required to identify the RANS grid point inside the body volume of the real propeller. Using this method, the body force will be able to be distributed in an arbitrary 3D shape such as a realistic propeller geometry.