A ship propeller design methodology of multi-objective optimization considering fluid–structure interaction

ABSTRACT This paper presents a multi-objective optimization methodology that applies the Non-dominated Sorting Genetic Algorithm-II(NSGA-II) to propeller design, and realizes Fluid-Structure Interaction (FSI) weak-coupling based on Panel Method (PM) and the Finite Element Method (FEM). The FSI iterative process and the convergent pressure coefficient distribution and pressure fluctuation of HSP (a propeller installed on a Japanese bulk freighter – Seiun-Maru) are numerical calculated. The FSI results turn out to have higher precision than those without FSI. The appropriate optimization parameters are chosen after studying five cases. The Sobol method, a global Sensitivity Analysis (SA) algorithm, is used to quantify the dependence of objectives and constraints on the input parameters. In the multi-objective optimization methodology, efficiency, unsteady force, and mass are chosen as optimum objectives under certain constraints. Effectiveness and robustness of the methodology are validated by running the program starting from four different random values, which all improve the objectives and converge to the similar results. The proposed multi-objective optimization methodology could be a promising tool for propeller design to help improve design efficiency and ability in the future.


Introduction
The increasing demand for ship propellers meeting synthetically requirements has been putting great pressure on propeller designers. These propellers must offer higher efficiency with lower vibrations, less noise, and lighter mass. To meet these seemingly mutually exclusive objectives, designers have been paying more attention to contemporary soft computing techniques and optimization algorithm. The contemporary soft computing techniques have been successfully applicated in many real-life fields, such as image identification (S. Zhang & Chau, 2009 September), time series forecasting (Chau & Wu, 2010;Wang, Chau, Xu, & Chen, 2015;Wu, Chau, & Li, 2009), Input Variable Selection (IVS) algorithms development (Taormina & Chau, 2015a), etc. The ability to estimate the hydrodynamic performance and structure response has been greatly improved according to the rapid development of computing techniques.

Fluid-structure interaction
The theories to predict hydrodynamic performance can be approximately divided into three categories: Momentum theory or Actuator Disc Method (ADM), potential flow theory and viscous flow method based on Computational Fluid Dynamics (CFD) tools. ADM, constituting a widely employed analysis/design tool both in its potential and CFD version (Bontempo & Manna, 2013, 2016, 2017Conway, 1998), focuses more on the global effects of the propeller after the hull. Potential flow theory, including Lifting Line Method (LLM) (Çelik & Güner, 2006;Coney, 1989;Mishkevich, 2006), Lifting Surface Method (LSM) (Kawakita & Hoshino, 1998;Kerwin & Lee, 1987;Yamasaki & Ikehata, 1992) and Panel Method (PM) (Hess, 1990;Karim, Suzuki, & Kai, 2004;Kerwin, 1987;Suciu & Morino, 1976), is a commonly utilized tool to analyze or design ship propeller. Although CFD tools based on RANS or LES have been widely applied, the computational cost is too high to be practical for use in propeller optimization, especially when the number of iterations is large. The Panel Method (PM), developed from the Boundary Element Method (BEM), has been widely used in propeller design, which represents a compromise between computational accuracy and cost. As for structure response prediction, FEM (Atkinson & Glover, 1988) is rapid developed and utilized in this paper for overcoming the fairly fundamental shortcomings of the beam theory and shell theory.
The ability to simultaneously estimate the fluidstructure interaction (FSI) effect can not only shorten design period, but also be important to hydrodynamic and structural design. Lin and Lin (1996) developed an iterative procedure, which coupled FEM and the noncavitating lifting surface method, but did not consider the fatigue strength of the propeller. Liu (2000) focused on investigating the physical nature of FSI problems by decomposing the system response into multiple frequency/wave number bands. Liu accomplished this by developing a multiple-scale Reproducing Kernel Particle Method (RKPM) and applying it to the 2D airfoils. Benaroya and Wei (2000) concentrated on accurately describing the response of structures to unsteady fluid loadings, including waves, currents and vortex shedding. Isaac and Iverson (2003) developed an automated FSI analysis program, which couples Fluent and ABAQUS. In our methodology, the FSI is obtained by iterating hydrodynamic loads predicted by potential PM and structure deformation analyzed by FEM automatically.

