Design Method for Low-Ice-Class Propellers Based on Multi-Objective Optimization

Abstract

The objective of this paper was to establish a comprehensive methodology for the optimized design of propellers for ice-class vessels, aiming to enhance hydrodynamic efficiency while ensuring structural integrity. This paper begins by introducing a novel approach for calculating blade stress, which takes into account both extreme ice loads and hydrodynamic loads, to be utilized in the propeller strength design process. Subsequently, a backpropagation (BP) neural network model was developed based on the data obtained from B-series propeller charts and integrated with a genetic algorithm to achieve a preliminary optimized design of the propeller’s hydrodynamic performance. To illustrate the application of this methodology, a case study of an ice-breaking tug propeller design is presented, detailing the optimization design process, including the preliminary, intermediate, and final design stages. The study also addresses key aspects such as geometric parameterization, the selection of optimization variables, the implementation of optimization algorithms, and the balance of multi-objective trade-offs. The proposed design approach can serve as a valuable reference for the practical engineering design of propellers for ice-class vessels, providing a systematic framework for achieving optimal performance in challenging operating conditions.

1. Introduction

Propellers, as essential propulsion devices for ships, play a crucial role in determining a vessel’s speed, economic efficiency, and maneuverability through their hydrodynamic performance. When navigating in ice-covered waters, propellers face severe ice loading conditions. Optimizing the design of ice-class propellers to balance performance and strength is vital for ensuring the safety and efficiency of ships operating in such environments.

In terms of hydrodynamic design, early prediction methods for conventional propeller performance primarily relied on potential flow theories, such as the lifting line theory and lifting surface theory. Although these methods are computationally efficient, they struggle to capture the effects of viscous flow. With the development of computational fluid dynamics (CFD), methods solving the Reynolds-averaged Navier–Stokes (RANS) equations have gradually become mainstream. Kerwin [1] were among the first to apply CFD methods to propeller hydrodynamic analysis, achieving promising results. Subsequently, Mishima [2] and Chen [3] employed RANS methods to numerically simulate propellers, revealing the complex flow characteristics in the propeller wake. To further improve prediction accuracy, Bosschers [4] proposed a coupled solution strategy based on the boundary element method and RANS method. In recent years, high-fidelity methods such as the large Eddy simulation (LES) and direct numerical simulation (DNS) have been applied to propeller hydrodynamics research [5], providing new means to deepen the understanding of the underlying mechanisms [6,7].

Regarding propeller structural strength prediction, early researchers primarily relied on empirical formulas and simplified theoretical methods, such as the beam theory and the finite element method (FEM). Søntvedt [8] was one of the first to analyze propeller blade strength using the beam theory. Sunnersjo [9] employed FEM to conduct strength analysis on highly skewed propeller blades, obtaining stress distribution characteristics. The choice of element type and mesh quality significantly influence the accuracy of the calculations. Using higher-order elements, appropriate mesh density, and layout are the key to obtaining reliable results. As computer performance has improved, high-precision numerical simulation methods have become widely used. Hu Zhikuan [10] performed strength analysis of propellers under ice loads using FEM based on empirical ice load formulas and five typical calculation conditions. He also established an ice propeller collision calculation model to study the variation in contact forces and the dynamic response of the propeller during the collision process. Young [11] was among the first to apply a one-dimensional fluid–structure interaction (FSI) method to analyze propeller blade loads and deformations. Ye [12] proposed a strength assessment method for ice-class propellers based on IACS URI3 and FEM, simplifying the modeling process through automatic mesh generation techniques and analyzing blade stress and deformation under different ice loading conditions, providing a reference for the strength design of ice-class propellers. Kang [13] developed a dynamic numerical simulation method for ice–propeller coupling based on ordinary state-based peridynamics and non-ordinary state-based peridynamics. This method enables the dynamic prediction of ice loads, sea ice fracture, and structural response.

Experimental testing, including model tests and full-scale measurements, is an essential means of studying the structural strength characteristics of propellers. Boswell [14] conducted static tests on highly skewed propellers and found that the stress distribution differed from that of non-skewed propellers, with stress concentrating in a narrow band on the leading edge from the blade root to the 0.9 R section. Zhao Bo [15] carried out static and dynamic stress testing on highly skewed propellers and compared the results with ANSYS calculations. Yang Xianghui [16] investigated blade stress and strain distributions under different advance speed and rotational speed conditions by attaching strain gauges to a propeller model. The dynamic test results revealed stress concentration in the region near the 0.55 R section on the suction side of the blade. In recent years, advanced testing techniques have emerged. Su [17] applied the digital image correlation (DIC) method to measure the surface strain of propeller blades, enabling the visualization of the full-field strain distribution. Ducoin [18] explored the use of fiber Bragg gratings for the structural health monitoring of marine propellers, providing new ideas for full-lifecycle strength assessment.

In terms of optimization design, balancing hydrodynamic performance and strength requirements is a key challenge. Traditional design methods, primarily driven by experience, struggle to achieve multiple design objectives simultaneously. To address this, Bose and Doucet [19,20] proposed a strength design method for ice-class propellers that comprehensively considers the influence of hydrodynamic loads and ice loads on propeller strength. They later demonstrated the application of this strength design method using the MV Ikaluk icebreaker as an example. Kader [21], in a research report from the Taylor Model Basin, described the design process of a U.S. Coast Guard cutter propeller, focusing on meeting hydrodynamic performance goals. The strength design was achieved by satisfying the minimum thickness requirements specified by the ABS classification society, without employing strength analysis methods to calculate blade stress. Liu [22] utilized their in-house-developed PROPELLA program to optimize the propeller design for a twin-screw ice-class vessel, with the design objective primarily focused on hydrodynamic performance. The strength design was realized using chart data. Teunis conducted optimization research on ice-class propellers and qualitatively discussed that the leading edge of ice-class propellers should be thicker to ensure sufficient strength during propeller–ice milling interactions.

This paper systematically elaborates on the design methods for ice-class ship propellers from five aspects: strength design methods, expression of hydrodynamic chart data, preliminary hydrodynamic performance design, preliminary strength performance design, and optimization design considering both hydrodynamics and strength. The article points out that compared to conventional propellers, ice-class propellers need to consider the additional impact of ice loads on strength. It proposes a preliminary hydrodynamic design method based on expressing chart data using neural network models. On this basis, preliminary strength design is carried out using load bending moment estimation and FEM analysis as the primary means. Finally, the paper provides a detailed introduction to the optimization design process for ice-class propellers, considering both hydrodynamics and strength, using a specific design case. From a methodological perspective, this paper systematically summarizes the key technologies for an ice-class propeller design, providing valuable references for guiding and conducting practical engineering design.

This paper systematically elaborates on the design methods for ice-class ship propellers from five aspects: strength design methods, expression of hydrodynamic chart data, preliminary hydrodynamic performance design, preliminary strength performance design, and optimization design considering both hydrodynamics and strength. The article points out that compared to conventional propellers, ice-class propellers need to consider the additional impact of ice loads on strength. It proposes a preliminary hydrodynamic design method based on expressing chart data using neural network models. On this basis, preliminary strength design is carried out using load bending moment estimation and FEM analysis as the primary means. Finally, the paper provides a detailed introduction to the optimization design process for ice-class propellers, considering both hydrodynamics and strength, using a specific design case. From a methodological perspective, this paper systematically summarizes the key technologies for ice-class propeller design, providing valuable references for guiding and conducting practical engineering design.

2. Propeller Strength and Hydrodynamic Design Methods

2.1. Propeller Strength Design Methods

For both conventional “open-water” propellers and ice-class propellers, stress conditions under extreme states are the primary factors considered in strength design. For ice-class propellers, not only do the extreme hydrodynamic loads need to be considered, but also two types of extreme ice loads: the extreme ice loads during normal operation under the ship’s design ice class and the extreme ice loads under off-design conditions.

BOSE N [19] and Carlton [23] considered different extreme states separately and related them using Equation (1):

(1)$1={\gamma }_{ice}{\displaystyle \frac{{\sigma }_{ice}}{{\sigma }_{aice}}}+{\left({\gamma }_{water}{\displaystyle \frac{{\sigma }_{water}}{{\sigma }_a}}\right)}^b+{\gamma }_{mean}{\displaystyle \frac{\left({\sigma }_{mean}+{\sigma }_{res}\right)}{{\sigma }_u}}$

where $\gamma$ is the safety factor, ${\sigma }_{ice}$ is the blade stress caused by ice loads, ${\sigma }_{aice}$ is the stress obtained from the blade material’s fatigue curve based on the number of ice-propeller contact occurrences within the blade’s life cycle, ${\sigma }_{water}$ is the blade stress caused by hydrodynamic loads, ${\sigma }_a$ is the fatigue stress caused by hydrodynamic loads, and ${\sigma }_{mean}$, ${\sigma }_{res}$, and ${\sigma }_u$ are the mean stress, residual stress, and material’s ultimate stress under steady extreme loads, respectively.

Typically, the fatigue term in the above equation can be combined with the number of load cycles using Miner’s rule. Madayag [24] introduced a more conservative stress ratio estimation method:

