Hull Form Optimization Study Based on Multiple Parametric Modiﬁcation Curves and Free Surface Reynolds-Averaged Navier–Stokes (RANS) Solver

Featured Application: Hull-form optimization for a containership based on free-surface RANS solver. Abstract: In this study, the hull form optimization process to minimize resistance of KCS (KRISO containership) at Fn = 0.26 is described. The bow hull form of KCS was modiﬁed by varying such design parameters as sectional area curve (SAC), section shape, bulb breadth, and bulb height using multiple parametric modiﬁcation curves devised by the authors. The resistance performances of modiﬁed hull forms were analysed by the viscous ﬂow Reynolds-Averaged Navier–Stokes (RANS) solver of WAVIS ver.2.2. With a view to saving computational time during iterative analyses in the optimization process, the sinkage and trim were set to the ﬁxed values which had been obtained for the original hull form with free condition. The validity of such constant sinkage/trim was then veriﬁed by conducting analysis for the optimal hull form with free condition. Optimization to minimize the cost function of the total resistance coefﬁcient of model C TM was performed by sequential quadratic programming (SQP), which is one of the gradient-based local optimization methods. Utilization of parallel computing led to the simultaneous calculation of the gradient, thereby speeding up the whole optimization process. At the design speed of 24 knots, the optimal hull yielded C TM reduction by 1.8%, which is extrapolated to 3.1% reduction of effective power P E in full scale.


Introduction
With a view to reducing the greenhouse gas (GHG) emission from newly built ships, the International Maritime Organization (IMO) established the Energy Efficiency Design Index (EEDI) regulation, which mandates stepwise reinforcement of GHG emission allowance every five years. As of 2025, the attained EEDI level should be reduced more than 30% from the baseline level set in 2013. Through this and other regulatory measures, IMO envisages halving the global maritime GHG emission in 2050 compared to that in 2008. Consequently, the development of energy efficient ships using every possible means has become a matter of grave concern for all shipbuilders.
Hull form optimization toward minimum resistance is one of relevant strategies for the development of energy efficient ships. However, the hull form optimization is a very complicated task because it involves the variation of sectional area curve (SAC), waterline shape, longitudinal section shape, profile shape, etc. under the design constraints such as principal particulars and the general arrangement of the ship. The skill and experience of designers significantly affect the hydrodynamic performance of hull form. In addition, the resistance performance evaluation requires time-consuming computational fluid dynamics (CFD) simulation. Thus, the simulation-based design (SBD) has long been chosen as the mesh. In particular, the mesh quality should be maintained for each modified hull form during optimization.
A literature survey indicates that the existing hull form optimization studies have been predominantly based on potential flow solver [1][2][3][4]6,7,[10][11][12][13][15][16][17][18][19][20][21][23][24][25]27]. There are a few notable studies based on free surface RANS solver; Campana et al. [8] minimized the total resistance of model R TM of naval vessel (DTMB5414) based on GA method. Tahara et al. [9] conducted multi-objective optimization study for the same hull to minimize both R TM and motion response amplitude operator (RAO). In order to save computational time, they used coarse grid for the iterative steps of optimization. Duy and Hino [14] presented optimization process of the transom part of KCS to minimize the pressure resistance coefficient C P using SQP method. Miao et al. [26] carried out the optimization of demihull shape and spacing of S60 catamaran to minimize total resistance at Froude numbers Fn = 0.40 and Fn = 0.45.
Recently, the present authors devised an advanced hull form modification method named as multiple parametric modification curves [43]. Compared to the conventional parametric modification function presented in Kim et al. [16] and Kim [42], this method has the advantages in terms of improved fairness and constraint satisfaction under wider range of hull form variation. The present study is aimed at combining such cutting-edge hull form modification method with free-surface RANS solver to present an example of hull form optimization study which is readily applicable in industry. Toward this end, the present study carried out an optimization study for forebody hull form of KCS (KRISO containership). By using the multiple parametric modification curves, SAC, section shapes, and bulb geometry were optimized using SQP method. The cost function was set to be the total resistance coefficient of model C TM , which was evaluated by means of WAVIS ver. 2.2 RANS solver. The C TM was then extrapolated to full scale effective power P E by the 2-D extrapolation of the ITTC 1978 performance prediction method.

