Design by Optimization on the Nozzle and the Stator Blades of a Rim-Driven Pumpjet

Abstract

The design of the stator and nozzle of a rim-driven pumpjet thruster (RDPJ) is addressed through a simulation-based design optimization approach built on a parametric description of the main geometrical characteristics of the system, a RANS solver with actuator disk model, and a genetic algorithm. As the propeller blades’ geometry is fixed, the rotor/stator (RDPJ-R/S) configuration is considered for the optimal design from a multi-objective optimization process aimed at minimizing the resistance keeping the cavitation inception index at the lowest possible value. Steady-state (moving reference frame plus mixing plane interface) and unsteady simulations (sliding meshes) with fully resolved rotor geometry were finally carried out on six selected optimal geometries to validate the optimization process and the performance improvements provided by the RDPJ configuration when compared with the original rim-driven thruster (RDT).

1. Introduction

In an era where environmental considerations are paramount, the ability to minimize underwater radiated noise while enhancing efficiency has led to a growing interest in unconventional propulsion systems (Su et al. [1]). In this context, the development of rim-driven thrusters has emerged as a ground-breaking advancement in maritime propulsion technology, reviving concepts that date back to the pioneering works of Kort [2] and Saunders [3], who first proposed driving a propeller via its blade tips (Grümmer [4]). Recent innovations in materials and manufacturing techniques, particularly the advent of brushless permanent magnet electric motors (Hughes et al. [5]), have enabled the practical application of this propulsive concept, which had remained at a purely conceptual stage until recent years. This type of thruster has several advantages among which, from a hydrodynamic efficiency point of view, the most interesting one is the minimization of the energy losses by effectively eliminating the gaps between the propeller tips and the nozzle inner surface. This reduces the formation of the propellers’ characteristic tip vortices like the leakage vortex of ducted propellers (Baltazar et al. [6]) and enhances the overall performance of the propulsor system. However, despite their promising advantages, the hydrodynamic design of rim-driven thrusters poses unique challenges, particularly due to the limited literature available on this novel propulsion mechanism. In fact, although visually this type of thruster is very reminiscent of traditional ducted propellers, the absence of tip clearance and the presence of a gap between the rim and duct where peculiar flow structures are formed (Yan et al. [7]), combined with the different load distribution on the blade (Hughes et al. [5]), make it nearly impossible to successfully use the methods developed for traditional propellers. Numerous studies (Liu and Vanierschot [8], Jiang et al. [9]) on the hydrodynamic efficiency of RDTs emphasize that, in the absence of clearance between the nozzle inner surface and the blade tip, the radial circulation—which is zero at the tip for conventional propellers—can have a finite value. Under specific conditions related to pitch, chord, and camber distributions, this circulation may be the maximum (Cao et al. [10]). In contrast, zero circulation is typical at the blade root, which increases the likelihood of stronger root vortices. Enhancing or maintaining a sufficiently high load in the outer part of the blade is advantageous for propulsive efficiency. Furthermore, sealing the gap between the blade and the nozzle inner surface greatly minimizes the risk of tip vortices (and avoids leakage), even under high load conditions, thus preventing the occurrence of noisy (cavitating) tip vortices. However, it is essential for the hydrodynamic design of the blade to account for the possibility that increased loading may result in excessively high suction pressures and, consequently, a higher risk of bubble cavitation. From a noise perspective, this can be as detrimental as cavitating tip vortices. Additionally, traditional design methods based on potential flow may prove too simplistic (Kerwin et al. [11]), particularly in the case of unconventional configurations and challenging operating conditions.

Early attempts to design rim-driven propellers, such as those by Sparenberg [12], Pashias and Turnock [13], and Mishkevich [14], employed modified lifting line theories to explore blade and nozzle loading distributions. Yakovlev et al. [15] proposed an even more simplified method using a hydrodynamic model consisting of 2D calculations for a lattice of hydrofoils and considering the pitch distribution as the only design parameter of the approach. Recent advancements in computational fluid dynamics, particularly the application of Reynolds-Averaged Navier–Stokes (RANS) calculations (Cao et al. [10], Song et al. [16]), have further enabled the accurate prediction of hydrodynamic performance for these propulsors (Dubas et al. [17], Liu et al. [18], Lu et al. [19]). In fact, conventional approaches may suffer from crude approximations (Kinnas et al. [20], Lea et al. [21]), particularly when the design space is expanded to investigate unusual configurations, like RDT, aimed at optimizing efficiency and cavitation inception. In contrast, the RANS calculations often employed nowadays allow for a precise examination of the unique characteristics of rim-driven thrusters, thus facilitating performance maximization. In fact, RANS calculations have been recently utilized for design purposes, since they more effectively address the unique characteristics of the propulsors and help to optimize their performance in terms of both efficiency and noise reduction.

There are few examples in the literature that involve the use of RANS in a design context, which, due to the nature of RANS calculations, naturally align with optimization paradigms (Liu et al. [22], Zhai et al. [23]). Gaggero et al. [24] and Gaggero [25] introduced a simulation-based design by optimization (SBDO) approach entirely based on RANS calculations for the design of both conventional and unconventional marine thrusters. This methodology is then also applied for designing rim-driven propellers in both accelerating and decelerating ducts [26]. In both scenarios, the optimal rim configurations demonstrated superior performance compared to conventional ducted propellers, achieving similar efficiency levels and significantly higher cavitation inception speeds.

In this work, the same methodology is used to investigate a rim-driven pumpjet thruster (RDPJ) in which, following the same idea that made the transition from ducted propellers to pumpjets, a stator is added downstream of the rotor blades to increase the overall thruster efficiency. Pumpjet propulsors, indeed, represent a unique type of ducted propeller that integrates an additional system of stator blades alongside the propeller blades and nozzle. These stator blades can be positioned either upstream (stator/rotor, or pre-swirling) or downstream (rotor/stator, or post-swirling) of the rotating stage. In recent years, they have garnered significant interest as highly efficient propulsive systems, particularly for high-speed underwater vehicles and submarines, due to their capability to deliver high thrust under heavy loading while significantly reducing radiated noise levels (McHugh [27], Kinnas et al. [28], Gaggero and Martinelli [29]).

The focus of this work is to progress towards the design of pumpjet propulsors in which the rotor stage is realized through the rim-driven concept. By applying a design-by-optimization methodology, the aim is to obtain the optimal design of the duct and the stator for an assigned rotor geometry, investigating in this first step of research the possibility of decoupling the rotor design from the rest (nozzle and stator) of the system. Compared to the simultaneous design of all elements of the RDPJ, this would result in a much lighter optimization problem (since, without the need to parametrically describe the rotor, the number of design parameters is sensibly reduced). Therefore, this allows one, in a preliminary design stage, to explore wider design spaces and include higher geometric variability in the search for optimal configurations that, for their unconventional nature, may benefit from unusual geometrical shapes.

