Investigating the Performance of Twin Marine Propellers in Different Ship Wake Fields Using an Unsteady Viscous and Inviscid Solver
https://doi.org/10.1007/s11804022002796

Abstract
In this study, the performance of a twinscrew propeller under the influence of the wake field of a fully appended ship was investigated using a coupled Reynoldsaveraged Navier—Stokes (RANS)/boundary element method (BEM) code. The unsteady BEM is an efficient approach to predicting propeller performance. By applying the timestepping method in the BEM solver, the trailing vortex sheet pattern of the propeller can be accurately captured at each time step. This is the main innovation of the coupled strategy. Furthermore, to ascertain the effect of the wake field of the ship with acceptable accuracy, a RANS solver was developed. A finite volume method was used to discretize the Navier—Stokes equations on fully unstructured grids. To simulate ship motions, the volume of the fluid method was applied to the RANS solver. The validation of each solver (BEM/RANS) was separately performed, and the results were compared with experimental data. Ultimately, the BEM and RANS solvers were coupled to estimate the performance of a twinscrew propeller, which was affected by the wake field of the fully appended hull. The proposed model was applied to a twinscrew oceanography research vessel. The results demonstrated that the presented model can estimate the thrust coefficient of a propeller with good accuracy as compared to an experimental selfpropulsion test. The wake sheet pattern of the propeller in open water (uniform flow) was also compared with the propeller in a real wake field.Article Highlights• The performance of a twinscrew propeller under the influence of the wake field of a fully appended ship was investigated using a coupled RANS/BEM aproach.• The trailing vortex sheet pattern of the propeller accurately captured at each time step.• Selfpropulsion experiment was performed numarically using this coupled solvers and results compared with physical wave tank data. 
1 Introduction
Propeller design in ship propulsion systems is usually confined to the analysis of steadystate conditions, deeply submerged propellers, or open water tests. However, when precious estimations (e.g., special maneuverability) are required, a sophisticated consideration should be applied to adequately determine propeller performance. In numerical studies, such a complex consideration can be addressed by Reynoldsaveraged NavierStokes (RANS) solvers, such as the finite volume method (FVM) or finite element method. In these methods, the propeller behind the ship can be simulated with many details, such as the influence of the velocity wake filed of a ship (and its appendages) on the propeller. Nonetheless, RANS solvers are timeconsuming methods, especially in simulating (two or more) moving objects that should be coupled with a movinggrid technique, e.g., overset method (Benek et al. 1983; Meakin and Suhs 1989; Henshaw and Schwendeman 2006).
Furthermore, inviscid approaches, such as the boundary element method (BEM)(Hess and Smith 1962, Hess 1972, Hess and Valarezo 1985), are efficient techniques for predicting the pressure force of lifting bodies. In recent decades, many scientists have studied and developed boundary element formulations to improve their accuracy. Politis (2004) used a threedimensional (3D) unsteady BEM to solve flows around a marine propeller. He used the timestepping method (TSM) to model the vortex sheet emanating from each blade trailing edge and showed that this approach could produce a very stable rollup wake pattern related to the prescribed wake shape (PWS) and wake relaxation method (WRM). Politis (2011) applied his solver to a wide range of working conditions of a propeller, such as a propeller in a straight/circular course, an azimuth propeller, and a ducted propeller in a circular course. He also showed that the burst starting of a propeller may be modeled with a good agreement through his method compared with experimental data. Subsequently, Politis (2016) introduced an artificial viscosity filtering technique to a domain that suppresses smallscale instabilities and leaves largescale vortices. These vortices were used to determine the rollup pattern of free wake vorticity. He applied his model to a number of flows around liftproducing configurations, such as a burst starting propeller. Finally, he concluded that this model could be used as an attractive alternative in the case of preliminary force predictions. Durante et al. (2013) compared different hydrodynamic models to analyze the load delivered by the marine propeller PROP 3714 in an axial nonuniform wake flow to find a reliable and fast model to balance the accuracy and computational time. They used the BEM model as a reference point for their calculations. Wang et al.(2016) described four different numerical methods to implement wake alignment algorithms for a marine propeller in steadystate flows. They applied all four methods to calculate the performance of a KP505 propeller. To show the wake geometry, a rectangular foil was simulated, and the wake shape was compared. In a similar effort to determine an efficient method to solve offdesign aspects related to ships and their propulsion system in the earlier design stages, Gaggero et al. (2019) investigated the capability of the BEM and blade element momentum theory (BEMT) propeller solver to predict propeller performance in a maneuvering condition. The comparison between the BEM and BEMT was performed behind the hull condition of a twinscrew ship. The results show that although the BEM is at least one order slower than the BEMT, it is more accurate and could predict propeller performance behind ships with acceptable accuracy.
Meanwhile, BEMs have been extensively improved, but designers still need the RANS solver. This fact motivated scientists to combine the efficiency of the BEM with the accuracy of the RANS solver to reduce the numerical effort without compromising the predictive accuracy. Greve, WöcknerKluwe et al. (2012) analyzed the dynamics of a singlescrew propeller ship in seaways by combining the BEM and RANS solvers. In their algorithm, the propeller force was computed in the BEM solver, while the inflow condition for the propeller was provided by the RANS method. Although the WRM was used to simulate propeller vortex sheets, the propeller thrust and ship resistance were the main investigated parameters in their work (no wake shape was shown). Rijpkema et al. (2013) numerically investigated the propeller–hull interaction of the wellknown singlescrew Kriso Container Ship (KCS) using a combination of the steadystate RANS solver for the ship hull flow and the unsteady potential method for propeller loading. Compared with experimental data, this approach predicted the propeller thrust forces to be within 2% to 3%, with a significant reduction in the computation time. Conversely, Queutey et al. (2013) tried to find an efficient way to simulate marine propellers behind a ship. They compared a full RANS solver with a steady coupled RANSBEM to predict the propulsion performance of a propelled ship. Code PROINS, which was developed by CNRINSEAN, was applied as the BEM computational model, while the unstructured finite volume solver ISISCFD was applied as the RANS approach. The STREAMLINE tanker and propeller were considered in the validation test case. Because of the massive computation time, the RANS simulation was performed only for the design speed, while the steady hybrid RANSBEM computations were performed for the speed range. They show that the CPU time for the steady RANSBEM coupling is divided by a factor of 50 as compared to the full RANS (with reasonable accuracy). Hence, it can be used for preliminary design purposes. Gaggero et al. (2014) analyzed the reliability of different numerical approaches (BEM and RANS solver) to predict the performance of the wellknown highly skewed marine propeller SeiunMaru, which operates in different wakes. In their BEM model, a frozen wake geometry was adopted as the result of the steady alignment with the propeller in a uniform inflow. They used StarCCM+ commercial software as a RANS solver. The results showed that the numerical code presents a good capability in predicting propeller characteristics and pressure distribution when the thrust identity approach is adopted. Rao and Yang (2017) presented a hybrid approach (panel method and RANS) to predict the selfpropulsion performance and effective wake field of an underwater vehicle. Numerical simulations were performed for the full appendage SUBOFF submarine. They used two different ways to evaluate the wake field and showed that the effective wake fractions are approximately 4% larger than the nominal. Conversely, Gaggero et al. (2017) presented a coupled BEM/RANS approach to simulate the selfpropulsion prediction of a KCS hull. The total hull resistance was computed by the RANS solver without a propeller, but they used the BEM method for propeller simulations, which was much faster than the finite volume moving mesh technique (e.g., overlap grids). Gaggero et al. (2019) presented a propeller solver based on the BEM to analyze the propeller performance in different flow conditions. The nominal wake of a twinscrew ship was first evaluated using the unsteady commercial RANS (StarCCM+) solver, and then it was used as an offline input to the propeller solver. The wake shape of the propeller was modeled by the PWS method. They also compared the results to a steadyblade element approach to determine a fast and less timeconsuming method to simulate propellers.
In the current study, a numerical solver was developed based on the coupled BEM/RANS approach to solve flows around twinscrew marines, considering the influence of the ship wake field. The main motivation of this study is to combine the advantages of both methods to reduce the computational time without compromising the predictive accuracy.
In the BEM solver, each propeller is supposed to consist of lifting bodies (blades). The TSM was used to simulate the freely moving unsteady trailing vortex sheet, which emanates from the trailing edge of each blade. This is the main novelty in combining the BEM and RANS solvers. This approach correctly simulates and visualizes the vortex bladewake interaction at each time step. The BEM solver was applied to a highly skewed controllable pitch marine propeller, and the results were compared with experimental data.
Meanwhile, in the RANS solver, two main equations must be simultaneously solved. The FVM was used to discretize the unsteady Navier–Stokes (NS) equation in arbitrary unstructured grids. The second equation is the scalar transport equation, which is solved to calculate the volume of the fraction. The highresolution compressive interface capturing scheme for arbitrary meshes (CICSAM) approach was applied to capture the fluid interface on the meshes with an arbitrary topology. Combining the two equations with rigid body motions makes it possible to simulate ships in calm water. The RANS solver was also validated by calculating the resistance on a highspeed catamaran and fully appended 50m oceanography research vessel (RV50) in calm water. The results were compared with a towing tank test and similar work. Finally, the coupled approach was applied to estimate the twin propeller performance behind the ship. The wake sheet pattern was also shown to be influenced by the ship wake field, which was consequently affected by ship appendages.
2 Governing equation
2.1 BEM method (inviscid solver)
As pointed out earlier, the propeller forces were calculated using the BEM. Considering incompressible, irrotational, and nonviscous flows, the BEM is based on the solution of a Laplacian equation:
$$ \nabla^{2} \phi(\boldsymbol{x}, t)=0 $$ (1) where ϕ(x, t) is the velocity potential in the discretized domain. ϕ(x, t) can be written as a linear combination of fundamental potential functions, such as sink/source, dipole, and vortex, which could be computed on each cell center. Based on Green's second identity, the 3D problem of Eq. (1) can be solved as a simple integral problem on the boundary surface of the computational domain:
$$ 2 \pi \phi(\boldsymbol{x}, t)=\int\limits_{S} \phi(\boldsymbol{x}, t) \frac{\partial \lambda}{\partial n} \mathrm{~d} S+\int\limits_{S} \lambda \frac{\partial \phi(\boldsymbol{x}, t)}{\partial n} \mathrm{~d} S $$ (2) where $\lambda=\frac{1}{r}$ is defined as the fundamental solution for the 3D potential flow. n and r are the unit normal vectors of the panel (pointing outward of the body) and the distance between the source point and field point, respectively. As a result, Eq. (2) can be written in the form of
$$ \begin{gathered} 2 \pi \phi(\boldsymbol{x}, t)\int\limits_{S_{W}+S_{B}} \phi(\boldsymbol{x}, t) \frac{\partial\left(\frac{1}{r}\right)}{\partial n} \mathrm{~d} S \\ =\int\limits_{S_{W}+S_{B}} \frac{1}{r} \frac{\partial \phi(\boldsymbol{x}, t)}{\partial n} \mathrm{~d} S \end{gathered} $$ (3) The surface cells in the discretized domain are divided into three main zones. Depending on the related zone, each cell is treated differently. We consider a threeblade propeller, which operates on an unbounded flow, as shown in Figure 1. Surface cells are divided into the following three main zones:
Zone 1 (blade/lifting cells (S_{Bl})): These cells are contributed to generating the lifting force. The wake sheet is generated from the leading edge of lifting cells.
Zone 2 (Kuttastrip cells (S_{k})): These cells are located in the first row and are adjacent to the trailing edges of the blade body.
Zone 3 (free wake sheet cells (S_{fw})): These cells include all cells on the wake sheet except Kuttastrip cells.
The wake sheet cells are all cells in the second and third zones (S_{W} = S_{k} ∪ S_{fw}).
The discretization of Eq. (3) is shown in Eq. (4) (Politis 2011):
$$ \begin{aligned} 2 \pi \phi(\boldsymbol{x}, t) &+\sum\limits_{S_{B}} \phi(\boldsymbol{x}, t) \int \frac{\boldsymbol{n} \cdot \boldsymbol{r}}{r^{3}} \mathrm{~d} S+\sum\limits_{S_{k}} \mu(\boldsymbol{x}, t) \int \frac{\boldsymbol{n} \cdot \boldsymbol{r}}{r^{3}} \mathrm{~d} S \\ &=\sum\limits_{S_{B}} \int \frac{\boldsymbol{n} \cdot \boldsymbol{q}}{r} \mathrm{~d} S\sum\limits_{S_{\mathrm{fr}}} \int \mu \frac{\boldsymbol{n} \cdot \boldsymbol{r}}{r^{3}} \mathrm{~d} S \end{aligned} $$ (4) The first and second terms on the righthand side are obtained from the previous time step; hence, they can be explicitly calculated (Najafi and Abbaspour 2017). The dipole intensity over body cells (ϕ) and Kuttastrip cells (μ) should be calculated by forming a system of linear equations. The pressuretype Kutta condition states that pressure should be continuous when the flow approaches the trailing edge from either the suction or pressure sides. In this study, the linear Morinotype Kutta condition was used at each time step:
$$ \phi_{\text {upper cells }}\phi_{\text {lower cells }}=\mu $$ (5) Calculating Eq. (4) for all body and Kutta cells leads to a linear system of equations in the form of (AX = B), which could be solved by common iterative methods. This equation also gives the dipole intensity over the body and Kutta cells.
As mentioned before, the TSM was used to investigate the wake sheet pattern in all simulations. In the TSM scheme, the wake sheet is developed through time, which implies that at each time step, a new row of nodes is added while the position of previous nodes is updated. It is thus necessary to calculate the velocity of wake sheet nodes at each time step. The velocity of wake sheet nodes is referred to as the perturbation velocity, which could be calculated by taking the gradient of Eq. (4) as follows (Zhu et al. 2002):
$$ \begin{aligned} \boldsymbol{v}_{m} &=\frac{1}{4 \pi} \int_{S_{B}}(\boldsymbol{n} \cdot \nabla \boldsymbol{\phi}) \frac{\boldsymbol{r}}{r^{3}} \mathrm{~d} S+\frac{1}{4 \pi} \int_{S_{B}}(\boldsymbol{n} \times \nabla \boldsymbol{\phi}) \times \frac{\boldsymbol{r}}{r^{3}} \mathrm{~d} S \\ &+\frac{1}{4 \pi} \int_{S_{W}}(\boldsymbol{n} \times \nabla \boldsymbol{\phi}) \times \frac{\boldsymbol{r}}{r^{3}} \mathrm{~d} S\frac{1}{4 \pi} \int_{L^{\prime}} \phi \frac{\boldsymbol{d} \boldsymbol{l} \times \boldsymbol{r}}{r^{3}} \mathrm{~d} S \end{aligned} $$ (6) The integral zones for the first three terms are shown in Figure 2. The fourth term is calculated over the lines bounding the free wake sheet. Moreover, it is a line integral over the bounding of the wake surface and is depicted in Figure 2 for a single blade (Katz and Plotkin 2001). If the last term is ignored, the starting vortex and wake rollup pattern will not be correctly formed.
After calculating the perturbation velocity for all wake sheet nodes, the new position of the wake sheet can be calculated by a simple linear formula as follows (Zhu et al. 2002):
$$ \boldsymbol{X}_{\text {new }}=\boldsymbol{X}_{\mathrm{old}}+\left(\boldsymbol{v}_{m}+\boldsymbol{V}_{\mathrm{rel}}\right) \mathrm{dt} $$ (7) where V_{rel} is the relative velocity of each cell:
$$ \boldsymbol{V}_{\mathrm{rel}}=\boldsymbol{V}(t)+\boldsymbol{\Omega}(t) \times \boldsymbol{r} $$ (8) When the potential velocity distribution is calculated over each body cell using the unsteady Bernoulli equation, the pressure can be computed by Najafi and Abbaspoor (2019).
$$ \frac{PP_{\infty}}{\rho}=\frac{\mathrm{d} \phi}{\mathrm{d} t}\frac{1}{2}\left\nabla \boldsymbol{\phi}\boldsymbol{V}_{\mathrm{rel}}\right^{2}+\frac{1}{2}\left\boldsymbol{V}_{\mathrm{rel}}\right^{2} $$ (9) The hydrodynamic force and moment over each blade are calculated by integrating the pressure (normal stress) and virtual friction (shear stress) forces over the body. More details can be found in Politis (2005).
2.2 RANS solver (viscous method)
The twophaseflow 3D NS equation was used as the governing equation for the viscous solver to simulate flows around the ship hull. In almost all marine problems, the Mach number is less than 0.3, so the NS equations are used in their incompressible form. Accordingly, the momentum and continuity equations can be simplified as follows:
$$ \begin{gathered} \rho\left(\frac{\partial u}{\partial t}+u \frac{\partial u}{\partial x}+v \frac{\partial u}{\partial y}+w \frac{\partial u}{\partial z}\right)=\frac{\partial p}{\partial x} \\ +\mu\left(\frac{\partial^{2} u}{\partial x^{2}}+\frac{\partial^{2} u}{\partial y^{2}}+\frac{\partial^{2} u}{\partial z^{2}}\right) \\ \rho\left(\frac{\partial v}{\partial t}+u \frac{\partial v}{\partial x}+v \frac{\partial v}{\partial y}+w \frac{\partial v}{\partial z}\right)=\frac{\partial p}{\partial y} \\ +\mu\left(\frac{\partial^{2} v}{\partial x^{2}}+\frac{\partial^{2} v}{\partial y^{2}}+\frac{\partial^{2} v}{\partial z^{2}}\right)+\rho g_{y} \\ \rho\left(\frac{\partial w}{\partial t}+u \frac{\partial w}{\partial x}+v \frac{\partial w}{\partial y}+w \frac{\partial w}{\partial z}\right)=\frac{\partial p}{\partial z} \\ +\mu\left(\frac{\partial^{2} w}{\partial x^{2}}+\frac{\partial^{2} w}{\partial y^{2}}+\frac{\partial^{2} w}{\partial z^{2}}\right)+\rho g_{w} \\ \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0 \end{gathered} $$ (10) In cases where there are twophase fluid simulations, one should solve the volume fraction transport equations:
$$ \frac{\partial \alpha}{\partial t}+\nabla \cdot(\alpha \boldsymbol{u})=0 $$ (11) In this equation, α is the fractional value demonstrating the amount of water in the cell:
$$ \alpha=\left\{\begin{array}{lc}1 & \text { Cell inside water } \\ 0 & \text { Cell inside air } \\ 0<\alpha<1 & \text { Transitional area }\end{array}\right. $$ (12) Generally, Eqs. (10) and (11) should be solved simultaneously. For this purpose, the fluid density and viscosity of each cell in the computational domain can be expressed as a function of the volume fraction α as follows:
$$ \begin{gathered} \rho=\rho_{1} \alpha+\rho_{2}(1\alpha) \\ \mu=\mu_{1} \alpha+\mu_{2}(1\alpha) \end{gathered} $$ (13) where subscripts 1 and 2 denote water and air, respectively.
2.2.1 Discretization of NS equations
In the present study, the FVM was applied to discretize the NS equations. Integrating these equations over each control volume can result in the following equation:
$$ \begin{array}{c} \frac{{\rm{d}}}{{{\rm{d}}t}}\int\limits_V {\mathit{\boldsymbol{u}}{\rm{d}}} V + \int\limits_A {\mathit{\boldsymbol{u}}(\mathit{\boldsymbol{u}} \cdot \mathit{\boldsymbol{n}}){\rm{d}}} A = \int\limits_A {v\nabla \mathit{\boldsymbol{u}} \cdot \mathit{\boldsymbol{n}}{\rm{d}}} A\\  \frac{1}{\rho }\int\limits_V {\nabla \mathit{\boldsymbol{P}}{\rm{d}}} V + \int\limits_V {\mathit{\boldsymbol{g}}{\rm{d}}} V \end{array} $$ (14) For the fluid velocity u_{i}, the diffusion term (the first term on the righthand side of (14)) was discretized using overrelaxed interpolation (Jasak 1996):
$$ \int\limits_A v \nabla \mathit{\boldsymbol{u}} \cdot \mathit{\boldsymbol{n}}{\rm{d}}A = \int\limits_V v \nabla \cdot \left({\nabla {\mathit{\boldsymbol{u}}_i}} \right){\rm{d}}V = \sum\limits_{{\rm{faces }}} {{v_f}} {\mathit{\boldsymbol{A}}_{\dot f}}{\left({\nabla {\mathit{\boldsymbol{u}}_i}} \right)_f} $$ (15) To discretize the convection term (the second term on the lefthand side of (14)), the fluid velocity at the control volume faces needed to compute fluxes. This condition can be evaluated as follows:
$$ \int\limits_A {{u_i}} \left({{\mathit{\boldsymbol{u}}_i} \cdot \mathit{\boldsymbol{n}}} \right){\rm{d}}A = \sum\limits_{{\rm{faces }}} {{u_f}} {F_f} $$ (16) To calculate u_{f}, the Gamma interpolation scheme was used, which was introduced by (Jasak 1996).
The Crank–Nicholson scheme was also used for the time discretization of diffusion and convection terms in Eq. (14).This scheme has secondorder accuracy and increases stability.
The pressure term (second term on the righthand side of (14)) is discretized as
$$ \int\limits_{A} P \boldsymbol{n} \cdot \mathrm{d} \boldsymbol{A}=\sum\limits_{\text {faces }} P_{f} A_{f} $$ (17) where A_{f} is the direction component of the face area vector.
Using the common linear interpolation between two neighboring control volume centers results in severe oscillations in the velocity field. This is of great importance, especially when there are two fluids with highdensity ratios, e. g., water and air. Here a piecewise linear interpolation was used, which was developed by Panahi et al. (2009).
2.2.2 Discretization of the volume of fluid equation
The finite volume discretization of the volume fraction transport Equation (11) is based on the integration over the control volumes and time step:
$$ \int\limits_{t}^{t+\delta t}\left(\int\limits_{V} \frac{\partial \alpha}{\partial t} \mathrm{~d} V\right) \mathrm{d} t+\int\limits_{t}^{t+\delta t}\left(\int\limits_{V} \nabla \cdot \alpha \boldsymbol{u} \mathrm{d} V\right) \mathrm{d} t=0 $$ (18) The first term is a simple integral form. By applying the Gauss theorem on the second term and assuming a small variation of F_{f}, we can derive
$$ \left(\alpha_{P}^{t+\delta t}\alpha_{P}^{t}\right) \frac{V_{P}}{\Delta t}=\sum\limits_{\text {faces }} \frac{1}{2}\left(\alpha_{f}^{t}+\alpha_{f}^{t+\delta t}\right) F_{f} \delta t $$ (19) where α_{f} must be approximated on the face of each control volume. A simple interpolation results in nonphysical values, leading to the use of a highorder composite interpolation. In this study, the compressive interface capturing scheme for arbitrary meshes (CICSAM) was applied (Ubbink and Issa 1999). CICSAM uses convection boundedness criteria (CBC) and ultimate quickest (UQ) interpolations by introducing a weighting factor γ_{f} (Gaskell and Lau 1988; Leonard 1991). This factor was evaluated according to the slope of the free surface and the direction of fluid motion (Ubbink and Issa 1999). Based on the normalized variable diagram (NVD), a normal face value was obtained (Leonard 1991):
$$ \tilde{\alpha}_{f}=\gamma_{f} \tilde{\alpha}_{f_{\mathrm{CBC}}}+\left(1\gamma_{f}\right) \tilde{\alpha}_{f_{\mathrm{UQ}}} $$ (20) NVD definitions are shown in Figure 3, where donor, acceptor, and upwind cells are defined according to the flow direction for each control volume face:
$$ \tilde{\alpha}_{f}=\frac{\alpha_{f}\alpha_{U}}{\alpha_{A}\alpha_{U}} $$ (21) Substituting Eq. (20) into (21) results in the estimation of α_{f}
In the volume of fluid (VOF) approach, the Courant – Friedrichs–Lewy (CFL) number plays an important role in the stability of a solution. Here, the Courant number was calculated according to Eq. (22) for each cell and was set to be less than 0.2 for all twophase flow test cases Panahi et al. (2009).
$$ \mathrm{CFL}=\sum\limits_{\text {faces }} \max \left(\frac{\text { Flux } \times \mathrm{d} t}{\text { Cell volume }}, 0\right) $$ (22) As mentioned before, a body attached mesh was used to simulate body motions. In this way, the relative velocity should be applied in all equations:
$$ \boldsymbol{u}=\boldsymbol{u}_{\text {fluid }}\boldsymbol{u}_{\text {body }} $$ (23) where u_{fluid} is the fluid velocity vector and u_{body} is the velocity of a moving body (which is equal to the mesh velocity).
2.2.3 Coupling of the pressure and velocity
The method used to couple the pressure and velocity fields is the modified fractional steps, which were first proposed by Kim and Choi (2000). We considered the NS equations for the velocity component u_{i} in its discretized form on a given control volume. The first step was to solve the first intermediate velocity (û_{i}) with the pressure gradient in the previous time step G_{i}(P^{n}):
$$ \frac{\hat{u}_{i}u_{i}^{n}}{\Delta t}=\frac{1}{2}\left[H\left(u_{i}^{n}\right)+H\left(\hat{u}_{i}\right)\right]\frac{1}{\rho} G_{i}\left(P^{n}\right)+K_{i} $$ (24) where
$$ \begin{array}{c} H\left({{u_i}} \right) = \int\limits_A v \frac{{\partial {u_i}}}{{\partial n}}{\rm{d}}A  \int\limits_A {u_i^n} \mathit{\boldsymbol{U}}_f^n \cdot \mathit{\boldsymbol{n}}{\rm{d}}A\\ {G_i}(P) = \int\limits_A P {n_i}{\rm{d}}A\\ {K_i} = \int\limits_V {{g_i}} {\rm{d}}V \end{array} $$ (25) Then, the second intermediate velocity is calculated as
$$ u_{i}^{*}=\hat{u}_{i}+\frac{\Delta t}{\rho} G_{i}\left(P^{n}\right) $$ (26) Using Eq. (26) results in the pressure Poisson equation, which gives a new pressure field P^{n+1}:
$$ \oint_{A} \frac{1}{\rho} \frac{\partial P^{n+1}}{\partial n} \mathrm{~d} A=\left.\frac{1}{\Delta t} \oint_{A} u_{i}^{*} \mathrm{~d} A \frac{\partial P^{n+1}}{\partial n}\right_{\text {Boundary }}=0 $$ (27) By writing Eq. (27) for all cells in the computational domain, a system of equations in the form of (AX = B) was obtained. Solving the system of equations, the velocity at a new time step was computed as follows:
$$ u_{i}^{n+1}=u_{i}^{*}+\frac{\Delta t}{\rho} G_{i}\left(P^{n+1}\right) $$ (28) Because this new pressure was used to calculate "modified" velocities, global mass conservation was automatically satisfied. After the calculation of all velocity components in an outer iteration, the control volume face velocity vector was calculated for the next time step from Eq. (29). It includes the effect of the pressure gradient on the calculation of the face velocity vector to overcome the checkerboard pressure in the nonstaggered (collocated) arrangement:
$$ \boldsymbol{U}_{f}^{n+1}=L I\left(\boldsymbol{u}_{i}^{*}\right)+\left(\frac{\Delta t}{\rho} \frac{\partial P^{n+1}}{\partial n}\right)_{f} \boldsymbol{n} $$ (29) A comprehensive inhouse computational program was developed based on the mentioned approaches. More details about the implemented algorithms and their validation can be found in Pourmostafa and Ghadimi (2020b).
3 Validation of solvers
Before coupling the inviscid/viscous solvers, they should be validated separately for related test cases. Although each solver was separately validated in the different studies of (Pourmostafa and Ghadimi 2020a; 2020b), here, according to the final test case (research vessel and its highly skewed propeller), two validations of the solvers are presented. A viscous solver was used to simulate the motion of a highspeed catamaran and a 50 m research vessel (RV50) and calculate their resistance forces. The calculated resistance forces were compared with experimental towing tank data and similar works. Meanwhile, an inviscid solver was applied to predict the performance of a highly skewed marine propeller, which was designed for the same research vessel. The BEM results were compared to open water test results. In particular, in all cases, a mesh independent analysis was performed, and the best mesh size was properly selected.
3.1 Resistance and motion of a highspeed catamaran (validation of the viscous solver)
To validate the viscous solver, the motion of a highspeed catamaran was investigated at three velocities (Vel=5, 16.5, 24 kn). The main dimensions of the highspeed catamaran are shown in Table 1. Due to hydrodynamic forces, the catamaran tended to exhibit heave motions and pitching angles, which can be modeled by the developed solver with two degrees of freedom.
Length (m) 12.3 Breadth (m) 4.6 Draught (m) 0.95 Mass (t) 17.850 Vertical center of gravity (m) 0.45 Longitudinal center of gravity (m) 3.81 Block coefficient 0.33 The computational domain was discretized by 180992 hexahedron structural grids as shown in Figure 4.
When the main hydrodynamic parameters of the catamaran become steady, the heave, trim, and resistance forces were compared with those in the similar work of Panahi et al. (2009), as shown in Figure 5.
The viscous solver could predict the hydrodynamic force and motions of the highspeed catamaran with good accuracy. The catamaran at Vel = 24 kn is illustrated in Figure 6.
3.2 motion of RV50 in calm water (validation of the viscous solver)
RV50 is a 50 m twin propeller research vessel, which was designed for a cruise velocity of 12.3 kn. In this study, the viscous solver was used to calculate the vessel resistance in still water. Table 2 shows the main characteristic of this ship.
Length between perpendiculars (L_{BP}) (m) 44.70 Breadth (m) 10.00 Draught at midship (m) 3.30 Draught at forward perpendicular (m) 3.30 Draught at aft perpendicular (m) 3.30 Displacement mass (t) 896 Block coefficient 0.592 3 Cruise and maximum speed (kn) 12.314.5 The computational domain is displayed in Figure 7. A total of 306222 unstructured grids were used to discretize the RV50 with full appendages, as shown in Figure 8.
The simulations were performed at seven different velocities, and the results are compared with experimental data of (RV50 Test Report 2011) in Figure 10. The model resistance of the research vessel was measured at the towing tank of Vienna Model Basin with a wooden model on a 1: 12 scale (Figure 9), and the total resistance for the full scale was calculated using International Towing Tank Conference (ITTC) relations. In particular, minor resistance components, such as air resistance and hull roughness effects, were also added to the numerical model based on the experimental relations advised by the technical committee of the (ITTCRecommended Procedures and Guidelines 2011). Figure 11 displays the free surface elevation near the vessel at V = 14.5 kn.
The wake field velocity in the propeller plane (in the absence of a propeller) is usually referred to as a nominal wake. The nominal wake of he vessel at its top speed is shown in Figure 12, which was calculated based on (Rao and Yang 2017) as follows:
$$ 1w=\frac{\int\limits_{0}^{2 \pi} \int_{r_{h}}^{R} \frac{R V_{x}(r, \theta)}{V_{m}} r \mathrm{~d} r \mathrm{~d} \theta}{\int_{0}^{2 \pi} \int_{r_{h}}^{R} r \mathrm{~d} r \mathrm{~d} \theta}=0.9188 $$ (30) To calculate the disk average wake fraction, one should integrate velocity in this field. In this way, infinitesimal radial elements shown in Figure 12 are defined on the propeller plane.
Having velocities on these points, the average velocity could be calculated using Eq. (30). In the present research, the computational wake was compared with an experimental value measured in the towing test using 5hole pitot tube probes equal to w = 0.088. Table 3 shows the comparison between the results.
Nominal wake (w) Error (%) Computational Experimental 0.081 2 0.088 7.73 3.3 Performance of the RV50 marine propeller (validation of the inviscid solver)
The inviscid solver was validated through different test cases, such as the DTMB4119 and KP505 marine propellers used in Pourmostafa and Ghadimi (2020a). In this section, the validation of the RV50 propeller is presented because this propeller was used for the simulation in the coupling section. It is a controllable pitch propeller designed for the research vessel. The main particulars of this propeller are shown in Table 4.
Number of blades 4 Propeller diameter (m) 2.1 Hub Diameter (m) 0.56 Expanded blade area (m^{2}) 1.95 Mean pitch (mm) 1 691.1 Design J = V/nD 0.623 The BEM inviscid solver was used to calculate the hydrodynamic parameters of RV50 propeller. Figure 13 shows the comparison between the computation results and open water test results.
The solver could predict the thrust coefficient with reasonable accuracy. Furthermore, as viscous effects were virtually applied in the solver, the torque coefficients were estimated with great errors. This condition results in the efficiency of predicting a greater value than that of the real condition. Figure 14 shows the rollup pattern of the RV50 propeller at a design advance ratio.
4 Coupling algorithm
After separately validating the solvers, they can be coupled to simulate the twin propeller behavior in the wake field of the RV50 hull. In the coupled solver, first, the fully appended hull was simulated in the NS solver until all the hydrodynamic variables (velocities and pressure field) became steady. Subsequently, the wake field behind the hull was used as the inlet velocity for the propeller in the BEM solver. To accomplish this task, the corresponding volumetric cell should be found (at each time step) for each surface panel for velocity matching. Figure 15 schematically shows the following algorithm.
5 Results and discussion
Based on the presented coupling algorithm, the simulation of the twin propellers of RV50 was simultaneously performed with both propellers, as shown in Figure 16. The propellers' direction of rotation was outward, which implies that the right propeller rotates clockwise, whereas the left propeller rotates counterclockwise.
The simulations were performed based on the design parameters (Table 5) on three different ship velocities. As shown in the table, the propeller rotational speed was calculated on its design advance ratio (J = 0.67 − 0.83).
RV50 speed (kn) Propeller rotational speed (r/min) Advanced ratio (J) Time step (s) 7.50 133.21 0.83 0.008 10.50 198.38 0.78 0.005 14.50 (Top) 319.33 0.67 0.003 At each velocity, the simulation was performed (based on the algorithm in Figure 15) until the propeller forces became steady. Figure 17 shows the time history of the left propeller force for the ship at a top speed (the right propeller force is almost the same, so its difference is not recognizable in this chart). Due to the completely transient nature of the BEM solver, the burst starting force appears in the chart in the early moments of the propeller rotation.
Because each lifting body was separately modeled (Figure 16) in this study, the forces of each blade can be shown in a similar way (Figure 18).
As shown in Figure 18, each propeller force has a harmonic behavior, and this property is similar for opposite blades. Blades No. 1 and 3 or blades No. 2 and 4 are in opposite directions (Figure 16). The frequency of this harmonic behavior is almost the same as the propeller revolution (propellers rotate 4.3 rounds per second; see Table 5). Thus, the shaft bracket imposes a harmonic behavior on the blade force.
In the real selfpropulsion test, at each specified velocity of the model, the propeller speed gradually increased until the net force on the carriage tended to zero. Table 5 shows the respective velocities for the ship and propeller extracted from the experiment. However, in the numerical selfpropulsion (coupled solution in the present research), the fixed propeller speed was applied to the inviscid model, and as a result, the forces of the rotating disk were designated to the RANS model until the vessel's surge velocity converged to a steady value. Accordingly, the time history of balancing forces and ship velocity for n_{p} = 319.33 r/min (Table 5) is shown. As depicted in Figure 19, the total thrust for the twin propellers reached the vessel drag after almost 66 s. Accordingly, as shown in Figure 20, the vessel speed converged at 15.55 kn.
The TSM was applied to capture the unsteady vortex sheet behind the propeller at each time step. Figure 21 shows the wake sheet of the propeller on three different time steps from the top view. Only the wake sheet pattern of blades No.2 and 8 (Figure 16) are shown for a good distinction.
Figure 22 shows the streamline of the left propeller from the top view. The 3D view of the propeller wake sheet pattern is also illustrated in Figure 23.
As a fully unstructured mesh was used, an asymmetrical wake field was produced, which led to different wake sheet patterns of the left and right propellers. Figure 24 shows the ship and propellers and the propeller wake sheet below the water surface.
By applying the solver on the presented test case at different velocities, based on Table 5, the thrust coefficient curve for the propeller behind the appended ship was achieved. Figure 25 compares the open water thrust coefficient with the fully appended ship results.
The results show that the propeller thrust in the open water test does not have much difference from the propeller thrust behind the RV50 hull. In fact, the average difference is $\frac{\mathrm{KT}_{\text {behind ship }}}{\mathrm{KT}_{\text {open water test }}} \cong 1.06$. Equation (31) defines the relative rotative efficiency η_{r} value and its range for rough estimations (Carlton 2018).
$$ 0.96<\eta_{r}=\frac{\text { propeller efficiency behind ship }\left(\eta_{B}\right)}{\text { propeller efficiency in open water test }\left(\eta_{o}\right)}<1.04 $$ (31) To get a good view of the subject, we consider the selfpropulsion test of the RV50, which was performed in Vienna Model Basin with the stock propeller No. 1090a4 as shown in Table 6 (RV50 Test Report 2011).
Number of blades 4 Propeller diameter (mm) 163.79 Hub Diameter (mm) 46.60 blade length at 0.7R(mm) 60.00 propeller pitch (mm) 180.01 The thrust coefficient of propeller No. 1090a4 for the open water and selfpropulsion tests is plotted in Figure 26. As shown, the thrust coefficients for the open water and selfpropulsion tests are almost the same. Therefore, the computational results for a highly skewed propeller seem logical.
This issue should be taken into consideration when the total propulsion power is estimated, especially in the concept design stage.
6 Conclusions
In the conceptual design of a ship propulsion system, there is tremendous demand for an efficient and accurate solver to predict the hydrodynamic performance of a propeller. In this research, a coupled viscous/inviscid solver was presented to calculate the hydrodynamic behavior of the twin propellers behind a fully appended ship. The inviscid solver is based on the potential flow, which is very efficient in calculating lifting body forces, such as propellers. An unsteady BEM was used to simulate the flows around the propeller. The TSM simulates the free wake rollup and its interaction at each time step. The BEM solver was validated by simulating a highly skewed propeller. The numerical results were compared with available data from the experimental open water test. The wake sheet pattern of the propeller when it encountered a uniform flow was also shown.
Moreover, an unsteady NS solver was used in combination with the VOF approach to capture the free surface. The moving mesh technique was also applied to simulate ship motions. To verify the RANS solver, the resistance of an oceanography research vessel in calm water at different velocities was investigated. The results were compared with a towing tank test results to validate the performance of the RANS solver.
Subsequently, the coupled RANS/BEM solver was used to simulate the twin propellers behind the RV50. In particular, the present study is not a real selfpropulsion simulation because the effect of a propeller in disturbing the pressure field around a ship transom was not considered. For this reason, the calculated nominal wake was very low. Indeed, this approach may be called a quasiselfpropulsion because a virtual disk that represents the real propeller was implemented into the RANS solver. However, there is an interaction between the two solvers, as the propeller inflow wake computed in the RANS solver was fed into the unsteady BEM solver.
The proposed solver could simulate propeller thrusts with good accuracy. The solver could also show the wake sheet pattern of a propeller affected by the wake field velocity of the ship.
Nomenclature ϕ(x, t): Velocity potential r Position vector between the source point and field point μ Doublet intensity of wake panel V(t): Linear velocity V_{rel} Relative velocity Re Reynolds number C_{D} Drag coefficient D Propeller diameter J Advanced coefficient K_{Q} Torque coefficient F Hydrodynamic force ρ Fluid density g Gravity ν_{m} Mean perturbation velocity n Unit normal vector of panel pointing outward of the body X_{i} Wake sheet node position Ω(t) Angular velocity P Pressure on panel P_{∞} Farfield pressure C_{p} Pressure coefficient n Propeller rotational speed V_{a} Freestream fluid velocity K_{T} Thrust coefficient M Hydrodynamic moment v Fluid kinematic viscosity 
Figure 3 Flow direction determines donor, acceptor, and upwind cells (Panahi et al. 2006)
Figure 9 Pictures taken from Vienna Model Basin (RV50 Test Report 2011)
Table 1 Dimensions of the highspeed catamaran (Panahi et al. 2009)
Length (m) 12.3 Breadth (m) 4.6 Draught (m) 0.95 Mass (t) 17.850 Vertical center of gravity (m) 0.45 Longitudinal center of gravity (m) 3.81 Block coefficient 0.33 Table 2 Main characteristics of the RV50
Length between perpendiculars (L_{BP}) (m) 44.70 Breadth (m) 10.00 Draught at midship (m) 3.30 Draught at forward perpendicular (m) 3.30 Draught at aft perpendicular (m) 3.30 Displacement mass (t) 896 Block coefficient 0.592 3 Cruise and maximum speed (kn) 12.314.5 Table 3 Comparison of the computational and experimental values of the nominal wake
Nominal wake (w) Error (%) Computational Experimental 0.081 2 0.088 7.73 Table 4 Main particulars of the RV50 propeller
Number of blades 4 Propeller diameter (m) 2.1 Hub Diameter (m) 0.56 Expanded blade area (m^{2}) 1.95 Mean pitch (mm) 1 691.1 Design J = V/nD 0.623 Table 5 Propeller conditions behind the ship
RV50 speed (kn) Propeller rotational speed (r/min) Advanced ratio (J) Time step (s) 7.50 133.21 0.83 0.008 10.50 198.38 0.78 0.005 14.50 (Top) 319.33 0.67 0.003 Table 6 Main particulars of the propeller No. 1090a4
Number of blades 4 Propeller diameter (mm) 163.79 Hub Diameter (mm) 46.60 blade length at 0.7R(mm) 60.00 propeller pitch (mm) 180.01 Nomenclature
ϕ(x, t): Velocity potential r Position vector between the source point and field point μ Doublet intensity of wake panel V(t): Linear velocity V_{rel} Relative velocity Re Reynolds number C_{D} Drag coefficient D Propeller diameter J Advanced coefficient K_{Q} Torque coefficient F Hydrodynamic force ρ Fluid density g Gravity ν_{m} Mean perturbation velocity n Unit normal vector of panel pointing outward of the body X_{i} Wake sheet node position Ω(t) Angular velocity P Pressure on panel P_{∞} Farfield pressure C_{p} Pressure coefficient n Propeller rotational speed V_{a} Freestream fluid velocity K_{T} Thrust coefficient M Hydrodynamic moment v Fluid kinematic viscosity 
Benek J, Steger J, Dougherty FC (1983) A flexible grid embedding technique with application to the Euler equations. 6th Computational Fluid Dynamics Conference, Danvers Carlton J (2018) Marine propellers and propulsion. ButterworthHeinemann Durante D, Dubbioso G, Testa C (2013) Simplified hydrodynamic models for the analysis of marine propellers in a wakefield. Journal of Hydrodynamics 25(6): 954–965 https://doi.org/10.1016/S10016058%2813%2960445X Gaggero S, Dubbioso G, Villa D, Muscari R, Viviani M (2019) Propeller modeling approaches for off—design operative conditions. Ocean Engineering 178: 283–305 https://doi.org/10.1016/j.oceaneng.2019.02.069 Gaggero S, Villa D, Viviani M (2014) An investigation on the discrepancies between RANSE and BEM approaches for the prediction of marine propeller unsteady performances in strongly nonhomogeneous wakes. International Conference on Offshore Mechanics and Arctic Engineering, American Society of Mechanical Engineers Gaggero S, Villa D, Viviani M (2017) An extensive analysis of numerical ship selfpropulsion prediction via a coupled BEM/RANS approach. Applied Ocean Research 66: 55–78 https://doi.org/10.1016/j.apor.2017.05.005 Gaskell P, Lau A (1988) Curvaturecompensated convective transport: SMART, a new boundednesspreserving transport algorithm. International Journal for Numerical Methods in Fluids 8(6): 617–641 https://doi.org/10.1002/fld.1650080602 Greve M, WöcknerKluwe K, AbdelMaksoud M, Rung T (2012) Viscousinviscid coupling methods for advanced marine propeller applications. International Journal of Rotating Machinery 2012: 743060 https://doi.org/10.1155/2012/743060 Henshaw WD, Schwendeman DW (2006) Moving overlapping grids with adaptive mesh refinement for highspeed reactive and nonreactive flow. Journal of Computational Physics 216(2): 744–779 Hess J, Smith A (1962) Calculation of nonlifting potential flow about arbitrary three dimensional bodies. J ship Res 8(4): 22–44 https://doi.org/10.5957/jsr.1964.8.4.22 Hess JL (1972) Calculation of potential flow about arbitrary threedimensional lifting bodies: Douglas Aircraft Co Long Beach CA Hess JL, Valarezo WO (1985) Calculation of steady flow about propellers using a surface panel method. Journal of Propulsion and Power 1(6): 470–476 https://doi.org/10.2514/3.22830 ITTCRecommended Procedures and Guidelines (2011) Performance Prediction Method 7.5020301.4: 9 Jasak H (1996) Error analysis and estimation for the finite volume method with applications to fluid flows. PhD thesis, University of London, London Katz J, Plotkin A (2001) Lowspeed aerodynamics. Cambridge University Press Kim D, Choi H (2000) A secondorder timeaccurate finite volume method for unsteady incompressible flow on hybrid unstructured grids. Journal of Computational Physics 162(2): 411–428 https://doi.org/10.1006/jcph.2000.6546 Leonard B (1991) The ULTIMATE conservative difference scheme applied to unsteady onedimensional advection. Computer Methods in Applied Mechanics and Engineering 88(1): 17–74 https://doi.org/10.1016/00457825%2891%2990232U Meakin R, Suhs N (1989) Unsteady aerodynamic simulation of multiple bodies in relative motion. 9th Computational Fluid Dynamics Conference Najafi S, Abbaspoor M (2019) Numerical investigation of flow pattern and hydrodynamic forces of submerged marine propellers using unsteady boundary element method. Proceedings of the Institution of Mechanical Engineers, Part M: Journal of Engineering for the Maritime Environment 233(1): 67–79 Najafi S, Abbaspour M (2017) Numerical study of propulsion performance in swimming fish using boundary element method. Journal of the Brazilian Society of Mechanical Sciences and Engineering 39(2): 443–455 https://doi.org/10.1007/s4043001606138 Panahi R, Jahanbakhsh E, Seif MS (2006) Development of a VoFfractional step solver for floating body motion simulation. Applied Ocean Research 28(3): 171–181 https://doi.org/10.1016/j.apor.2006.08.004 Panahi R, Jahanbakhsh E, Seif MS (2009) Towards simulation of 3D nonlinear highspeed vessels motion. Ocean Engineering 36(3–4): 256–265 https://doi.org/10.1016/j.oceaneng.2008.11.005 Politis G (2016) Unsteady wake rollup modeling using a mollifier based filtering technique. Dev. Appl. Ocean. Eng 5(1): 1–28 Politis GK (2004) Simulation of unsteady motion of a propeller in a fluid including free wake modeling. Engineering Analysis with Boundary Elements 28(6): 633–653 https://doi.org/10.1016/j.enganabound.2003.10.004 Politis GK (2005) Unsteady rollup modeling for wakeadapted propellers using a timestepping method. Journal of Ship Research 49(3): 216–231 https://doi.org/10.5957/jsr.2005.49.3.216 Politis GK (2011) Application of a BEM time stepping algorithm in understanding complex unsteady propulsion hydrodynamic phenomena. Ocean Engineering 38(4): 699–711 https://doi.org/10.1016/j.oceaneng.2011.01.001 Pourmostafa M, Ghadimi P (2020a) Applying boundary element method to simulate a highskew Controllable Pitch Propeller with different hub diameters for preliminary design purposes. Cogent Engineering 7(1): 1805857 https://doi.org/10.1080/23311916.2020.1805857 Pourmostafa M, Ghadimi P (2020b) Unsteady 2D and 3D NavierStokes Solver with application of multigrid scheme to pressure poisson fractional step on arbitrary unstructured grids in various applications with emphasis on ship motion. Mathematical Problems in Engineering 2020: 1–28 https://doi.org/10.1155/2020/3248958 Queutey P, Deng GB, Guilmineau E, Salvatore F (2013) A comparison between full RANSE and coupled RANSEBEM approaches in ship propulsion performance prediction. International Conference on Offshore Mechanics and Arctic Engineering, OMAE 2013–10366 Rao, ZQ, Yang CJ (2017) Numerical prediction of effective wake field for a submarine based on a hybrid approach and an RBF interpolation. Journal of Hydrodynamics 29(4): 691–701 https://doi.org/10.1016/S10016058%2816%29607813 Rijpkema D, Starke B, Bosschers J (2013) Numerical simulation of propellerhull interaction and determination of the effective wake field using a hybrid RANSBEM approach. 3rd International Symposium on Marine Propulsors, Launceston, Australia RV50 Test Report (2011) Hydrodynamic model tests. Test Report No. 2512/01, Model No. 2512, Vienna Model Basin (www. sva. at): 43 Ubbink O, Issa R (1999) A method for capturing sharp fluid interfaces on arbitrary meshes. Journal of Computational Physics 153(1): 26–50 https://doi.org/10.1006/jcph.1999.6276 Wang Y, AbdelMaksoud M, Song B (2016) Convergence of different wake alignment methods in a panel code for steadystate flows. Journal of Marine Science and Technology 21(4): 567–578 https://doi.org/10.1007/s0077301603750 Zhu Q, Wolfgang M, Yue D, Triantafyllou M (2002) Threedimensional flow structures and vorticity control in fishlike swimming. Journal of Fluid Mechanics 468: 1–28 https://doi.org/10.1017/S002211200200143X