MULTI-DIMENSIONAL SURROGATE BASED AFT FORM OPTIMIZATION OF SHIPS USING HIGH FIDELITY SOLVERS

Summary Surrogate (metamodel) based optimization has numerous potential applications in the field of naval architecture. It is aimed here to establish a methodology for the aft form optimization for minimum viscous resistance, thus the present study is focused on the aft form where the viscous effects become dominant. It is necessary to solve this problem within acceptable time span from practical naval architectural point of view which requires metamodeling techniques currently under investigation. Accordingly, the present paper investigates the metamodeling ability of the Kriging interpolation and attempts to explore its capabilities and limitations in the aft form optimization from viscous resistance point of view. As metamodeling techniques become more widely used, their constraints are more apparent. Especially in highly nonlinear design spaces, the effect of dimensionality should be taken into consideration. Taking all those factors into account, the present paper is to examine the capabilities of Kriging and to establish the learning performance in terms of RMS error, correlation coefficient and required number of training points according to selected optimization algorithm for multidimensional ship design problem. The results show that, at least 5% reduction in viscous pressure drag can be attained by the present optimization methodology.


Introduction
There is a growing need to reduce the fuel consumption of the ships by hull form optimization for economic considerations and to meet the International Maritime Organisation's (IMO) regulation related with Energy Efficiency Design Index (EEDI) [1] which urges ship-owners to reduce CO2 emissions. Hull form design is a complex problem which requires solutions under the influence of many constraints that affect each other positively or negatively, hence there are many studies which show the complexity of this issue in the literature [2][3][4][5]. Particularly aft form design has its own difficulties, as the flow in this region of the hull is dominated by viscous effects on the one hand and there are several objectives in this region to be minimized or optimized such as viscous resistance, flow uniformity, propeller design, cavitation, noise and vibration on the other. For this reason, many constraints should be studied (see for example; [6]) and meantime it is important to choose the correct objective functions to optimize the aft form [7]. In addition, there are a couple of studies which focused on the aft form design to increase the propulsive efficiency of the hull [8][9][10][11].
Difficulties with the experimentation such as required time, financial matters, scaling issues on the one hand and the improvements in computing power and numerical techniques, on the other, lead the engineers to include CFD methods in their optimization [12][13][14][15][16][17][18]. Recently, with the practical advantages of CFD, surrogate modelswhich are able to replicate the CFD with a reasonable error and without requiring too much computing timeare employed in engineering optimization studies. They can replace computationally expensive simulation processes in ship design problems and costly model experiments. Today, no unique metamodel technique has developed on which a consensus is reached among naval architects. Every technique has its own strengths and weaknesses. Thus many research groups are trying either to develop tuning parameters of the methods that is currently in use or introduce new metamodels [19][20][21][22]. In the literature there is variety of different studies that uses metamodeling from chemical industry [23] to aerospace [24] or bio-inspired robots [25]. One can also find several comprehensive review papers that describe different types of metamodeling and their advantages/disadvantages in every aspect [26][27][28][29][30]. Thus, metamodeling in engineering is widely used in design optimization studies and/or in design space exploration.
In the present study, Kriging interpolation is chosen for aft form optimization problem. Kriging has many advantages; such as being a proper approach to multi-dimensional problems and its prediction capability as compared to other metamodeling techniques [31][32] and being able to produce error estimates. On the other hand, Kriging technique is still developing by many research groups from all around the world [33][34][35][36][67][68]. Those development studies are rather concentrate on the sampling types [37], on the tuning parameter exploration such as variogram adaptation [38] or on adding some descriptive new information into the algorithm such as gradient/hessian enhanced types [39][40][41][42][43]. Thus, based on the above discussions, Kriging metamodeling is chosen as a viable tool in the present hull form optimization study.
The aim of this study is to develop an algorithm to optimize the vessel's aft form by using metamodels. Accordingly, the study presents a new design procedure which accomplishes the optimization of the aft form for minimum viscous resistance via surrogate modeling on which a very limited number of studies are recorded in the literature. For example, a different technique can be cited which employs an adjoint based optimization method, Rung [6]. Another study in the literature uses generic hull shapes for modifications for larger re-design region and uses viscous flow analysis, [44][45]. The study given in [44], used a similar approach as we used here but focuses on a global optimization procedure by employing Lackenby's form transformation technique whereas the present approach uses a point-based hull form transformation. Thus, the most important output of the study was able to find and eliminate the problem of local flow separation, with the advantage of focusing on a local design patch.
Design of experiments (DoE) required for the training of the metamodel is achieved here by a viscous flow solver. In this context, 15-20% of the ship's length is selected in the aft part as the optimization zone and the wetted surface geometry of the prescribed zone is defined by a limited number of control points. Subsequently, variant hull forms obtained by the variation of the control points were subjected to viscous flow analysis. Variations (transformations) of the hull form are obtained by Akima surface generation method introduced in [46]. The dimensionality reduction in optimization problems is a different 87 research topic, and some important examples can be found in [47][48]. In this study, the dimension of the design space is 6 and DoE are developed by using these 6 parameters. The optimization study is then carried out for minimum viscous pressure resistance. The form obtained from optimization procedure is then tested by viscous solvers. Although, 5% reduction seems to be a limited gain because of strict constraints, it shows the promising ability of the present methodology in aft form optimization of ships.
The remainder of the paper is structured as follows: The next section provides theory followed by design of experiments and then results and our conclusions with future work comments will be given.

