# Viscous-Inviscid Coupling Methods for Advanced Marine Propeller Applications

Authors: Martin Greve; Katja Wöckner-Kluwe; Moustafa Abdel-Maksoud; Thomas Rung
Journal: International Journal of Rotating Machinery
 (2012)

---

## Abstract

The paper reports the development of coupling strategies between an inviscid direct panel method and a viscous RANS method and their application to complex propeller ows. The work is motivated by the prohibitive computational cost associated to unsteady viscous flow simulations using geometrically resolved propellers to analyse the dynamics of ships in seaways. The present effort aims to combine the advantages of the two baseline methods in order to reduce the numerical effort without compromising the predictive accuracy. Accordingly, the viscous method is used to calculate the global flow field, while the inviscid method predicts the forces acting on the propeller. The corresponding reaction forces are employed as body forces to mimic the propeller influence on the viscous flow field. Examples included refer to simple verification cases for an isolated propeller blade, open-water validation simulations for a complete propeller, and more challenging investigations of a manoeuvring vessel in seaways. Reported results reveal a fair predictive agreement between the coupled approach and fully viscous simulations and display the efficiency of the coupled approach.


---

## Body

## 1. Introduction

The design of propellers is usually confined to the analysis of open-water, deeply submerged, and steady-state conditions. When attention is directed to the dynamics of ships operating in harsh conditions, for example, offshore supply vessels, a more sophisticated investigation might be necessary to adequately determine the propeller loads. Offshore supply vessels are frequently employed with pod propulsors and feature a pronounced influence of the ship, the pod housing, and the free surface on the blade loads when the ship is exposed to heavy seas. Such complex phenomena can only be addressed by (at least) Reynolds-averaged Navier-Stokes (RANS) simulation methods for turbulent flows, which are fairly expensive due to the wide range of involved time scales. Hence, only few examples associated to propeller operations behind a ship (e.g., Carrica et al. [1]) or off-design conditions (e.g., Sileo and Steen [2]) were reported for viscous flow simulations. Furthermore, multiple domain approaches using sliding or overlapping grid techniques (see, e.g., Muscari and Di Mascio [3]) are necessary to simulate the relative motions between propeller and the hull. On the contrary, methods based on the potential-flow theory are not able to predict the frictional contributions to the wake flow and cannot therefore be used for the simulation of propeller flows in behind conditions.A viable solution to reduce the computational effort without compromising the predictive accuracy is to couple a potential-flow method for the analysis of the propeller with a RANS method for the simulation of the flow around the vessel in seaways.The basic coupling principle is fairly mature, thus only few examples are mentioned here. First approaches to model propeller effects by means of a force field in an inviscid flow have been described as early as 1972 by Sparenberg [4, 5], who aimed to find the optimum actuator disc for a given thrust. Schetz and Favin [6, 7] reported attempts to model propeller forces in an axisymmetric flow field employing a one-way coupling to the Navier-Stokes equations in a stream-function formulation. A more general two-way coupling was firstly reported by Stern et al. in 1988 [8]. Computations were performed in steady-state based on circumferential averages of the propeller forces. The propeller forces were determined by a vortex-lattice lifting-surface method and adjusted to the flow field upstreaming the propeller in an iterative procedure.Several two-way coupling approaches with different levels of complexity have recently been developed. Typically the ship and the appendages are resolved (i.e., geometrically modeled) by the RANS mesh, and the effect of the propeller is computed externally to limit the computational and mesh generation effort. In most cases, steady potential-flow methods or even more simple assumptions, for example, prescribed force parameterisations [9] are used to determine the propeller influence.An example for such an approach was presented by Phillips et al. [10]. The authors describe a coupling procedure using the blade element momentum theory (BEMT) to model the propeller forces. These forces are applied in a RANS domain with a ship hull and a rudder included. An averaged RANS-calculated nominal wake fraction evaluated at prescribed radial positions serves as input to the BEMT. The calculated thrust and torque values obtained from the BEMT are transformed into radius-averaged momentum source terms acting in a cylindrical region at the propeller position inside the RANS domain. With these body forces, the steady RANS simulation is iterated to convergence. The authors employ the method to compute the forces acting on a rudder with fair accuracy.Another possibility to model the propeller forces in a RANS computation is described by Müller et al. [11]. The authors compare the performance of a coefficient-based force model with a geometrically resolved propeller. The propeller influence is not computed from an inviscid numerical analysis but modeled based on established force coefficients evaluated in different sections of the blade. They studied the influence of the propeller representation on the velocity distribution in the aft ship region as well as the associated computational and grid-generation efforts. Good accuracy was observed for the velocity distributions in conjunction with a significantly smaller computational effort for the simulations with the propeller model due to the smaller and nonsliding grids. In the approach of Müller, a simplified velocity distribution located at a short-distance upstream, the propeller serves as input for the coefficient model. The thrust and torque contributions are determined for prescribed radial blade sections. The propeller model thus considers the variation of the thrust and torque in radial and circumferential direction but not the actual blade positions.The present paper reports the development and application of a coupling procedure between a RANS method and a boundary element method (BEM) with respect to pronounced unsteady behavior of the propeller thrust and torque due to large amplitude ship motions. To achieve this the time-accurate propeller loads computed by the inviscid method are applied to the RANS method without any circumferential or radial averaging. Additionally, no rotating meshes are required to simulate the propeller and significantly smaller grids for the RANS simulation can be used—a feature which is of particular relevance for unsteady simulations involving multiple time scales.The remainder of the paper is structured as follows. The second and third sections are devoted to the viscous and inviscid baseline methods. The fourth section outlines the rationale of the coupling strategy, while Section5 refers to grid studies for both involved solvers and to a verification study based on a comparison of the propeller-induced velocities.Subsequently, a validation of the coupling method for a free-running propeller and the application for an investigation of propeller-hull interaction effects in seaways is described. The objective of this study is to compute the speed losses and the additional resistance for a maneuvering vessel in waves. The paper ends with an outline of the conclusions and future aims in Section 7.


## 2. Viscous Method