Propeller optimization
Optimization algorithms have been applied to propeller design by many researchers. Benini (2003) developed a multi-objective design program based on the evolutionary algorithm to optimize a B-screw propeller, while the open water performance is calculated using regression formulas. Han, Bark, Larsson, and Regnstrom (2006) chose the Dynamic Hill Climbing method as the optimization mechanism to optimize both the efficiency and induced pressure fluctuation, while the structural response is not taken into consideration. Zhao and Yang (2010) optimized pitch and camber distributions, obtaining the surface pressure distributions more uniformly without reducing the propeller efficiency. Cai, Ma, Chen, Qian, and Zhang (2014) applied an improved particle swarm algorithm to optimize and analyze propeller skew distribution. Klesa (2014) optimized the circulation distribution including the influence of viscosity. Meanwhile, optimization algorithm has been improved by Chau. He introduced multi-sub-swarm particle swarm optimization (MSSPSO) (J. ) and multi-objective fully informed particle swarm (MOFIPS) (Taormina & Chau, 2015b) to the standard particle swarm optimization (PSO) and got better performance.
After reviewing previous similar work, we find out that most of the propeller optimization researches are focused on single object optimization and few combine FSI in the optimization program. To improve design efficiency and ability, we propose and develop a multi-objective optimization design methodology using NSGA-II as the optimization algorithms and realizing FSI based on PM and FEM. A careful selection of design parameters, objectives and constraints is made by practical experience and global sensitivity analysis. The results turn out that the objectives are improved. Robustness of the methodology is validated by running the methodology four times starting from different random values, which all converge to the similar results. The design process is shown in Figure 1.

Panel method
The panel method is used to predict the steady and unsteady hydrodynamic performance of the propeller. The surface of propeller and the wake vortex are divided into a series of small units, and each unit is approximated by a hyperboloidal quadrilateral panel. Some assumptions are made for practicality. 1: the geometry of wake surface changes linearly. The wake surface shed to downstream from trailing edge along the direction of the mean camber surface and its pitch linearly changes to average value of the geometrical pitch distribution of the propeller with respect to the rotating angle; and 2: perturbation potential ϕ, potential jump ϕ and − → V 0 · − → n Q are all assumed to be constants within each panel. Then, the integral equation can be written as liner equations: Where N and N W stand for the numbers of the panels on the propeller blade surface and the wake surface respectively, and δ is the Kronecker function. B ij , C ij , W ij are influence functions: where S j and S l are the surface panels on the propeller surface and the wake surface. R ij and R il are the distances from the boundary point on the ith panel to integration points on S j and S l . All these influence functions are calculated by the analytical formulation proposed by Morino and Kuo (1974). The Kutta condition is applied to achieve convergence of the iterative linear Eq. (1). To avoid some unexpected results drawn by direct numerical derivation, the Yanagizawa and Kikuchi (1982) proposed a method to determine the velocity on the propeller surface, and then the pressure on the propeller surface is calculated using the Bernoulli equation. By integrating the pressure on all propeller surfaces we can obtain the thrust and torque of propeller: where n xi ,n yi ,n zi are components of unit normal vector of ith panel, x i ,y i ,z i are coordinates of ith panel centroid, S i is the area of ith panel and p i is the pressure of ith panel.
Then the thrust and torque coefficients and efficiency can be calculated as Meanwhile, to overcome the panel method's shortcoming of not taking viscosity into account, we added the viscous elements by combining the Prandtl-Schlichting law with the frictional resistance formula. The formula is The Reynolds number of the blade section is generally defined as But we determined the Reynolds number by expression Nc k=1 v k l k is used here to replace the conventional v r0 c r in the definition expression of Reynolds number. This is considered be able to include the effects of the geometry of the blade section and the non-uniform velocity.

Finite element method
As previously mentioned, FEM is widely used to predict the structural response of propeller (Young, 2007;Young & Liu, 2007). The blades are assumed to be made of homogeneous, isotropic and linear elastic material. In the blade-fixed coordinate system, the dynamic equation of motion can be written as follows: are the matrixes of structural mass, damping, and stiffness, respectively.  It should be noted that the Coriolis force and damping are assumed to be negligible due to the small deformation and static analysis.

Weak coupling of FSI
The interaction of fluid-structure is considered to be weak coupling in this paper. The calculation process is: Firstly, the fluid loading of the undeformed blade is calculated by PM. Then, the deformation of the blade is calculated according to the fluid loading by FEM. Secondly, the fluid loading of the deformed blade is again calculated, and then the deformation is predicted based on the last deformed blade according the difference of last two fluid-loadings. The process, which takes the 'feedback-effect' into account, is repeated until the deformation and fluid loading is converged. The details are shown in Figure 2. In Figure 2, Q is the set of endpoints of panels on the blade; ε 1 and ε 2 are the convergence indexes. F i,Q and U i,Q mean the fluid loading and deformation at the nodes of the i th iteration. During the iteration, 3D linear interpolation is utilized to achieve data exchange between PM and FEM.

