Research on bulbous bow optimization based on the improved PSO algorithm

In order to reduce the total resistance of a hull, an optimization framework for the bulbous bow optimization was presented. The total resistance in calm water was selected as the objective function, and the overset mesh technique was used for mesh generation. RANS method was used to calculate the total resistance of the hull. In order to improve the efficiency and smoothness of the geometric reconstruction, the arbitrary shape deformation (ASD) technique was introduced to change the shape of the bulbous bow. To improve the global search ability of the particle swarm optimization (PSO) algorithm, an improved particle swarm optimization (IPSO) algorithm was proposed to set up the optimization model. After a series of optimization analyses, the optimal hull form was found. It can be concluded that the simulation based design framework built in this paper is a promising method for bulbous bow optimization.


Introduction
In recent years, optimization methods have been widely used in the design of aircraft, ships, machinery and other engineering structures. From the optimization of design parameters to the selection of structures, more and more design problems have been solved. Optimization methods can be divided into three categories: unconstrained optimization, constrained optimization and intelligent optimization. Intelligent optimization is a heuristic optimization algorithm with high speed and strong practical applicability. It includes the genetic algorithm (GA), ant colony optimization (ACO) algorithm, particle swarm optimization (PSO) algorithm and so on. As a genetic algorithm, the PSO algorithm (Kennedy and Eberhart, 1995) is also an optimization algorithm based on iteration. It is easy to be implemented and has few parameters to adjust such that it can be applied to different kinds of optimization problems. In a new multi-scheme, a multi-object group decision-making method proposed by Hou et al. (2012), the individual evaluation value was adjusted by the PSO algorithm, and the best shiptype was obtained. Wang and Liu (2009) built an evalu-ation model for collision risk factors and found the optimal collision rate by the PSO algorithm. With the development of the PSO algorithm, a series of new algorithms have emerged at a historic moment improving optimization efficiency and the performance of solutions. Hou et al. (2016) used the PSO algorithm with changed acceleration coefficients in a ship hull optimization. Liu et al. (2009) introduced the method of dynamically changing inertia weight in the velocity evolutionary equation and developed a modified PSO algorithm to study the parameter identification of a creep constitutive model of rock. By changing the method of initialization, using particle updates and adding variation factors, Luo et al. (2009) used an improved PSO algorithm in a wind power system. Li (2012) calculated the inertia weight of the PSO algorithm by self-adaptation, proposed an improved PSO algorithm to evaluate the objective function and obtained the best hull form.
Although the PSO algorithm is a relatively mature method, it easily falls into the local extreme. In the later stages, the convergence rate is slow, and the accuracy is poor (Li et al., 2014). Therefore, in order to improve the global search ability and local improvement of the PSO algorithm, an improved particle swarm optimization (IPSO) algorithm with random weight and the crossbreed algorithm is proposed. First, the PSO and IPSO algorithms are used to evaluate each of the four functions in order to verify the applicability of the IPSO algorithm. Second, an overset mesh technique is used for the mesh generation and the Reynolds-averaged Navier-Stokes (RANS) method is used to calculate the total resistance of the hull. The arbitrary shape deformation (ASD) technique is introduced to modify the bulbous bow, and the optimization mathematical model is built with the IPSO algorithm. After the completion of the optimization, the best hull form is obtained.

Governing equation
The whole computational domain uses the continuity equation and the RANS equation as the control equations: and are the fluid density, static pressure, fluid viscosity, Reynolds stresses, and body forces per unit volume, respectively.

Turbulence model
The RNG k-ε model is derived by the mathematical method of renormalization with instantaneous N-S equation. The transport equations can be written as: where k is the turbulence kinetic energy; ε is the turbulence dissipation rate; the quantities a k and a ε are the inverse effective Prandtl numbers for k and ε, respectively; μ eff is the effective dynamic viscosity; G k is the generation of turbulent kinetic energy by the mean velocity gradients; G b is the generation of turbulent kinetic energy by buoyancy; Y M represents the contribution of the fluctuating dilatation in compressible turbulence; and C 1ε , C 3ε and C 2ε are empirical constants.