Surrogate Modelling: The Kriging Interpolation
Kriging was originally developed in geostatistics by Danie G. Krige [49]. The mathematical basis for Krige's idea was then developed by Matheron [50]. Following Matheron, Kriging metamodels were applied to the input/output data sets of deterministic simulation models. Kim et al. [51] show that Kriging is a superior method when applied to nonlinear, multimodal problems. Kriging technique appears to be better than the other surrogate modelling methods such as Radial Basis Function (RBF) and Artificial Neural Network (ANN) as shown by Peri [63].
Kriging is mostly presented as a way of 'modelling the function as a realization of a stochastic process'. The reason is, if we want to make a guess at a point x in the domain before sampling any point, the prediction will have an uncertainty. Modelling this uncertainty resembles a random variable Y(x) that is normally distributed with the mean,  , and the variance, 2  , [52]. This means that the function has a value  which varies between ±σ.
This implies that closer data points will tend to have nearly the same function values () i y x , () j y x , where x1 and x2 are data points. In other words, this interpolation technique assumes correlation between closer observed data as can be statistically modelled in the following expression: Expression (1) has an intuitive property that if ij x = x then the correlation is 1, and as ij − →  xx the correlation tends to 0 for each d (number of dimensions). The interpolation parameters l  are related with correlation change by moving in the th coordinate direction and if this parameter tends to have greater values, function values change rapidly even if the distances between the data points are small. The p makes the model smooth when it is near 2. On the one hand, the additional parameters in (1) make the model more complex, but on the other hand enhance the accuracy and capability of the prediction. The uncertainty of the function values at the n (number of sample points), by using a random vector, can be written as: (2) and covariance matrix equal to; where R vector is a nxn matrix which is represented in (1)  -as: To simplify the likelihood estimation, we maximize (4) by taking its natural logarithm and derivative; (5) and end up with the optimal values of  and 2  as functions of R; Then substituting (6) and (7) into the (5) gives the concentrated log-likelihood function ignoring the constant termas: p is specified as 2as described aboveand  is searched according to the procedure given in Jones [52]. Now, the same procedure is applied to a pseudo observation which means a new point added to the previous data set for a new estimation. The estimator will be derived by evaluating the quality of the new estimation which means again maximizing the likelihood function with the model parameters obtained so far. By adding (n+1) th observation to the data as The correlation matrix for the augmented data R ; * Covariance is a measure of the correlation between two or more sets of random variables which one can derive correlation from covariance as; ( , ) cov( , ) /  (10) and by using (10) in the concentrated likelihood function and solving the equationagain by maximizing the relation by taking derivatives and setting it to zerofor the new * y which gives us the standard Kriging predictor (see [52][53][54]) as: The metamodel is tested with respect to Root Mean Square Error (RMSE) and correlation coefficient (r).
RMSE in (12) and the correlation coefficient in (13) provide a quantitative measure of model accuracy and also is useful to have a numerical understanding of the quality of the surrogate model.