Objective Ship
The KCS, which is a 3600 TEU containership developed by KRISO, is selected as the objective ship type in this study. The principal particulars of KCS are listed in Table 1. The design speed V S of KCS is 24.0 knots, which corresponds to the Froude number Fn = 0.26. The body plan and side view of the initial hull are presented in Figure 1.

Hull Form Variation and Design Parameters
Hull form variation is important in the entire hull form SBD process because it affects many significant aspects such as the variety of hull forms, the correlation between design variables and resulting hull form, and the fairness of hull form.
In this study, the hull form variation was carried out by varying the SAC, the sections shape, and the bulb shape of the original hull form of KCS. In order to maintain the length overall (LOA), the bulb length was kept unchanged. The hull form variation can be expressed in terms of parametric modification function f n as follows: Here, the subscript n denotes variation of different kind of design variable; n = 1 for the SAC variation, n = 2 for the section shape variation, n = 3 for the bulb height variation, n = 4 for the bulb breadth variation. As pointed out earlier, the fairness of the new hull form H new is automatically guaranteed as long as the fairness of old hull form H old and the modification function f n are maintained. The initial hull form of KCS and the coordinate system are illustrated in Figure 2. Each parametric modification function affects hull surface coordinates differently; in other words, the SAC variation changes only X coordinates, the section shape variation involves only Y coordinates, and the bulb height and breadth variation correspond to the Z and Y coordinate translation of bulb part. These can be expressed as follows: Each coordinate translation defines the incremental change of coordinate, which can then be regarded as the design parameters to be optimized. The changes in hull form are illustrated in connection with the corresponding design parameters in Figure 3. Table 2 describes the design parameters in detail.  The SAC modification function f 1 (X) was consisted of two B-splines defined by three design parameters: X SAC C , ∆X SAC 0C , and ∆X SAC C1 . The section shape modification function f 2 (X) was a calculated by combining two waterline modification curves and section shape modification function which are defined by six design variables X C _Z Section D , ∆Y 0C _Z Section D , ∆Y C1 _Z Section D , X C _Z Section L , ∆Y 0C _Z Section L , and ∆Y C1 _Z Section L . As pointed out by Park et al. [43], the utilization of two waterline modification functions defined at two different drafts enabled one to achieve the hull form variation with higher degree of freedom in terms of design variables as well as higher accuracy in constraint satisfaction. This is the central idea of the multiple parametric modification curves method. Regarding the bulb height modification function f 3 (X, Z) and breadth modification function f 4 (X, Z), each function was given one design variable.
Among the eleven design variables listed in Table 2, the values of ∆X SAC 0C for the SAC variation and ∆Y 0C _Z Section L for section shape variation were adjusted to maintain the displacement of the modified hull unchanged. Thus, the number of independent design variables used in this study was nine. The change in the displacement was kept below 0.015% of the original hull displacement.

Governing Equations
The CFD solver employed in this study was WAVIS ver. 2.2 [44,45] which was developed in KRISO. The governing equations for 3-D unsteady incompressible viscous turbulent flow are the Reynolds-Averaged Navier-Stokes (RANS) equations, which are expressed in Cartesian coordinates as follows: Here, U i = (U, V, W) is the average velocity component in the x, y, z direction. u i = (u , v , w ) is the fluctuating velocity component. −u i u j is the Reynolds stress defined as k is the turbulent kinetic energy, ν t is the eddy viscosity calculated by Realizable k − turbulence model. The WAVIS ver. 2.2 is the FVM (Finite Volume Method) based RANS solver which is widespread among Korean shipyards. Discretization of convective term and diffusion term were carried out using the 3rd order MUSCL method and the central difference, respectively. The Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm was applied to the velocity-pressure coupling. A level-set method was employed to compute the evolution of free surface.

