HULL SHAPE OPTIMIZATION OF SMALL UNDERWATER VEHICLE BASED ON KRIGING-BASED RESPONSE SURFACE METHOD AND MULTI-OBJECTIVE OPTIMIZATION ALGORITHM

Summary Small underwater vehicles have unique advantages in ocean exploration. The resistance and volume of a vehicle are key factors affecting its operation time underwater. This paper aims to develop an effective method to obtain the optimal hull shape of a small underwater vehicle using Kriging-based response surface method (RSM) and multi-objective optimization algorithm. Firstly, the hydrodynamic performance of a small underwater vehicle is numerically investigated using computational fluid dynamics (CFD) method and the value range of related design variables is determined. The mesh convergence is verified to ensure the accuracy of the calculation results. Then, by means of the Latin hypercube sampling (LHS) design of simulation, the Kriging-based RSM model is developed according to the relation between each design variable of the vehicle and the output parameters applied to the vehicle. Based on the Kriging-based RSM model, the optimal hull shape of the vehicle is determined by using Screening and MOGA. As results, the vehicle resistance reduces and volume increases obviously.


Introduction
The ocean is an important asset for the sustainable development of human society [1], and investment in marine surveying and naval equipment manufacturing is increasing in various countries. Underwater vehicles are suitable for surveying, photographing, and collecting information in deep water areas. Therefore, their application fields are gradually expending, including underwater observation, underwater reconnaissance, aquaculture, underwater military and so on. Small underwater vehicles have unique advantages in both civilian and military applications due to their outstanding concealment [2,3]. Xiaodong Xing, Haixia Gong, Xiujun Xu response surface method and multi-objective optimization algorithm 112 Due to the small size of underwater vehicles, they cannot carry excess energy supply [4]. As a result, it is necessary for small underwater vehicles to overcome as little resistance as possible when moving underwater. Optimizing the hull shape of small underwater vehicles can effectively reduce their resistance from water. As an important tool for underwater vehicles analysis [5][6][7][8], computational fluid dynamics (CFD) has been introduced into the underwater vehicle resistance analysis process because of its efficient, accurate and fast calculations. Rans equations and three different turbulence including RKE, SSTKO and RSM were used for the description of turbulent flow in [9]. The authors studied the influence of bow on three ship types. The results showed that the experimental and numerical results are similar and the total resistance decreased 7% using bulbous bow. Ignacio et al. [10] used Navier-Stokes equations with Favre's average and k-ε model in SolidWorks Flow Simulation. The design variables of an autonomous underwater vehicle (AUV) of torpedo body type was optimized and the expected low resistance was verified through the experiment. The authors of [11] used RANS equations and SST k-ω model. And a parallel multidisciplinary optimization design method was proposed, which used a design of experiment (DOE) of Latin hypercube sampling (LHS) design combined with a radial basis function model to improve the shape of underwater vehicles. Li et al. [12] used RANS equations and SST k-ω model in ANSYS Fluent. They proposed a simplified shape optimization (SSO) strategy for blendedwing-body under-water gliders (BBWUGs) by using a Kriging-based constrained global optimization (KCGO) method to optimize the aerofoil shape. And the results showed that the lift-resistance ratio was improved. Zhao et al. [13] used the realisable k-ε two-layer model to describe the influence of turbulence in STAR-CCM+. Lackenby method was used to modified the main hull of the trawler and NSGA-II algorithm and Sobol + Tsearch algorithm are used to optimize the trawler. The reliability of optimization results was verified by comparing the experimental and numerical results. Wang et al. [14] used RANS equations and SST k-ω model in ANSYS CFX. The free-form deformation (FFD) parametric method and the surrogate-based optimization (SBO) algorithm were proposed to optimize the hydrofoil shape. Finally, the lift-resistance ratio of the optimized hydrofoil was improved by more than 90%. Vasudev et al. [15] used RANS equations and SST k-ω model in CFD solver and used genetic algorithm (GA) and CFD solver to optimize the hull shape at a single velocity. Liu et al. [16] used RANS equations and SST k-ω model in OpenFOAM. The multidisciplinary design optimization (MDO) method is applied to the multiple objective MDO design of the heavierthan-water underwater vehicle (HUV). They also combined the non-dominated sorting genetic algorithms (NSGA-II) and the Kriging model to search the optimal performance of the HUV. Wu et al. [17] used RANS equations in OpenFOAM. An optimization method which combines the lifting line theory and the continuous concomitant method was invented to solve the conflict between shape design and propeller design. The numerical results showed that the AUV resistance was reduced about 17% at rated velocity. In addition to the optimization of traditional hull shapes, new types of structures have also been designed through bionics. For example, an innovative design underwater glider was proposed in [18] which used SST k-ω model in SolidWorks Flow Simulation and mimicked a whale shark. The results showed that the whale shark-like glider has less resistance, higher lift-resistance ratio and significant gliding capability in practice.
For multivariate problems, it is difficult to accurately establish the connection between input and output parameters in practical engineering. Response surface quadratic model is widely used for analysis and improvement because of its low process complexity and high fitting accuracy. Sun et al. [19] optimized the space utilization and range of AUV by using response surface method (RSM) model and demonstrated the effectiveness of using RSM for shape optimization based on experiments. Liu et al. [20] investigated the strength and stability of carbon fibre composite AUV shell with reinforcement ribs based on finite element analysis (FEA) and RSM. Kriging is an exact interpolation method that can replace the multiple quadratic regression equations in RSM because of its powerful and flexible properties [21][22][23]. Kriging-based RSM has higher accuracy and reliability. To obtain the optimal solution, the regression equations analysis of the RSM model is necessary. Guan et al. [24] developed a optimization procedure which integrating a parametric geometric model, a numerical analysis of the total resistance, a RSM model and an evolutionary optimization algorithm to illustrate the effectiveness of the optimization by using a mapping strip. Guan et al. [25] also used the particle swarm optimization (PSO) algorithm and sequential quadratic programming (SQP) algorithm to solve the approximation model established by the second-order RSM and verified the effectiveness of the automatic optimization design method with a 7m unmanned surface vehicle (USV).
In this paper, a small underwater vehicle is designed with the Myring curve as a reference to adsorb the target. Therefore, the hull shape of the small underwater vehicle is controlled by a parametric Myring curve with six design variables. In view of the vehicle's special hull shape, it is necessary to constrain the range of each design variable by CFD so that the optimized hull shape can be found rapidly and accurately. The convergence of the mesh is studied to ensure the correctness of the calculation. A technique integrating Krigingbased RSM and multi-objective optimization algorithm is proposed to optimize the vehicle's design variables to reduce the resistance and increase the volume when moving underwater. The calculation results show that the optimization method is applied to optimize the hull shape of the vehicle effective.