The coupling approach introduced in the present paper uses the finite-volume Navier-Stokes procedureFreSCo+ as viscous baseline method. The tool is a spin-off of the FreSCo suite [12], a joint development of Hamburg University of Technology (TUHH), the Hamburg ship model basin (HSVA), and the Maritime Research Institute Netherlands (MARIN), which was developed under the aegis of the FP6 EU-initiative VIRTUE [13]. The simulation package is supplemented by an adjoint flow solver AD-FreSCo+, dedicated to shape optimization [14], an overset-grid technique OVER-FreSCo+ [15] and a dynamic grid-adaptation approach [16]. The algorithm solves the Navier-Stokes equations for an incompressible fluid composed of two (or more) immiscible phases (i.e., air and water), namely,
(1)ρDUiDt=-∂p∂xi+∂∂xj[(μ+μt)(∂Ui∂xj+∂Uj∂xi)]+fi,∂ρ∂t+∂(ρUi)∂xi=0,
with
(2)ρ=cρA+(1-c)ρW,μ=μA+(1-c)μW,
where Ui, xi, and fi represent the Cartesian coordinates of the velocity vector, the spatial coordinates and the body-force vector, and t, p, ρ, μ, μt denote to the time, pressure, fluid density, molecular and turbulent dynamic viscosity. The property c refers to the local volumetric air concentration (a.k.a. mixture fraction), the suffix A refers to air, and the suffix W refers to water.The segregated algorithm is based on the strong conservation form of the momentum equations and employs a cell-centered, colocated storage arrangement for all transport properties. Structured and unstructured grids with arbitrary polyhedral cells or hanging nodes can be used. The implicit numerical approximation is of second-order accuracy in space and time. Integrals are approximated using the conventional mid-point rule. The solution is iterated to convergence using a SIMPLE-type pressure-correction scheme. Various turbulence-closure models using statistical (RANS) or scale-resolving (LES and DES) approaches are available. Two-phase flows are modeled by interface-capturing methods based upon the level-set or volume-of-fluid (VOF) technique. To obtain the solution for the linear equation systems, methods offered by the PETSc library are employed. The algorithm is parallelized using a domain-decomposition technique based on a single program multiple data (SPMD) message-passing model. Interprocessor communication employs the MPI communications protocol. Load balancing is achieved using the ParMETIS partitioning software.


## 3. Inviscid Method

The potential-flow methodpanMARE [17] is developed at the Hamburg University of Technology. It uses a BEM to solve Laplace's equation: ∇Φ*=0, describing a potential flow field. The total velocity potential Φ*=Φ∞+Φ consists of the potential due to the inflow velocity Φ∞ and the perturbation potential Φ due to the presence of the body. The latter is modeled using a linear combination of sources and doublets distributed on the surface of all physical boundaries present in the simulated fluid domain. When applying Green's third identity to Laplace's equation, the velocity potential at a point x in the fluid domain becomes
(3)ϕ(x)=-14Π(∑j=1NB‍σj∫Aj‍1rj(x)ds∑j=1NB‍σj∫Aj‍1rj(x)+∑j=1NB+NW‍μj∫Aj‍nj∇1rj(x)ds),
where NB denotes the number of quadrilateral panels on the body surface and NW the number of panels on the trailing wake surface of the propeller blade. The source strength of panel j is σj, and the dipole strength is μj. Aj is the panels area, and nj is its normal vector pointing into the flow. The distance from the evaluation point x to the panel j is rj(x).The Neumann boundary condition requires the velocity normal to a boundary surface to be equal to zero (∇Φ*·n=0). The source strength is calculated based on the inflow velocity:
(4)σ=-nv∞.Further, the dipole strength of the trailing wake panels is defined by the linear Kutta condition. It is applied on the panels directly connected to the trailing edge and the dipole strength of the first panel in the trailing wake of the blade (μTE) becomes
(5)μTE=μupper-μlower,
with the dipole strengths of the directly attached panels of the body surface on the suction (“upper”) and pressure side (“lower”) of the propeller blade. The panels in the propeller wake sheet are deformed to align with the local flow.Using the above-described boundary conditions the equation system with the unknown dipole strengths can be solved. The panel pressure is determined by evaluating Bernoulli's equation at the collocation pointx:
(6)-p(x)ρ=∇Φ(x)2+∇Φ∞∇Φ(x)+∂μ(x)∂t.To take into account viscous effects, an empirical friction-correction formulation is implemented. For calculating the forces acting on a propeller behind a ship arbitrary distributions of inflow velocities (e.g., nonaxisymmetric) can be used as inflow condition.


## 4. Coupling Strategy