Computational Domain and Boundary Conditions
The WAVIS solver is capable of automatic grid generation based on hull surface mesh and computational domain, thereby saving efforts required for grid generation. Details on the automatic grid generation can be found in Kim et al. [44,45]. Figure 4 displays the computational domain and the boundary conditions. The computational domain was defined for −1.5 ≤ x/LBP ≤ 1.5, 0 ≤ y/LBP ≤ 1.0 and −1.0 ≤ z/LBP ≤ 0.023. Undisturbed free surface was set at z = 0. The FP (front perpendicular) and AP (after perpendicular) of the model were located at x/LBP = 0.5 and x/LBP = −0.5, respectively. y = 0 corresponded to the centerplane around which the symmetry of the hull form allowed half domain simulation. No-slip condition was applied on the hull surface with the standard wall function being utilized to compute turbulence quantities at y + = 50.

Validation and Verification of CFD Results
Grid convergence test was performed for three-grid systems (coarse, medium, fine) generated for the original KCS hull form subject to 2-DOF motion. As had been suggested by Kim et al. [44], the finer grid was prepared by refining hull surface mesh of the coarser grid by the ratio of √ 2. The numerical simulations were conducted for the design speed of KCS at Fn = 0.26, which corresponded to model Reynolds number Rn = 1.40 × 10 7 .
The numerical results from the grid convergence test are listed in Table 3. The row denoted as "Exp" corresponds to the model test results provided in the Tokyo 2015 Workshop on CFD in Ship Hydrodynamics [46]. In Table 3, it appears that the coarse grid yielded error in C TM of 3.1%. For medium and fine grid, the errors in the C TM were found to be 2.2% and 2.3%, respectively. The differences in sinkage in model scale between the experimental result and numerical results are about 2 mm for coarse and fine grids and 1.5 mm for medium grid. In case of the trim angle, the largest deviation from experimental result was obtained for the fine grid; the differences in the drafts at bow and stern are found to be 1 mm in model scale. In summary, all three grids showed reasonable agreements with experimental results both in terms of the total resistance coefficient C TM and the running attitude. In this study, 24 CPU cores (Xeon ® Gold 6252 CPU @ 2.10GHz, Intel Corporation) in the cluster computing system were allocated for the numerical simulations. In order to secure convergence, the simulations were time marched until 100 s, which corresponded to 20,000 time steps (∆t = 0.005 s). As listed in Table 3, the computational times for each grid were three, five, and ten hours. In consideration of the accuracy of the results and the computational time, the medium grid was selected for the simulations hereinafter.

CFD Validation with Respect to Running Attitude
The computational time for the medium grid selected was five hours, which was not considered to be short enough to carry out iterative simulations for the SQP optimization. Consequently, further reduction in computational time without sacrificing the accuracy was necessary. In order to shorten the initial converging period before the running attitude is settled from the force and moment balance, the simulation started with the running attitude being fixed to the values of original hull. This assumes that the running attitude would be essentially independent on the hull form variation. Table 4 and Figure 5 compare the results obtained from various running attitude conditions for the original KCS hull. The first row denoted as "Free" represents the simulation with the running attitude being freely set by the force and moment balance. The second row ("Fixed") means that the running attitude is fixed at the static equilibrium, i.e., zero sinkage and trim throughout the simulation. For the third row ("Fixed with static sinkage & trim"), the sinkage and trim values are first fixed to the values of "Free" condition and then simulation is conducted all along. In Table 4, it is observed that the fixed with static sinkage & trim case yielded very close value of C TM , which is only 0.2% different from that of free condition in a considerably shorter computational time of two hours. In Figure 5, "Free" and "Fixed with static sinkage & trim" conditions give rise to almost identical wave profiles, which support the close match of the resistance coefficient. On the other hand, the "Fixed" condition fails to give a close match to the "Free" solution, which differs by 3.2%. This is because the running attitude was kept fixed to the incorrect values all along the simulation. The convergence history in Figure 5 indicates that the "Fixed with static sinkage & trim" case converged much faster at around 6000~7000 iterations. This is in support of the significant saving of computational time, approximately 60%, enabled by this treatment.  Although the present running attitude treatment is shown to be effective for the original KCS hull, this may not necessarily be the case for other modified hulls. The comparison between the results using two different running attitude conditions for the optimal hull is listed in Table 5. As a matter of fact, these results are taken from the final result of optimal hull, which will be given in Section 6.2. Even though the underlying assumption of independence of running attitude on hull form variation may be seemingly unreasonable, the excellent agreement between the results in Table 5 is just remarkable. The result "Fixed with static sinkage & trim", which was calculated with "incorrect" running attitude of original hull, yielded the C TM value which was only 0.2% different from that obtained using correct force and moment balance. The close match between sinkage and trim for two conditions is indeed supportive of the underlying assumption. In summary, the present treatment of predetermined running attitude is shown to provide a useful means to speed up the iterative hull form optimization process involving the time-consuming RANS simulations.

