Abstract

Propulsor crashback is an off-design operating condition where a propulsor rotates in the reverse direction to yield negative thrust. Crashback is characterized by the interaction of the free stream with the reverse flow generated by propulsor rotation. This causes a highly unsteady vortex ring which leads to flow separation and unsteady forces and moments on the blades. Large eddy simulation (LES) is performed for marine propulsors in crashback for various configurations and advance ratios and validated against experiments. The predictive capability of LES as a tool for propulsor crashback is demonstrated on an open propulsor, open propulsor with a submarine hull, and ducted propulsor with and without stator blades. LES is in good agreement with experiments for the mean and RMS levels, and spectra of the unsteady loads on the propulsors.

1. Introduction

The advent of supercomputers has enabled the routine use of computational fluid dynamics for the development of engineering devices. Various methodologies exist for the numerical simulations of fluid flows, each with their inherent approximations and models. Reynolds-Averaged Navier-Stokes (RANS) is the most commonly used tool for engineering applications. RANS equations with a turbulence model are based on time average or ensemble average, and they model the Reynolds average of the turbulent fluctuations (Reynolds stresses). In Direct numerical simulation (DNS), all the scales of motion are resolved accurately in space and time, and no model is used. However, the computational cost associated with DNS has limited its use to study turbulent flows in simple geometries [1]. Large eddy simulation (LES) is an unsteady, three-dimensional simulation methodology intermediate between RANS and DNS. There is increasing interest in applying LES to complex turbulent flows because of more reliability over RANS in massively separated flows and relaxed computational requirements over DNS.

RANS typically uses upwind biased schemes which provide numerical dissipation and make the solution procedure robust. Upwind schemes compromise on accuracy since the numerical dissipation often overwhelms the effect of the subgrid model. On the other hand, standard nondissipative central-difference schemes have numerical instabilities related to aliasing errors. Mahesh et al. [2] developed a nondissipative but robust finite volume method for LES by satisfying discrete kinetic energy conservation. This numerical method provides the capability to perform LES of complex turbulent flows. LES is performed on marine propulsors in the crashback mode of operation in this paper.

The “four-quadrant” performance test evaluates propeller performance in the four operational quadrants of the vessel velocity ()-propulsor rotational speed () plane as shown in Figure 1(a). The design operating condition is forward (, ) while the off-design conditions are crashback (, ), backing (, ), and crashahead (, ). Crashback is the operating condition where the vessel moves in the forward direction while the propulsor rotates in the reverse direction to decelerate the forward moving vessel.

Flow around a propulsor in crashback is characterized by large scale unsteadiness and separation. Low frequency, high amplitude off-axis forces and moments produced by the unsteadiness are transmitted to the vessel, inhibiting the ability to maneuver efficiently. The crashback condition is dominated by the interaction of the free stream flow with the strong reverse flow from reverse propulsor rotation as shown in Figure 1(b). This interaction forms a highly unsteady vortex ring around the propulsor. Crashback displays unsteady massively separated flow over the propulsor blades. LES is therefore a suitable simulation methodology for crashback.

Traditionally, crashback has been studied experimentally. Four-quadrant evaluations for a series of propellers were conducted by Hecker and Remmers [6] using open water (OW) experiments. Thrust and torque coefficients were measured over a wide range of advance ratios. Jiang et al. [7] studied the structure of the unsteady vortex ring using particle image velocimetry (PIV) measurements for propeller P4381 in water tunnel (WT). Jessup et al. [4] presented more detailed measurements of flow around the same propeller using PIV and laser Doppler velocimetry (LDV). Bridges et al.'s experiment [5] documented the effect of hull on an open propulsor. The effect of a duct was studied experimentally by Jessup et al. [8] and Donnelly et al. [9].

Simulations of the flow around open propulsors have been performed using unsteady Reynolds-Averaged Navier-Stokes equations (RANS) [10, 11]. These studies showed that RANS yielded good results for forward and backing modes (attached flow regime) but produced significant discrepancies in crashback and crashahead modes (separated flow regime). Vyšsohld and Mahesh [12] performed LES of flow around P4381 in crashback at an advance ratio of . They showed that LES yields good agreement for mean and RMS of unsteady loads as well as spectra. Jang and Mahesh [13] investigated the effect of the duct in crashback for the same propeller. Chang et al. [3] performed LES at other advance ratios, and , and used the LES-generated surface forces to predict the structural response of the blades. That demonstrated the fluid-structure interaction (FSI) capability of LES when coupled with a finite-element structural solver. Jang and Mahesh [14] introduced two quantities for pressure contributions to forces to understand the origin of thrust and side force and used conditional averaging to study extreme amplitude events. Verma et al. [15] performed LES of the flow past a propulsor in crashback attached to an upstream submarine hull.