Following a brief description of the addressed problem, the model-scale rim-driven thruster used as reference for the design is presented in Section 2. Section 3 discusses the methodology, hypothesis, design workflow, and computational meshes, implemented using StarCCM+ [30] and in-house developed code. The optimization process is conducted at a model scale (D = 0.23 m), leveraging extensive data derived from previous work (Gaggero et al. [31]) verified against towing tank measurements and observations of inception in the cavitation tunnel. Section 4 summarizes the results of the optimization activities, which include different geometries on the Pareto front coming out from the multi-objective optimization driven by efficiency and cavitation behavior.

2. Problem Description and Reference Geometry

In recent decades, the design of both conventional and unconventional marine propellers has evolved significantly. Traditional design tools primarily focused on maximizing efficiency. Recently, instead, propeller-related issues, such as pressure pulses and radiated noise, have become crucial constraints to meet comfort standards and stricter environmental regulations. Key challenges for designers include preventing or managing cavitation during operations, creating highly efficient propellers for various operating conditions, and minimizing on-board vibrations. Traditional design tools often struggle with these issues. At the same time, the propeller design process has transformed due to new simulation tools that can address these modern requirements. Advancements in numerical algorithms and hardware have enhanced the robustness and accuracy of simulations, allowing for methodologies that can reduce or even eliminate the need for experimental testing. These simulations can provide essential data early in the design process to meet additional criteria and objectives related to increasing technical and commercial demands. As a result, simulation-based design (SBD) frameworks have become standard in propeller design, facilitating efficient connections among tools and methods for multi-disciplinary analyses. These advanced simulations aim to replicate relevant phenomena with high accuracy. Achieving the best fidelity at the lowest cost is critical but challenging. While precise design analysis methods require significant computational resources, they are essential for ensuring flexibility and validity, especially for novel concepts and unconventional configurations that challenge the limits of computational efficient, low- and mid-fidelity, methods. Evaluating design performances or comparing trade-off alternatives is vital for continuous improvements in a competitive market. This requires exploring larger design spaces and more options, which necessitate the use of high-fidelity simulations and complex geometrical representations to enable substantial performance enhancements. In this framework, SBD evolved into simulation-based design by optimization (SBDO), where simulation tools, optimization algorithms, and parametric geometry descriptions are seamlessly integrated to manage numerous decision parameters and conflicting design objectives and constraints. Navigating high-dimensional design spaces is essential, as exploring unconventional configurations beyond traditional design rules is often the simplest route to performance improvements. However, relying on low- or medium-fidelity computational approaches limits the design space exploration, particularly for unconventional designs. It is crucial to avoid design optima that are skewed by insufficient accuracy and the limitations of low- and medium-fidelity flow solvers. Consequently, RANS calculations are favored over the previously used, more computationally efficient lifting surface or BEM analyses.

In this context, the goal of the proposed work is to apply the SBDO paradigm to improve the performances of rim-driven thrusters combined with the philosophy of pumpjets by adding stator blades to the propulsors. Since one of the crucial objectives of the design is the reduction in the risk of cavitation, this work focuses on the optimization of the performance of thrusters equipped with divergent nozzles (divergent rim-driven, DRD), leading to rim-driven pumpjets (RDPJs) of which an example geometry is shown in Figure 1a.

This thruster, in addition to the rotor, rim, and duct, is fitted with stator blades downstream with the aim of converting the rotational energy of the rotor wake into useful thrust for the system. Note that for structural reasons and in order to avoid strange flow behavior on the stator blade approaching $r/R=0$, a hydrodynamic ogive (hub) is fitted at the stator blade’s root.

The complete design of an RDPJ is a relatively complex problem due to the number of components involved and their interactions. In particular, to apply the design approach through optimization it is necessary to parametrize the geometry of all the components involved and their connections. Furthermore, for the optimization process, to produce interesting results, it is necessary to leave as many as possible parameters free to vary. For example, for the rotor blade alone there are distributions of pitch, chord, skew, rake, and thickness, as well as the sectional shape of each hydrofoil describing the blade. By adding the stator, duct, and hub parameters we obtain an optimization problem with a high number of parameters, which requires the analysis of a large number of geometries to converge. Moreover, as suggested in Lugaresi et al. [32], for a rim-driven configuration it is advisable to include the gap geometry between the rim and the internal nozzle surface directly in the design phase to consider its effect on the performance of the thruster. This is because, as shown in Cao et al. [33], there are strong interactions between the pressure fields generated by the rotor and duct, also through the gap. Considering that each geometry is analyzed using an RANS simulation, the calculation time to obtain a final geometry can be high, if not prohibitive from a preliminary design perspective.

For this reason, it makes sense to ask whether it is possible to design the RDPJ components distinctly from each other, obtaining less complex optimization problems (with fewer variables) and consequently becoming more convenient from a computational point of view. In the case of an RDPJ, the element that increases the complexity of the design is the rotor as it requires a high number of parameters to be described and to be optimized.

In this work, then, the possibility of optimizing the duct and the stator shapes keeping the rotor geometry fixed is investigated. It is important to note that even if the rotor is not the object of the optimization, its effect on the flow must be included in the simulations as it is fundamental to correctly align and design the duct and stator under the correct local velocity field modified by action of the rotor itself. In the proposed approach, the rotor geometry is not geometrically resolved in the simulations (Figure 1b) but its effects on the flowfield are simulated using an actuator disk. This simplified approach should allow the stator, duct, and hub to be correctly optimized while guaranteeing at the same time three important advantages:

• It reduces the number of cells in the calculation grid (i.e., the solution time), which instead should be refined in the areas surrounding the complex geometry of the rotor blades;
• It reduces the realization time of the grid, which does not have to be generated around the complex geometry of the rotor blades;
• It does not require the exact definition of the rotor geometry, which may be unknown. Only the open water curve and the torque and thrust distributions, from previous experience or stock propellers, are the required inputs of the model.

Clearly, decoupling the rotor and the stator design is a major limitation of the design since some of the mutual interactions among rotor, stator, and nozzle are neglected. In the case of the design of a stator/rotor device, neglecting the geometry of the rotor would have been impossible, since the risk of cavitation of the rotor is much higher and much more heavily influenced by the modified inflow from the stator blades in front of it. This would have required a dedicated and simultaneous rotor geometry design to maximize the performance gains. Instead, given the type of pumpjet under consideration for this application (rotor/stator) and the aim of this work, this choice can be accepted in a preliminary design phase considering the following:

• This initial investigation is aimed at mainly discussing the role of nozzle performance plus the influence of the thrust added by the stator. For this, an actuator disk model that modifies the momentum accordingly to the required total thrust of the propulsor is sufficient to account for local flowfield modifications to the nozzle and the stator themselves;
• The design of a device working in non-homogeneous inflow (propeller blades in the ship wake or, as in this case of a rotor/stator configuration, stator blades under the action of periodic flow fluctuations) is usually performed assuming the spatial average of the incoming flow as inflow. In this respect, the actuator disk loaded with radial distributions of axial and tangential forces inherently provides the averaged flow to the stator behind it, in a sort of “design for the averaged flow”;
• In order to investigate the importance of reducing the resistance of the nozzle (also thanks to the thrust contribution of the optimally designed stator blades) or increasing the static pressure at the propeller plane as key performance indicators for cavitation inception, which are contrasting objectives, the choice of not modeling the rotor, as already mentioned, grants the computational efficiency needed to address the problem in reasonable computational time.