Volume of fluid method
The free surface is captured by the volume of fluid (VOF) method (Hirt and Nichols, 1981). This is a surface tracking method fixed under the Euler grid. It simulates the multiphase flow model by solving the momentum equation and the volume fraction of one or more fluids. Within each control volume, the sum of the volume fractions of all the phases is 1. As to Phase q, its equation is: , q=1 represents the air phase, q=2 represents the water phase, a q is the volume fraction of q-th phase, and a q =0.5 is defined as the interface of air and water.

Mesh generation
With further research in the field of ships, the motion of the ship must be taken into account in the calculation of the hull resistance. This is rather difficult in simulations of the motion with traditional structured grids or unstructured grids. The overset mesh technique is a new type of the mesh generation method. It can generate high-quality grids and may be more accurate for solving large amplitude motions of the ship. It has been widely used thus far. In the overset mesh method, grids are divided for each part of the model and then nested into the background grids. First, the hold points and interpolated points need to be marked. Second, overlapping units can be removed by digging the hole. Then, interpolation is carried out to complete data exchange in the interface. Finally, the entire flow field is calculated. Two acceptor cells are shown by dashed lines, one in the background mesh, and the other in the overset mesh (CDadapco, 2014) as shown in Fig. 1.
The fluxes through the cell face between the last active cell and the acceptor cell are approximated in the same way as between two active cells. However, whenever the variable value at the acceptor cell centroid is referenced, the weighted variable values at the donor cells are substituted: where α i are the interpolation weighting factors, φ i are the values of the dependent variable φ at donor cells N i and sub- script i runs over all donor nodes of an interpolation element (denoted by the green triangles in Fig. 1). In this way, the algebraic equation for the cell "C" in Fig. 1 involves three neighbor cells from the same mesh (N 1 to N 3 ) and three cells from the overlapping mesh (N 4 to N 6 ). The coefficient matrix of each solved equation (for both the segregated and the coupled solution method) is updated accordingly to ensure that the equations can be solved up to the round-off level of the residuals. There are many interpolation methods, and the linear interpolation uses shape functions to connect the center of the variable grids. The interpolation unit is transferred from one interpolation element to another from the center of the mesh. Although this method may have low computational efficiency, it is more accurate for resistance prediction. In this paper, the overset mesh is used to divide the computational domain with the linear interpolation, as shown in Fig. 2.

Boundary condition
According to the requirements of the overset mesh, the whole model needs two individual blocks, i.e. the background block and the overlapping block.
(1) Background block: The front, top and bottom boundaries are selected as the velocity inlets. The back boundary is selected as the pressure outlet, and both sides of the tank are selected as the symmetric planes.
(2) Overlapping block: The right side of the cuboid is set as the symmetric plane, and the rest of the surface is set as the overset region. The hull is set as the wall.

Calculation process
By taking David Taylor Model Basin (DTMB) model 5512 as an example, the total resistance of the hull in calm water is calculated using the RANS method.
(1) Pre-calculation: The geometric model is established, the grids are divided, the mesh quality is checked, and boundary conditions are set.
(2) Calculation process: The hull is fixed, and the speed is added along the x-direction at the front boundary. The continuity equation (1) and the RANS equation (2) are selected as the control equations of the whole computational domain. The RNG k-ε model is selected to close the equations. The VOF method is used to simulate the free surface. The SIMPLE algorithm is selected to solve the equations, and to start the iteration, the ship speed is taken as the initial velocity. The methodology is illustrated in Fig. 3.