The objective of the present paper is to (1) summarize the use of LES at the University of Minnesota to predict unsteady loads on marine propulsors in crashback and (2) demonstrate the utility of LES in isolating the physical effects of the hull, duct, and stator. The paper is organized as follows. Section 2 introduces the governing equations and numerical method for crashback simulations. LES for open propulsor is performed in Section 3.1. The effect of a hull is investigated in Section 3.2. Simulations are performed for a ducted propulsor without stator blades in Section 3.3 and with stator blades in Section 3.4.

2. Governing Equations

Since all stationary parts are axisymmetric, the open propulsor in Section 3.1, the open propulsor with an axi-symmetric hull in Section 3.2, and the ducted propulsor without stator blades in Section 3.3 can be solved in a rotating frame of reference. On the other hand, the ducted propulsor with stator blades in Section 3.4 requires a moving grid technique, called the sliding interface method.

2.1. Navier-Stokes Equations in a Rotating Reference Frame

The simulations are performed in a noninertial frame of reference that rotates with the propulsor. The incompressible Navier-Stokes equations in the rotating frame of reference are written in terms of the absolute velocity. In LES, large scales are directly solved from the spatially filtered N-S equations, while small scales are accounted for by modeling the subgrid stress. The filtered N-S equations in the rotating frame of reference are as follows: where is the absolute velocity in the fixed frame, is the pressure, are coordinates in the rotating frame, is the angular velocity of the rotating frame, is the kinematic viscosity, is the permutation symbol, the overbar denotes the spatial filter, and is the subgrid stress. Note that (1) has two additional terms which account for the system rotation ( and ). The dynamic Smagorinsky model proposed by Germano et al. [16] and modified by Lilly [17] is used to model the subgrid stress in the paper.

Equation (1) is solved using a numerical method developed by Mahesh et al. [2] for incompressible flows on unstructured grids. The algorithm is derived to be robust without numerical dissipation. It is a finite-volume approach which stores the Cartesian velocities and the pressure at the centroids of the cells, and the face normal velocities are stored independently at the centroids of the faces. A predictor-corrector approach is used. The predicted velocities at the control volume centroids are first obtained and then interpolated to obtain the face normal velocities. The predicted face normal velocity is projected so that the continuity equation is discretely satisfied. This yields a Poisson equation for pressure which is solved iteratively using a algebraic multigrid method. The pressure field is used to update the Cartesian velocities using a least-squared approach for minimizing the conservation error. Implicit time advancement is performed using the Crank-Nicholson scheme. The algorithm has been validated for a variety of problems over a range of Reynolds numbers [2].

2.2. ALE Formulation for a Sliding Interface Method

The sliding interface method is used to handle relatively rotating grid systems. Relatively rotating grids do not have to overlap since the movement of the grid is a rigid rotation. The only intersection is an interface between the grids, which is called the sliding interface. Data exchange between the grids can be performed only at the sliding interface so that no remeshing, deformation, or hole cutting is required. The sliding interface method has been developed on multiblock structured grids [18, 19] and on unstructured grids [20, 21], and is popular to handle relatively rotationg grid systems and on unstructured grids.

The sliding interface method developed for this study is capable of being applied on arbitrarily shaped unstructured grids on massively parallel computing platforms. Figure 2 shows a schematic for the sliding interface method on an unstructured grid. Sliding elements of a sub-domain are created by extrusion to the adjacent sub-domain. However, sliding elements of the subdomain and boundary cells of the adjacent sub-domain do not match due to the relative rotation between the sub-domains. As shown in Figure 2, a sliding element (black dashed) is extruded from a boundary control volume (black shaded) in the left sub-domain. Since the grids of the two sub-domains do not align, the control volume containing the sliding element is not known without a search process. The control volume is named the host element and shaded in blue.

The governing equations for the sliding interface method are expressed in the Arbitrary Eulerian-Lagrangian (ALE) formulation [22] as follows: where is the grid velocity, which is given by in this study due to the rigid rotation of the grid. Since there is no mesh deformation, the geometric conservation law is automatically satisfied so that no spurious mass sources or sinks are created.

3. Results