Reference Geometry

The reference geometry, used both as a performance reference to compare the results, is shown in Figure 2.

It is a rim-driven thruster with a decelerating duct (DRD), made up of the six-blade propeller (z6-4097) obtained in Gaggero [26] and associated with the divergent nozzle of the ducted propeller reference presented in Gaggero et al. [31]. This geometry was designed to avoid the inception of cavitation as much as possible and, in fact, the purpose of the divergent duct is to increase the pressure upstream of the propeller. Figure 3 shows the “inception area”, i.e., the portion of the rotor blades subjected to a pressure lower than the vapor tension at the design functioning point, which is fixed at an advance coefficient close to 1 and at a cavitation index based on the propeller rate of revolution equal to 1.625. Reduction of the maximum blade suction, as better specified in Section 3.3 and of the extension of the pressure region under the vapor tension on the blade surface are the simplified design criteria adopted for the design of this reference thruster. These are needed to address cavitation in a computationally efficient manner to keep, when using RANS, the calculation times within a limit compatible with the large number of simulations to be carried out from an optimization perspective. Thanks to the optimization process, the geometry of the reference rim-driven thruster is designed specifically to work inside the nozzle, which was itself modified by the optimization to achieve the highest possible performance from their mutual interactions. The resulting open water efficiency of this DRD thruster at the design point is ${\eta }_0=0.585$. Cavitation inception is limited to a very small portion of the blade at its leading edge.

3. Methodology

As mentioned in the previous sections, the aim of this work is to apply the SBDO paradigm to optimize the stator, duct, and hub geometries of an RDPJ while keeping the rotor geometry unchanged and described using an actuator disk model. To achieve this, there are three fundamental steps:

• Parametrization of the geometry;
• Set-up of the simulation;
• Set-up of the optimization process.

The parametric representation of the duct, the stator blades, and the hub utilizes B-Spline curves to directly describe the geometrical shape or some parameter distributions like pitch or camber (Section 3.1). In particular, the classical design table serves as the inherent parametric description of the stator blades while in the duct and hub cases B-Splines are used to describe directly their sectional hydrofoil shape. The flexibility and geometric variability offered by B-Splines are unmatched by other parametric methods. The relatively small number of control points required for the parametrization allows for manageable design spaces while ensuring that the design variability remains linked to the expected hydrodynamic shapes.

The simulations to predict the performance of the considered geometry are based on the solution of the RANS equations for an incompressible and turbulent flow. This approach, described in detail in Section 3.2, is commonly accepted today in various engineering problems involving fluid flows as a good compromise between accuracy and computational cost.

Ultimately, the entire process, ranging from geometry pre-processing to post-processing RANS results and monitoring constraints and objectives, has been organized within the modeFRONTIER environment [34] using a genetic optimization algorithm and the workflow presented in Section 3.3.

3.1. Parametric Description of Geometries

The parametric description of the geometry has been realized using the approach already devised for conventional propellers (Bertetta et al. [35]). In the stator blade case, B-Spline curves are used to describe the radial distributions of the canonical quantities adopted in the blade design table. Their control points turn into the free parameters of the optimization process. In this particular case, the chord and the maximum thickness are variable between different configurations but in each geometry are kept constant along the radius. In all geometries, the sectional profile is that of conventional NACA hydrofoils (NACA66 thickness, a = 0.8 camber line), and being less influential on stator performance, rake and skew were not considered in the design. The maximum sectional camber and pitch radial distributions are, each, described by a mean of three points B-Spline curve. An example of such parametrization (e.g., the pitch distribution) is given in Figure 4a. In this way, the stator geometry is described by 8 free parameters.

The hub and nozzle geometries are parametrized through the definition of their two-dimensional sections. The three-dimensional geometries are then generated, exploiting the axial symmetry, through a 360 degree revolution around their axis. Considering its small influence on the overall RDPJ performance, only two parameters were used for the hub. In particular, its dimensionless shape is kept constant in all geometries and only the chord and the maximum thickness are free to vary. Ad hoc control criteria have been forced for the minimum values of these parameters to avoid the generation of geometries with problematic intersections between stator blades and hub.

The duct section is defined using two B-Spline curves, one for the back, defined by 4 control points, and one for the face, defined instead by 6 control points. The two additional points on the face side were needed to ensure a cylindrical internal nozzle surface to accommodate the rim locally. Since the leading- and trailing-edge points of the nozzle sectional profile are shared between back and face sides, and taking into account the vertical alignment of the points at the leading edge to guarantee $C^1$ continuity of the curve plus the horizontal alignment of the central points on the face, the 20 degrees of freedom of these two splines are reduced to only 10 parameters to fully define the shape of the section. A generic duct section geometry with the corresponding control points is shown, as an example, in Figure 4b. Finally, to describe the RDPJ geometry, the number of stator blades (6 to 12) is added as a last parameter; in this way, a 21-dimension optimization problem is obtained.

3.2. RANS Analyses for SDH Optimization

Calculations were performed using RANS for a single-phase, incompressible fluid. The solver StarCCM+ [30] was employed on an unstructured mesh composed of polyhedra, utilizing the Finite Volume method and second-order discretization schemes for all physical quantities. To simplify the analysis, only one blade passage, with local mesh refinements focused on the actuator disk, stator blade, and duct wake (Figure 5a), was examined by leveraging the periodicity of the problem. The computational domain, shown in Figure 5b, features a cross-sectional area at the propeller plane that is one hundred times larger than the area of the propeller disk ($D_{domain}=10D_{prop}$). The inlet is positioned four propeller diameters in front of the propeller, the outlet is located six diameters behind, and uniform inflow is assumed as the initial condition. The $k-\omega$ SST turbulence model is utilized to calculate the turbulent viscosity, with turbulence intensity set at 1% and the ratio of turbulent to molecular viscosity set at 10. The actuator disk model used to simulate the rotor’s effect is the body force propeller method implemented in StarCCM+ [30].

The computational grid was designed according to established guidelines (Lugaresi et al. [32]) and the reference mesh (about 1.3 million cell, $h_1/h_i=1.6$, depending on the specific geometry) was selected based on the verification analysis shown in Figure 6.

The verification study, carried out using the methodology developed by Eça and Hoekstra [36] on a representative geometry, considers the total force (resistance) of the fixed elements of the propulsors (nozzle, stator blade, and hub) as the variable to assess the grid convergence and the uncertainties associated with the numerical calculations. Six simulations were conducted by varying the number of cells from 1.2 to 5.8 million while keeping the mesh topology (distribution of refinements) unchanged. Prism layer cells (total thickness, number of layers, grading) remained constant regardless of the grid. Data are collected as a function of the ratio $h_1/h_i$, which corresponds to the cubic root of the ratio between the number of cells of the finest grid under investigation (₁) and the number of cells of the current grid (_i). The results of the analysis show second-order convergence to the extrapolated value (i.e., for “infinitely dense mesh”) of the resistance ($R_{t_0}=65.41N$). The uncertainty on the chosen mesh is relatively high (about 7.5%). However, the error with respect to the extrapolated value, is lower than 2%. These values are considered acceptable, especially considering the fact that increasing the number of cells, even significantly, does not lead to a relevant reduction in the uncertainty and in the error with respect to the extrapolated value. On the contrary, reducing the calculation time of a single simulation in an optimization problem allows a greater variety of cases to be explored in the same amount of time, making this choice of grid the necessary compromise for a usable design tool.