Design of Experiments
The design of experiments was accomplished by a two-step process which contains a hull form transformation procedure and a sensitivity analysis to properly locate the control points defining the surface patch of the optimization region.
Main particulars of the hull form in model scale are given in Table 1. The geometry of the proposed hull form can be seen in Figure 1. The hull form used in this study is not a wellknown benchmark form, but a form with problematic flow characteristics around the aft body.

Hull form variations and Geometric transformation
For creating the data set of 280 variant hull form, cubic interpolation technique of Akima [55] is employed. There are some other interpolation techniques (see [46] and [56]). But for our case we focused on cubic interpolation techniques; among them are natural, monotonic and Akima's surface polynomials. Natural cubic interpolation is a piecewise curve on each data point and second derivative is assumed to be zero at the beginning and the end of the curve. As a result, the curve oscillates around the boundary points. The monotonic cubic interpolation uses monotonicity of the points by modifying the tangent values of the data to prevent oscillation. Akima method is a continuously differentiable interpolation, built from piecewise third order polynomials and applicable to successive intervals of the given points. The slope of the curve is locally determined at each given point, by the coordinates of five points centered on the studied point. In consequence, this spline type creates a smooth and natural curve between the waypoints and always passes directly through them. In Figure 2, the comparison of various cubic interpolation techniques is made. Akima surface polynomials are obviously superior to other cubic interpolation methods and moreover, according to our experience, it provides smooth natural curves/surfaces and it is computationally efficient and reduces oscillatory effects [46]. In conclusion, Akima surface representation is chosen to obtain the variants of the aft form.

RANS Computations
The governing equations are the continuity equation; and the Navier-Stokes equations approximated by the RANS equations for the steady, threedimensional, incompressible flow: P demonstrates the mean pressure, ρ the density and ν the kinematic viscosity of the fluid where the velocity U can be decomposed into mean velocity i U component and fluctuating velocity part as in (16), The k-ε turbulence model is applied in order to simulate the turbulent flow around the hull. This turbulence model is applicable when there are no high pressure changes along the hull and quite economical in terms of CPU time [58], compared to, for example, the SST turbulence model, which increases the required CPU time by nearly 25% [59]. During the analyses, Reynolds stress tensor is calculated as in (17) while Cμ is an empirical constant (Cμ =0.09). k is the turbulent kinetic energy and ε is the turbulent dissipation rate. Also, transport equations are solved for k (18) and ε (19): where comprehensive explanations can be found in [60].
The solver uses a finite volume method which discretizes the governing equations. A second order convection scheme was used for the momentum equations and a first order temporal discretization was used. The flow equations were solved in a segregated manner. The continuity and momentum equations were linked with a predictor-corrector approach. The pressure field is solved by using SIMPLE algorithm which is based on pressure-velocity coupling. All governing equations are discretized using a cell based finite volume method and the advection terms are discretized with a first-order upwind interpolation scheme.
The computational domain was discretized by three dimensional finite volume cells and appropriate mesh structure was created around the hull using hexahedral elements and trimmer meshing was applied. Systematic studies were then performed in order to carry out a grid sensitivity study and to predict the numerical uncertainties. For this purpose, three different types of mesh refinements were tested given in Table 2. A mesh refinement factor, was chosen as 2 . In order to estimate uncertainties of numerical model, a verification study should be carry out on grid structure. The Grid Convergence Index (GCI) method can be employ for monotonically convergent cases. In this method, , represent the solution of key parameter at corresponding grid type. If the convergence condition ( ) is between 0 and 1, this means that the solution is monotonically convergent [61]. In the present study, the convergence condition is calculated as 0.7 from Table 2. The discretization uncertainty of numerical model is evaluated as % 2.28. As a result, the fine mesh configuration is selected which shows the best compromise between the accuracy and run-time. Detailed information of GCI procedure can be found in [64].
Subsequently, double-body, fully turbulent viscous flow computations were performed for all of the variant hull forms. The usage of double-body model is preferred, since the present study particularly aims at obtaining optimal forms for minimum viscous resistance. Thus, the investigation of wave resistance is excluded from the meta-modeling and from the optimization procedure to see the net viscous pressure effect of the variant hull forms.
The computation time was about an hour approximately for one case via using computational fluid dynamics software CD-Adapco's Star-CCM + with Intel Xeon i7 2.4GHz CPU with 64 GB of RAM. Three residuals monitored for convergence: continuity, forces acting on the hull and velocity values. The convergence criteria is selected as 0.0001 for all these residuals for all hull variant analysis.