Large eddy simulations are performed under the crashback condition at negative advance ratios and Reynolds number . The advance ratio and Reynolds number are defined as where is the free-stream velocity, is the diameter of the propeller disk, , and is the kinematic viscosity. The thrust is defined as the axial component of the force, and torque is the axial component of moment of the force. and denote horizontal and vertical components of the force whose vector sum yields side-force . Nondimensional thrust coefficient , torque coefficient , and side force coefficient are defined as where a rotational velocity is used as the reference velocity, and is used as the reference area in this normalization. Hereafter, denotes mean value, and denotes standard deviation.

3.1. Open Propulsor
3.1.1. Problem Description

The simulations are performed for marine propeller P4381, which is five bladed, right handed with variable pitch, and has no skew and no rake. The propeller has been used in various experiments [4, 7, 8] and computations [3, 1012, 14]. Details of the blade and hub geometry are given in Jessup et al. [4]. The present simulation is compared to Jessup et al.’s [4] water tunnel (WT) experiment. Statistics of unsteady loads are also compared to open water (OW) experiments conducted at David Taylor Model Basin.

The computational domain is a cylinder with the diameter of and the length of where is the diameter of the propeller disk. The computational domain is based on the WT geometry. Free-stream velocity boundary conditions are specified at the inlet and the lateral boundaries. Convective boundary conditions are prescribed at the exit. Since velocities in the governing equations are prescribed in the fixed frame of reference, boundary conditions on solid walls are also prescribed in the fixed frame. A schematic of the computational domain and boundary conditions is shown in Figure 3(a).

The computational grid is shown in Figure 3(b). The whole grid system is composed of two components: a cylindrical subdomain including blades and hub surface, and the outer region of the sub-domain. The cylindrical sub-domain is filled with tetrahedral elements to match the complex geometry of the propeller, and the outer region consists of hexahedral elements. Prism layers are extruded from the blade surface in order to improve resolution near the wall at the first height of 0.0017. Recall that propulsor crashback is inherently a massively separated flow regime. Hence, resolution of the near-wall flow is not as strict a requirement as it is for attached flows (such as forward/backward condition). The computational grid consists of 19.3 million control volumes.

3.1.2. Time History of Unsteady Loads

A time history of in Figure 4(a) shows large fluctuations and low frequency oscillations in the blade loading. Table 1 compares statistics of the unsteady loads to experimental results at three advance ratios (, , and ). Good agreement is observed, and the computed mean values of and are between results from WT and OW experiments. The agreement with experimental data is within the scatter observed between the experiments. Similarly, standard deviations and side-force magnitude show good agreement with the WT experiment. The power spectral density also shows similar agreement. Overall, power spectral density in Figure 4(b) shows good agreement except for the high-frequency region () where experimental uncertainty is large. According to Jessup (private communication), the differences at high-frequency region could occur due to blade bending and vibrations, or other shaft-related resonances affecting the measured data. The most noticeable peak is found at , the blade passage frequency. Figure 4(c) shows the presence of a highly unsteady vortex ring with an isosurface of constant low pressure in the LES results.

3.1.3. Flow Field

The computed results are averaged in time over a period of 135.7 revolutions. Since time-averaged flow field in the WT experiment was measured in the plane of the fixed frame, prior to comparison the time-averaged field from the LES is further averaged over the circumferential direction. Profiles from circumferentially averaged axial velocity are extracted from five -locations and compared to experimental data in Figure 5. The profiles show good agreement at all locations.

In order to describe the distribution of thrust and side force on the blade surface, a blade surface is divided into ten constant-radius sections on both pressure and suctions sides, mean thrust and side force are computed in each radial section and compared in Figure 6. on the suction side linearly increases as the radius increases except for the blade tip. This exception is due to the smaller area in the section at the blade tip. The suction side generates more and than the pressure side. on the suction side does not show radial dependency, unlike thrust.

Circumferentially averaged streamlines at three advance ratios are compared in Figure 7. Note that the vortex ring core is located right above the blade top at . This behavior may be explained by stronger reverse flow induced by the faster propeller rotation at .5 overcoming the free stream momentum so that the vortex ring moves closer to the propeller. On the other hand, the reverse flow at is weaker due to the slower rotation. The vortex ring moves downstream at high negative .

3.2. Open Propulsor with Hull
3.2.1. Problem Description