For fair comparisons of the geometry performance under the action of the actuator disk throughout the optimization process, the only constraint is that the thrust obtained from each tested geometry is equal to the net thrust of the reference configuration. This translates into the condition:

$$T_{vd}+T_{duct}+T_{stator}+T_{hub}:=T_{net}=T_{ref}=T_{rotor}+T_{refDuct}$$

(1)

where $T_{vd}$ is the actuator disk thrust, $T_{duct}$ is the duct thrust, $T_{stator}$ is the stator thrust, $T_{hub}$ is the hub thrust, $T_{net}$ is the net thrust of the entire system (defined as the sum of the thrust of all elements of the RDPJ), and $T_{ref}$, $T_{rotor}$, and $T_{refDuct}$ are, respectively, the net thrust, the rim-driven rotor thrust, and the duct resistance of the reference configuration at design working point. Note that for consistency, all forces are expressed as thrust, but typically $T_{duct}$ and $T_{hub}$ are resistance and therefore they have, using this notation, negative values. The actuator disk model is driven by thrust, i.e., it inserts a momentum force, into the specified volume, equivalent to the thrust imposed by Equation (1). Here, $T_{ref}$ is computed for the reference geometry, and remains the same for all simulations, while $T_{duct}$, $T_{stator}$, and $T_{hub}$ are computed every iteration by integrating the pressure field and the viscous forces on the actual surfaces subsequently built by the optimization tool.

For a more consistent inclusion of the action of the rotor and a more accurate prediction of the performance of the RDPJ, especially the effect on the stator, it is necessary that, not only are the forces in equilibrium in an integral manner, but also that the flow from the actuator disk is as similar as possible to that of the real rotor the actuator disk is replacing. For this reason, the radial distributions of torque and thrust were extracted from a preliminary simulation of the reference rotor and then used to feed the actuator disk. As the geometry of the duct and stator varies, however, the overall resistance of the propulsion system changes. For example, if the total resistance decreases, to maintain the same net thrust it is necessary that the rotor thrust also decreases. Since the rotor geometry is fixed and not geometrically adapted to provide the required thrust, this translates into reducing its revolution rate. Conversely, when configurations with higher resistance than the reference case are tested by the optimizer, the rotor will have to increase the revolution rate. The calculation of the equilibrium point in terms of torque and revolution rate of the equivalent “propeller” as the geometry changes and the thrust varies is possible since the actuator disk is driven by the open water diagram of the reference rotor.

In particular, the $K_t$, $K_q$, and ${\eta }_0^{*}$ curves as a function of $J^{*}$ are passed to the actuator disk model. The superscript * indicates that the velocity used for the performance coefficients’ calculation of the rotor is not the free-stream (as for usual open water) but is the mean velocity extracted in correspondence of a sampling disk, which has a diameter equal to that of rotor D and is located at $10\%D$ upstream from the rotor disk. So, $J^{*}$ represents the advance coefficient calculated as:

$$J^{*}={\displaystyle \frac{U_{10\%D}}{n_{ref}D}}$$

(2)

where $U_{10\%D}$ is the averaged x-component of the velocity on this extraction plane. From the reference simulation, varying the free-stream velocity, it is possible to compute the open water performance as:

$$K_t={\displaystyle \frac{T_{rotor}}{\rho n_{ref}^2{D}^4}}$$

(3)

$$K_q={\displaystyle \frac{Q_{rotor}}{\rho n_{ref}^2{D}^5}}$$

(4)

$${\eta }_0^{*}={\displaystyle \frac{J^{*}}{2\pi }}{\displaystyle \frac{K_t}{K_q}}$$

(5)

where $T_{rotor}$, $Q_{rotor}$, and $n_{ref}$ are, respectively, the rim-driven rotor thrust, torque, and the revolution rate of the reference configuration. At the end, it is possible to obtain the * open water curves shown in Figure 7, which are the values of $K_t$, $K_q$, and ${\eta }_0^{*}$ as a function of this modified advance coefficient $J^{*}$, and use them to feed the actuator disk model.

This procedure allows us to decouple the reference rotor from its duct. In a traditional open water diagram obtained using the advance coefficient calculated with the free-stream velocity, the duct effect (i.e., local change in the velocity field, then further modification of the flow seen by the rotor additional to the usual self-induced velocities) is inherently included in the diagram. To isolate, as much as possible, the behavior of the rotor, its performances have to be computed considering the actual rotor inflow velocity. In this way, it is possible to capture at least some of the influence that different ducts may have when their geometry varies during the optimization. However, it must be kept in mind that when using this simple body force propeller model, the effect of the duct on the rotor is seen simply as a change in the flow rate; no information about the angle of attack of the blades is considered. To overcome this limit, one could use a more sophisticated actuator disk model based on blade element method (BEM) that considers, for the momentum source calibration, locally the effective angle of attack resulting from the modified nozzle action. This, for sure, could improve the accuracy of the interaction since it can more accurately take into account how the radial distributions of rotor blade load (and then induced downstream velocities) are affected by the different inflow given by the modification of the nozzle. An approach of this type would be mandatory to perform an analysis of a stator/rotor pumpjet configuration, in which the variation in angle of attack induced by the stator stage on the propeller blades represents the main effect and must therefore be captured. On the other hand, on a rotor/stator layout and, in particular, since the design of the stator, as already mentioned, is addressed in the “averaged inflow condition”, small local modifications of the inflow should not significantly alter the solution. Given the focus of the design process, thus, only a basic actuator disk model based only on the 1D conservation of momentum has been employed in the proposed optimization process.

3.3. Optimization Strategy, Constraints, and Objectives

Beyond the parametric description of the geometry, which requires customization based on the reference geometry, the proposed optimization-based design process is basically the same used for the rim-driven optimization in Gaggero [26]. As mentioned in the introduction, the unique aspects of this configuration and the necessity to accurately account for blade/nozzle interactions underscore the need for a viscous solver. In this study, the computational set-up was significantly simplified with respect to the previous works by using an actuator disk in place of the propeller, as the focus was on how nozzle shape and downstream stator blades influence the internal pressure field and the performance of the thruster. Three objectives were prioritized in the design: minimizing the (risk of) cavitation on the stator and rotor blades while simultaneously minimizing drag of the entire RDPJ. The minimization of the nozzle resistance (also through an increase in the thrust provided by the stator) is considered as a design objective since any change in the force paid by the duct corresponds to a different, eventually more favorable, loading of the rotor. In the end, moreover, a reduction in drag can be correlated to a change (increase) in the total propulsive efficiency of the system. Here, the problem is that the rotor geometry is not resolved so the cavitation on its blade is hard to assess. Cavitation on the rotor blades depends mainly on three factors: the geometry of the rotor, the revolution rate of the rotor, and the pressure increase generated by the duct on the rotor plane. As already mentioned, the analysis conducted aims at optimizing the elements of the RDPJ by decoupling them from the rotor. The first two factors, therefore, cannot be considered and only the third can be included in the optimization process. For these reasons, the “delay in cavitation” properties of the duct is estimated thanks to the average value of the pressure immediately upstream of the rotor, and this is a quantity that must be maximized. Cavitation on the stator is, instead, controlled by maximizing the pressure on the blades.