SQP Algorithm
The mathematical formulation of the optimization problem is expressed as subject to the equality and inequality constraints h j (x) = 0, j = 1, · · · , p (10) where f i (x) is the objective function, K is the number of objective functions, p is the number of equality constraints, q is the number of inequality constraints, and x = (x 1 , x 2 , · · · , x N ) ⊆ S is a solution or design variable. The search space S is defined as an N-dimensional rectangle in R N (domain of variables defined by their lower and upper bounds) The constraints define the feasible area, which means that if the design variables vector x is in agreement with both constraints h j (x) and g j (x), it belongs in the feasible area. In this study, the sequential quadratic programming (SQP) method was used to obtain the numerical solution of constrained nonlinear optimization problem. ∇ f (x) is the gradient of f (x) and it can be expressed as follows: When the gradient cannot be evaluated directly with respect to the design variable x i , the gradient of objective function f can be evaluated by finite differencing as follows: in which i is the step size of design variable x i . This method, based on the iterative formulation and solution of quadratic programming subproblems, obtains subproblems using a quadratic approximation of the Lagrangian and by linearizing the constraints. The following quadratic programming problem is solved to obtain the search direction vector d: H is an approximation to the Hessian matrix of the Lagrange function. When the design variables x k i at k-th iteration step are determined, the design variables x k+1 i at k + 1-th iteration step are evaluated as follows: Here, α k represents the step length at k-th iteration step. During the iterations for optimization, the principal particulars such as LBP, B, T were kept unchanged, and the inequality constraints were set in terms of maximum increment in design variables. As pointed out earlier, the displacement of ship was allowed to vary within 0.015%. Table 6 lists the relevant parameters of hull form optimization such as the initial value, the range, and the step size for calculation of gradient for nine design variables introduced in Section 3. For instance, the initial value of X SAC C for the SAC variation was set at 16.297 ST, where the prismatic coefficient C P becomes equal to 0.5. Here, all the quantities related to the design variable defined in the X direction were expressed in terms of station (ST), which is 1/20 of LBP. Similarly, the quantities defined in the Y and Z direction are normalized with respect to the breadth B and the draft T, respectively. Regarding the section shape variation, the reference Z positions were placed at Z D = T (at design draft) and Z L = T/3. The initial values of the design variables X C _Z Section D and X C _Z Section L were set at X = 17.215 ST and X = 15.852 ST, respectively. These positions were obtained at the X coordinates where the waterlines at Z D and Z L cross the quarter-breadth position Y = B/4. Figure 6 illustrates how these initial values were determined in terms of the SAC and the waterline shapes of the original KCS hull form. Table 6. Initial value and range of design variables for SQP optimization.