The objective of this section is to evaluate the ability of LES to predict the interaction of an upstream hull with the propulsor during crashback. Bridges et al.'s experiment [5] on an open propulsor with an upstream submarine hull in the Large Cavitation Channel (LCC) noted that the side force increased dramatically below an advance ratio of compared to Jessup’s experiment [8] without the hull in the WT. A standard axisymmetric hull (DTMB Model 5495–3) is used for the hull geometry. Details of the hull geometry are given in [5]. In the simulations, half of the hull body is used and stabilizing fins are ignored. The size of the computational domain is exactly the same as the open propulsor and boundary conditions are similarly specified. The computational grid is created using a similar strategy as in Section 3.1.1. The unstructured grids for both the cases (with and without hull) are shown in Figure 8.

3.2.2. Time History of Unsteady Loads

Simulations are performed for the propulsor with and without hull at and . is lower than the critical advance ratio of , and is higher than the critical advance ratio. Table 2 shows that computed side-force magnitude are in good agreement with the experimental results for propulsor with and without hull at both advance ratios. The side-force magnitude at is slightly lower than experiment with hull, but importantly the LES predicts the experimentally observed increase in side force with hull at . On the other hand, the hull does not significantly change the side-force magnitude at .

The power spectral density is shown for horizontal force with and without hull in Figure 9 at . The most important peak is at which has also been observed in [3, 12, 14]. This corresponds to the passage of the blades of the five-bladed propulsor. A low frequency, high amplitude modulation of the side-force has important ramifications for the maneuverability of the vehicle. For the spectra, there is a peak at a lower frequency of with the hull which agrees very well with that observed by Bridges et al. [5] at . The comparison with experiments is summarized in Table 3.

3.2.3. Flow Field

Time-averaged statistics of the flow field are computed over 170 propulsor rotations for the propulsor with hull and 172 rotations without hull at . Circumferentially averaged streamlines in Figure 10 reveal a recirculation zone upstream of the blades in the presence of the hull. No such recirculation zone appears near the shaft without hull. In the presence of the hull, the vortex ring appears to be more compact and its center is closer to the blades compared to the more elongated and stretched out vortex ring without the hull. Note that the vortex ring is much closer to the tip of the blade when the hull is present. Table 4 shows good agreement for the distance of the center of the vortex ring relative to the propulsor center between LES and experiment [5].

In order to compare side-force distributions with and without hull, side-force magnitudes are computed in ten radial sections in Figure 11. Note that higher side force magnitude is observed with hull than without hull. For both cases, the suction side generates more side-force than the pressure side. Most of the side force is generated from the blade root area in both cases.

3.3. Ducted Propulsor without Stator Blades
3.3.1. Problem Description

To investigate the effect of a duct on crashback performance, a neutrally loaded duct was added to fit around propeller 4381 in Jessup et al. [8]. The duct was designed to not yield additional propulsor loading at the design advance ratio of . The designed ducted propulsor also had 13 stator blades upstream of rotor blades. Simulations are performed at the same Reynolds number and the same advance ratio as in the open propulsor of Section 3.1. The flow around the open propulsor has been solved in a rotating frame of reference. However, the ducted propulsor with stator blades cannot be solved in the rotating frame since the stator blades are not axisymmetric. In this section, stator blades are ignored as shown in Figure 12(a) to simulate the ducted propulsor in the rotating reference frame and to isolate the effect of the duct and the stator blades. The stator blades are considered using the sliding interface method in Section 3.4.

The size of the computational domain is exactly the same as the open propulsor and boundary conditions are similarly specified. The computational grid is created using a similar strategy as in Section 3.1.1. The computational grid in the propeller neighborhood is shown in Figure 12(b). The total number of the computational c.v.'s is approximately 19 million.

3.3.2. Time History of Unsteady Loads

Statistics of , , and are computed and compared with experimental results of Jessup et al. [8] and Donnelly et al. [9]. Since duct forces were not measured in Jessup et al.’s [8] experiment, statistical values on the blade surface are compared with the experimental data. Donnelly et al. [9] measured side-force magnitude on the duct surface. Table 5 compares computed means and standard deviations of unsteady loads to the experimental results. The computed results show good agreement with the experimental results. Recall that and from the LES are located between the WT and the OW results in the open propulsor case. Even though the computational and are larger than the WT experiment, the geometric assumptions are likely a source of error, nevertheless, the LES results are expected to be close to the experimental results. Also, on the blade and duct surfaces are quite close to the WT data. The effect of the duct for the thrust is relatively small and that for the torque is negligible. However, the side force shows totally different behavior. on the duct is about 4 times larger than on the blade. This indicates that the duct surface is the primary contributor to side force.

3.3.3. Flow Field