Design Space and Sensitivity Analysis
Firstly, sensitivity analysis was carried out to determine the location of the control points defining the surface patch of the optimization region. We sought for the points of which their small variations end up with greater resistance changes. The main idea in the selection procedure of control points via sensitivity analysis is based on the approach that the search is made to figure out the control points having significant effects on viscous resistance when small geometric variations are applied transverse-wise. Then the location of control points is rearranged in a way that they concentrated around the points having greater changes (sensitivity) in resistance due to small changes in transverse-wise variation.
From the surface view of the hull it can be seen that the two lines (1-3-6 and 2-4-5) are the two sensitive directions along which the geometry is changed rapidly. 6 control points were selected for this purpose, depicted in Figure 3. 12 variant hull forms were then created by allowing ± 10 % changes horizontally (by the width) to the original hull form (Form0). For a producible hull surface, we confine the transverse-wise changes by ± 10 % of the original offsets. The results of the sensitivity analyses of the control points with respect to viscous resistance by RANS computations can be seen in Table 3. As a descriptive example; the hull variant y-1_0.9 means that the first control point moves 10% of the original half-breadth inwards and y-1_1.1 means moving again the first point 10% outwards of the hull. As can be seen from Table 3, the considerable differences in viscous resistance are due to the changes in the 3 rd and 4 th control points. One can see in Figure 4 how effectively the changes in 3 rd and 4 th points result in the flow pattern on the hull surface. Accordingly, the distribution of control points was re-arranged in a way that they are positioned by taking into account of the effectiveness and importance of the 3 rd and 4 th points. The new distribution and location of the control points are given in Figure 5.   From Table 4, initial and revised control point distribution can be compared numerically according to the aft extremity of the hull form.
As a next step; using the 6 control points, a set of design of experiments were created by randomly varying their half-breadths ±10%, which contains 280 number of hull variants nondimensionalized with respect to Form0 given in Figure 6. This figure summarizes the results of the design-of-experiments-study obtained by the computational evaluation of the variant hull forms.

Training of the Metamodel
The training output of Kriging is given in Table 5. From Table 5, it is obvious that increasing the number of training points affects the correlation coefficient greatly. According to the general understanding when the correlation coefficient is greater than the 0.8 the metamodel is assumed to be satisfactory [54]. From the Table 5, 140 training data points almost doubled the correlation coefficient which implies that there should at least be around 170 training input data to receive a satisfactory prediction from the Kriging metamodel in its present form. It is concluded from the table that, metamodel is very sensitive to the additional new data at first which means the surrogate model need more data to represent the relation between the control points variation and the corresponding resistance results to find a satisfactory minimum for the output. As we know from the literature every engineering problem has a saturation point where no more data can increase the approximation accuracy [54]. Thus it is concluded to continue the metamodeling with 240 data points for training.