Optimization Process
In this paper, SolidWorks and ANSYS are used to optimize the hull shape of a small underwater vehicle. The optimization process is shown in Fig. 2. First, the head and tail parts of the vehicle is built using Myring curve which is controlled by equation. Rotation and stretching are adopted to create the vehicle entity. Changing five input design variables to realize the model reconstruction. Then, the vehicle model is imported into ANSYS Workbench. In geometry module, basin is created, the faces of velocity input and pressure output are determined. The structured mesh is generated in ICEM module and its convergence is verified by studying the first layer height above the vehicle and the number of grid elements. Finally, analytical settings are adjusted and the vehicle resistance is calculated in Fluent module.
After creating these modules, the values of five design variables is changed respectively to calculate and obtain the constraints of each design variable. Then, set the range, name and unit of the input and output variables in the parameter module. The Latin hypercube sampling design is set to obtain a certain number of sample values. Finally, the Kriging-based response surface is established and the goal is to minimize resistance and maximize volume. In optimization, Screening and MOGA are selected to optimize the output variables.

Overall Design
A new small underwater vehicle designed in this paper aims to adsorb the bottom of the ship which sailing on the water. After the vehicle is equipped with the sub-robot, it can corrode and damage the ship, to reduce the combat readiness of enemy ships and further consume the maintenance cost of ships. The overall design of small underwater vehicle is shown in Fig. 1. The body is divided into four parts: adsorption device, sealed pressureresistant bin, expansion module and propeller. The vehicle is driven by propeller, and the magnet controlled by sensor ejects and adsorbs the target when approaching it.
The design parameters of small underwater vehicle are shown in Table 1. Where d is the maximum diameter of the vehicle, L is the total length of the vehicle, vm is the maximum navigation velocity of the vehicle, tm is the minimum navigation time of the vehicle and hm is the maximum working depth of the vehicle.