The present approach models the propeller by a set of equivalent body forces embedded in the viscous flow. The body forces are derived from the propeller loads computed by the inviscid method which in turn requires a prescribed velocity distribution upstream the propeller. The strategy is based on the exchange of two properties, that is, the transfer of inviscid propeller-induced body forces to the viscous solver and the transfer of the wake velocities to the inviscid solver. As regards the implementation, two details of the coupling procedure related to the (a) spatial and (b) temporal transfer of information are of specific interest. The former issue is associated to mapping techniques between two grids of substantially different nature. For the velocity transfer a plane which is typically located 0.5R–1.0R upstream of the propeller—where R denotes the propeller radius—is defined. Inside this plane, the velocities computed by the RANS method are mapped onto a predefined set of locations in a circular domain (specified by corresponding angular and radial positions) as depicted in Figure 1.Figure 1
Exemplary distribution of locations used to transfer velocities from the viscous to the inviscid method.The spatial transfer of body forces from the potential flow solver to the viscous solver follows a different route. Since one force per BEM panel is transferred, the forces of the pressure- and suction-side of the blade are transferred separately. These forces are then distributed on the adjacent fluid cell centers of the RANS mesh as volumetric body forces (Figure2). Although the grids of the viscous and the inviscid method are significantly different, in particular with respect to the grid density, the mapping procedure maintains force conservation. It utilizes a cell-search algorithm [15] that connects each panel to the closest cell center and the corresponding neighboring cells in a conservative way.Figure 2
Illustration of RANS mesh cells containing propeller blade body forces (approximately 2500 cells are involved).For the temporal transfer of information, three different techniques are developed. The first refers to an explicit strategy, where one data transfer per time-step is performed, while in the implicit mode, the respective data is transferred several times per time-step as illustrated in the flowchart in Figure3.Figure 3
Flowchart of the employed coupling algorithm.At the beginning of each coupling cycle, the viscous velocity distribution upstream of the propeller plane is transferred to the inviscid method. In the explicit mode, these viscous wake velocities are based on the informations of the previous time-step. The wake field computed by the viscous solver inheres the influence of the propeller by means of the induced velocities. As the inviscid solver requires the effective wake as inflow condition, the propeller induced velocities have to be subtracted from the viscous wake. Hence, in each coupling cycle the inviscid method starts with a calculation of these induced velocities to obtain the effective wake. The algorithm then proceeds with the solution of the equation system insidepanMARE to determine the propeller forces. The force distribution is transferred to the viscous method, and in the explicit mode the time-step is finalized by the solution of the governing equations in FreSCo+.In the implicit mode, multiple coupling cycles are performed in each time-step. Here, the inviscid solution is embedded into the outer iteration process of the segregated viscous solver. For each time-step, a sequence of outer iterations is performed in line with the description for the explicit mode. An updated time-accurate viscous velocity distribution is transferred to the BEM in every outer iteration ofFreSCo+. The effective and induced wake velocities employed by panMARE are, therefore, also determined in a time-accurate manner. The solution advances in time after a converged solution is achieved in both solvers.The third mode refers to a semiexplicit technique, where the global coupling procedure follows the explicit mode, but the propeller induced velocities are iterated to convergence in each coupling cycle. Using the inviscid solution of the previous time-step as an initial guess for the induced velocities, the latter are updated for a sequence of new inviscid solutions. The respective effective wake field is thus composed from the former viscous wake field and the present induced velocities.


## 5. Verification

This section is devoted to a comparison between the implicit, explicit, and semi-implicit coupling approach as described in Section4. Emphasis is given to the attainable computational efficiency by means of the time-step size and the upstream position of the velocity transfer plane. The results are primarily analyzed in terms of the induced velocities computed for an exemplary one-blade propeller. Comparison is made between the induced velocities derived from the velocity field in the RANS-domain and the induced velocities computed by panMARE. Both are obtained from one coupled simulation for a deeply submerged propeller. For this case, two grid studies are conducted beforehand, where the first refers to open water calculations using only the BEM and the second to different RANS-grids in a coupled open-water simulation. Finally, results for the corresponding four-blade propeller and a comparison between a coupled and a viscous simulation with a geometrically modelled propeller are presented. The investigated propeller and the applied coordinate system can be found in Figure 4. More details about the test case can be found in Section 6.1.Figure 4
Geometry of investigated propeller.


### 5.1. Grid Study

#### 5.1.1. BEM Grid

Inside the BEM, the blade geometry is discretised withn×m panels, where n refers to the number of panels in radial direction and m to the number in circumferential direction of the blade section. The hub geometry is not modeled. Five different grids are used, and the deviation of the calculated kt value to the value for the finest grid is shown in Table 1.Table 1
Deviation ofkt to value of finest grid for the BEM grid study.
J
8
×
18
16
×
36
24
×
54
32
×
72
40
×
90
0.5
+4.7%
+2.0%
+0.8%
+0.2%
0.351
0.7
+2.2%
+0.8%
+0.0%
+0.0%
0.290
0.9
−1.5%
−0.8%
−0.6%
−0.2%
0.220With increasing number of panels in the radial direction, the shape of the trailing wake surfaces especially the blades tip vortex shape can be determined more accurately.The grid with 16×36 panels is chosen for further investigations.


#### 5.1.2. RANSE Grid

With the above-selected BEM-grid a grid study for the RANS mesh in a coupled simulation is conducted. The box-shaped grid used for the open-water coupling process does not inherit any geometry. Different refinement areas are generated inside the domain, and the grid consists of unstructured hexahedral volume cells. One refinement box encloses the propeller area and two others the propeller wake area. The plane, where the velocity field is passed to the BEM, is also within the refinement boxes. Due to the discretisation, the propeller forces are distributed on a different number of grid cells resulting in dissimilar propeller force distribution and thus also dissimilar propeller inflow. Therefor, the calculated thrust and torque values differ between the grids. The grids generated and the predictedkt values for the advance ratio J=0.5 are given in Table 2.Table 2
Deviation of blade thrust to value of finest grid for the RANS grid study (J=0.5).
grid
27 k
65 k
152 k
300 k
662 k
Thrust (N)
71.94
71.92
72.48
72.56
72.55
Deviation (%)
−0.9
−0.9
−0.1
0.0
—The value range for the computed thrust is only about1% between the coarsest and the finest grid. This result is encouraging because a coarse grid is sufficient to capture the relevant effects related to transferring the velocities to the BEM and distributing the forces in the RANS grid. Of course, when the focus lies on determining the transient forces on any kind of structure (e.g., a ship or a thruster) the grid density will have to be investigated further.For the following calculations, the 300 k cells grid is chosen because the computed thrust values match perfectly with the finest grid, and the calculation time is still short.


### 5.2. Investigation of Induced Velocities