Initial Value Range
Step Size, As addressed above, the most significant concern in the implementation of SBD hull form optimization is the computational time required for the repeated RANS simulations. Here, a practical consideration on this is given. Since the SQP is a gradient-based optimization algorithm, the finite differencing to calculate gradient demands (2N + 1) times RANS simulations for one iteration with N design variables. If the iteration number K is needed to obtain the converged solution to the SQP optimization, then at least (2N + 1) × K simulation times will be required. If we assume the computational time for each simulation to be H hours, then the total computation time will become (2N + 1) × K × H hours. As discussed in Section 4, the value of H can be decreased from 5 to 2 by adopting "Fixed with static sinkage & trim" running attitude.
The total computational time could be further reduced by means of parallel computing. If two simulations required for the calculation of the gradient can be performed simultaneously, then the computational time will be reduced to (1 + 1) × K × H hours. However, the number of CPU cores allocated should be increased by 2N (N = 9) times in return for speeding up the process. For the present study, this means that in total 2N × 24 = 2 × 9 × 24 = 432 cores should be used for optimization process. Table 7 shows how the total computational time can be shortened by adopting predetermined running attitude and parallel computing for the case of K = 10. It appears that the total time will be greatly shortened from 950 h to 40 h. In summary, the hull form optimization based on RANS simulation has become a practically feasible process at the expense of a reasonable number of computational resources. It is worth mentioning that the total computational time in Table 7 is an exemplary estimate and the actual time can vary depending on the conditions of CPUs and the interface between them.

Hull Form Optimization Results
Based on the methodology described in the previous sections, the hull form optimization of KCS was conducted to minimize the total resistance coefficient of model C TM at design draft and design speed Fn = 0.260. During the SQP optimization process, the running attitude was fixed to the values for the original hull. By using 432 CPU cores, the gradients with respect to nine design variables were calculated simultaneously. After 10 iteration steps, the C TM was converged to minimum value and total computational time was 44 h including the original hull form calculation. Table 8 lists the convergence history of the C TM , the pressure resistance coefficient C P and the frictional resistance coefficient C F . The iteration 0 denotes the original hull form calculation and the iteration 10 corresponds to optimal hull form, which yielded 1.8% reduction in C TM . The subscript 0 denotes the original value, so C TM /C TM0 represents the ratio of C TM between modified hull and original hull. The time histories of C TM /C TM0 and C P /C P0 are plotted with respect to the iteration number in Figure 7. As seen, the reduction in C TM is closely associated with the reduction in C P and occurs mainly until the 4th iteration. The amount of reduction in C TM in the initial phase until the 4th iteration appears to be about 70% of total improvement. Table 8. Convergence history of resistance coefficient with respect to the iteration number.

Iteration
C P × 10 3 C P /C P0 (%) C FM × 10 3 C TM × 10 3 C TM /C TM0 (%)  Body plans of three hull forms at iteration step 0 (original), 4, and 10 (optimal) during the optimization process are presented in Figure 8. The comparison between original and iteration 4 in Figure 8a demonstrates how the hull form is modified in the initial optimization period. It appears that the volume in the entrance part ahead of 17 ST increased whilst the volume in the shoulder part after 17 ST decreased. In order to keep the displacement unchanged, the amount of volume increase in some part should be balanced by the same amount of volume increase in the other part of the hull. A closer inspection reveals that the volume increase in the entrance part is distributed mainly around the bulbous bow near Z L rather than at the design draft Z D . Figure 8b represents the hull form modification during the latter period of optimization. It is observed that the trend of initial hull form modification is still preserved except that the volume near Z D in the entrance part has been slightly diminished.  Table 9 compares the hydrostatics of the optimal hull form with that of original hull form. Firstly, the displacement and the block coefficient C B are essentially unchanged. This is because the displacement of the modified hull forms was allowed to vary within 0.015% of the original hull displacement, which is equivalent to |∆C B | ≤ 0.0001. The wetted surface area has been changed only slightly by 0.021%. The longitudinal center of buoyancy LCB, however, is increased noticeably, which implies that more volume is displaced to the frontal part. This is consistent with the trend of hull form variation described earlier, i.e., the volume increase at the bulbous bow in the entrance part ahead of 17 ST and the volume decrease in the shoulder part.  Figure 9 presents the SAC, the body plan, the waterplane, and the profile of optimal hull form in comparison with those of original hull form. The forward movement of volume for the optimal hull is also discernible from the SACs as well as the body plans. On the other hand, the waterplane and the profile shapes scarcely exhibit the difference between original hull and optimal hull. This indicates that the hull form modification was not significant in terms of volume distribution at the design draft Z D and the bulb height. It is worth mentioning that the hydrodynamic performance of the optimal hull described in the previous section is not the final one. This is because the running attitude of the optimal hull was assumed to be the same as that of the original hull. Therefore, it needs to be recalculated under the free condition with the correct force and moment balance being accounted for. Table 10 contains such results for the optimal hull and further compares with the original hull. Although the C TM value in Table 10 is slightly different from the value in Table 8, the amount of C TM reduction of 1.8% remains unchanged. −0.0020 −0.174 Figure 10 compares the wave profiles on the hull between the original hull form and the optimal hull form. It is noteworthy that the optimal hull form yielded generally milder wave profile except for the slightly higher 1st peak. This suggests the reduction in the wavemaking resistance for the optimal hull form, which is also consistent with the 6.5% reduction in the pressure resistance coefficient C P .