FSI of HSP blade
The weak coupling effect of FSI is carefully studied in this subsection. HSP is a propeller installed on a bulk freighter called Seiun-Maru in Japan, and the ship wake field is shown in Figure 3. The parameters of HSP and the experimental status are listed in Table 1. There are 20 panels in both chord wise and radial direction, which are shown in Figure 4. As the appropriate value of ε 1 and ε 2 are unknown now, the iteration numbers are temporary set as 5. The dimensionless F i , max and U i , max during the iteration are given in Figure 5.
The FSI process is evaluated in Figure 5. The thrust coefficient (K T ), the first order of thrust coefficient ( K T ), torque coefficient (K Q ), efficiency (η), first order (f 1 ) and second order (f 2 ) of natural frequency during iteration are shown in Figure 6. The results are nondimensionalized by the values when Iter = 0. From   Figure 6, we can see the final results are not quite different from those without FSI because of the small deformation caused by metal material, but still not the exact same. While during iteration, K T and K Q first decrease obviously and then increase a little cased by the vibration swing back and forth. The final results are a litter smaller than those without FSI because of deformation. The same trends of K T and K Q result in small change of η. Meanwhile, during the iteration, K T and natural frequency (f 1, f 2 ) almost keep unchanged. From Figure 5 and Figure 6, ε 1 and ε 2 is appropriate to set as 0.005* F 0,max and 0.005* U 0,max hereinafter. The pressure coefficient distribution at different rotation angles (θ = 0, 90,180, Figure 7) and the pressure coefficient fluctuation at different chord position (x/c = 0.25, 0.6, 0.8, Figure 8) are calculated and compared with the experimental data and the numerical results of Hoshino (1998). The results of this paper considered FSI, while Hoshino didn't. From Figure 7 and Figure 8 can see that the unsteady hydrodynamics, no matter in the pressure coefficient distribution or pressure fluctuation, considering FSI turn out to have higher precision than Hoshino which is without FSI. Meanwhile, the uncertainties of numerical results may be raised by the linear wake assumption and the low-order potential with the equivalent-plate viscosity correction of the panel method.

Optimization parameters study
As an improved version of the Nondominated Sorting Genetic Algorithm (NSGA, Srinvas & Deb, 1994), NSGA-II was proposed by Deb, Agrawal, Pratap, and   Meyarivan (2000). It has alleviated all three of the main criticisms of NSGA: high computational complexity of non-dominated sorting, lack of elitism and the need for specifying the sharing parameter. NSGA-II is chosen as the optimization algorithm in this paper. The optimization parameters (population size, number of generations, crossover probability and distribution index) are evaluated using the following optimization    Table 2. In order to construct a suitable multi-objective optimization function, we provided minimize torque as the optimization objective function, which can effectively accelerate the convergence and save the computing resources without changing the nature of the optimization problem. Thrust and torque are set as the constraints to make sure the optimization results can satisfy the matching conditions of the enginepropeller. Five cases of different optimization parameters are carried out. The η/η 0 calculation histories in all cases are shown in Figure 9. The optimization solution here means the weight coefficient of each objective is equal. Table 3 lists the optimization solution of the cases. Although design parameters R k /R k0 , f /f 0 and c/c 0 of Case3-1, Case3-2 and Case3-3 are different, the results of optimization objectives are almost the same. This may because the values of D/D 0 and P/P 0 are similar, which are the most important influence on the objectives learned from global sensitivity analysis in Sec.4.2. Compared Case3-3, Case3-4 and Case3-5, decrease crossover probability and distribution index may get worse results. Hereinafter, considering calculation accuracy and cost, the values of design parameters are set to be the same as Case3-3.