The subsequent analysis is based on the induced velocities at a fixed evaluation radiusr=0.7R. The propeller mostly operates at an advance ratio of J=0.5. When the blade is at the topmost position, the values of the induced velocities are depicted for all angular positions θ between 0° and 360°. At the evaluation radius, the leading edge of the blade is located at approximately 340° and the trailing edge at approximately 20° as indicated by vertical lines in Figure 5. The induced velocities are evaluated in a plane at a distance dx=0.5R in front of the propeller. The employed time-step refers to an angular rotation of 3.78°. Results are displayed by means of axial, radial, and tangential velocity components. The values of the induced velocities are nondimensionalised by the inflow velocity.Figure 5
Normalised induced velocities in cylindrical coordinates for a single-blade propeller (explicit coupling:r/R=0.7, J=0.5, dt·n=3.78°, and dx=0.5R).According to the coordinate system, negative axial-induced velocities denote for an acceleration of the flow. In Figure5, the axial velocity values computed by the BEM show a steep increase until they reach a maximum in the region directly affected by the blade's trailing edge at approximately 30°. After this position, the induced axial velocity decreases gradually until reaching a minimum near the position of the leading edge. The induced radial velocity is directed inwards and shows a small variation over the circumferential direction. The pronounced maximum of the tangentially induced velocities is directed in the direction of rotation and is found at approximately one-third of the chord length. For the following analysis, only the axial component of the induced velocity vector is used.When comparing the different coupling algorithms in Figures6 and 7, only the results for the upper half of the angular positions are displayed for the sake of more clearance. It is obvious that the results of the different coupling modes are quite similar. Dissimilarities can be found with respect to the steep increase of induced axial velocities between the leading and the trailing edge. The implicit algorithm shows a slightly earlier and more accentuated increase of the velocities. This can also be found for the semi-implicit algorithm but less distinctly visible. The finding is in line with the small inherent time lag of nonimplicit couplings. Both baseline methods display a good agreement of the predicted maxima, while more pronounced minima are predicted by the inviscid BEM. Apart from this, the values do not differ significantly between the coupling approaches. For this quite simple case, the explicit approach seems to be sufficient to capture the relevant effects. Henceforth, this variant of the coupling algorithm is used for the following analysis and the simulations in Section 6, where due to the pulling pod arrangement no significant ship wakefield can be observed.Figure 6
Comparison of the normalised induced axial velocities obtained from the explicit and the implicit coupling approach (r/R=0.7, J=0.5, dt·n=3.78°, and dx=0.5R).Figure 7
Comparison of the normalised induced axial velocities obtained from the explicit and the semi-implicit coupling approach (r/R=0.7, J=0.5, dt·n=3.78°, and dx=0.5R).In Figure8, a comparison between the normalised induced velocities in the axial direction for three different time-step sizes is given. The medium time-step was also applied in the previous investigations. Again attention is restricted to the upper half of the angular positions. It can be observed that the steepness of the increase in front of the maximum value is enhanced with decreasing time-step size. Because the overall agreement is fair, the medium time-step is chosen for further investigations.Figure 8
Comparison of the normalised induced axial velocities for different time-step sizes (explicit coupling:r/R=0.7, J=0.5, and dx=0.5R).When analysing the induced velocities for different upstream distances between the blade and the transfer plane of the velocity field, a clear dependency becomes obvious. The upstream influence of the propeller is locally limited and the transfer plane should be as close to the propeller as possible (Figure9). For the larger distances the correlation between the solvers is still sufficient, but the total value of induced axial velocity decreases significantly and reduces to about one-third when evaluating at a distance of dx=1R. This implies that the induced velocities should be determined as close to the propeller position as possible. However, when the distance is reduced to dx≤0.5R, the plane collides with the leading edge of the blade inside the BEM. The smallest distance shown in Figure 9 is thus employed for the other investigations presented here. The induced velocities computed by the BEM are slightly smaller than the ones computed by RANS, which results in an overpredicted inflow velocity to the propeller. This also accounts for the induced velocities at radii other than r/R=0.7 and depends on the local distance to the blade and the hydrodynamic angle of attack of the profile section (see Figures 9 and 10).Figure 9
Comparison of the normalised induced axial velocities for different axial distances between the velocity transfer plane and the propeller (explicit coupling:r/R=0.7, J=0.5, and dt·n=3.78°).Figure 10
Comparison of the normalised induced axial velocities for different propeller advance ratios explicit coupling: (r/R=0.7, dt·n=3.78°, and dx=0.5R).An investigation of the induced velocities for different propeller loadings can be found in Figure10. The characteristics of the curves are similar, but the absolute difference between the viscous and inviscid computed induced velocities is larger for lower advance ratios. Under consideration of the difficulties in simulating highly loaded propellers with inviscid methods, this result still seems to be sufficient for the purpose of the coupling algorithm.Figure11 shows the three components of the induced velocities for a four-blade propeller at J = 0.5. Also the correlation between the solvers for the radial and tangential components is very good and illustrates that the propeller faces equivalent working conditions in the involved baseline methods. Finally, the induced velocities computed during a coupled simulation and a RANS simulation with geometrically modelled propeller are compared for J=0.5 in Figure 12. Therein, the values predicted by the coupling approach are derived from the viscous computed velocity field. Although the grid size differs significantly between the simulations, the characteristics of the curves are quite similar, but the peak values of the RANS simulation cannot be fully captured by the coupled simulation.Figure 11
Comparison of normalised induced velocities for a four-bladed propeller (explicit coupling:r/R=0.7, J=0.5, dt·n=3.78°, dx=0.5R).Figure 12
Comparison between the normalised induced velocities obtained from a coupled (explicit mode) and a RANS simulation with a geometrically modeled propeller (r/R=0.7, J=0.5, dt·n=3.78°, and dx=0.5R).


## 5.1. Grid Study

### 5.1.1. BEM Grid

Inside the BEM, the blade geometry is discretised withn×m panels, where n refers to the number of panels in radial direction and m to the number in circumferential direction of the blade section. The hub geometry is not modeled. Five different grids are used, and the deviation of the calculated kt value to the value for the finest grid is shown in Table 1.Table 1
Deviation ofkt to value of finest grid for the BEM grid study.
J
8
×
18
16
×
36
24
×
54
32
×
72
40
×
90
0.5
+4.7%
+2.0%
+0.8%
+0.2%
0.351
0.7
+2.2%
+0.8%
+0.0%
+0.0%
0.290
0.9
−1.5%
−0.8%
−0.6%
−0.2%
0.220With increasing number of panels in the radial direction, the shape of the trailing wake surfaces especially the blades tip vortex shape can be determined more accurately.The grid with 16×36 panels is chosen for further investigations.


### 5.1.2. RANSE Grid