The main difference compared to previous works is that in this case the net thrust of the system is not an output of the computation but is forced into the simulation by modifying the thrust provided by the actuator disk to compensate for the change in the overall resistance of the system due to the variation in the nozzle and stator geometries. In this case, geometries with high resistance that would be acceptable, are avoided by the optimizer thanks to the objective of minimizing resistance. This approach is preferable with respect to imposing the working point of the rotor and discarding the geometries that do not fall within a thrust tolerance, for mainly two reasons. The first and most fundamental reason is that highly performing duct and stator geometries, which would therefore produce higher thrusts than the reference geometry, would be discarded. In this case, one could think of limiting the thrust only in one direction, by imposing that it satisfies a minimum equal to the reference geometry thrust and then requesting the maximization of the overall thrust of the system as an objective. This approach, which is in principle correct for obtaining the best configuration, has the flaw of not preserving the thrust, which is typically the main parameter required by a thruster. To ideally keep all the designed thrusters usable on a hypothetical reference ship, the first approach was therefore chosen.

A simple flowchart outlining the optimization process is presented in Figure 8.

Since geometrically unfeasible designs are not a concern, given that the geometric constraints are directly embedded within the ranges of variation of the free parameters, all cases must be evaluated in terms of performance to inform the optimization process. After selecting the free parameters (point 1); generating the 3D model of the stator blade, duct, and hub (point 2); and arranging the mesh (point 3), the solution of the RANS equations is obtained (point 3). The RANS simulation is set up to run for a minimum 400 iterations that, in the light of the convergence trends observed for the reference geometry, are sufficient to complete the numerical transient and to preliminarily assess with a certain accuracy the propulsive performance of the configuration. After it, a stopping criterion on the asymptotic value of resistance is enabled and a limit of 800 iterations is set to avoid bad non-converging geometries that waste computational time. From the solution, the values of resistance and the pressures on the stator blade and on the rotor disk are computed and passed to the optimizer (point 5), which select a new set of parameters.

The RDPJ propeller design approach can thus be summarized as a bound-constrained optimization problem. Identifying any candidate as the set of free parameters that describe its geometry $\mathit{p}=[p_1,p_2,\dots ,p_i,\dots ,p_{21}]$, the problem is stated as follows:

$$\{\begin{array}{c}minimize\left(R_t\left(\mathit{p}\right)\right)\hfill \\ minimize\left(cavitation\right(\mathit{p}\left)\right)\hfill \\ maximize\left(P_{up}\left(\mathit{p}\right)\right)\hfill \\ \mathrm{subject}\mathrm{to}:min\left(p_i\right)\le p_i\le max\left(p_i\right)\hfill \end{array}$$

(6)

where

$$R_t=-\left(T_{duct}+T_{stator}+T_{hub}\right)$$

(7)

and $P_{up}$ is the mean pressure extracted at 10% D upstream the rotor from a radial position of $r/R=0.5$ to $r/R=1$. The choice of the sampling plane for the pressure was made by analyzing the pressure variation upstream of the rotor in different duct configurations. A typical pressure trend for $P_{up}$ is shown in Figure 9.

We note that the pressure reaches its maximum near the leading edge of the nozzle (stagnation point) and then drops and reaches a minimum near the actuator disk due to its action (suction in front of it) on the flow. Inside the disk, the pressure then rises due to the volume force imposed by the body forces. Given this regular trend of the mean pressure from nozzle leading edge to the actuator disk, in a comparative process like optimization any extraction position would have been appropriate to estimate how the modification of the nozzle shape influences the internal pressure field and then delays or anticipates the cavitation inception. The sampling at a position 10% D upstream of the propeller disk, where the pressure was observed to reach its minimum, was then selected.

While predicting the total resistance of the system is relatively straightforward, cavitation occurrence on the stator presents challenges, particularly from a computational perspective when using RANS. In theory, cavitation, especially leading-edge sheet phenomena, can be accurately predicted using multiphase approaches with a homogeneous mixture assumption, and there are many successful examples in the literature. In the progressive process that will lead to a final RDPJ design (and the best-balanced design tool from a computational and an accuracy point of view), multiphase calculations are planned. These will be of particular importance when the shape of the rotor is considered, because the rotor is expected to have the most severe issues in terms of cavitation.

For this particular application, however, considering the relatively low advance coefficient and the consequently small forces from stator blades, it is expected that stator geometries well-aligned with the streamlines from the rotor (and then with very low risk of cavitation occurrence) will be identified by the optimization process. Consequently, no particular cavitation problems are expected.

For these reasons, mostly confirmed by the final results of the optimization in which only few geometries exhibit marginal cavitation on the stator, instead of minimizing cavitation as the area of the stator blade covered by vapor, the objective focused on the risk of cavitation inception, which translates to maximizing pressure (or reducing suction peaks) on the stator blades. Cavitation inception is therefore assessed by monitoring suction peaks at the stator blade leading edge—indicative of sheet cavitation risk—on both the suction and pressure sides, as well as the pressure value at mid-chord to account for bubble cavitation risk. For this analysis, only computationally efficient steady, non-cavitating calculations are required. Moreover, also the extent of the low-pressure (lower than vapor tension) region on the stator (the “inception area”) is considered, and in particular is minimized. Thus, minimizing cavitation can be viewed as maximizing the pressure data extracted, at convergence, from the stator surface while minimizing the low-pressure area. To achieve this, the stator blade surface is divided into four zones where pressure is collected and processed separately after each case. In a design focused on avoiding cavitation, the goal is to minimize the highest suction peak across the entire stator blade. Collecting data from various locations, however, allows for a more detailed characterization of the phenomena occurring on the stator, which could be valuable for post-processing when selecting the “optimal” design among the Pareto configurations.

To prevent “dirty” data from affecting pressure peak estimations (e.g., data from warped, skewed, or low-quality cells, particularly at the leading edge), the highest suction values collected during the design process are averaged for each zone. The reliance on non-cavitating analyses to characterize RDPJ performance implies that the thrust is produced under non-cavitating conditions, as thrust breakdown from severe cavitation could make unfeasible interesting geometries. Since the design aims to minimize cavitation, such occurrences should be avoided.

Ultimately, the optimization process described in Figure 8 and Equation (6) consists of a constrained multi-objective optimization problem with ten design objectives (minimize total resistance, maximize the averaged pressure inside the nozzle and in front of the rotor disk, maximize the pressure and minimize the low-pressure areas on four different stator blade locations) and twenty-one free variables. This is managed using a genetic algorithm within the modeFRONTIER environment. The initial population for each design run comprises 210 configurations generated using a Uniform Latin Hypercube sequence. Twenty generations are evolved, employing elitism, mutation, and direct crossover among candidates, resulting in a total of 4200 analyzed geometries.

4. Results

The simulation-based design by optimization was used to design the optimal combination of duct, stator, and hub given a fixed rotor geometry, satisfying the problem statement expressed in Equation (6). All geometries obtained from the optimization process are shown in the $R_t$ (total resistance) vs. $P_{up}$ (the radial mean pressure measured at 10%D upstream the rotor disk) diagram shown in Figure 10.

The color-code of each configuration shown in the diagram represents the case ID and, thus, the progress of the optimization process toward convergence. We can see that while initially the exemplars turn out to be equally distributed within the space of solutions, advancing with generations results in increasingly clustered results. We also note that the results do not accumulate along the entire Pareto frontier but tend to remain in a well-defined area. The clustering of the solutions in the “central” part of the Pareto is the usual behavior of a multi-objective algorithm that concentrates the investigation, predilecting solutions (and pushing the choice of subsequent generations) closer to the “utopia” optimal point (i.e., the ideal configuration that minimizes all the contrasting objectives of the design). In this case, the clustering zone is not exactly symmetrical because resistance and over-pressure are not the only objectives of the optimization. In fact, low resistance and over-pressure correspond to faster flow and, consequently, a higher probability of cavitation on the stator. The goal of minimizing cavitation on the stator, in this case, comes into play by “shifting” the clustering to higher-resistance solutions.

The clustering of configurations on the Pareto is also a qualitative indicator of optimization convergence and helps determine the number of generations at which it is possible to stop the process. Unfortunately, assessing the convergence of a multi-objective optimization, especially when the number of objectives is high, is never easy. In this case, testing 20 generations is a practical choice since no substantial improvements were recognized during the last few generations. Moreover, it is a limit of the computational efficiency of the process that, regardless the simplifications introduced, requires several days (weeks) of calculations.

The same set of data is analyzed, in Figure 11, in terms of influence of “global” geometric parameters of the duct, such as its chord, the angle of attack, or the maximum sectional camber, on the key performance indicators (resistance and over-pressure in front of the disk) selected as objectives of the design. This analysis, even if limited, establishes a sort of sensitivity of the objectives to the most representative geometrical parameters of the design. As expected, both drag and over-pressure are observed to decrease at low angles of attack and camber. On the other hand, it can be observed that the chord has a substantially negligible influence on the over-pressure, while, as it increases, it increases the resistance, probably due to viscous effects on larger wetted areas.

In Figure 12, moreover, some ducts’ sectional geometries picked from the Pareto frontier are shown and compared each other to sustain the exploration capabilities of the optimization process when a sufficiently wide space of the design parameters is used. Although the aim of the optimization is a divergent duct RD, which theoretically is the obvious choice to increase the static pressure inside the nozzle and postpone the risk of cavitation, some ducts’ sections are closer to a convergent configuration. This happens because the divergent-duct condition has not been forced geometrically but, as might be expected, is achieved naturally as the configurations evolve toward the objectives. In order to facilitate the convergence, the ranges of parameters were chosen to explore preferably divergent nozzle geometries; however, slightly convergent configurations were not ruled out a priori. As a consequence, looking at Figure 12 we note that, as expected, such geometries exist as feasible cases of the optimization and that they are correctly positioned in the low resistance and low over-pressure part of the Pareto diagram.

Once only the geometries with zero cavitation area on the stator were considered, and based on the comparison with the performance of the reference rim-driven geometry in terms of total resistance and pressure upstream the rotor (the red square in the Pareto diagram), six different configurations were selected with the following criteria:

• The geometry with equal resistance of the reference but with maximum pressure upstream of the rotor: ID 2961 (blue star);
• The geometry with equal pressure downstream of the rotor compared to the reference but minimum resistance: ID 2932 (purple star);
• A geometry on the Pareto boundary with lower resistance and higher upstream pressure than the reference geometry: ID 4314 (green star);
• A geometry with resistance slightly higher than that of the reference configuration but pressure upstream of the rotor significantly higher: ID 968 (brown star);
• Two geometries with much lower resistance and slightly lower upstream pressure than the reference: IDs 3112 and 2426 (yellow and orange stars).

These geometries, which were all chosen on the Pareto front of geometries that avoid any risk of cavitation on the stator, all represent “optimal” and not-dominated configurations among which the optimization process, for as set, it is not able to choose. For each of these solutions, indeed, there is not the possibility to improve one of the objective without worsening the other.

Only a detailed critical analysis conducted retrospectively, then, can help understanding which geometries should be preferred. As expected, the two objectives imposed are conflicting. In particular, we can observe that the geometries with lower resistance are those that also have a lower pressure upstream of the rotor. To have low resistance, the duct must be aligned as much as possible with the flow, avoiding strongly divergent geometries and therefore preventing large increases in pressure. In contrast, geometries that present higher pressures on the rotor disk are associated with nozzles with more strongly divergent geometries and therefore higher resistance.

In addition to the duct effect, it must be considered that all the geometries analyzed, given how the optimization problem was set up, give the same net thrust. For this reason, since the total resistance varies from one geometry to another, the total thrust of the rotor changes as well. Because the rotor geometry is fixed through the open water curves that feed the actuator disk, this implies that the propeller is forced to work at a different number of revolutions. Two contrasting effects then emerge on the determination of cavitation on the rotor. On one hand, the effect of the divergent duct, which by increasing the pressure on the propeller disk should reduce the cavitation (geometries on the right of the Pareto diagram). On the other hand, a nozzle that creates less over-pressure also produces less resistance and, consequently, requires less thrust from the rotor, which is then unloaded and free to work at a lower rps to guarantee the required net thrust. Operating at a lower number of revolutions implies an unloaded blade and therefore a less accentuated low-pressure on the suction side, with a consequent reduction in the risk of cavitation. It is therefore necessary to investigate which of the two effects is more effective for the objectives of the design: whether it is better to have a rotor working at lower rps and a duct that produces little over-pressure, or whether it is better to accept a higher number of revolutions because the suction generated by the blade is compensated by the increase in pressure provided by the more divergent nozzle.

To investigate this aspect, the six selected geometries were further investigated using unsteady RANS analysis with the resolved rotor geometry and truly rotating mesh via sliding interfaces. For the comparison to be meaningful, the unsteady simulations had to provide the same net thrust obtained for each analyzed geometry from the actuator-disk simulations, i.e., the net thrust of the reference geometry. It was therefore necessary to determine the correct propeller revolution rate of the fully resolved RANS simulations to feed the rotor in the unsteady simulations for each geometry. In order to do so while limiting the calculation times as much as possible, steady RANS analyses were carried out with a resolved rotor by exploiting both the periodicity of the domain and the equations written in a rotating reference frame (MRF). Starting from the revolution rate estimated during the optimization process by the actuator disk model, the number of revolutions at equilibrium was determined iteratively for each geometry using the following algorithm, which usually achieved convergence in fewer than three iterations:

• A steady RANS simulation with the rotor geometrically resolved is carried out by imposing the number of revolutions equal to the value $n_0$ obtained from the actuator disk simulation during the optimization;
• The value of the total net thrust $T_{ne{t}_0}$ is calculated and compared with that imposed during optimization $T_{net}$;
• The next value of the revolution rate is obtained by increasing or decreasing $n_0$ by a fixed value of $0.1$ rps; in particular

  $$n_1=\{\begin{array}{c}{n}_0+0.1\mathrm{if}:T_{ne{t}_i}<T_{net}\hfill \\ n_0-0.1\mathrm{if}:T_{ne{t}_i}>T_{net}\hfill \end{array}$$

  (8)
• A steady RANS simulation with the rotor geometrically resolved is carried out by imposing the number of revolutions equal to the value $n_1$ and the value of the total net thrust $T_{ne{t}_1}$ is calculated;
• The next value of the revolution rate $n_2$ is calculated by linear interpolation or extrapolation depending on the value of the obtained thrusts; in particular

  $$n_2=n_1-{\displaystyle \frac{n_1-n_0}{T_{ne{t}_1}-T_{ne{t}_0}}}(T_{ne{t}_1}-T_{net})+k$$

  (9)

  where $k=0$ if $T_{net}$ results between the values of $T_{ne{t}_1}$ and $T_{ne{t}_0}$, $k=-0.1$ if $T_{net}$ results greater than both $T_{ne{t}_1}$ and $T_{ne{t}_0}$, and $k=0.1$ if $T_{net}$ results smaller than both $T_{ne{t}_1}$ and $T_{ne{t}_0}$.
• The previous step is repeated using the closest values to $T_{ref}$ among those calculated for the interpolation, until the obtained $T_{re{f}_i}$ differs from $T_{ref}$ by less than $1\%$.

Once obtained, the revolution rate at equilibrium was used to feed the unsteady simulations. Since the thrust obtained from the unsteady analyses using the rate of revolution from the MRF-based iterative process resulted equal to the required thrust of the propulsor within a tolerance always lower than 1%, no additional iterations using the unsteady calculations were necessary.

The comparison between actuator disk data from the optimization process and fully unsteady analyses provides a posteriori validation of the design assumption and permits discussing correlations between forces as well as the importance of rotor unloading (i.e., reduction in nozzle resistance) or pressure increase on the cavitating performance of the rotor through the direct evaluation of the pressure field (and then, the cavitation risk) on the rotor blades.

The correlation diagrams of resistances and thrusts obtained with the simplified simulation set-up adopted for the optimization (actuator disk) and the unsteady RANS are shown in Figure 13 and Figure 14.

From the correlations shown in Figure 13, it can be observed that the set-up based on the actuator disk model works well to predict the total resistance of the system. The forces developed by the individual elements are also predicted quite accurately when using the actuator disk. The stator resistance is the only quantity that presents a fairly significant difference between the simulation with the actuator disk and the unsteady RANS. However, we note from Figure 14a that this difference is an approximately constant offset for the different cases, which worsens the performance of the stator in the actuator disk simulations. This is not a major problem for several reasons. Firstly, the difference in absolute terms is small compared to the total resistance of the system that, anyhow, is still well predicted. Furthermore, the fact that the offset is approximately constant implies that in relative terms the best cases are correctly identified and that their absolute performances are simply actually underestimated compared to those predicted with the more accurate simulation set-up.

Good correlations among forces mean good correlations among pressure fields, as it can be observed in Figure 15 where the pressure distributions on the stator blades from actuator disk analyses are compared to the unsteady RANS. With the exception of local differences associable to the interaction with the resolved blade in the case of unsteady RANS, since the pressure field is not averaged over an entire rotor revolution but taken from an instantaneous relative position between rotor and stator, the most relevant features of the pressure fields are very similarly predicted. This is representative of the very similarly predicted functioning of the stator, as already discussed when forces were addressed, and represents a further proof of the equivalence. It also supports the usability of different modeling strategies (the actuator disk, in present case) of the rotor action on the flow for design purposes.

The predicted risk of cavitation, identified by the portion of the blade surface at a pressure lower than the vapor tension from the resolved rotor in the case of the unsteady RANS analyses is shown in Figure 16.

Comparing the results of the six selected geometries, we note that configuration IDs 3112, 2932, and 4314 are not affected at all by cavitation. Geometry IDs 2961 and 968, instead, present a risk of cavitation both on at the leading edge of the rotor blade and on the nose of the duct. The risk of cavitation occurrence on the duct indicates strong suction associated with a highly divergent nozzle. This characteristic, as already mentioned and known, causes an increase in resistance of the duct and, consequently, requires a higher thrust from the rotor to deliver the requested net thrust. One could therefore expect that the final consequence is that the rotor is forced to work at a higher number of revolutions to satisfy the thrust request and that this condition increases the suction on the back of the blade. This consequently generates the area at risk of cavitation. In reality, this is not the case as it is observed that the rate of revolution at equilibrium for the two geometries considered is lower than in the others. This can be explained by considering the angle of attack of the incident flow. In fact, since the nozzle is very divergent, the pressure immediately before the rotor increases and the speed decreases, as seen in Figure 17. This causes an increase in the angle of attack seen by the rotor itself, and it is precisely this increase in angle of attack that is the cause of the strong suction on the leading edge of the blade, which generates the area at risk of cavitation and, at the same time, provides the required thrust at a lower rate of revolution.

The situation is different for ID 2426. In this case, the resistance of the duct is very low, indicating a well-aligned and not particularly divergent duct. Consequently, the velocity is not slowed upstream of the rotor and the angle of attack does not increase accordingly. The fact that the angle of attack remains low implies that the rotor must rotate faster to guarantee the required thrust. This increase in the revolution rate may cause higher suction, especially in the areas near the tip where the effect of a higher number of revolutions is mostly transformed into an increase in incident speed seen from the blade. In this case, the slight cavitation on the nozzle is due to the fact that the increase in revolution rate causes a greater suction of the rotor and, consequently, a slight acceleration of the fluid on the duct nose with a consequent increase in suction.

Based on the unsteady RANS analyses of the optimized geometries, the risk of cavitation appears to be influenced by two factors: on the one hand the increase in the number of revolutions in cases where, due to the high incident velocity, the angle of attack of the rotor is too low to produce the required thrust; on the other hand, by the large increase in the angle of attack given by the highly divergent ducts. Both are consequences mainly of the nozzle shape, which changes the inflow and the loading of the rotor. This discussion suggests that accounting for cavitation on the rotor is possible in a consistent way only by including its geometry in the optimization and analysis process. However, the developed tool is able to select stator, nozzle, and hub configurations that somehow optimize the efficiency of the system, with the total resistance directly related to the efficiency obtained with the unsteady RANS simulations. As expected, the lower the total resistance, the higher the efficiency. Compared to the initial configuration, ID 2426, 3112, 2932, and 4314, which all have total resistance of nozzle, hub, and stator lower or very close to the reference geometry, have higher efficiency while maintaining (ID 2426) or improving (ID 3112, ID 2932, ID 4314) the cavitation performances of the rotor. The nozzle, as seen in the correlation diagrams of Figure 13 and Figure 14, is the major contributor to the total resistance, with the hub almost negligible, and any increase in the thrust delivered by the modification of the stator is more than balanced by a significant increase in the resistance of the corresponding nozzle.

A reasonable design strategy is, then, to look primarily for solutions with low nozzle resistance. For example, ID 2426 presents a duct resistance (which also translates into a total resistance) equal to approximately half compared to the reference configuration, at the price of a reduction in pressure upstream of the rotor of only 1000 Pa (approximately 5%), which is not sufficient to cause a significant increase in cavitation. ID 3112, which based on the combined results from the optimization and the unsteady RANS can be considered the best choice, is characterized itself by a total resistance equal to approximately half compared to the reference geometry. The resulting efficiency, which in any case is that given by a rotor working in off-design conditions due to the changes in loading, incoming flow, and rate of revolution, is equal to 0.5912, which is about 1% higher than the reference rim-driven geometry. The cavitation risk is completely nullified, showing the effectiveness of the design approach since simultaneous improvements in cavitation risk and efficiency were achieved. ID 2932 and ID 4314 nullify the risk of cavitation as well, but their calculated efficiencies using the unsteady RANS are practically identical to that of the reference geometry (0.58 and 0.585, respectively) as a consequence of total resistances of nozzle, hub, and stator only marginally lower to that of the reference case. ID 2916 and ID 968, which were selected to investigate the importance of substantial increases in the internal pressure as a possible counteraction to cavitation inception, instead, pay the substantial increment in resistance associated to very diverging nozzle geometries: the overloaded rotor blades, working in very off-design conditions at a substantially higher angle of attack, show a risk of cavitation similar to that of the reference propeller at efficiencies (0.567 and 0.558) much lower. Based on these unsteady, fully resolved, RANS analyses, which allowed us to assess the efficiency of the propulsors and the risk of cavitation of the rotor, it seems that under the given simplifications and assumptions for the design of nozzle and stator of RDPJs, reducing the drag of the nozzle is overall more effective than increasing the static pressure at the propeller disk. A reduction in resistance allows for an unloading of the rotor blade, with effects on local pressure distributions and functioning points favorable for cavitation postponement. This effect seems dominant with respect to the relatively small increases in the static pressure inside the nozzle achievable by exaggerating the divergent nature of the duct and, consequently, the cost in terms of additional resistance to be recovered by overloading the rotor.

5. Conclusions

In this work, exploiting an SBDO paradigm, the optimization of the stator, duct, and hub was conducted for a rim-driven pumpjet-type thruster. The proposed test case is a tentative one to decouple the design of the fixed parts of the propulsor (nozzle, hub, stator) from the rotor by using simplified models and criteria for controlling cavitation and efficiency of the system. Moreover, the proposed analyses allow us to study, after detailed unsteady RANS calculations of the optimal geometries identified by the design process, which features mostly influence the performance of the design.

The proposed approach makes use of an actuator disk to replace the action of the rotor blades on the flowfield. This choice, which validity for the design only of nozzle and stator blades has been assessed through detailed comparisons of predicted forces and pressure fields with the results from fully resolved RANS analyses, was dictated by several reasons. Primarily, there is the need of a computationally efficient tool well-suited for an optimization-based design process. Since the focus was on stator blades and the nozzle of a rotor/stator configuration, the momentum sources of the actuator disk were more than appropriate to model the average influence of the rotor wake on the stator blades as well as of the resulting streamlines for the definition of the optimal nozzle shape. Avoiding the discretization of real rotor blades allowed us to save computational time in terms of both reduction in the grid size and use of more efficient steady solvers in place of moving reference frames and iterative coupling between rotating and fixed parts of the computational domain. This permits analyses in less than 3 min for each geometry compared to about 8 min in the case of a resolved rotor with a moving reference frame, and a few hours in the case of fully unsteady analyses. A simple scaling of the momentum source strength, moreover, permitted us to easily satisfy the design constraint of fixed delivered total net thrust by the propulsion system in place of an iterative adjustment of the rate of revolution of the given rotor geometry (used instead for the verification by unsteady RANS), or of a modification of the rotor geometry (which is the problem to be addressed when rotor, stator, and nozzle are designed simultaneously) to realize the required thrust.

The criteria of the design were identified by the risk of cavitation of the stator blades (which was monitored through the pressure directly sampled on their surface), the efficiency of the propulsor system, and the risk of cavitation of the rotor blades. Without modeling the rotor, propulsive efficiency and risk of cavitation of the rotor were replaced by resistance of the fixed parts of the propulsors (nozzle, hub, stator), to be minimized, and static pressure inside the nozzle, to be maximized. In this respect, the results of the optimization activity, enhanced with unsteady RANS analyses of a set of optimal geometries, also permitted us to discuss the opportunity to unload the rotor by reducing the resistance of the nozzle, or to increase the internal pressure by exacerbating the divergent action of the duct, for cavitation avoidance and efficiency increase.

Collected results evidenced the importance of the reduction in resistance of the nozzle, from which both efficiency and cavitation gained advantage. At fixed rotor geometry (potentially in non-optimal functioning conditions as a consequence of inflow modifications), any reduction in resistance corresponded to an increase in efficiency but also to a reduction in the risk of cavitation thanks to the significant unloading of the rotor blades. The same was not observed with the increase in static pressure, which was never sufficient to counteract the excessive blade overloading due to the sudden increase in nozzle resistance. Given this evidence, it was possible to identify a certain number of configurations with performance, both in terms of efficiency and cavitation avoidance, improved with respect to the reference rim-driven thruster, proving in the end the effectiveness of the proposed design approach at least on a comparative level. In particular, it is the case of ID 3112 which, thanks to a substantial reduction in the nozzle resistance, permitted (unsteady RANS analyses) simultaneously an increase in efficiency (about 1%) and the reduction in the risk of cavitation, that is completely nullified, or of ID 2426 and ID 2932 that nullified the cavitation as well without worsening the efficiency of the reference test case.

The successful application of this very preliminary optimization-based design process to rim-driven pumpjets, for now limited to nozzle and stator geometries, opens interesting possibilities for even more reliable procedures that would encompass the inclusion of the actual rotor geometry and the use of more sophisticated numerical models for the assessment of the key performance indicators of the design. Expected steps are the inclusion of the actual rotor by using a blade element method, which may improve the interaction with the stator (and provide quantitative information on the performance of the rotor itself) thanks to a local characterization of the flowfield based on the real inflow due to the nozzle, or by the use of a fully parametrized rotor model. The latter, in particular, would allow for a combined rotor/stator design, optimizing the rotor shape for the current inflow and required thrust, plus a direct characterization of the risk of cavitation not simply based on correlations between expected blade loading or static pressure inside the nozzle. From the point of view of numerical solvers, but only with the aid of surrogate models given the required computational effort, the application of multiphase flow models could allow for a more robust estimation of the cavitation (and of its consequence in terms of thrust and efficiency) on rotor, stator, and nozzle, replacing the simplistic assumption of cavitation inception based on pressure distributions.