Adsorption Device Design
Considering the influence of the adsorption device on the vehicle resistance, the part of adsorption is hidden in the device. Since the contact surface of the adsorbed target is plane, the magnet surface is selected as plane. The untriggered and triggered states of the annular adsorption device is shown in Fig. 3. The proximity sensor located in the head of the vehicle makes the electromagnet lose its magnetic force so that the magnet is ejected by the spring and adsorb the target when approaching it.
Claw adsorption device is shown in Fig. 4. The adsorption type and working principle are like the annular adsorption device. The claw adsorption device has better hydrodynamic performance compared with the annular adsorption device. Its magnets are evenly distributed on the four claws which are controlled by the electromagnet. Four claws are required to be opened until all magnets are coplanar to achieve the best adsorption effect. At the same cost, the efficiency of claw adsorption device is lower than the annular adsorption device. Therefore, the annular adsorption device is selected in the study.

Sealed Pressure-Resistant Bin Design
There are some electronic components in sealed pressure-resistant bin, including control circuit board, wires, electronic governor and signal receiver. The sealing ring is selected to seal and prevent the electronic components from being eroded by water. PMMA is selected to make sealed pressure-resistant bin whose geometric model is shown in Fig. 5. The design parameters of sealed pressure-resistant bin are shown in Table 2. Where ls is the length of bin, ds is the internal diameter of bin, ts is the thickness of bin and E is the modulus of elasticity.  The maximum working depth of the vehicle is shown in Table.1. To protect the internal devices, the critical pressure of the vehicle is designed to be greater than the interaction pressure underwater. The critical pressure, Pc, is calculated by Eq 1 which is referenced in [ where ρ is the density of water and g is gravity. P is calculated at about 0.62 MPa. It can be seen from the results that Pc is greater than P.

Expansion Module Design
As shown in Fig. 6, there is a cavity with diameter of 80mm and length of 100mm in the expansion module. The expansion module is used to store a sub-robot which is ejected from the vehicle and adsorbs the ship to corrode and damage the target. The length of the sub-robot is 65mm, and the maximum diameter is 30mm. A door opening mechanism is designed on the module surface. The length of the door is 85mm and the width is 50mm. Considering the space size of the module, the door opening mechanism is simplified. A torsion spring is installed at the connection between the door and the module. Choosing the electromagnet likes adsorption device to control the state of the door. The door is opened by the torsion spring and the sub-robot is ejected when the electromagnetic switch triggered.

Hull Shape Design
The hull shape of underwater vehicle is usually designed to be spherical, ellipsoidal or cylindrical. The head and tail of underwater vehicle are a combination of these shapes. We need to install a functional mechanism in the vehicle head and a propeller in the tail. Therefore, the vehicle's space utilization and hydrodynamic performance should be considered at the same time in the design. In particular, the improvement of hydrodynamic performance is the focus of our attention.
In general, the hydrodynamic performance of the underwater vehicle depends on its hull shape. Myring and NURBS are often used in shape design [26,27]. In the unmanned vehicles, the rotary Myring curve is most used, such as Remus [28]. In this paper, Myring curve is used to build the profile of small underwater vehicle.
As shown in Fig. 7, it is a schematic diagram of Myring curve with six design variables. The head curve equation of the vehicle is: Xiaodong Xing, Haixia Gong, Xiujun Xu response surface method and multi-objective optimization algorithm 118 where d is the maximum diameter of the vehicle, a is the length of the head and n is the dimensionless variable, which controls the profile of the head. Especially, the larger n is, the more circular the head is. The equation of the vehicle tail curve is given as follows: where b is the length of the middle, c is the length of the tail and θ is the tail semi-angle, which controls the profile of the tail. Similarly, the larger the θ is, the rounder the tail is.

CFD Analysis
As shown in Fig. 9, the initial model of the vehicle is established in SolidWorks, and a is 70 mm, n is 1, b is 480 mm, c is 100 mm, θ is 50°, d is 90 mm. The volume of the vehicle is 0.00358 m 3 by measuring.

Controlling Equation and Computation Domain
The commercial software ANSYS Fluent 2021 R1 is selected as the CFD solver. The underwater environment is regarded as an incompressible fluid and is a turbulent field. The simulation takes Reynolds-averaged Navier-stokes (RANS) equations as the controlling equation. The author of [15] compares the errors of resistance coefficients of different turbulence models, taking SST k-ω turbulence model has least errors. The RANS have an advantage of being computationally efficient. Its average mass and momentum transfer equations are as follows [29]: where ui is the time-averaged value of the velocity component, t is the time, p is the timeaveraged value of the pressure, ρ is the fluid density, μ is the dynamic viscosity coefficient, ij ρuu  is the Reynolds stress term, and Sj is the generalized source term of the momentum equation.
Computation domain is established in Fig. 10. Considering the development of the downstream wake. The domain with length of 6L, width of 10d and height of 10d is built, where turbulence can be considered. The vehicle is placed at 2L to the left boundary and 4L to the right boundary. The left wall is set as velocity inlet with a fixed horizontal velocity value and the right wall is set as pressure outlet.

Grid Convergence Study for Resistance
A structured mesh is generated for CFD simulation using ANSYS ICEM. The detailed grid elements of the vehicle are shown in Fig. 11. The grids density is high around the vehicle. The y + for the first layer of the elements from the vehicle is around 1. This y + is suggested for the use of the SST k-ω turbulence model [30]. The kinematic viscosity, vk, can be calculated by where μ is the dynamic viscosity of the fluid. The Reynolds number, Re, is calculated as where v is the velocity of the vehicle. The skin friction coefficient on the vehicle is Wall shear stress, τw, is given by Friction velocity, Uτ, is calculated by And the first layer height above the vehicle, y, is calculated by In all calculation, a value of y=0.01 mm is taken to satisfy y + =1.

Fig. 11 Mesh details
The numerical results can be accurate if they are not affected by the number of grid elements. As shown in Table 3, there are three different mesh resolutions to evaluate the influence of the number of grid elements on the vehicle resistance. The fine mesh with 5.0 million cells, the medium mesh with 2.7 million cells and the coarse mesh with 0.9 million cells are calculated at v=2 m/s. It can be seen from the table that the results of the fine mesh and the medium mesh are given approximately. Therefore, the vehicle resistance will not change a lot by increasing the number of grid elements. Considering the efficiency of calculation, the mesh with 2.7 million cells is selected for the following calculation.

Numerical Calculation
The vehicle subjected resistance is normally divided into pressure resistance, viscous resistance and wave making resistance [31]. The wave making resistance is so small that it can be ignored when the sailing depth is more than one-third of the vehicle length. Therefore, only the pressure resistance and viscous resistance are considered in this paper.
In simulation, the vehicle is analysed for resistance at the depth of 60 m and v=0.5 m/s, 1 m/s, 1.5 m/s and 2 m/s respectively. The resistance conditions at different velocities are shown in Table 4. Where Ft is the total resistance of the vehicle in the water, Fp is the pressure resistance and Fv is the viscous resistance. As shown in Fig. 12, it can be seen that Fp is larger than Fv and Ft is equal to the sum of Fv and Fp when the vehicle navigates underwater.
The main reason why Fp is much larger than Fv is that the vehicle head lacks a streamlined part for the design of the sensor and the adsorption module. The curve shows that Ft, Fv and Fp all increase with the raise of the vehicle's velocity and Ft is proportional to the quadratic of the velocity. Ft and Fp change significantly with the increase of velocity and Fv is little affected by velocity. The static pressure distribution image is shown in Fig. 13.  To investigate the effect of each design variable on the resistance, the variables a, n, b, c, and θ are used as variables to perform CFD of the vehicle at 2 m/s respectively. The relationship between Ft and design variables of the vehicle is shown in Fig. 14. It can be seen from the image that each variable affects Ft. Ft decreases rapidly when a increases between 20 mm and 35 mm. Ft shows a tendency to increase when a increases between 35 mm and 70 mm. Therefore, the range of a is taken from 25 mm to 55 mm. Ft increases when n increases between 0.2 and 1.6. Therefore, the range of n is taken from 0.2 to 1.2. Ft grows linearly when b increases. Obviously, the smaller b is, the Ft smaller is. However, most electronic components are stored in the middle part of the vehicle and b should be taken to satisfy the constraint of L shown in Table 1. Therefore, b is taken from 400 mm to 550 mm. Ft decreases rapidly when c increases between 30 mm and 50 mm. Then, Ft fluctuates when c increases between 50 mm and 130 mm. The minimum Ft is located in c=80 mm. Therefore, the range of c is taken from 50 mm to 120 mm. The relationship between Ft and θ is similar to the relationship between Ft and c. The range of θ is taken from 10° to 50° to obtain an accurate θ under this trend.

Optimization
In the previous section, the discrete analysis on different design variables of the vehicle is performed to estimate the resistance of the vehicle and determines the range of design variables. In this section, Workbench 2021 R1 is applied to connect the Fluent module and Response Surface Optimization module to build an automated modelling analysis. A Krigingbased response surface method is used to build a mathematical model and the optimal hull shape is found by Screening and MOGA used in [20,32]. The Screening optimization method uses a simple approach based on sampling and sorting. It supports multiple objectives and constraints as well as all types of input parameters. The MOGA method is a variant of the popular NSGA-II based on controlled elitism concepts. It supports multiple objectives and constraints and aims at finding the global optimum.

Kriging-Based Response Surface Method
Response surface method (RSM) is a statistical method to solve multivariate problems by using reasonable experimental design methods and obtaining certain data through experiments. It uses a multiple quadratic regression equation to fit the functional relationship between factors and response values and to seek the optimal variables.
There are various experimental designs which can perform response surface analysis such as Latin hypercube sampling (LHS) design, central composite design (CCD) and boxbehnken design (BBD) [33,34]. In this paper, the design variables are range-based and LHS design is chosen as an experimental design method.
LHS design is used to generate random values of each design variable and to make each variable appear only once in any row or column of the experimental design matrix. Therefore, the randomized LHS design method is used to design a more homogeneous and complete area which is covered with a minimum number of sample points. LHS design method is suitable for the Kriging-based RSM which is an interpolation method. CCD and BBD design methods are valid for traditional polynomial based RSM. Xiaodong Xing, Haixia Gong, Xiujun Xu response surface method and multi-objective optimization algorithm 126 The Kriging method is a spatial interpolation method based on random positions. Its theoretical approach was proposed in 1970 and was later applied to the engineering discipline. The Kriging model has high prediction accuracy and good nonlinear fitting ability and can replace the second-order polynomial as a meta model to obtain higher accuracy of the response surface. The Kriging model expression is [35]: where X is the input variable set, y(X) is the predicted response value, ( ) , F β X is the regression model, β is the regression coefficient and Z(X) is the measure of the deviations from the regression model. The covariance of Z(X): where R(xi,xj) is correlation function between any two values in the input variables, and its expression is: where θk which is used to estimate the correlation relationship between the test points is the parameter of the correlation function. The distance between two points is changed by changing the value of θk, which changes the correlation between two points and can generate a smoother surface.
In Eq 13, is estimated using Eq 16: ( ) (16) where f is a column vector, R is the correlation matrix of training data and y is a row vector.
The output values of final design points in the Kriging-based RSM is: where r(x) is the correlation vector between the untried x and the sampled points. Five design variables which control the vehicle's hull shape are analysed in previous section. Considering the problem of variables interaction occurs in the stage of model reconstruction, four design variables a, n, b, c, and θ are introduced to Kriging-based RSM.

Response Surface Optimization
The Kriging-based response surface of the vehicle hull shape is established and optimized. The developed response surface model is used to analyse any two input variables effects on two output variables by three-dimensional response surface. The relationship between output variables and input design variables is shown in Table 5 by automated modelling analysis. The generated Kriging model is verified with ten sample points different from the set of designs used for generation. The goodness of fit of response surface is shown in Fig. 15. It can be seen from the fitness image that the verification points for V are distributed on the midline and the verification points for Ft are distributed compactly on both sides of the midline. The root mean square error of V is 2.749E-06 and the root mean square error of Ft is 0.0263. The response surface established in ANSYS Workbench is shown in Fig.  16. It can be seen from the image that b and c have a more apparent effect on V, a and n have a more apparent effect on Ft. These conclusions are consistent with ANSYS Workbench sensitivity analysis results, as shown in Fig.17.
We determine the optimization strategy and establish a multi-objective optimization model.
In this paper, Screening and MOGA are used to optimize the hull shape of the vehicle. The objective functions are to maximize V and to minimize Ft. Fig.18 shows the obtained Pareto front for V and Ft. Each red point represents a Pareto optimal solution and all red points constitute the Pareto front. Ft plays more important role in optimization. Therefore, the objective importance of Ft is set higher than V in decision support process. By this decision technique, five candidate points (C1, C2, C3, C4 and C5) are selected from Fig.18 and the values of the objective functions are shown in Table 6 and Table 7. The candidate points are calculated by CFD to verify the effectiveness of optimization results. As the results shown in Table 6 and Table 7, the maximum errors of V and Ft are 0.11% and 1.54% using Screening, and they are 0.19% and 1.48% using MOGA, which is acceptable in the optimization. As shown in Table 8, compared with the initial hull shape of the vehicle, V increased 1.40% and Ft decreased 8.98% using Screening and V increased 3.07% and Ft 9.52% decreased using MOGA. Both optimization methods can obviously improve the hull shape of the vehicle and the hull shape of the vehicle in MOGA is superior to in Screening. Therefore, the design variables are round to a=45 mm, n=0.8, b=506 mm, c=96 mm and θ=49° based on the principle of proximity. The generated optimal model is shown in Fig. 19.  Fig. 15 The goodness of fit of response surface

Conclusion
A hull shape optimization design of a small underwater vehicle which is used to adsorb targets is carried out with regard to minimizing the resistance and maximizing the volume in this paper. Firstly, Myring curve with six design variables is used to control the vehicle profile. Then, the value range of design variables and resistance are calculated by CFD. What's more, a Kriging-based RSM model which can predict the resistance and volume is established to optimize the hull shape of the vehicle. Multi-objective optimization algorithm is proposed to find the optimal hull shape: 1. The small underwater vehicle is parametrically modelled using the Myring type as a reference and the profile of the vehicle is controlled by six different variables. This Myring curve can provide a reference for the shape design and optimization of underwater vehicles in the future. 2. The vehicle requires separate functional modules in the head, so the range of each variable based on experience may not meet the design criteria. Therefore, it is important for the optimization that the influence of each variable on the resistance is analysed by using the control variable method to obtain the range of each variable. 3. Through solving the established Kriging-based RSM model using MOGA, the hull shape of the vehicle is optimized obviously compared with initial model. The Kriging-based RSM model proposed in this paper solves the multivariate consistency constraint and decoupling problems in the overall optimization of small underwater vehicle, and provides a theoretical basis and technical support for the design optimization of similar underwater vehicles.