With the above-selected BEM-grid a grid study for the RANS mesh in a coupled simulation is conducted. The box-shaped grid used for the open-water coupling process does not inherit any geometry. Different refinement areas are generated inside the domain, and the grid consists of unstructured hexahedral volume cells. One refinement box encloses the propeller area and two others the propeller wake area. The plane, where the velocity field is passed to the BEM, is also within the refinement boxes. Due to the discretisation, the propeller forces are distributed on a different number of grid cells resulting in dissimilar propeller force distribution and thus also dissimilar propeller inflow. Therefor, the calculated thrust and torque values differ between the grids. The grids generated and the predictedkt values for the advance ratio J=0.5 are given in Table 2.Table 2
Deviation of blade thrust to value of finest grid for the RANS grid study (J=0.5).
grid
27 k
65 k
152 k
300 k
662 k
Thrust (N)
71.94
71.92
72.48
72.56
72.55
Deviation (%)
−0.9
−0.9
−0.1
0.0
—The value range for the computed thrust is only about1% between the coarsest and the finest grid. This result is encouraging because a coarse grid is sufficient to capture the relevant effects related to transferring the velocities to the BEM and distributing the forces in the RANS grid. Of course, when the focus lies on determining the transient forces on any kind of structure (e.g., a ship or a thruster) the grid density will have to be investigated further.For the following calculations, the 300 k cells grid is chosen because the computed thrust values match perfectly with the finest grid, and the calculation time is still short.


## 5.1.1. BEM Grid

Inside the BEM, the blade geometry is discretised withn×m panels, where n refers to the number of panels in radial direction and m to the number in circumferential direction of the blade section. The hub geometry is not modeled. Five different grids are used, and the deviation of the calculated kt value to the value for the finest grid is shown in Table 1.Table 1
Deviation ofkt to value of finest grid for the BEM grid study.
J
8
×
18
16
×
36
24
×
54
32
×
72
40
×
90
0.5
+4.7%
+2.0%
+0.8%
+0.2%
0.351
0.7
+2.2%
+0.8%
+0.0%
+0.0%
0.290
0.9
−1.5%
−0.8%
−0.6%
−0.2%
0.220With increasing number of panels in the radial direction, the shape of the trailing wake surfaces especially the blades tip vortex shape can be determined more accurately.The grid with 16×36 panels is chosen for further investigations.


## 5.1.2. RANSE Grid

With the above-selected BEM-grid a grid study for the RANS mesh in a coupled simulation is conducted. The box-shaped grid used for the open-water coupling process does not inherit any geometry. Different refinement areas are generated inside the domain, and the grid consists of unstructured hexahedral volume cells. One refinement box encloses the propeller area and two others the propeller wake area. The plane, where the velocity field is passed to the BEM, is also within the refinement boxes. Due to the discretisation, the propeller forces are distributed on a different number of grid cells resulting in dissimilar propeller force distribution and thus also dissimilar propeller inflow. Therefor, the calculated thrust and torque values differ between the grids. The grids generated and the predictedkt values for the advance ratio J=0.5 are given in Table 2.Table 2
Deviation of blade thrust to value of finest grid for the RANS grid study (J=0.5).
grid
27 k
65 k
152 k
300 k
662 k
Thrust (N)
71.94
71.92
72.48
72.56
72.55
Deviation (%)
−0.9
−0.9
−0.1
0.0
—The value range for the computed thrust is only about1% between the coarsest and the finest grid. This result is encouraging because a coarse grid is sufficient to capture the relevant effects related to transferring the velocities to the BEM and distributing the forces in the RANS grid. Of course, when the focus lies on determining the transient forces on any kind of structure (e.g., a ship or a thruster) the grid density will have to be investigated further.For the following calculations, the 300 k cells grid is chosen because the computed thrust values match perfectly with the finest grid, and the calculation time is still short.


## 5.2. Investigation of Induced Velocities