Geometry regeneration
Based on the B-spline technique, the arbitrary shape deformation (ASD) technique is a method to change the shape of different models using the commercial software (Sculptor). It is necessary for the ASD volume to be set up outside the body with many control points and connections. When the control points are moved, the shape of the relevant areas can be changed. A new surface can be created with movements of the control points, for which third-order continuity can be maintained. This can insure the smoothness of the new model, even under the conditions of large deformation. This direct deformation method makes the optimization of a complex geometry possible. The specific steps are as follows.
(1) The ASD volume is built around the CAD model with different control points and connections. The deformation volume can be finer or coarser, depending on the de-  ZHANG Sheng-long et al. China Ocean Eng., 2017, Vol. 31, No. 4, P. 487-494 489 sired control of the shape change.
(2) Control points are changed by the movements of the positions and directions.
(3) The new model is created according to the model generation algorithm with modified control points. Take DTMB5512 as an example. First, the ASD volume is established, and the surface control points 1-3 of the bulbous bow are taken as design variables, as shown in Fig. 4. Second, the positions and directions from points 1 to 3 are changed. Finally, the new surface is generated.

Particle swarm optimization (PSO) algorithm
The particle swarm optimization (PSO) algorithm is a recently developed evolutionary algorithm. Through iteration, the PSO algorithm can find the optimal solution starting from a random solution. Because it lacks crossovers and mutations, it is much easier to use than the genetic algorithm. It can find the global optimal solution following the optimal value of the current result and evaluate the quality of the solution by fitness. It also has a strong processing capability and good robustness.
In d-dimensional space, the i-th particle's velocity and position can be written as and X i = , respectively. pbest denotes the local best solution, and gbest denotes the global optimal solution. In each iteration, the particle updates itself by tracking pbest and gbest. After finding these two optimal values, the velocity and position of each particle is updated according to the following equations: where ω is the inertia weight coefficient; c 1 and c 2 are positive-valued acceleration coefficients called the cognitive parameter and social parameter, respectively; r 1 and r 2 are random variables between 0 and 1. Generally, the number of particles takes from 20 to 40, the acceleration constants are equal to 2, and the range of ω is based on the specification of the problem.

Improved particle swarm optimization (IPSO) algorithm
The inertia weight coefficient ω is a parameter that influences the trade-off between global and local searches within the solution domain. In order to improve the precision of the PSO algorithm, adjusting ω is very important.
Large ω values may avoid local optimal solutions, and small ones can effectively accelerate convergence speed. ω can be calculated by the constant method, linear decrement method, self-adaptation method and so on. In this paper, ω follows a random distribution, so that it can overcome the instability caused by the linear decline of the two aspects.
The formula of ω can be written as: where ψ is the mean value of the random weight; σ is the variance of the random weight; N (0, 1) is the random variable of the standard normal distribution; ψ min is the minimum value of the random weight; ψ max is the maximum value of the random weight; rand (0, 1) is the random variable between 0 and 1. According to the crossbreeding probability, a number of particles are put into the pool in each of the iterations. Particles in the pool are hybridized randomly by twos, and the same numbers of child-particles are generated. Parentparticles are then replaced by child-particles. The positions of child-particles are determined by those of parent-particles: where ξ is a random variable between 0 and 1, and parent 1 (x) and parent 2 (x) are the positions of the parentparticles. The velocities of child-particles can be calculated using the following formula: where parent 1 (v) and parent 2 (v) are the velocities of the parent-particles. The IPSO algorithm is developed with MATLAB (R2012b) and can be expressed as follows: (1) Initialize the positions and velocities of a set of particles randomly.
(2) Evaluate the fitness value of each particle. Store the position and fitness value of each particle in pbest. Store the lowest of all the pbests values in gbest. (3) Use Eqs. (7)-(8) to update the velocities and positions of the particles.
(5) Compare the position and fitness value of each particle to find the best position, and update gbest.
(6) Drawing on the idea of the crossover algorithm, the positions and velocities of child-particles are calculated according to Eqs. (11)-(12).