Optimal Hull
In order to facilitate the comparison between the wave generation in a wider space, Figure 11a plots the contours of wave elevation for two hull forms. It is observed that the area of wave trough region, colored blue, near the shoulder is diminished for the optimal hull form. This is also consistent with the milder wave profile for the optimal hull form observed in Figure 10. Figure 11b provides a comparison in different perspectives-hull pressure distribution. Except at the confined region near the bulbous bow, the optimal hull form exhibits relatively wider spacing between pressure contours in the entire forebody than the original hull form. This also leads to the contraction of pressure minima near the midship. These are closely associated with the diminished shoulder wave observed in Figures 10 and 11a. In summary, the free surface RANS simulation results indicate that the wavemaking performance of the optimal hull has been notably improved over the original hull form.

Extrapolation to Full Scale Performance
The full-scale hydrodynamic performance of the optimal hull form can be predicted by the 2-D extrapolation of ITTC 1978 performance prediction method. Firstly, the residuary resistance coefficient C R needs to be calculated from the C TM obtained by the RANS solver and the frictional resistance coefficient in model scale using C R = C TM − C FM and C FM = 0.075 (log 10 Rn M −2) 2 . Then the extrapolation to full scale is carried out as follows: Here, the model-ship correlation allowance C A = −0.0002 was used. The air resistance coefficient is given as C AA = 0.001 A TS S S using the transverse area over waterline A TS = 817 m 2 . S S and S BK = 123 m 2 denote the wetted surface area and the bilge keel area of full scale, respectively. Table 11 summarizes the relevant quantities involved in the full-scale extrapolation procedure for the original hull and the optimal hull. As seen, the present hull form optimization study achieved a 3.1% reduction in the effective power at the design speed of 24.0 knots (Fn = 0.260). Table 11. Comparison of full-scale performance prediction between original and optimal hull forms at the design speed V S = 24.0 knots (Fn = 0.260 ).

Hull Form
Wetted Surface Area (m 3 ) C TM × 10 3 C FM × 10 3 C R × 10 3 C FS × 10 3 C TS × 10 3 R TS (kN) P E (kW) Finally, Figure 12 displays the predicted values of P E of the optimal hull as a function of speed ranging from 20 to 26 knots in comparison to those of the original hull. At each speed, the percentage reduction in P E is marked, which increases with increasing speed. Since the present hull form optimization is aimed at the reduction of wavemaking resistance, the performance improvement becomes larger at higher speeds when the wavemaking resistance becomes more dominant. In conclusion, the optimal hull form yielded performance improvement over the original hull form in a wider speed range, including the design speed.

Conclusions
The present study provides the results of hull form optimization process for KCS (KRISO Containership) based on the multiple parametric modification curves and free surface RANS solver. Optimization was performed using SQP algorithm to minimize the resistance at design speed Fn = 0.260. Hull form was modified by varying in total of nine design variables. By adopting predetermined running attitude and parallel computing with 432 CPU cores, it was possible to shorten the overall computational time down to 44 h. At the design speed of 24 knots, the optimal hull yielded C TM reduction by 1.8%, which was extrapolated to 3.1% reduction of effective power P E in full scale. The present study is noteworthy in that it combines effective hull form variation technique with accurate performance prediction by RANS solver, thereby implementing the simulation-based design (SBD) for hull form in a practically meaningful way.