The computed results are averaged in time over a period of 93.8 revolutions. Figure 13 compares the circumferentially averaged flow field for the ducted propulsor (left figure) with that for the open propulsor (right figure). The vortex ring core is located outside the duct in Figure 13(a). The pressure difference between pressure and suction sides are much higher for the ducted propulsor. A high pressure region is located near the tip gap of the duct, which can cause tip leakage flow.

To investigate where the majority of the side force is generated on the duct, contours of RMS pressure are plotted on the duct surface in Figure 14. The duct surface is divided into inner surface and outer surface and unfolded on the plane. The slanted lines in Figure 14(a) represent locations of the blade tips. Note that the levels of RMS pressure on the inner surface are much higher than those on the outer surface, especially in the vicinity of blade tips. This behavior indicates that the side force on duct mostly originates from blade-duct interaction.

3.4. Ducted Propulsor with Stator Blades
3.4.1. Problem Description

In the previous section, the stator blades were ignored when computing the ducted propulsor in the rotating frame of reference. Even though the stator blades are designed to contribute no additional loading at the design advance ratio, their performance in crashback is unknown. Thus, the stator blades are considered with the sliding interface method introduced in Section 2.2. The computed rotor-stator problem with a relatively rotating grid system is shown in Figure 15(a).

The size of the computational domain is exactly the same as the open and ducted propulsors. Since velocities in the governing equations are specified in the fixed frame of reference similar to the rotating reference frame method, boundary conditions are similarly specified. The computational grid is created using a similar strategy as in Section 3.1.1, but hexahedral elements are used near the sliding interface to obtain better accuracy from the interpolation at the nonmatching boundary. Figure 15(b) shows the computational grid for the ducted propulsor with the stator blades. In the figure, the red-colored line represents the sliding interface. The total number of the computational c.v.'s is approximately 27 million.

3.4.2. Preliminary Results

A preliminary simulation with the sliding interface method is performed for the ducted propulsor with stator blades. The simulation is performed at the same Reynolds number and advance ratio of . It is performed for 34.7 revs, and statistics of the unsteady loads are computed over 18.9 revs. Time history of of is shown in Figure 16. Table 6 compares computed means and standard deviations of unsteady loads with the experimental results. The computed results show reasonable agreement with the experimental results. Since the stator blades are attached on the duct surface, the duct force can be considered as the sum of loads on duct and stator blades. computed on the duct and stator surfaces are in favorable agreement. The effect of the duct on thrust is relatively small, and that on the torque is significantly large. Also, the duct still is the major contributor to side force.

Instantaneous flow fields are shown in Figure 17. The highly unsteady flow interaction is observed, and pressure contours are smoothly connected at the sliding interface in Figure 17(a). The pressure distribution is plotted in Figure 17(b). High-pressure regions are observed at the leading edge of the stator blades near the inner surface of the duct.

4. Summary

The predictive capability of LES as a tool for highly complex flow over marine propulsors is demonstrated. LES has been performed at an off-design condition of marine propulsors, called crashback, which is well known as one of the most challenging propeller operations to analyze. The interaction of the free stream with reverse flow due to reverse propulsor rotation causes a highly unsteady vortex ring which leads to flow separation and unsteady forces and moments on the blades. LES has been performed for an open propulsor, open propulsor on a hull, and ducted propulsor with and without stator blades in crashback. The LES uses a discretely kinetic energy conserving algorithm developed by Mahesh et al. [2] for unstructured grids.

For the open propulsor, the simulations are validated against experiments performed by Jessup et al. [4], and the results show good agreement with the experiments. LES is able to reproduce the experimentally observed increase of side force for an open propulsor with hull below an advance ratio of . The mean, RMS, and spectra of unsteady loads are in good agreement with the experiments of Jessup et al. [4] and Bridges et al. [5]. For a ducted propulsor, most of the thrust is still due to the blades, but the blade-duct interaction causes side force to be generated from the duct surface. A sliding interface method is developed to enable simulations of rotor-stator propulsors on parallel computing platforms. Preliminary results from the application of the sliding interface method to a ducted propulsor with stator blades show good agreement with experiment.

Acknowledgments

This work was supported by the Office of Naval Research under ONR Grant N00014-05-1-0003 with Dr. Ki-Han Kim as program manager. Computing resources were provided by the Arctic Region Supercomputing Center (ARSC) of the High Performance Computing Modernization Program (HPCMP) and the Minnesota Supercomputing Institute (MSI). The authors are grateful to Dr. Stuart Jessup, Dr. Martin Donnelly, Dr. Christopher Chesnakas, and their colleagues at NSWCCD for providing them with experimental data. They would also like to express their gratitude to Dr. Peter Chang for useful discussions.