(2)${\displaystyle \frac{N}{N_r}}={\left({\displaystyle \frac{{\sigma }_r}{\sigma }}\right)}^b$

where $N$ and $N_r$ are the cycle numbers associated with stress $\sigma$ and maximum reference stress ${\sigma }_r$, respectively. For nickel–aluminum bronze, a common propeller material, $b$ is approximately 6.3.

In a preliminary strength design, the cantilever beam theory is often used to calculate blade stress, especially when the blade skew is small, often yielding satisfactory results. The stress concentrates at the thickest section location:

(3)$\sigma ={\displaystyle \frac{M}{Z}}$

where $Z$ is the section modulus and $M$ is the bending moment. For hydrodynamic loads, the in-plane and axial bending moments are small, and $M$ usually refers to the out-of-plane bending moment. During ice–propeller interactions, the axial bending moment is significant and cannot be ignored, forming shear stress at the same section location:

(4)${\tau }_{\mathrm{m}\mathrm{a}\mathrm{x}}={\displaystyle \frac{5.09Q_{sice}}{c{t}^2}}$

where $Q_{sice}$ is the axial bending moment, and c and t are the section chord length and thickness, respectively.

There are two commonly used equivalent stress evaluation methods: the root mean-square calculation method proposed by Osgood [25] and the energy-based calculation method proposed by Timoshenko [26]. This paper adopts the more conservative latter method to calculate the equivalent stress:

(5)${\sigma }_{eq}=\sqrt{{\sigma }_{bend}^2+3{\tau }_{\mathrm{m}\mathrm{a}\mathrm{x}}^2}$

Assuming the maximum steady bending moment, $M_{mean}$ occurs in the bollard pull condition and has a proportionality coefficient $k_w$ with the maximum unsteady bending moment $M_{wake}$, and the residual stress ${\sigma }_{res}$ has a proportionality coefficient $k_{res}$ with the yield stress ${\sigma }_p$, while introducing a correction factor $k_c$ to consider the influence of blade rake on the centrifugal bending moment, combining Equations (1) and (3–5) yields:

(6)$Z={\gamma }_{ice}{\displaystyle \frac{\sqrt{{\left(M_{ice}\right)}^2+3{\left(\frac{5.09Q_{sice}Z}{c{t}^2}\right)}^2}}{{\sigma }_{aice}}}+Z{\left({\gamma }_{wabe}{\displaystyle \frac{k_w{M}_{mean}}{Z{\sigma }_a}}\right)}^b+{\gamma }_{mean}\left({\displaystyle \frac{k_c{M}_{mean}}{{\sigma }_u}}+{\displaystyle \frac{k_{res}Z{\sigma }_p}{{\sigma }_u}}\right)$

The blade section modulus $Z$ can be approximated by $Z=0.15Ec{t}^2$, where the coefficient E is 1.0 or 1.25, depending on whether the section has a “washback” or not. From Equation (6), it is apparent that under the condition of known section chord length, once the bending moments formed by the hydrodynamic and ice loads are determined, the blade section thickness t can be solved. Thus, the problem of designing the blade section thickness is transformed into solving the blade load bending moments.

2.2. Expression of Propeller Hydrodynamic Chart Data

Chart data summarize the experimental rules of propeller performance, helping designers grasp the key factors affecting propeller hydrodynamic performance and assisting in selecting appropriate main parameters, playing a vital role in optimization design. Over the years, multiple series of charts have been developed, such as the B-series, AU-series, Gawn-series, Ka-series, and JD-CPP series. However, the use of chart data involves cumbersome steps such as storage, retrieval, and iterative calculations. With the improvement of computer performance, the chart design method has been extended, and methods using recursive polynomials to express the nonlinear relationship of chart data have emerged. The application of polynomials simplifies the use of charts, but issues such as regression coefficient solving and dimensional curse limit the scalability of chart data. Neural networks, with their strong nonlinear mapping ability and higher approximation accuracy, have overcome the limitations of polynomial methods to a certain extent and can be used to express massive chart data. Numerous studies have shown that BP neural networks are widely used in chart data expression due to their advantages in network structure and mapping ability [27].

BP neural networks are widely used multilayer feedforward neural network models trained using the error backpropagation algorithm. A BP network consists of an input layer, hidden layers, and an output layer, with full connections between layers and no connections within layers. Information is transmitted and processed through activation functions. During the learning process, data are propagated backward layer by layer, and network weights are corrected forward layer by layer from the output layer through the hidden layers in the direction of error reduction. Figure 1 shows the structure of a three-layer BP network with the hyperbolic tangent sigmoid as the activation function, where X, B, and Y represent the neurons in the input layer, hidden layer, and output layer, respectively, and $\omega \cdot v$ represents the connection weights. The input ${\alpha }_h$ of the $h$-th neuron in the hidden layer and the input ${\beta }_j$ of the $j$-th neuron in the output layer can be calculated using Equation (7):

(7)${\alpha }_h=\sum _{i=1}^d{\upsilon }_{ih}{X}_i{\beta }_j=\sum _{h=1}^q{\omega }_{hj}{B}_h$

where $d$ and $q$ represent the number of neurons in the input layer and hidden layer, respectively, $X_i$ is the output value of the i-th neuron in the input layer, and $B_h$ is the output value of the h-th neuron in the hidden layer.

3. Preliminary Design of Propeller Hydrodynamic Performance

The optimization design of marine propellers is a multi-objective problem involving hydrodynamics, cavitation noise, and strength. This section focuses on the hydrodynamic design, with the core being solving the blade load distribution, i.e., the design of pitch and camber distributions. According to the literature and existing conclusions, once the design task is determined, the approximate range of propeller blade pitch can be determined [24]. Among the various design methods and data, the series diagram data formed based on experimental results can assist in evaluating the final scheme. Therefore, the preliminary design research of propeller hydrodynamic performance carried out in this section mainly revolves around the diagrams. Through diagram analysis, a better understanding of the influence of design parameters on propeller performance can be obtained, providing an important reference for optimization design.

3.1. Establishment of a Low-Fidelity Neural Network Proxy Model for Ice-Going Vessel Propeller Optimization Design

Based on the open-water performance data of 14 propeller types in the four-quadrant B-series, this section constructs a neural network model for a propeller optimization design suitable for the complex maneuvering conditions of ice-going vessels. The experimental data of the B-series propellers used in this study originate from the systematic propeller performance tests conducted by Van Lammeren W P A [28] at the Netherlands Ship Model Basin (NSMB) in Wageningen, and the parameters are shown in Table 1. These tests comprehensively measured the open-water performance of B-series propellers under various operating conditions, providing a high-quality original dataset for training the neural network models in this research.

Through systematic collation and analysis, training samples are generated and divided into training set, test set, and validation set in the ratio of 70%, 15%, and 15%, in order to obtain a high-precision and robust propeller performance prediction model, providing strong support for the optimization design of ice-going vessel propellers. In addition, the specific physical meanings of the parameters are as follows. Number of blades (Z): represents the quantity of blades on a propeller, which plays a crucial role in determining its performance characteristics. Expanded area ratio (EAR): defined as the ratio of the total blade area to the propeller disc area, this parameter reflects the size and aspect ratio of the blades relative to the overall propeller dimensions. Pitch angle (β): describes the angle formed between the blade section chord line and a plane perpendicular to the propeller shaft axis, varying along the blade radius and influencing the hydrodynamic properties of each section. Pitch-to-diameter ratio (P/D): expresses the relationship between the axial distance a propeller would theoretically advance during one complete revolution in an ideal fluid and its diameter, providing insight into the propeller’s geometrical pitch and its effect on propulsive efficiency.

When navigating through ice-infested waters, vessels often need to perform reversing maneuvers to extricate themselves from challenging situations or evade obstacles. The focal point of this research’s network model investigation was to predict the performance of propellers under reverse conditions, specifically in the fourth quadrant. We typically employ thrust coefficient (${CT}^{\mathrm{*}}$) and torque coefficient (${CQ}^{\mathrm{*}}$) to characterize propeller performance. During the network model training, the number of blades $Z$, expanded area ratio $EAR$, advance angle $\beta$, and pitch ratio $P/D$ in the diagram were first used as input variables, and the thrust coefficient ${CT}^{\mathrm{*}}$ and torque coefficient ${CQ}^{\mathrm{*}}$ were used as output variables for preliminary training. Details are shown in Table 2.

(8)$$\begin{gathered} C{T}^{*}={\displaystyle \frac{T}{{\rho \ast n^2\ast D}^4}} \\ C{Q}^{*}={\displaystyle \frac{Q}{{\rho \ast n^2\ast D}^5}} \end{gathered}$$

where $T$ is the thrust produced by the propeller, $Q$ is the torque absorbed by the propeller, $\rho$ is the density of the fluid, $n$ is the rotational speed of the propeller, and $D$ is the diameter of the propeller.

The number of neurons in the hidden layer of the network model was taken as 9, 10, 11, 12, 13, 14, and 15, respectively, and the approximation accuracy of the network was measured by the correlation coefficient $R^2$. To reduce the influence of random errors, each network configuration is trained multiple times, and the average value of the correlation coefficients was taken as the training result, as shown by the black fold line in Figure 2. After analyzing the results, it was found that when the number of hidden layers is 11, the approximation accuracy of the network is relatively high, with $R^2$ around 0.98. However, if the network model is to be used for propeller design, the prediction accuracy of the network needs to be further improved.