Optimization and verification
To verify the applicability of the IPSO algorithm, four functions shown in Eqs. (13)-(16) are studied. The PSO and IPSO algorithms are used to find the minimum value of each of the four functions. The parameters of the two algorithms are shown in Table 1. ; (13) (14) ; (15) Table 2 provides a comparison of the optimization results with the two algorithms. Fig. 5 shows the iterative pro-cess. The optimization results of the IPSO algorithm are much closer to the theoretical results with a faster convergence rate in the initial optimization. The PSO algorithm gets trapped in a local optimum (Figs. 5a, 5c and 5d), especially for multi-modal functions f 3 (x) and f 4 (x). Although both algorithms produce almost the same results in Fig. 5b, the IPSO algorithm has a faster convergence speed. It can be concluded that the IPSO algorithm has better performance and faster convergence speed for optimization, especially for multi-optimization problems.

Hull form optimization
According to Section 5.4, the IPSO algorithm has a better global search ability and higher efficiency than the PSO algorithm. Therefore, in this section, the IPSO algorithm is used to find the optimal bulbous bow in calm water.

Geometry and optimization strategy
(1) Geometry model The modified area of DTMB5512 is shown in Fig. 6. The major characteristics are shown in Table 3.  PSO  40  2  2  0.5  -----IPSO  40  2 2 -0.8 0.5 0.2 0.9 0.2 Notes: N is the number of particles, P c is the crossover probability, and S p is the pool ratio. ZHANG Sheng-long et al. China Ocean Eng., 2017, Vol. 31, No. 4, P. 487-494 491 (2) Optimization objective The total resistance coefficient C cT in calm water has a minimum at Fr=0.28. C cT can be obtained from: where F is the total resistance of the hull, ρ is the fluid density, and v hull is the hull speed.
(3) Design variables The design variables are a 11 , a 12 , a 21 , and a 22 . The control parameters are shown in Table 4. The control positions are shown in Fig. 4.
, where ∆ opt is the variable of the new hull, and ∆ org is the variable of the parent hull.
(5) Optimization processes The optimization flow chart used in this paper is shown in Fig. 7, and it includes five steps: ① According to the characteristics of the parent hull, use the ASD technique to modify the bulbous bow. Then use the CAD technique to build a series of new hulls and computational domains. Change the design variables with the IPSO algorithm.
② Calculate the displacement of each new hull and compare it with the constraints. If the conditions are met, continue the calculation. Otherwise, return to Step ①.
③ Divide the computational domains according to the method given in Section 3.
④ Calculate the total resistance coefficients in calm water according to the method described in Section 4, and save the results.
⑤ Repeat Steps ①-④ until the stopping criterion is met. Then the new hull can then be found.

Comparison of the results and discussion
After the completion of the optimization, a new hull with smaller total resistance is obtained. Table 5 shows the comparison of the optimization results. It can be seen from Table 5 that the total resistance decreases 5.49% for the optimized hull-B. Fig. 8 shows the comparison of the total resistance coefficients C cT for the optimized hull-B and parent hull with the change of Fr in calm water (Gui et al., 2001a(Gui et al., , 2001bLongo and Stern, 2005). Fig. 8 shows that the total resistance decreases at all speeds and that it decreases more in the design speed and high speed. Fig. 9 presents body-plans of optimized hull-B and the parent hull, which show that the hull lines are smooth and slightly concave. Fig. 10 compares the wave profile for the optimized hull-B and the parent hull at y/L=0.105 (h represents the depth of the water). It can be shown that the amplitude of the waves has been reduced, which indicates a reduction in the total resistance for the optimized hull-B. Figs. 11 and 12 present a comparison of wave patterns and static pressure for the optimized hull-B and the parent hull, re-    spectively. Optimized hull-B has a smaller splash than the parent hull near the bow, as shown in Fig. 11. The pressure distribution of the bulbous bow undergoes a significant change, as shown in Fig. 12.

Conclusion
In this paper, the total resistance of the hull in calm water was selected as the objective function, and the overset mesh technique was used for the mesh generation. The ASD technique was used to change the shape of the bulbous bow. Under the displacement and constraints, the IPSO algorithm was proposed to evaluate the objective function. The obtained hull lines were smooth and close to the parent hull, and the total resistance of the obtained hull form decreased 5.49% compared with the parent hull. The results show that the established platform is feasible for the hull form optimization and the IPSO algorithm converges faster and performs better than the PSO algorithm.