The subsequent analysis is based on the induced velocities at a fixed evaluation radiusr=0.7R. The propeller mostly operates at an advance ratio of J=0.5. When the blade is at the topmost position, the values of the induced velocities are depicted for all angular positions θ between 0° and 360°. At the evaluation radius, the leading edge of the blade is located at approximately 340° and the trailing edge at approximately 20° as indicated by vertical lines in Figure 5. The induced velocities are evaluated in a plane at a distance dx=0.5R in front of the propeller. The employed time-step refers to an angular rotation of 3.78°. Results are displayed by means of axial, radial, and tangential velocity components. The values of the induced velocities are nondimensionalised by the inflow velocity.Figure 5
Normalised induced velocities in cylindrical coordinates for a single-blade propeller (explicit coupling:r/R=0.7, J=0.5, dt·n=3.78°, and dx=0.5R).According to the coordinate system, negative axial-induced velocities denote for an acceleration of the flow. In Figure5, the axial velocity values computed by the BEM show a steep increase until they reach a maximum in the region directly affected by the blade's trailing edge at approximately 30°. After this position, the induced axial velocity decreases gradually until reaching a minimum near the position of the leading edge. The induced radial velocity is directed inwards and shows a small variation over the circumferential direction. The pronounced maximum of the tangentially induced velocities is directed in the direction of rotation and is found at approximately one-third of the chord length. For the following analysis, only the axial component of the induced velocity vector is used.When comparing the different coupling algorithms in Figures6 and 7, only the results for the upper half of the angular positions are displayed for the sake of more clearance. It is obvious that the results of the different coupling modes are quite similar. Dissimilarities can be found with respect to the steep increase of induced axial velocities between the leading and the trailing edge. The implicit algorithm shows a slightly earlier and more accentuated increase of the velocities. This can also be found for the semi-implicit algorithm but less distinctly visible. The finding is in line with the small inherent time lag of nonimplicit couplings. Both baseline methods display a good agreement of the predicted maxima, while more pronounced minima are predicted by the inviscid BEM. Apart from this, the values do not differ significantly between the coupling approaches. For this quite simple case, the explicit approach seems to be sufficient to capture the relevant effects. Henceforth, this variant of the coupling algorithm is used for the following analysis and the simulations in Section 6, where due to the pulling pod arrangement no significant ship wakefield can be observed.Figure 6
Comparison of the normalised induced axial velocities obtained from the explicit and the implicit coupling approach (r/R=0.7, J=0.5, dt·n=3.78°, and dx=0.5R).Figure 7
Comparison of the normalised induced axial velocities obtained from the explicit and the semi-implicit coupling approach (r/R=0.7, J=0.5, dt·n=3.78°, and dx=0.5R).In Figure8, a comparison between the normalised induced velocities in the axial direction for three different time-step sizes is given. The medium time-step was also applied in the previous investigations. Again attention is restricted to the upper half of the angular positions. It can be observed that the steepness of the increase in front of the maximum value is enhanced with decreasing time-step size. Because the overall agreement is fair, the medium time-step is chosen for further investigations.Figure 8
Comparison of the normalised induced axial velocities for different time-step sizes (explicit coupling:r/R=0.7, J=0.5, and dx=0.5R).When analysing the induced velocities for different upstream distances between the blade and the transfer plane of the velocity field, a clear dependency becomes obvious. The upstream influence of the propeller is locally limited and the transfer plane should be as close to the propeller as possible (Figure9). For the larger distances the correlation between the solvers is still sufficient, but the total value of induced axial velocity decreases significantly and reduces to about one-third when evaluating at a distance of dx=1R. This implies that the induced velocities should be determined as close to the propeller position as possible. However, when the distance is reduced to dx≤0.5R, the plane collides with the leading edge of the blade inside the BEM. The smallest distance shown in Figure 9 is thus employed for the other investigations presented here. The induced velocities computed by the BEM are slightly smaller than the ones computed by RANS, which results in an overpredicted inflow velocity to the propeller. This also accounts for the induced velocities at radii other than r/R=0.7 and depends on the local distance to the blade and the hydrodynamic angle of attack of the profile section (see Figures 9 and 10).Figure 9
Comparison of the normalised induced axial velocities for different axial distances between the velocity transfer plane and the propeller (explicit coupling:r/R=0.7, J=0.5, and dt·n=3.78°).Figure 10
Comparison of the normalised induced axial velocities for different propeller advance ratios explicit coupling: (r/R=0.7, dt·n=3.78°, and dx=0.5R).An investigation of the induced velocities for different propeller loadings can be found in Figure10. The characteristics of the curves are similar, but the absolute difference between the viscous and inviscid computed induced velocities is larger for lower advance ratios. Under consideration of the difficulties in simulating highly loaded propellers with inviscid methods, this result still seems to be sufficient for the purpose of the coupling algorithm.Figure11 shows the three components of the induced velocities for a four-blade propeller at J = 0.5. Also the correlation between the solvers for the radial and tangential components is very good and illustrates that the propeller faces equivalent working conditions in the involved baseline methods. Finally, the induced velocities computed during a coupled simulation and a RANS simulation with geometrically modelled propeller are compared for J=0.5 in Figure 12. Therein, the values predicted by the coupling approach are derived from the viscous computed velocity field. Although the grid size differs significantly between the simulations, the characteristics of the curves are quite similar, but the peak values of the RANS simulation cannot be fully captured by the coupled simulation.Figure 11
Comparison of normalised induced velocities for a four-bladed propeller (explicit coupling:r/R=0.7, J=0.5, dt·n=3.78°, dx=0.5R).Figure 12
Comparison between the normalised induced velocities obtained from a coupled (explicit mode) and a RANS simulation with a geometrically modeled propeller (r/R=0.7, J=0.5, dt·n=3.78°, and dx=0.5R).


## 6. Validation and Application

### 6.1. Propeller in Open-Water Condition

The first application example refers to open-water calculations of a propeller, where model-scale experiments are performed by MARINTEK under the aegis of PROPSEAS (Kozlowska et al. [18]). The investigated propeller has a diameter of D = 0.2 m and operates with n = 18 rps. Different advance coefficients are achieved by a variation of the inflow velocity. The numerical grid of the propeller used in the BEM consists of 576 panels per blade (cf. Figure 13(a)). Figure 13(b) displays the propeller regime of the viscous grid. About 2.500 of 15.000 control volumes in the propeller disc are subjected to body forces. The Reynolds number based on the chord length at r/R=0.7, and the resulting velocity at the propeller blade reads Re=0.5·106. The flow is considered fully turbulent, and the standard k-ω turbulence model is applied.Numerical grid and resulting pressure distribution for the BEM simulation (a). FV grid in the propeller disc (b). Only cells at the actual blade positions contain body forces (colored).
(a)
(b)The resulting open-water diagram of the test case can be found in Figure14. Results of the coupling procedure are shown with dashed lines, while the results obtained from the BEM-solver panMARE are shown with dotted lines. Comparison is made to the measurements of MARINTEK indicated by the solid curves and to RANS simulations with a geometrically modelled propeller (marked with closed circles). A fair agreement is seen for the predicted thrust values, while larger deviations occur for the torque values. A possible reason for this deterioration can be the empirical determination of the surface-friction losses inside panMARE, which have a larger influence on the torque than on the thrust. Furthermore, small deviations between the coupled results and panMARE can be obtained. This might be influenced by different resolutions of the pressure peak at the leading edge between the involved methods.Figure 14
Comparison of predicted and measured open-water results for the investigated propeller.


### 6.2. Investigation of a Self-Propelled Ship in Regular Waves