Optimization model
The optimization model is given in this subsection. The objectives are to maximize propeller efficiency η, minimize the mass W and first order of unsteady thrust coefficient K T . To avoid changing the design operation greatly, thrust and torque are limited in a certain range. Meanwhile, deformation, safety redundancy, cavitation performance and avoiding blade-hull resonance (evaluated by dimensionless parameters F 1 , F 2 , S 1 and S 2 ) are also chosen as the constraints. The design parameter, constraints and objectives are listed in Table 4. All the parameters with subscript 0 are the values of HSP.
We take natural frequency into consideration in order to set the foundation for future research to study the effect of interaction among the hull, shaft and propeller. The first and second orders of natural frequency of the hull are defined as f H−1st and f H−2nd , and those of blade are Figure 9. η/η0 calculation history in the five cases. (a) η/η 0 calculation history in Case3-1 (b) η/η 0 calculation history in Case3-2 (c) η/η 0 calculation history in Case3-3 (d) η/η 0 calculation history in Case3-4 (e) η/η 0 calculation history in  f B−1st and f B−2nd . We define four parameters to evaluate: We define the safety redundancy SR as: Where n is rotational speed, σ b is permissible tensile stress, and σ MAX_N is the maximal Mises stress of the blade when the propeller is at maximum speed. Meanwhile, σ /ξ of decompression coefficient and cavitation number is defined as: Where p v is vapor pressure and it's 174 (kgf/m 2 ) at 15°C of water. And p 0 is calculated as Where p a is atmospheric pressure (10330 kgf/m 2 ), γ is the weight density of water and h s is submersible depth of propeller hub (h s = 1.75D in this paper).

Global sensitivity analysis
In this subsection, a simulation study is carried out to quantify how constraints and objectives in this model depend on the input factors. The Sobol method (Saltelli & Bolado, 1998;Saltelli, Tarantola, & Chan, 1999) is utilized as the global sensitivity analysis method, and the design matrix is generated by the Optimal Latin Hypercube technique. The main liner effects of several outputs are listed in Table 5. Figure 10 is the Pareto Charts including the quadratic effects of inputs. In these charts, for convenience, the parameters (for example D) stand for the nondimensional parameters (D/D 0 ). The black solid boxes stand for the positive effects and the white hollow ones stand for the negative effects.
From these charts, the η is affected by several parameters, of which the interactions effects are included. P, D, P-D are the most three significant positive factors affecting thrust, while are the negative ones for K T and the σ MAX . The W is only affected by t, D and D-t. P has the most significant effect on transformation, followed by D-P and D, which have almost the same effect coefficients.

Automated optimization
In this subsection, the automated optimization methodology is carried out four times, along with design parameters starting from different random values in the range. HSP is chosen as the parent form. And the optimization parameters are set the same as Case 3-3 in Sec.3.2. Figure 11 shows the Pareto front solutions and the optimization solution of the four cases. We make the weight coefficient of different objectives equal. The three objectives are plotted with respect to the most important input factor according to the sensitivity analysis. The results of both Pareto solution and optimization solution in the four cases turn out to be almost the same. Table 6 lists the optimization results, we can see that the four cases get the almost the same objectives, although design parameters are different, especially in R k /R k0 , f/f 0 and c/c 0 . Case 4-2 have lower value of K T / K T0, sacrificing the value of η/ η 0 compared with  the other three cases. The differences of values of W/ W 0 seem very slight in different cases.

Optimization comparison with/without FSI
In this subsection, two cases without FSI are carried out. After we get the results only optimize the hydrodynamics, we calculated the structure response. Case4-5 is the Pareto-front solution in the hydrodynamic optimization and Case 4-6 is the superior solution in all the feasible results which can satisfy the constraints listed in Table 4. Compare Table 7 with Table 6, we can draw the conclusion that if the optimization process only considers the aspect of hydrodynamics, the results may not satisfy the constraint condition of structure response.
Meanwhile, choosing optimized result from optimization results without FSI not only adds work but also has worse results than those with FSI.

Conclusions
In this paper, we propose and develop a multi-objective optimization design methodology that applies NSGA-II to propeller design, and that takes FSI into optimization process. By comparing with the experimental data of pressure coefficient distribution and pressure fluctuation, we validate the necessity to consider FSI as it has higher precision and helps improve design ability and efficiency. In the optimization design methodology, we maximize propeller efficiency, minimize blade mass  and also minimize the unsteady thrust coefficient, which have all been improved. The global SA helps to quantify the dependence of the objectives and the constraints on design factors. Only one case, HSP, is studied in this paper, which seems to limit the application of the methodology. Future work we will try to apply the methodology to more propellers design work, like composite propeller and propellers for different purposes. We will also focus on taking hull-shaft-propeller interaction into consideration, which can greatly enrich the methodology.

Disclosure Statement
No potential conflict of interest was reported by the authors.

Funding
Support for this research is provided by NSFC No. 51079158 (National Natural Science Foundation of China).