Considering that the “open-water” performance data of the four quadrants have a certain periodicity, the sine value $sin\beta$ and cosine value $cos\beta$ of the advance angle $\beta$ are used as input variables for network training, and the training results are plotted in Figure 2, where Mid is the number of neurons in the hidden layer of the network. By comparison, it was found that after increasing the input variables, the approximation accuracy of the neural network has been significantly improved. When the number of hidden layers is 14, $R^2$ is about 0.997. The standard “open-water” data and predicted data of propellers with four blades, an expanded area ratio of 0.7, and pitch ratios of 0.5, 0.6, 0.8, 1.0, 1.2, and 1.4 are plotted in Figure 3a. By comparison, it was found that the network prediction results of the thrust coefficient are generally good, but there are large errors in the prediction results of the torque coefficient, which will seriously affect the optimization design of the propeller.

By analyzing the prediction results, it can be found that the locations where the neural network has large prediction errors for the torque coefficient are mainly near the advance angles β equal to 90 degrees and 270 degrees. These two states are the critical points of emergency braking of the propeller, where the torque coefficient fluctuates greatly, leading to large forecast errors. This paper adopts two methods to improve the neural network model: one is to establish separate forecast models for the thrust coefficient and torque coefficient; the other is to divide the data into four parts according to the advance angle for training, with each part roughly centered on the critical state advance angle, and the specific division method is shown in Table 3.

By comprehensively applying the above two improvement methods, the network is re-trained, and the torque coefficient training results of area D are plotted in Figure 4. It can be seen from Figure 4 that the prediction accuracy of the re-trained neural network model is greatly improved, with an overall correlation coefficient $R^2$ greater than 0.998. The standard data and the prediction data of the new network model are re-plotted in Figure 3b, showing that the prediction accuracy of the network model has been greatly improved and can be used for propeller performance prediction and optimization design. In summary, through adjusting and training the network model, the numbers of neurons in the input layer, hidden layer, and output layer of the BP network are finally determined to be 6, 14, and 1, respectively.

3.2. Design Example

After establishing the approximate model for predicting the hydrodynamic performance of propellers, intelligent optimization algorithms (such as particle swarm optimization, genetic algorithms, simulated annealing, etc.) can be utilized to optimize the hydrodynamic performance of propellers. This paper employs a genetic algorithm to optimize propeller parameters, avoiding cumbersome steps such as interpolation and iterative calculations when using charts. The optimization process is shown in Figure 5. More detailed optimization and genetic algorithm execution processes can be found in the literature [29] and will not be repeated here.

In the hydrodynamic design of propellers with the single objective of improving propulsive efficiency, under a limited power constraint, the thrust of the propeller needs to be increased or the torque reduced. The propeller optimization design problem can be described with the following mathematical model:

• (1). Variables: $X=\left[Z,D_{\delta },t,P,EAR,\cdots \right]$, where $Z,D_{\delta },t,P,EAR$ represent the number of propeller blades, relative diameter, thickness, pitch, and expanded area ratio, respectively. The ellipsis indicates other parameters involved in the optimization process, such as camber, skew, or their combination.

• (2). Objective function: $f_2=\left[T,Q,\cdots \right]$, where $T$ and $Q$ represent thrust and torque, respectively. The ellipsis may represent other design objectives such as strength, noise, etc.

• (3). Constraint set: $\mathrm{R}\left[Z_1\le Z\le Z_2; D_1/D_{ref}\le D_{\delta }=D/D_{ref}\le D_2/D_{ref},\cdots \right]$, where the number of blades $Z$ is a positive integer, $D_{ref}$ is the reference diameter and set it to 0.5 m, the diameter $D$ is a real number, and the ellipsis represents other constraints. The constraints need to be set based on the specific design task.

Take the design task of the propulsion unit for a company’s Global 3.0 outboard motor as an example: single propeller, main engine power 2.8 KW, speed 1400 RPM, diameter 0.3 m, forward speed 6 knots, thrust greater than 410 N, reverse speed 0 knots, thrust greater than 450 N, forward power reserve 5%.

Based on the design task, the design objectives and constraints are set as follows:

• Under the premise of meeting the design index requirements, the design objectives are to increase forward thrust and reverse thrust.

• Set the design diameter $D_{\delta }=$ 0.6 and number of blades Z = 3. The main engine power for this task is very small and a three-bladed propeller can meet the design requirements.

• After conversion, considering the signs in different quadrants, set the forward design advance angle to 11.33°, thrust coefficient $CT{1}^{\mathrm{*}}\ge 0.04706$, torque coefficient $CQ{1}^{\mathrm{*}}\le 0.00694$; reverse design advance angle to 180°, thrust coefficient $CT{2}^{\mathrm{*}}\le -0.05373$, torque coefficient $CQ{2}^{\mathrm{*}}\ge -0.00752$.

• Pitch ratio $0.4\le P/D\le 1.4$, expanded area ratio $0.3\le EAR\le 0.85$. The design range can be set based on experience.

The optimization design of propellers is a multi-objective optimization problem. The solution set of the optimization design can be visualized through a Pareto chart, and the designer can analyze it based on the specific problem to determine the final optimization design scheme.

Based on the neural network formed in the previous section and using a genetic algorithm for optimization design, the Pareto chart for this design example is obtained. In Figure 6, the horizontal coordinate is the negative value of the forward thrust coefficient, and the vertical coordinate is the reverse thrust coefficient (the fitness function in the genetic algorithm is to minimize the objective function). It can be seen from Figure 6 that a relatively uniform solution set is obtained using the text method, covering a variety of optimization design schemes with forward thrust coefficients ranging from 0 to 0.2 and reverse thrust coefficients ranging from 0 to −0.25.

Based on constraint c, the final optimization scheme is selected from the Pareto chart. Two schemes are screened out and marked in Figure 6. The average value of the two schemes is taken as the final design scheme, determining the final pitch ratio as 0.61 and the expanded area ratio as 0.55.

Before delivering the propeller design, the hydrodynamic performance was further verified using the CFD software STAR-CCM+ (version 2020.2.1 Build 15.04.010), and the calculation model is shown in the Figure 7. Extensive studies have shown that the CFD results are in good agreement with experimental measurements, demonstrating the reliability and maturity of this method. In this study, unstructured polyhedral meshes were employed, with local refinements applied to the blade surfaces and wake regions. The non-dimensional wall distance (y+) of the first layer of cells was appropriately set.

A mesh convergence study was conducted with the total number of grids ranging from 3.32 million to 16.77 million. When the grid count increased to 7.38 million, the changes in the thrust coefficient (CT) and torque coefficient (CQ) were less than 2%, indicating grid independence. As the grid size further increased to 16.77 million, the relative error decreased to 1.76%. The specific values of CT and CQ for each grid size are presented in Table 4. Ultimately, a medium-density mesh with a total grid count of 9.40 million was selected for the propeller model simulations, striking a balance between computational accuracy and efficiency.

The CFD simulations revealed that at a forward speed of 6 knots, the designed propeller generated a thrust of 420 N and a torque of 15.5 N m, with a power consumption of 2272 W. In the astern condition at 0 knots, the propeller produced a reverse thrust of −475 N and a torque of −18.9 N·m, consuming 2770 W of power. These results confirm that the propeller design meets the specified requirements, and the propeller is shown in Figure 8.

4. Preliminary Design of Propeller Strength Performance

The ice loads on propeller blades lead to in-plane, out-of-plane, and spindle bending moments. Under the combined action of various bending moments, an accurate stress distribution and concentration on the blades can only be obtained through detailed stress analysis methods. However, before accurately evaluating the blade strength using the FEM method, a preliminary design of the blade strength is required to obtain an initial blade thickness distribution. This section introduces the preliminary design method for blade strength using a design example.

4.1. Load Moment Estimation Methods

• (1). Hydrodynamic Load Moments

There are various methods to calculate the hydrodynamic loads on propellers. The maximum thrust of the blades in bollard pull condition can be accurately calculated using the panel method, estimated using neural network models, or directly estimated using the fitting formula by VAN [30]:

(9)$\left\{\begin{array}{c}{K}_T=0.485{\left({\displaystyle \frac{P}{D}}\right)}^2-0.2532{\displaystyle \frac{P}{D}}+0.1064\\ K_Q=0.0709{\left({\displaystyle \frac{P}{D}}\right)}^2-0.042{\displaystyle \frac{P}{D}}+0.0166\end{array}\right.$

Subsequently, the thrust $F_T$ and torque $F_Q$ of a single blade can be calculated. Once the load amplitude is determined, the maximum steady bending moment $M_{mean}$ can be calculated using the following equation:

(10)$M=F_T\left(a-r_0\right)cos\theta -{\displaystyle \frac{F_Q}{b}}\left(b-r_0\right)sin\theta$

where $\theta$ is the geometric pitch angle, and $a$ and $b$ are the moment arms from the propeller shaft to the load center, which are assumed to be equal in this analysis.

• (2). Ice Load Moments

For the strength design of propellers in ice-infested waters, two types of extreme ice loads need to be considered: design operating conditions and non-design conditions. For the former, refer to the IACS rules. However, the rules do not provide a calculation method for the out-of-plane bending moment due to ice loads. Therefore, an estimation formula needs to be established. Assuming the ice load center is at the 0.95 R section, the out-of-plane bending moment can be calculated using the following equation:

(11)$\left\{\begin{array}{c}{M}_{ice}={\displaystyle \frac{D}{2}}\left(0.95-r_h\right)F_{ice} \mathrm{f}\mathrm{o}\mathrm{r} \mathrm{r}\mathrm{o}\mathrm{o}\mathrm{t} \mathrm{s}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\\ M_{ice}={\displaystyle \frac{D}{2}}\left(0.95-0.6\right)F_{ice} \mathrm{f}\mathrm{o}\mathrm{r} 0.6 \mathrm{R} \mathrm{s}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\end{array}\right.$

where $F_{ice}$=$\mathrm{m}\mathrm{a}\mathrm{x}\left(F_b,F_f\right)$. Although the rules provide different calculation formulas for the maximum backward bending load $F_b$ and the maximum forward bending load $F_f$ on the blades, this section assumes that the maximum ice load can occur in any direction during the entire lifecycle of the blades. Regarding the spindle bending moment $Q_{sice}$, the IACS rules provide a detailed calculation method:

(12)$Q_{sice}=0.25F_{ice}{c}_{0.7}$

where $c_{0.7}$ is the chord length at the 0.7 R section.

For extreme loads under non-design conditions, the ice load amplitude is closely related to the number of load cycles (number of ice–propeller contacts). This paper adopts the ice load cycle calculation method proposed in the IACS rules:

(13)$N_{ice}=k_1{k}_2{k}_3{k}_4{N}_{class}n$

where $N_{class}$ is the reference number of impacts per propeller revolution for different ice classes (Table 5); $k_1$ is the propeller location factor, which is 1 for propellers on the centerline and 1.35 for propellers on the wings; $k_2$ is the propeller type factor, which is 1 for open propellers and 1.1 for ducted propellers; $k_3$ is the propulsor type factor, which is 1 for fixed propulsors an 1.2 for azimuthing propulsors; $k_4$ is the immersion factor, which can be determined by the following equation, where $h_0$ is the propeller shaft immersion, D is the diameter, and $H_{ice}$ is the ice thickness parameter:

(14)$k_4=\left\{\begin{array}{ll}0.8-f& f<0\\ 0.8-0.4f& 0\le f\le 1\\ 0.8-0.2f& 1<f\le 2.5\\ 0.1& 2.5<f\end{array}\right., f={\displaystyle \frac{h_0-H_{ice}}{D/2}}-1$

Afterward, the extreme ice loads under non-design conditions can be determined using fitting formulas, and the ice load moments can be determined using Equation (11).

4.2. Design Example

The preliminary design of blade strength does not require complex numerical processes. After preparing the necessary input parameters, Equation (6) can be programmatically solved. The design process is shown in Figure 9. This subsection illustrates the preliminary design process of blade strength using the R-Class propeller discussed in the previous chapters as a design example.

The R-class propeller is a fixed-pitch propeller mounted on the wings without a duct. The propeller diameter D is 4.12 m, the rotational speed n is 180 RPM, the number of blades Z is 4, and the material is nickel–aluminum bronze. The expanded area ratio EAR is 0.708 and the rake is 262.8 mm. The pitches at blade roots 0.6 R, 0.7 R, 0.8 R, and blade tip are 2.76 m, 3.2 m, 3.27 m, 3.3 m, and 3.24 m, respectively. The corresponding chord lengths are 1408 mm, 1783 mm, 1720 mm, and 1589 mm. The length between perpendiculars LBP of the Des Groseilliers ship is 78.8 m, the freeboard height T is 7.17 m, the blade root thickness is 280 mm, and the thickness of the 0.6 R section is 180 mm.

The specific design steps are as follows:

• (1). Determine the class parameters

  • (a). Determine the sea ice thickness parameter according to the rules. For the PC3 class R-class propeller, $H_{ice}=3.0 \mathrm{m}$.

  • (b). Calculate the propeller shaft immersion (assuming the blade bottom is above the ship bottom), $h_0=T-D/2=5.11 \mathrm{m}$.

  • (c). Determine the blade location and propeller type parameters, $k_1=1.35,k_2=1$.

  • (d). Determine the propulsor type parameter and determine the immersion factor according to Equation (12), $k_3=1,k_4=0.79$.

  • (e). Determine the number of load cycles during the blade lifecycle using Equation (11), $N_{total}=k_1{k}_2{k}_3{k}_4{N}_{class}n/Z=12\times {10}^6$.

• (2). Determine the hydrodynamic loads

  • (a). Calculate the thrust and torque coefficients of the blades in bollard pull condition. $k_T=0.4,k_Q=0.045$ (panel method); $k_T=0.39,k_Q=0.047$ (estimation formula). It can be seen that the difference between the two is small. For convenience, the estimation formula is used in the subsequent calculations.

  • (b). Calculate the blade thrust and torque using Equation (8) (seawater density is taken as 1025 kg/m³). $F_T=\mathrm{257,556}N,F_T=\mathrm{128,963}N$.

  • (c). Determine appropriate safety factors,${\gamma }_{wake}=1.5,{\gamma }_{mean}=1.5,{\gamma }_w=0.15,{\gamma }_c=1.05,{\gamma }_{res}=0.25$.

  • (d). Determine the cantilever beam length.

    Blade root section: $r_0=0.3\times R=0.618 \mathrm{m},a=b=0.75\times R=1.545 \mathrm{m}$. 0.6 R section: $r_0=0.6\times R=1.236 \mathrm{m},a=b=0.8\times R=1.648 \mathrm{m}$.

  • (e). Calculate the geometric pitch angle.

    Blade root section: $\theta ={\mathrm{t}\mathrm{a}\mathrm{n}}^{-1}\left(P/2\pi r\right)=0.618 \mathrm{r}\mathrm{a}\mathrm{d}$. 0.6 R section: $\theta ={\mathit{tan}}^{-1}(P/2\pi r)=0.391 \mathrm{rad}$. $\theta ={\mathrm{t}\mathrm{a}\mathrm{n}}^{-1}\left(P/2\pi r\right)=0.391 \mathrm{r}\mathrm{a}\mathrm{d}$.

  • (f). Calculate the hydrodynamic load moments using Equation (9).

    Blade root section: $M_{mean}=\mathrm{239,433} \mathrm{N}\mathrm{m}$. 0.6 R section: $M_{mean}=\mathrm{82,791} \mathrm{N}\mathrm{m}$.

• (3). Determine the material properties

  • (a). Determine the exponent b used for ice load calculations according to Equation (2) and the material fatigue curve. For the nickel–aluminum bronze material in this example, $b=6.3$.

  • (b). Look up and obtain the alternating stress limit of the material. For the nickel–aluminum bronze material in this example, ${\sigma }_a=130 \mathrm{M}\mathrm{P}\mathrm{a}$.

• (4). Calculate the blade thickness without ice loads

  • (a). Set b = 1, ignore the ice load term (first term on the right side of the equal sign) in Equation (6), and solve for the section modulus. Blade root section: $Z=\mathrm{3,937,343} {\mathrm{m}\mathrm{m}}^3$. 0.6 R section: $Z=\mathrm{1,361,464} {\mathrm{m}\mathrm{m}}^3$.

  • (b). Calculate the section thickness from the section modulus according to the section shape.

    Blade root section: t = 176 mm. 0.6 R section: t = 92 mm.

• (5). Calculate the blade thickness under design ice load conditions

  • (a). The ice load under design conditions, $F_{ice}=1,678,357 \mathrm{N}$.

  • (b). Calculate the out-of-plane bending moment $M_{ice}$ and spindle bending moment $Q_{sice}$ using Equation (11): Blade root section: $M_{ice}=1,063,394 \mathrm{Nm}.$ 0.6 R section: $M_{ice}=425,357 \mathrm{Nm}$. $Q_{sice}=721,540 \mathrm{Nm}$.

  • (c). Set $b=6.3$, solve Equation (6) to calculate the section modulus: Blade root section: $Z=\mathrm{9,750,976} \mathrm{m}{\mathrm{m}}^3$. 0.6 R section: $Z=\mathrm{6,878,462} \mathrm{m}{\mathrm{m}}^3$.

  • (d). Calculate the section thickness from the section modulus according to the section shape: Blade root section: $t=213\mathrm{mm}$. 0.6 R section: $t=142 \mathrm{mm}$.

• (6). Calculate the blade thickness under non-design ice load conditions

  • (a). Calculate the non-design condition ice load (approximately 50.7% higher than the design condition) based on the number of collisions and the extreme load curve, $F_{ice}=2,529,283 \mathrm{N}$.

  • (b). Calculate the out-of-plane bending moment $M_{ice}$ and spindle bending moment $Q_{sice}$ using Equation (11): Blade root section: $M_{ice}=1,602,535 \mathrm{Nm}$. 0.6 R section: $M_{ice}=641,014 \mathrm{Nm}$. $Q_{sice}=1,087,362 \mathrm{Nm}$.

  • (c). Set $b=6.3$, solve Equation (6) to calculate the section modulus: Blade root section: $Z=\mathrm{14,352,463} \mathrm{m}{\mathrm{m}}^3$. 0.6 R section: $Z=\mathrm{10,243,095} \mathrm{m}{\mathrm{m}}^3$.

  • (d). Calculate the section thickness from the section modulus according to the section shape: Blade root section: $t=260 \mathrm{mm}$. 0.6 R section: $t=175\mathrm{mm}$.

Table 6 presents the preliminary design results of blade strength. It can be seen that ice loads have a significant impact on blade strength. Without considering the sea ice loading conditions, the blade thickness required to withstand hydrodynamic loads is only 63% of the actual blade thickness. For the two extreme ice loads, the blade root thickness obtained from the design condition extreme load is about 76% of the actual blade root thickness, while the blade root thickness calculated from the non-design condition is about 93% of the actual blade root thickness. The designed thickness of the 0.6 R section differs from the actual section thickness by approximately 2.7%.

Through this design example, it can be found that in the preliminary design process of blade strength, the non-design condition strength has the highest requirement for section thickness. The extreme load estimation formula can accurately predict the extreme loads that blades may encounter under non-design conditions during their lifecycle, and thereby conduct a preliminary design of section strength.

5. Optimal Design of Ice-Class Propellers Considering Hydrodynamic and Strength Properties

The previous two subsections introduced the preliminary design methods for propeller hydrodynamics and strength. Based on these research foundations, this section uses the design of the propeller for the Biscayne (WTGB-104), a Bay-class U.S. Coast Guard icebreaking tugboat, as an example to introduce the optimal design method for ice-class propellers considering hydrodynamics and strength. It should be particularly noted that in specific applications, the constraints and design objectives of different design tasks often vary, and designers need to flexibly adjust the design process according to different requirements to achieve the design purpose.

5.1. Design Task

Multi-condition operation is a typical characteristic of ice-going vessels. Ice-class vessels have at least two typical operating conditions: ice and ice-free. Therefore, this section uses tugboats with large variations in operating conditions as an example to introduce the propeller design process. The operating conditions of icebreaking tugboats involve “bollard pull” conditions, low-speed towing conditions, and free-running conditions. Any of these conditions can be used as the design point to absorb full power. However, considering that fixed-pitch propellers will only absorb full power at one speed and one rotational speed, compromises need to be made in other operating conditions. The design of conventional propellers often selects the free-running condition as the design point, but considering that tugboats operate mostly in low-speed towing conditions, the design process needs to be appropriately adjusted to use the towing condition as the design point (the optimization objective function is mainly based on hydrodynamics, while strength optimization aims to ensure structural safety).

Table 7 provides the parameter information for the corresponding design task. Additionally, the U.S. Coast Guard specified that the ahead vessel’s towing speed is 5 knots, and the main engine can achieve a full-power output at any speed between 245 and 305 RPM. The Coast Guard also requires that the astern bollard pull is not less than 75% of the ahead bollard pull.

The optimization design process for ice-class propellers is shown in Figure 10. After preliminarily determining the blade geometry using the design methods introduced in the previous two subsections, the load distribution of the blades can be designed and optimized. The hydrodynamic and strength performance are then checked in turn to form the final scheme. It should be noted that the design in this paper uses the IK series airfoils, which have excellent performance and are widely used in low ice-class propellers.

5.2. Preliminary Design Process

In the initial design stage of propeller hydrodynamics, the influence of sea ice was not considered, and the preliminary design of propeller hydrodynamic performance was directly carried out through the established neural network surrogate model. There are two main advantages of using neural networks for preliminary design: First, neural networks have strong generalization ability, and the hydrodynamic performance prediction of propellers with arbitrary expanded area ratios and pitch ratios within a given range can be realized without re-performing Fourier analysis and compiling design programs. Second, the data source of neural networks is experimental data, which can verify the design results of the circulation theory (under a given design task, the design results of blade pitch often differ little, so the rationality of the theoretical design results can be judged by comparing the pitch ratio given by the surrogate model with the theoretically designed pitch ratio).

The primary design goal is to obtain maximum thrust at the design point, that is, to achieve maximum thrust at a speed of 5 knots and maximum main engine power. As discussed above, to calculate the hydrodynamic performance of the propeller, the advance speed (advance coefficient or advance angle) must first be determined. Currently, only the speed is determined, so the propeller speed needs to be assumed.

Combining the main engine speed data given in Table 4 and considering the fact that the main engine speed in actual towing is higher than the bollard pull condition speed, it can be assumed that the main engine speed range is 250~310 RPM. Combining the wake fraction given in the literature, the corresponding advance coefficient (converted to advance angle when using the neural network model) and torque coefficient can be interpolated, as shown in Table 8.

Using the established neural network model, the pitch ratio and thrust of the propeller absorbing full power at different expanded area ratios can be calculated. As shown in Figure 10, as the expanded area ratio increases, the pitch of the propeller absorbing full power at the same speed shows a decreasing trend. In addition, Figure 11 also shows the variation curve of blade pitch ratio with expanded area ratio in the bollard condition. The final design scheme needs to be selected on this curve. It can be seen that as the expanded area ratio increases, the blade pitch ratio gradually approaches 0.7.

Then, the preliminary design scheme is determined by the ratio of astern thrust to the ahead thrust. In fact, only the four-quadrant hydrodynamic data of propellers with specific pitch ratios and expanded area ratios are available. To evaluate the hydrodynamic performance in astern conditions, separate design programs usually need to be developed using methods such as Fourier series to extend the applicability of the charts. Strom et al. conducted related research and developed a design program. While the neural network model established in this paper shows potential for predicting the astern performance of a given propeller without the need for additional Fourier series methods to develop design programs, it is important to note that the model’s generalization ability may be limited due to the relatively small and dated propeller geometry dataset used for training. The pitch distribution, for example, is sometimes given as an average or assumed to be constant, whereas in reality, propeller geometries often have non-constant pitch distributions along the blade span. Further validation and testing with a more diverse and modern propeller geometry dataset would be necessary to conclusively demonstrate the model’s strong generalization capability.

Figure 11 shows the variation in the ahead and astern bollard pull thrusts of the propeller with the expanded area ratio at a pitch ratio of 1.0. In Figure 12, it can be seen that as the expanded area ratio increases, both the astern and ahead thrust show an increasing trend, but the thrust ratio first increases and then decreases, reaching a maximum at an expanded area ratio of 0.7. Therefore, the expanded area ratio of the preliminary propeller design scheme is determined to be 0.7, and the thrust ratio is about 0.76, meeting the design requirements. After determining the propeller expanded area ratio of 0.7, the pitch ratio can be determined as 0.73 (red circle in the figure, and represents the intersection of both the reverse thrust requirements and the forward hydrodynamic performance requirements) from the purple curve in Figure 11.

Through the above discussion, the preliminary design scheme for propeller hydrodynamics can be determined, that is, a B-series propeller with an expanded area ratio of 0.7 and a pitch ratio of 0.73. This provides input parameters for the preliminary design of blade thickness, as shown in Table 9.

According to the strength design method, the preliminary design of blade strength is carried out, as shown in Table 10. The table also gives the initial thickness of the B-series propeller and the design thickness of Kader [21] for this case as a comparison:

From Table 9, it can be seen that the preliminary design results of the section thickness in this paper and Kader’s design results are much higher than the initial thickness of the B-series propeller. This also shows that the blade strength enhancement brought about by ice loads is the most typical feature that distinguishes ice-class propellers from conventional propellers. Compared with Kader’s design results, the blade root section designed in this paper is thicker (about 16% in the design ice condition and about 40% in the non-design ice condition). This may be caused by two reasons. First, Kader used earlier ABS regulations for blade strength verification. As mentioned in the introduction of this paper, the early regulations only made simple provisions on blade thickness and did not specify the design ice load and the main shaft bending moment caused by ice load, nor consider non-design conditions, resulting in smaller section-thickness design results than those in this paper. Second, Kader may have modified the blade strength based on the actual operating environment and period of the tugboat. By considering the ice conditions in the operating area of the Biscayne vessel, the section thickness may have been appropriately reduced.

The final design of blade strength needs to be evaluated in combination with the finite element model. Here, the section thickness under the design ice condition is used as the preliminary design result. At this point, the preliminary design of the blade considering hydrodynamics and strength is completed. After the preliminary design of hydrodynamics and strength, an initial geometric model based on the B-series propeller diagram is formed, as shown in Figure 13a. Its main parameters are as follows: pitch ratio of 0.73, expanded area ratio of 0.7, blade root thickness of 191 mm, 0.6 R section thickness of 129 mm, and blade tip thickness of 40 mm. This model can be used for the next step of the intermediate design process. Based on the propeller blade structure, standard ice loads were applied in accordance with relevant regulations. The analysis results indicate that the maximum stress experienced by the propeller blades is approximately 50 MPa, which is within the acceptable limits specified by the structural strength design requirements.

5.3. Intermediate Design Process

The intermediate design process for propellers utilizes the lifting surface method proposed by Kerwin to determine the camber and pitch distributions of the blade sections. In this study, the lifting surface is set on the initial cambered surface. For detailed theory and numerical calculation procedures related to the lifting surface, such as boundary condition handling and discretization of the blade and wake regions, please refer to the research by Zhou [31]. The same calculation program is used here. In previous lifting surface designs, the results of the lifting line design were usually used as the initial approximation for iterative design. In this study, based on the propeller formed in the preliminary design chart, its parameters are adjusted to form the initial approximate propeller in the lifting surface program.

According to the characteristics of the airfoil profiles and pitch distributions used in ice-class propellers, this section replaces the airfoil profiles of the initial propeller model with the IK82 profile and modifies the pitch distribution form of the propeller. Moreover, it was found that the circulation distribution of ice-class propellers under design conditions is closer to the optimal circulation distribution form (i.e., in the equation, a=0.1, b=0.9). Therefore, the optimal circulation distribution form is also adopted in the design process.

Before conducting the lifting surface design, it is necessary to initialize the rake and skew distributions of the blades to determine the shape of the reference surface. For rake, as mentioned earlier, it often needs to be determined based on the consideration of the tip clearance and the stern profile. This paper does not discuss too much about the rake distribution and directly adopts Kader’s [21] design result, i.e., the rake angle is 8 degrees. Regarding the blade skew, its design purpose is to avoid different radial sections simultaneously rotating into high wake regions, causing large fluctuations in unsteady thrust. Therefore, the design principle of blade skew is to ensure that the skew distribution line has a large intersection angle with the section normal wake phase distribution curve based on the harmonic analysis of the wake field.

Ignoring the influence of radial wake, the normal wake $w_N$ of the blade section can be calculated from the axial wake $w_X$ and tangential wake $w_T$:

(15)$w_N=w_X\mathit{cos}{\beta }_p-w_T\mathit{sin}{\beta }_p$

Considering the symmetry of the wake field of single-screw ships and rewriting Equation (14) into the form of superposition of various harmonic components

(16)$w_N={\overline{w}}_{0N}+{\sum }_{n=1}^{\mathrm{\infty }}{C}_{nN}\mathit{cos}(n\theta +{\psi }_{nN})$

When $n\theta +{\psi }_{nN}=0$, the value of the normal wake $w_N$ is maximum, thus determining the direction angle position of the nth order maximum value of the normal wake. However, in the actual design process, due to the need to ensure the smoothness of the blade geometry, it is difficult to have large intersection angles at different radial positions. Therefore, the blade skew distribution is determined in conjunction with the intersection angle between the skew distribution line and the section normal wake phase distribution curve. In addition, before designing the skew, it is necessary to initialize the pitch distribution of the blades, keeping the average pitch ratio unchanged at 0.73.

For the convenience of display, combined with the data provided by Kader, the wake field distribution results of the tugboat Biscayne are redrawn in Figure 14. It can be seen from Figure 13 that the axial and tangential wake distributions behind the ship have symmetry, so the derivation of Equation (15) is reasonable. For the axial wake, the high-wake region is within about ±15° behind the ship, while for the tangential wake, the high-wake region ranges from about 15° to 135°, and 225° to 345°. The Fourier transform method is used to perform a harmonic analysis on the wake field to obtain the amplitude and phase angle of the wake components.

Figure 15 shows the results of the harmonic analysis of the tugboat Biscayne’s wake field. In fact, the results of the harmonic analysis of wake components are mainly used to assist in the selection of the number of blades. Although the number of blades is given in the design example of this section, it can be seen that the selection of a four-blade propeller is consistent with the conclusion about the relationship between the wake field and the number of blades in Zhou [31]; that is, the amplitudes of the KZ+1- and KZ-1-order wake harmonics are small, and the phase of the KZ order harmonic does not change much along the radial direction.

Figure 16 compares the phase angles of the wake field. Zhou [31] and Ye [32] proposed from the perspective of reducing blade frequency vibration that the blade skew distribution should have a large intersection angle with the Z-order wake phase angle. From this perspective, the skew distribution of propeller no. 8 is more favorable. It can be seen from Figure 14 that the fourth-order wake phase angle has the same bias as the skew angles of propellers No. 2 and No. 5, reducing the intersection angle between them, while propeller no. 8 will increase the intersection angle between the two (the fourth-order wake phase angle and the skew distribution of propeller no. 8 are bolded in the figure, and the skew of propellers no. 2 and no. 5 falls within the interval). At the same time, the design of the blade skew should also be combined with the ship type, operating tasks, and requirements for further analysis to determine the final skew distribution. Therefore, this section refers to the conclusion in Ye [32] and selects the skew distribution of propeller no. 8 to carry out the intermediate design, and compares it with the non-skewed propeller in the final design to determine the ultimate skew distribution scheme.

At this point, all the initial blade geometry parameters for the lifting surface design have been determined. The lifting surface program described in the literature can be used to design the blade cambered surface and determine the pitch and camber distributions of the blade sections. Figure 17 shows the results of the lifting surface cambered surface design. The dashed lines in the figure are the original design results. Since the Kerwin method sets the circulation at the hub to be non-zero, the blade pitch and circulation have inflection points at the blade root. Here, the camber and pitch distributions are smoothed, as shown by the solid lines in the figure.

From Figure 17a, it can be seen that the camber of the blade sections is mainly distributed between 0.015 and 0.03, which is consistent with the statistical results of the camber distribution of ice-class propeller sections, indicating that the section camber distribution after the lifting surface design is relatively reasonable. Similarly, analyzing the pitch distribution shown in Figure 17b, it can be found that the pitch distribution of the lifting surface design is similar to that of propellers no. 6, no. 10, and KP505. The pitch shows an increasing and then decreasing trend from the blade root to the blade tip, reaching a maximum value around 0.6 r. Although there is a difference from the pitch distribution of the initialized propeller no. 2, it also indicates that the pitch design results are relatively reasonable. The main reason for this difference is that the influence of sea ice action on the blade geometry parameters was not considered in the hydrodynamic design. In the subsequent final design process, the conclusions related to the projected area of the tip sections can be applied to the selection of blade parameters to consider the influence of sea ice action.

5.4. Finalization of the Design Process

The performance of the propellers designed based on the circulation theory is often not optimal. On one hand, the preliminary and intermediate design processes described above take the towing condition as the design point without considering the influence of the free-running condition in non-towing states. On the other hand, the idealized vortex system model in the lifting surface theory introduces errors in the calculation of propeller hydrodynamic performance. To evaluate the propeller performance more accurately, tools with higher precision are needed. This section utilizes the panel method calculation program to conduct the final design of the propeller.

Apart from efficient propeller performance evaluation, the key issues to be addressed in the final (optimization) design of propellers include the following: parametric representation of the propeller, selection of optimization variables and algorithms, and setting of optimization objectives. Due to space limitations and research focus, this section provides a brief description of these issues.

B-spline curves offer flexible control over curves and can obtain smooth radial distributions of blade geometric parameters with fewer control points. The selected control points can be used as design variables in the optimization process. Therefore, this paper employs B-spline curves for the parametric representation of propeller geometry. Referring to the optimization case in ATKINSON [29], this section preliminarily selects the camber values, chord lengths, pitch values, and skew values at 0.25 R, 0.5 R, 0.85 R, and 1.0 R as the design variables (rake values are not optimized; thickness values are optimized after hydrodynamic optimization), totaling 16 design variables. To enhance the approximation modeling and optimization design effects, the selection of design variables and their variation spaces is based on the following six points:

• The geometric parameters of high-performance propellers vary smoothly along the radial direction. Abrupt changes of geometric parameters at certain radial positions are not allowed.

• The chord length and camber at the blade tip are minimal and do not appear as design variables in the optimization process.

• The maximum chord length of the blade occurs at the outer radius, between 0.6 R and 0.7 R.

• According to the statistical and research results in Chapter 5, the maximum pitch of the blade will not occur after 0.6 R. At the same time, the rapid decrease in pitch in the outer radius is beneficial for reducing ice loads.

• The camber values of blade sections are all positive, and negative camber is not allowed.

• The skew distribution type is consistent with the intermediate design results, adopting a balanced skew distribution. The midpoint trajectory of the chord length is located on the right and left sides of the blade reference line in the inner and outer radius regions of the blade, respectively.

To improve the optimization design effect, the design space should be as large as possible during optimization. Considering the optimization efficiency and the above six parameter settings, the optimization design space for the design variables is selected, as shown in Table 11.

After determining the optimization variables and their design space, the optimization algorithm needs to be selected based on the optimization objectives. The final design of propellers is a multi-objective optimization problem. This section sets two optimization objectives: improving efficiency under towing conditions and improving efficiency under free-running conditions. The constraint conditions are the power and speed limitation conditions described in the preliminary design. To solve this multi-objective optimization problem, the optimization algorithm used is consistent with that in the hydrodynamic preliminary design process. For detailed calculation procedures such as genetic algorithm fitness calculation, selection, crossover, and mutation parameter settings, please refer to ATKINSON [29].

It is worth mentioning how to determine the weight factors of the two efficiencies. Assuming that the sailing distances of the Biscayne vessel under towing and free-running conditions are 1:1 in its life cycle, and the main engine operates at a maximum continuous power during the operation process, and the ratio of operation time under the two conditions is inversely proportional to the speed, $T_{tow}:T_{free}=V_{free}:V_{tow}$, with a value of approximately 3. Assuming that the power received by the propeller under both conditions is PD and the open water efficiencies are ${\eta }_{tow}$ and ${\eta }_{free}$, respectively, the effective work W_E of the propulsor can be expressed as

(17)$W_E=P_D{T}_{tow}{\eta }_{tow}+P_D{T}_{free}{\eta }_{free}$

Substituting the speed ratio into Equation (16) yields

(18)$W_E=P_D{T}_{free}(3{\eta }_{tow}+{\eta }_{free})$

The weight factors of the two efficiencies can be determined from Equation (18).

Figure 18 presents the results of 500 optimization designs, with the green curve representing the trend line. It can be seen that the propulsive efficiency of the optimized propeller under towing conditions is mainly between 0.265 and 0.285, while the efficiency under free-running conditions is roughly between 0.545 and 0.575. The schemes that meet the design requirements for towing conditions are within the red line range shown in Figure 17b. Based on the parameter distribution analysis results and to ensure the smoothness of the blade geometry, the average value near the trend line is selected as the optimization scheme.

After completing the hydrodynamic optimization, strength optimization methods are used to optimize the blade strength. After strength verification, the final optimized propeller is formed, and its blade geometric parameters are shown in Table 12.

Through preliminary design, intermediate design, and final design, the ultimate optimized scheme is formed. Comparing the geometry of the designed propeller with that of the original propeller of the Biscayne vessel, it was found that the distributions of chord length, pitch, thickness, and camber are similar, but there are certain differences in the skew distribution. The original propeller adopted a zero-skew distribution form, while the propeller designed in this paper incorporates skew based on the wake field analysis. This difference may be determined by the vessel’s operating tasks. For ice-class towing vessels with thick hulls, the resonance caused by an unsteady propeller thrust is not particularly considered. However, from the perspective of reducing resonance, this paper designs blade skew. After strength verification and hydrodynamic evaluation, it can achieve good propulsion effects. Regarding the open water performance of the propeller, the design performance has been improved by 0.6%.

It is important to emphasize that the variation of propeller three-dimensional geometric parameters in practice is quite complex and requires comprehensive consideration of the actual operating environment and vessel system. It is necessary to deeply explore the relationship between geometric parameters and propeller performance, and to achieve optimized propeller design based on practical constraints.

6. Conclusions

This paper proposes a systematic optimization design process and methodology for ice-class propellers, which can provide beneficial reference and guidance for the design of propellers for ice-going vessels. In practical applications, designers need to flexibly adjust the design process according to specific mission requirements to achieve design goals that balance efficiency and reliability.

(1) The design of ice-class propellers needs to consider factors for both ice loads and hydrodynamic loads. Ice loads mainly include ice loads under normal operating conditions and extreme ice loads under off-design conditions, while the hydrodynamic loads mainly come from extreme conditions. The design needs to comprehensively consider the impact of these loads on blade strength.

(2) Neural networks are a suitable method for expressing propeller hydrodynamic chart data. Approximate models of propeller hydrodynamic performance can be established based on BP neural networks to assist in the preliminary design and optimization design of propellers.

(3) In the preliminary design stage of propellers, the influence of design parameters on performance can be evaluated based on neural network proxy models. At the same time, it is necessary to determine appropriate design points according to the actual operating conditions of ice-class vessels, such as towing and ice-breaking.

(4) In propeller strength design, it is necessary to estimate the bending moments caused by hydrodynamic loads and ice loads, and to preliminarily design the blade thickness based on the fatigue characteristics of the material, providing initial values for subsequent finite element analysis.

(5) The final design of propellers needs to be carried out based on multi-objective optimization considering both hydrodynamics and strength. The optimization process involves key issues such as the selection of design variables, the application of optimization algorithms, and the balance of multi-objective weights.

Author Contributions

Conceptualization, C.W. (Chao Wang) and C.W. (Chunhui Wang); Methodology, C.G. and K.H.; Validation, C.G.; Resources, K.W.; Writing—original draft, C.G.; Writing—review & editing, C.G.; Funding acquisition, C.W. (Chao Wang) and C.W. (Chunhui Wang) All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data is contained within the article.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Three-layer BP neural network.

Figure 2. Network model accuracy comparison.

Figure 3. Propeller hydrodynamic performance results predicted by the neural network. (a) Forecast results before network optimization. (b) Forecast results after network optimization.

Figure 4. Verification of the network forecast accuracy.

Figure 5. Process of preliminary design of hydrodynamic performance.

Figure 6. Pareto chart for hydrodynamic design.

Figure 7. CFD calculation model. (a) Model blade. (b) Watershed grid division.

Figure 8. Propeller physical model.

Figure 9. Propeller strength preliminary design process.

Figure 10. Optimal design process for ice-class propellers (considering hydrodynamic and strength properties).

Figure 11. Predicted propeller parameters based on the neural network model.

Figure 12. Predicted astern and ahead thrusts and their ratio.

Figure 13. (a) Propeller preliminary design model. (b) Blade stress diagram.

Figure 14. Wake field distribution of the Biscayne tugboat.

Figure 15. Wake field coordination analysis results.

Figure 16. Wake field coordination analysis results.

Figure 17. Blade camber and pitch distribution. (a) Camber distribution of the blade. (b) Pitch distribution of the blade.

Figure 18. Hydrodynamic optimization design results. (a) Efficiency under two scenarios. (b) Torque coefficients under two scenarios.

Table 1

B Series propeller data sheet.

Serial Number	B-Series	P/D	EAR	Z
1	4-100	1	1	4
2	4-85	1	0.85	4
3	4-70	1	0.7	4
4	4-55	1	0.55	4
5	4-40	1	0.4	4
6	3-65	1	0.65	3
7	5-75	1	0.75	5
8	6-80	1	0.8	6
9	7-85	1	0.85	7
10	4-70	0.5	0.7	4
11	4-70	0.6	0.7	4
12	4-70	0.8	0.7	4
13	4-70	1.2	0.7	4
14	4-70	1.4	0.7	4

Table 2

Forecast model input and output parameter.

Input Data	Output Data
Number of blades, Z	Thrust coefficient $C{T}^{*}$
Blade area ratio, EAR
Timing, $\beta$	Torque coefficient $C{Q}^{*}$
Pitch ratio, P/D

Table 3

Neural network training area and boundary.

Area	Training Boundary	Output Boundary
A	−50 to 60	−45 to 50
B	30 to 150	50 to 140
C	130 to 240	140 to 230
D	220 to 330	230 to 315

Table 4

Mesh convergence.

Grid Quantity	3.32 Million	5.35 Million	7.38 Million	9.40 Million	11.86 Million	14.32 Million	16.77 Million
CT	0.0927	0.0932	0.0934	0.0935	0.0935	0.0936	0.0936
CQ	0.0113	0.0114	0.0115	0.0115	0.0115	0.0115	0.0115

Table 5

Reference number of impacts per propeller revolution for each ice class.

Ice Class	PC1	PC2	PC3	PC4	PC5	PC6	PC7
$N_{class}$	$21\times {10}^6$	$17\times {10}^6$	$15\times {10}^6$	$13\times {10}^6$	$11\times {10}^6$	$9\times {10}^6$	$6\times {10}^6$

Table 6

Preliminary design table of blade strength.

Blade Thickness	Thickness	Hydrodynamic Load	Design Ice Condition	Non-Design Ice Condition
Blade Root Section	280 mm	176 mm	213 mm	260 mm
0.6 R Section	180 mm	92 mm	142 mm	175 mm

Table 7

Parameter information for propeller design task.

Ice Class	Length Between Perpendiculars	Draft	Propeller Diameter	Number of Blades
PC6	43 m	4 m	2.59 m	4
Shaft power	Full power free-running speed	Full power bollard pull speed	Rotation direction	Hub diameter ratio
1864.25 kW	305 RPM	245 RPM	Right-hande	0.25

Table 8

Propeller initial design working condition table.

n	250	260	270	280	290	300	310
$w_t$	0.245	0.245	0.248	0.246	0.244	0.242	0.240
J	0.180	0.173	0.166	0.160	0.155	0.150	0.146
K q	0.034	0.031	0.027	0.024	0.022	0.020	0.018

Table 9

Main parameters of blades in the preliminary design.

Parameter	Expanded Area Ratio	Skew Angle	Pitch Ratio	Blade Root Pitch	Blade Root Section Chord Length	0.6 R Section Chord Length
Value	0.7	15°	0.73	1890 mm	805 mm	991 mm

Table 10

Preliminary design table of strength for tugboat propeller blades.

Blade Thickness	B-Series Propeller	Kader Design	Design Ice Condition	Non-Design Ice Condition
Blade Root Section	89 mm	135 mm	157 mm	191 mm
0.6 R Section	51 mm	78 mm	105 mm	129 mm

Table 11

Design variables and value ranges.

Design Variable	Initial Value	Lower Bound	Upper Bound
0.25 R camber–chord ratio	0.010	0.009	0.011
0.5 R camber–chord ratio	0.029	0.026	0.032
0.9 R camber–chord ratio	0.024	0.022	0.027
1.0 R camber–chord ratio	0.012	0.011	0.013
0.25 R pitch–diameter ratio	0.620	0.558	0.682
0.5 R pitch–diameter ratio	0.761	0.684	0.837
0.9 R pitch–diameter ratio	0.694	0.624	0.763
1.0 R pitch–diameter ratio	0.575	0.517	0.632
0.25 R skew angle	−7.261	−6.535	−7.988
0.5 R skew angle	3.732	3.358	4.105
0.9 R skew angle	14.111	12.700	15.522
1.0 R skew angle	15.320	13.788	16.852
0.25 R chord–diameter ratio	0.357	0.321	0.393
0.5 R chord–diameter ratio	0.430	0.387	0.473
0.9 R chord–diameter ratio	0.308	0.278	0.339

Table 12

Propeller geometry table.

	Chord Length (mm)	Pitch (mm)	Camber (mm)	Thickness (mm)	Rake (mm)	Skew (°)
0.25	999.9	1444.7	9.3	150	45.6	−6.6
0.3	1040.5	1566	15.2	143	54.6	−4.0
0.4	1092.1	1751.8	24.1	128.5	72.8	0.4
0.5	1109.2	1871.5	28.7	114.2	90.9	4.0
0.6	1098.6	1941.1	29.8	100	109.3	6.9
0.7	1059.7	1974.9	28.2	85.7	127.4	9.6
0.8	981.0	1972.5	24.6	71.5	145.6	12.0
0.9	819.4	1909.6	18.1	57.2	163.7	14.3
1.0	0	1415.2	0	42.7	182.1	16.6

References

1. Kerwin, J.E.; Lee, C.S. Prediction of Steady and Unsteady Marine Propeller Performance by Numerical Lifting-Surface Theory; Society of Naval Architects and Marine Engineers: Jersey City, NY, USA, 1978.

2. Mishima, S.; Kinnas, S.A. Application of a numerical optimization technique to the design of cavitating propellers in nonuniform flow. J. Ship Res.; 1997; 41, pp. 93-107. [DOI: https://dx.doi.org/10.5957/jsr.1997.41.2.93]

3. Chen, B.; Stern, F. Computational fluid dynamics of four-quadrant marine-propulsor flow. J. Ship Res.; 1999; 43, pp. 218-228. [DOI: https://dx.doi.org/10.5957/jsr.1999.43.4.218]

4. Bosschers, J.; Vaz, G.N.V.B.; Starke, A.R.; van Wijngaarden, E. Computational Analysis of Propeller Sheet Cavitation and Propeller-Ship Interaction. Proceedings of the RINA Conference “MARINE CFD2008”; Southampton, UK, 26–27 March 2008.

5. Song, K.; Gong, J.; Ge, S.; Sun, C.; Yang, H.; Zhang, W. Ship energy-saving device: Influence of the scale effect and hull roughness on the resistance reduction effect of an interceptor. Ocean. Eng.; 2023; 287, 115875. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2023.115875]

6. Muscari, R.; Di Mascio, A.; Verzicco, R. Modeling of vortex dynamics in the wake of a marine propeller. Comput. Fluids; 2013; 73, pp. 65-79. [DOI: https://dx.doi.org/10.1016/j.compfluid.2012.12.003]

7. Kumar, P.; Mahesh, K. Large eddy simulation of propeller wake instabilities. J. Fluid Mech.; 2017; 814, pp. 361-396. [DOI: https://dx.doi.org/10.1017/jfm.2017.20]

8. Søntvedt, T. Propeller blade stresses, application of finite element methods. Comput. Struct.; 1974; 4, pp. 193-204. [DOI: https://dx.doi.org/10.1016/0045-7949(74)90082-0]

9. Sunnersjo, C.S. Stress Analysis of Highly Skewed Propeller Blades Using the Finite Element Method. Proceedings of the SNAME Propeller and Shafting Symposium; Virginia Beach, VA, USA, 15–16 May 1984.

10. Hu, Z.K.; Gui, H.B.; Xia, P.P.; Wang, J.X. Finite element analysis of marine propeller strength under ice load. Ship Eng.; 2013; 35, 4.

11. Young, Y.L. Fluid—Structure interaction analysis of flexible composite marine propellers. J. Fluids Struct.; 2008; 24, pp. 799-818. [DOI: https://dx.doi.org/10.1016/j.jfluidstructs.2007.12.010]

12. Ye, L.Y.; Guo, C.Y.; Wang, C.; Wang, C.H.; Chang, X. Strength assessment method of ice-class propeller under the design ice load condition. Int. J. Nav. Archit. Ocean Eng.; 2019; 11, pp. 542-552. [DOI: https://dx.doi.org/10.1016/j.ijnaoe.2018.09.008]

13. Kang, H.; Kaiqiang, W.; Chunyu, G.; Chao, W.; Liyu, Y.; Zeping, W. Mechanism analysis of propeller-ice contact and rapid prediction of ice loads. Cold Reg. Sci. Technol.; 2023; 216, 103987. [DOI: https://dx.doi.org/10.1016/j.coldregions.2023.103987]

14. Boswell, R.J. Static Stress Measurements on a Highly-Skewed Propeller Blade; Naval Ship Research and Development Center: Washington, DC, USA, 1969; 3247.

15. Zhao, B. Strength Study of Large Skew Propellers; Huazhong University of Science and Technology: Wuhan, China, 2003.

16. Yang, X.; Cheng, E. Stress Test Study on Blades of Large Skew Propellers. Ship Ocean Eng.; 2004; 1, pp. 6-9.

17. Su, Z.; Pan, J.; Zhang, S.; Wu, S.; Yu, Q.; Zhang, D. Characterizing dynamic deformation of marine propeller blades with stroboscopic stereo digital image correlation. Mech. Syst. Signal Process.; 2022; 162, 108072. [DOI: https://dx.doi.org/10.1016/j.ymssp.2021.108072]

18. Ducoin, A.; Barber, R.B.; Wildy, S.J.; Codrington, J.D.; Baker, A. Experimental evaluation of the use of embedded fiber Bragg gratings to measure steady and unsteady flow-induced marine propeller blade deformation. Ocean Eng.; 2023; 281, 114889. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2023.114889]

19. Bose, N.; Veitch, B.; Doucet, J.M. A design approach for ice class propellers. Trans. Soc. Nav. Archit. Mar. Eng.; 1998; 106, pp. 185-211.

20. Doucet, J.M.; Bose, N.; Veitch, B. Selection of Scantlings for an Ice Class Propeller: A Case Study for the MV IKALUK. Proceedings of the 18th International Conference on Offshore Mechanics and Arctic Engineering; St. John’s, NL, Canada, 11–16 July 1999.

21. Kader, R.D. The Design of a Propeller for a US Coast Guard Icebreaker Tugboat; Taylor Naval Ship Research and Development Center: Bethesda, MD, USA, 1975.

22. Liu, P.; Islam, M.F.; Searle, S.; MacNeil, A.; Prior, A. Design and Optimization of an Ice Class Propeller Under Shallow Water, Semi-Tunnel Hull and Heavy Load Condition. Proceedings of the International Conference on Offshore Mechanics and Arctic Engineering; Hamburg, Germany, 4–9 June 2006; pp. 675-683.

23. Carlton, J. Marine Propellers and Propulsion; Butterworth-Heinemann: Oxford, UK, 2018.

24. Madayag, A.F. Metal Fatigue: Theory and Design; John Wiley and Sons: Hoboken, NJ, USA, 1969.

25. Osgood, C.C. Fatigue Design: International Series on the Strength and Fracture of Materials and Structures; Elsevier: Amsterdam, The Netherlands, 2013.

26. Timoshenko, S.; Maccullough, G.H. Elements of Strength of Materials; D. Van Nostrand Company: London, UK, 1935.

27. Wang, C.; Han, K.; Sun, C.; Guo, C. Optimization Design and Parameter Analysis of Marine Propellers. J. Huazhong Univ. Sci. Technol. (Nat. Sci. Ed.); 2020; 48, 6.

28. Van Lammeren, W.P.A.; van Manen, J.V.; Oosterveld, M.W.C. The Wageningen B-screw series. Soc. Nav. Archit. Mar. Eng.; 1969; 77, 43.

29. Atkinson, P. Prediction of marine propeller distortion and stresses using a superparameteric thick-shell finite-element model. Trans. Rina; 1973; 115, pp. 359-375.

30. Van, M.; Van, O.; Lewis, E.V. Resistance, Propulsion, and Vibration. Principles of Naval Architecture; 2nd ed. Society of Naval Architects and Marine Engineers: Jersey City, NJ, USA, 1988; Volume 2.

31. Zhou, J. Research on Design Methods of Large Skew Propellers; Harbin Engineering University: Harbin, China, 2011.

32. Ye, L. Wake-Adapted Design of Marine Propellers; Harbin Engineering University: Harbin, China, 2015.