As a second application of the coupling procedure a ship-propeller configuration in model scale is investigated. The ship is self-propelled and operates in calm water and in regular head waves with a wave length ofλ=Lpp and a wave amplitude of ζ^=0.0124·Lpp. When maneuvering in regular waves, the vessel is free to move in the vertical plane (surge, pitch, and heave). Aim of this test case is to investigate the ships behavior in head waves with special attention to the added resistance. The geometry of the vessel is depicted in Figure 15. The picture illustrates the hull shape of the vessel as well as the geometry of the pod housing and outlines the main dimensions of the vessel.Figure 15
Geometry of the investigated vessel and the aft body of the ship including the pod housing.The employed computational grid consists of 1.6 million control volumes and is depicted in Figure16.Figure 16
Computational grid of the investigated vessel.About 80% of the cells are located around the ship and in the free-surface region. The vessel refers to an offshore supply vessel equipped with two pods in pulling configuration. A symmetry condition is applied at the center plane and the propeller is not geometrically resolved in the viscous domain. This simplifies the grid generation inside the RANS method significantly because a relatively coarse and cartesian grid can be used in the propeller plane as displayed in Figure17. The panel grid employed in the inviscid method to model the propeller, which is the same as in Section 6.1, has 576 panels per blade.Figure 17
Details of the computational grid in the aft ship region with the local refinements in the propeller plane.The propeller operates at a constant number of revolutions ofn = 9 rps. In order to be able to simulate many wave periods, a relatively large timestep corresponding to an angular rotation of 22.5° of the propeller is chosen. The speed of the vessel is initialized with an estimated forward speed of Fn=0.22 in both cases. The evolution in time of the ship speed is depicted in Figure 18. Therein, the solid line displays the results of the calm water case, while the dashed line illustrates the time series of the forward velocity in waves. In calm water, a steady velocity corresponding to Fn=0.214 is reached at the end of the simulation. The predicted target speed of the vessel in waves is considerably smaller. In this case, the added resistance leads to a velocity reduction of almost 12.5% compared to the calm water case. Furthermore, larger fluctuations in the target velocity are observed. Thus, this result illustrates how even moderate wave heights can lead to a significant reduction of forward speed.Figure 18
Evolution of the forward speed of the vessel over time. The solid line refers to the calm water case, while the dashed line displays the speed losses in waves.The corresponding resistance values for both cases are plotted in Figure19 for a time interval of 10 wave periods.Figure 19
Time history of the computed resistance coefficient(RT/(ρn2D4)).Again, the calm water results are illustrated with the solid line, while the dashed line displays the ship resistance in waves. Large fluctuations are observed during the simulation in head waves. In total, the increase in the time averaged resistance value can be estimated to 23% compared to the calm water resistance.This additional resistance in head waves is mainly influenced by the pitch motion of the vessel. The time history of the pitch angle for a corresponding time frame of the simulation is depicted in Figure20. Mind that the vessel trims down by the stern in calm water, which results in a negative averaged pitch angle in head waves.Figure 20
Pitch angle (θλ/ζ^) of the vessel during the last 10 wave periods. The solid line displays the dynamic trim of the vessel in calm water, while the dashed line depicts the pitch motions due to head waves.In addition to the previous results, Figure21 displays the time history of the computed thrust values. Due to the smaller forward velocity in waves, the propeller operates at a decreased advance coefficient, which leads to larger values of the thrust coefficient.Figure 21
Comparison of the computed thrust coefficient(T/(ρn2D4)) during the last part of the simulation.In Figure22, the computed free surface elevation is shown for the time instant, when the wave crest is at the midship section. At that time the propulsor is effected by a trough and orbital velocities in the axial direction are pointing forward. The mean value of the axial velocities in the RANS solver for a slice of the wakefield between r/R=0.5 and r/R=0.7 is 8.5% smaller than for the calm water case. When the ships midsection lies in the area of a trough and the propeller therefor beneath a crest, the axial velocities are doubled. These large values are induced by the actual pitch motion of the ship. When looking at inner or outer radii, the differences are smaller as they are superposed either by the ships boundary layer or by the stagnation pressure of the hub.Figure 22
Isolines for free surface elevation of the ship in head waves.In the Figures23 and 24, a longitudinal and a transverse slice of the velocity field in the propeller region for the same time instant of the simulation is given. Also the grid with its refinement regions can be seen. The imprints of the four single blades are visible and also the imprints from the last timestep which are less accentuated. Due to the above-mentioned time-step size of 22.5°, they are in between the actual imprints. Although the grid is coarse a distinct representation of the single blades is feasible. In Figure 24 the position of the propeller in front of the pod housing and the slip stream of the propeller is clearly visible.Figure 23
Normalised velocities in a transverse cut aft of the propeller region.Figure 24
Normalised velocities in a longitudinal cut through the pod housing.The simulation included the ship model encountering 80 wave periods accounting for a real time of approximately 70 seconds. The simulation time was about two weeks on two eight-core computers.Although no validation data exists for this case, the presented results demonstrate the capabilities of the coupling procedure for the investigation of ship-propeller interaction.


## 6.1. Propeller in Open-Water Condition

The first application example refers to open-water calculations of a propeller, where model-scale experiments are performed by MARINTEK under the aegis of PROPSEAS (Kozlowska et al. [18]). The investigated propeller has a diameter of D = 0.2 m and operates with n = 18 rps. Different advance coefficients are achieved by a variation of the inflow velocity. The numerical grid of the propeller used in the BEM consists of 576 panels per blade (cf. Figure 13(a)). Figure 13(b) displays the propeller regime of the viscous grid. About 2.500 of 15.000 control volumes in the propeller disc are subjected to body forces. The Reynolds number based on the chord length at r/R=0.7, and the resulting velocity at the propeller blade reads Re=0.5·106. The flow is considered fully turbulent, and the standard k-ω turbulence model is applied.Numerical grid and resulting pressure distribution for the BEM simulation (a). FV grid in the propeller disc (b). Only cells at the actual blade positions contain body forces (colored).
(a)
(b)The resulting open-water diagram of the test case can be found in Figure14. Results of the coupling procedure are shown with dashed lines, while the results obtained from the BEM-solver panMARE are shown with dotted lines. Comparison is made to the measurements of MARINTEK indicated by the solid curves and to RANS simulations with a geometrically modelled propeller (marked with closed circles). A fair agreement is seen for the predicted thrust values, while larger deviations occur for the torque values. A possible reason for this deterioration can be the empirical determination of the surface-friction losses inside panMARE, which have a larger influence on the torque than on the thrust. Furthermore, small deviations between the coupled results and panMARE can be obtained. This might be influenced by different resolutions of the pressure peak at the leading edge between the involved methods.Figure 14
Comparison of predicted and measured open-water results for the investigated propeller.