Optimization and Results
Training the metamodel properly makes the procedure ready to go for the optimization stage. We restrict the transverse-wise variations within ±10% of the original offsets. The optimization process for minimum viscous resistance is performed by genetic algorithm (GA). GA solver, utilized from the Matlab library in the present study, is a commonly used optimization solver based on genetic algorithms, [65,66]. The basic components common to almost all genetic algorithms are; a fitness function for optimization, a population of chromosomes, a selection of which chromosomes will reproduce, crossover to produce next generation of chromosomes and random mutation of chromosomes in new generation [62].
In the genetic algorithm, the evolution begins with a population of randomly generated individuals. In each generation, the fitness of every individual in the population is evaluated; multiple individuals are selected from the current population according to fitness and modified to form a new population. Then, new population is used in the next iteration stage of the algorithm. The algorithm terminates by achieving the maximum number of generations or by exceeding a satisfactory fitness level for the population.
The present optimization study started with a set of data, x1 to x6, which correspond to multidimensional design input of each variant hull form, with observed responses, R (presently the viscous resistance). The lower bound is 0.9 and the upper bound 1.1 for x1 to x6 that we restrict the variation of points ±10% by the width of the hull model. Then the metamodeling stage started to find an expression for a predicted value at a new point (x1 * , … ,x6 * ). The basis function of the Kriging is given in (1). In this equation Y values are the observed data (viscous pressure results) which corresponds to the related set of control points x1 to x6 (half-breadth values of the selected design space) data set points. In the present GA optimization; the population type is taken as double vector, the population size is selected as 50. The other main parameters are the number of generations, taken as 100, and the tournament size, taken as 5, [54].
The optimization results correspond to various levels of training data are given in Table  6. The table particularly shows how the level of the number of training data contributes to the output of the optimization procedure. The percent gains in Table 6 were obtained as pointed out from the optimization procedure which depends on the metamodel. But post-CFD analyses show that the reduction attained by the optimal form in viscous pressure is 5% only instead of 12 % given in Table 6. Note that the optimal control points (half-breadths) in the table are normalized with respect to the original values (Form0) in Table 7. Sectional views of the optimal form and the initial form (Form0) are given in Figure 7 for comparison. The gain in viscous pressure resistance could not exceed 5%, though in Figure 8 the separation problem on the hull surface is improved in considerable amount. Besides the differences in the predictions of the flow solver and the surrogate model on the resistance values, there is a good agreement between the surrogate model and the flow solver on the flow characteristics. In addition, from the viscous resistance results, we can say that the optimal aft is able to make the flow relatively better streamlined as compared to Form0. One can conclude that, on the other hand, the present metamodel should further be studied to make the training process of the algorithm by fewer data.

Conclusions
The present work attempts to introduce an optimization procedure to optimize aft body of a ship for minimum viscous resistance by employing metamodeling of the viscous solver. Accordingly, the approach adopted here points out that the metamodel approach presented here can drastically reduce the computation time and is able to arrive at sound results from viscous resistance point of view.
According to established methodology, it can be concluded that Kriging is a convenient interpolation model for multidimensional hull form optimization problem. It can be seen from the results that the separation problem is reduced in considerable amount and 5 % decrease in viscous pressure resistance is achieved in the end. Although this figure (5%) seems to be relatively small one, the present study demonstrates the capability of optimizing the aft form of ships for minimum viscous resistance by means of metamodeling techniquesemploying Kriging interpolation presently.
On the other hand, it is understood from the results that further improvement in the prediction performance of the present metamodel using Kriging interpolation is necessary such as examining the effect of using different sampling type or adding some descriptive information such as gradient or hessian in the interpolation algorithm. The inclusion of the flow uniformity in the propeller disk in the optimization procedure may be taken as a future work as well as more powerful metamodeling techniques which are able to give higher sensitivity with less data.