## 6.2. Investigation of a Self-Propelled Ship in Regular Waves

As a second application of the coupling procedure a ship-propeller configuration in model scale is investigated. The ship is self-propelled and operates in calm water and in regular head waves with a wave length ofλ=Lpp and a wave amplitude of ζ^=0.0124·Lpp. When maneuvering in regular waves, the vessel is free to move in the vertical plane (surge, pitch, and heave). Aim of this test case is to investigate the ships behavior in head waves with special attention to the added resistance. The geometry of the vessel is depicted in Figure 15. The picture illustrates the hull shape of the vessel as well as the geometry of the pod housing and outlines the main dimensions of the vessel.Figure 15
Geometry of the investigated vessel and the aft body of the ship including the pod housing.The employed computational grid consists of 1.6 million control volumes and is depicted in Figure16.Figure 16
Computational grid of the investigated vessel.About 80% of the cells are located around the ship and in the free-surface region. The vessel refers to an offshore supply vessel equipped with two pods in pulling configuration. A symmetry condition is applied at the center plane and the propeller is not geometrically resolved in the viscous domain. This simplifies the grid generation inside the RANS method significantly because a relatively coarse and cartesian grid can be used in the propeller plane as displayed in Figure17. The panel grid employed in the inviscid method to model the propeller, which is the same as in Section 6.1, has 576 panels per blade.Figure 17
Details of the computational grid in the aft ship region with the local refinements in the propeller plane.The propeller operates at a constant number of revolutions ofn = 9 rps. In order to be able to simulate many wave periods, a relatively large timestep corresponding to an angular rotation of 22.5° of the propeller is chosen. The speed of the vessel is initialized with an estimated forward speed of Fn=0.22 in both cases. The evolution in time of the ship speed is depicted in Figure 18. Therein, the solid line displays the results of the calm water case, while the dashed line illustrates the time series of the forward velocity in waves. In calm water, a steady velocity corresponding to Fn=0.214 is reached at the end of the simulation. The predicted target speed of the vessel in waves is considerably smaller. In this case, the added resistance leads to a velocity reduction of almost 12.5% compared to the calm water case. Furthermore, larger fluctuations in the target velocity are observed. Thus, this result illustrates how even moderate wave heights can lead to a significant reduction of forward speed.Figure 18
Evolution of the forward speed of the vessel over time. The solid line refers to the calm water case, while the dashed line displays the speed losses in waves.The corresponding resistance values for both cases are plotted in Figure19 for a time interval of 10 wave periods.Figure 19
Time history of the computed resistance coefficient(RT/(ρn2D4)).Again, the calm water results are illustrated with the solid line, while the dashed line displays the ship resistance in waves. Large fluctuations are observed during the simulation in head waves. In total, the increase in the time averaged resistance value can be estimated to 23% compared to the calm water resistance.This additional resistance in head waves is mainly influenced by the pitch motion of the vessel. The time history of the pitch angle for a corresponding time frame of the simulation is depicted in Figure20. Mind that the vessel trims down by the stern in calm water, which results in a negative averaged pitch angle in head waves.Figure 20
Pitch angle (θλ/ζ^) of the vessel during the last 10 wave periods. The solid line displays the dynamic trim of the vessel in calm water, while the dashed line depicts the pitch motions due to head waves.In addition to the previous results, Figure21 displays the time history of the computed thrust values. Due to the smaller forward velocity in waves, the propeller operates at a decreased advance coefficient, which leads to larger values of the thrust coefficient.Figure 21
Comparison of the computed thrust coefficient(T/(ρn2D4)) during the last part of the simulation.In Figure22, the computed free surface elevation is shown for the time instant, when the wave crest is at the midship section. At that time the propulsor is effected by a trough and orbital velocities in the axial direction are pointing forward. The mean value of the axial velocities in the RANS solver for a slice of the wakefield between r/R=0.5 and r/R=0.7 is 8.5% smaller than for the calm water case. When the ships midsection lies in the area of a trough and the propeller therefor beneath a crest, the axial velocities are doubled. These large values are induced by the actual pitch motion of the ship. When looking at inner or outer radii, the differences are smaller as they are superposed either by the ships boundary layer or by the stagnation pressure of the hub.Figure 22
Isolines for free surface elevation of the ship in head waves.In the Figures23 and 24, a longitudinal and a transverse slice of the velocity field in the propeller region for the same time instant of the simulation is given. Also the grid with its refinement regions can be seen. The imprints of the four single blades are visible and also the imprints from the last timestep which are less accentuated. Due to the above-mentioned time-step size of 22.5°, they are in between the actual imprints. Although the grid is coarse a distinct representation of the single blades is feasible. In Figure 24 the position of the propeller in front of the pod housing and the slip stream of the propeller is clearly visible.Figure 23
Normalised velocities in a transverse cut aft of the propeller region.Figure 24
Normalised velocities in a longitudinal cut through the pod housing.The simulation included the ship model encountering 80 wave periods accounting for a real time of approximately 70 seconds. The simulation time was about two weeks on two eight-core computers.Although no validation data exists for this case, the presented results demonstrate the capabilities of the coupling procedure for the investigation of ship-propeller interaction.


## 7. Conclusions

The paper reports the development of a close-coupling algorithm between an inviscid boundary element method dedicated to propeller flows and a RANS method. The procedure uses the propeller forces computed by the inviscid method to mimic the propeller inside the viscous-flow solver, while the inflow condition for the propeller code is provided by the RANS method.A critical aspect in terms of the accuracy of the procedure is the determination of the induced velocity field. Here, the presented verification study reveals a sufficient accuracy, which assures equivalent working conditions for the propeller in both solvers. Encouraging results are obtained for the two application examples. Furthermore, the computational effort is reduced significantly compared to simulations with a geometrically resolved propeller. The grid generation for the RANS solver is simplified as coarse and simple hexahedral grids can be used in the propeller plane.Moreover, the coupling algorithm is robust and shows a favorable convergence behavior.



---

*Source: 743060-2012-07-11.xml*