1   Introduction

As temperatures in Arctic regions continue to rise under the effects of global warming, the coverage of the Arctic ice cap is decreasing. This is having an important impact on the safe navigation of Arctic ships. The extended opening period of Arctic shipping routes presents an opportunity for commercial shipping. Icebreakers are important vessels in the navigation of polar routes. By optimizing the design of icebreakers, it is possible to improve navigation performance in ice areas, reduce energy consumption and fuel costs, and promote environmental protection. Hence, there has been extensive research on the optimization of ship resistance in open water and polar navigation.

The shape of the stem is an important factor affecting the icebreaking load. Kim and Lee (2010) proposed a hull form for an icebreaking vessel that incorporated a low bulb bow to optimize the capability of the vessel for the operational conditions of the route, the thickness of flat ice, and the overlap degree of the ice. A study of the resistance performance of icebreaking cargo vessels based on the fluid–structure interaction method indicated that a large waterline angle is beneficial for reducing the navigation resistance in high-density floating-ice areas (Kim et al., 2014). The bow shape has a significant influence on the icebreaking pattern. The icebreaking load exhibits cyclical behavior, with crack propagation in one cycle occurring at the same time as crack initiation in the next cycle (Zhang et al., 2021b). Anzai et al. (2022) used a multi-objective genetic algorithm to optimize the performance of icebreakers in terms of reducing the ice resistance under level ice and the wave resistance in open waters, and obtained a ship form with lower resistance and a large flare bow. The icebreaking force at the waterline position of the bow is second only to the total resistance load in ships with non-typical icebreaking bow shapes, such as high stems and small waterline angles (Hisette and Myland, 2020). In ships with non-typical icebreaking bow shapes, the bow induces approximately 80% of the total resistance, with the stern accounting for the remaining 20%. The contribution of the bow increases with increasing ice thickness (Myland and Hisette, 2021; Myland et al., 2021). Backward-moving crushed ice is known to affect the performance of the propeller (Zhou et al., 2024). Sun and Huang (2023) reported that the icebreaking resistance in the stem area accounts for more than 60% of the total resistance when navigating through ice. Based on analysis of the ice resistance of three different bow forms, Ren et al. (2023) recommended that their characteristics be combined to produce an effective icebreaking design. Variations in the contact length on the frame ice have been observed to cause a non-monotonic upward trend in the scale parameter of the stress amplitude distribution (Han et al., 2024b).

Numerical methods and empirical formulas can be used to predict the icebreaking load in optimization problems. Computational fluid dynamics combined with the discrete element method (CFD-DEM) has been extensively used to simulate the ship–ice–fluid interaction of ships in recent years (Chen et al., 2024). An optimization algorithm combining CFD-DEM with Sobol sampling and the T-Search method was used to optimize the performance of propellers in icebreaking regions, effectively reducing the combined self-propelled power by 9.71% (Liu et al., 2024). To reduce the computation time and resources during the optimization process, the ice resistance is often estimated using empirical formulas. For instance, the Lindqvist formula considers multiple variables including the resistance in flat ice areas and the ship width, depth, and stem angle. Wang et al. (2021) proposed optimal ranges for the bow characteristic parameters based on the Lindqvist formula, and optimized the bow line shape considering the resistance performance of open water areas via the Sobol algorithm. To optimize the bow hull form, Lu et al. (2022) used the Rankine source method and the Lindqvist formula to calculate the water and ice resistance, respectively. The overall weight-averaged resistance based on the open-water and ice-covered-water sections was reduced by 5.42% using version 2 of the non-dominated sorting genetic algorithm (NSGA-Ⅱ). Wang et al. (2021) used Sobol's method to perform a global sensitivity analysis of the icebreaking performance, and found that there is greater sensitivity to the stem angle and waterline angle in the Lindqvist and Keinonen models. The ice resistance on the ship hull during navigation in level ice using the Lindqvist formula is consistent with that given by DEM simulations (Hu et al., 2021). Therefore, the Lindqvist model is better for rapid predictions of the ice load on the bow of a ship. Data-driven alternatives have been explored as a means of reducing the computational costs and improving the optimization efficiency (Walker et al., 2024). Wang and Shan (2006) explored cutting-edge metamodel-based techniques, with special attention on their ability to facilitate design optimization. They considered multiple facets of metamodeling, such as model approximation, the exploration of design spaces, and the formulation of problems, and tackled a wide array of optimization challenges. Aerospace design evaluation relies heavily on computationally demanding simulations, driving the need to leverage efficient surrogate-based methods for optimization (Forrester and Keane, 2009). Ang et al. (2015) presented an optimization study of a bulbous bow using the NAPA software, and discussed potential advances in parallel computing and machine intelligence. The use of simulation models in engineering design is computationally expensive. Surrogate-based optimization, such as shape-preserving response prediction, offers a promising solution (Leifsson and Koziel, 2016). Optimizing simulation-based or data-driven systems is challenging, whereas surrogate models are efficient. Kim and Boukouvala (2020) compared five subset selection regression techniques with nonlinear interpolating functions, and reported that subset selection-based regression performs well in low-dimensional tasks, whereas interpolation excels for higher numbers of dimensions. In any optimization process, it is important to consider uncertainties and environmental constraints (Han et al., 2024a).

Ship-type optimization is typically a multi-objective problem. Kim et al. (2006) conducted a series of model experiments to assess the navigational capabilities of different bow shapes in various ice conditions, including open water, level ice, and brash ice. Their innovative design of the bulbous bow significantly improved the icebreaking performance over that of traditional bulbous bows designed for open-water use, while having minimal effect on resistance in open water. The final design of a small-waterplane-area twin hull (SWATH) ship with a high seakeeping function was selected according to the median compromise of the Pareto front, using the heave motion and resistance as objectives (Renaud et al., 2022). To tackle the low efficiency of single-objective optimization, Wang et al. (2023) developed a sampling-based optimization system for the total resistance. Wei et al. (2024) employed the NSGA-Ⅱ algorithm to search for the maximum pseudo-expected improvement matrix value and identify the optimal resistance and nominal wake performance. The current approach leverages hybrid CFD–data-driven surrogate models to optimize key performance indicators such as drag or lift. A high degree of customization complicates the reuse of data-driven surrogate models in optimization with different parametrization schema (Walker et al., 2024). One study integrated the Rankine source method to calculate the wave resistance coefficient with the NSGA to determine the optimal Series 60 ship hull form (Wang et al., 2022). Zhang et al. (2024) used T-spline technology to represent ship hulls, enabling a unified parameter domain with a reduced number of control points to effectively depict the entire hull surface, and applied the multi-island genetic algorithm to find the stem shape with the minimum wave-making resistance calculated by the Rankine source method. Myland and Ehlers (2016) evaluated the contribution of the icebreaking force to the resistance in ice for different bow shapes and waterline angles. Lu et al. (2022) optimized a polar ship by varying the ship speed and the hull form angle of the bow, while Kondratenko et al. (2023) considered the freeboard, stem angle, and block coefficient in their optimization of an Arctic ship. Fasse et al. (2024) performed a surrogate model-based experimental optimization of a propeller with parameterized pitch laws, combining multi-objective optimization with a surrogate construction of Gaussian processes (Fasse et al., 2024). Kim et al. (2024) used CFD to identify potential hull forms with characteristic curves and predicted their performance using a deep neural network.

In addition to optimizing the ice resistance, the total resistance of the ship model and the effective power of the real ship in open water can be optimized to improve the icebreaking performance (Yu et al., 2010; Wang et al., 2018). The energy efficiency design index can be used as the objective of the optimization function and the main dimensions of the ship can be designated as design variables (Su et al., 2018; Tang et al., 2020). Duan et al. (2017) used multi-objective optimization to minimize the propulsion power of a polar ship in the ice zone for three different bow forms. Chen et al. (2025) evaluated the changes in fuel consumption produced by optimizing the total resistance. This is a crucial indicator of a vessel's operational efficiency, and is determined as a function of the main engine's power. For vessels operating in icy conditions, it is crucial to withstand ice–structure interaction loads, although this requirement may reduce a ship's economic feasibility. Wang (2013) conducted scantling, weight, and cost optimization for a platform supply vessel of polar class 2; MATLAB was used to facilitate weight optimization, and the RulesCalc tool and finite element modeling were used to verify the design compliance and strength, respectively. The significant increase in structural weight required to enhance durability can have a negative impact on costs. Balancing these factors is essential for effective and cost-efficient vessel design (Pedersen et al., 2014). Ruiz-Capel et al. (2023) designed the hull structure for the bow region of an ice-class ship according to the Finnish–Swedish Ice Class Rules, and compared the ship's bow weight obtained through direct calculations with that derived from current regulations. Studying the navigation area of polar ships helps optimize their resistance. The historic ice data and ship routes are considered in optimization studies of Lu et al. (2022) and Chen et al. (2025).

This study investigates the effect of the stem shape on the icebreaking load of an icebreaker using a fully parameterized modeling method and optimization scheme. The optimized resistance performance for an icebreaker is thus obtained. The stem of the icebreaker is modeled in terms of the waterline angle and buttock angle at the upper ice waterline considering the ship's position and weight factors. The cyclic process of contact compression bending failure is adopted to calculate the icebreaking loads and the secondary fracture of sea ice is considered (this is a common ice destruction mode in the interaction between icebreakers and flat ice) (Huang et al., 2018; Huang et al., 2016). The momentum generated during the secondary fracture process verifies the impact of this phenomenon on the speed of an icebreaker. The proposed optimization framework determines the appropriate design variables to produce the associated Pareto front through a trade-off among the multi-objective solutions. A differential evolution (DE) algorithm is then used to obtain the Pareto solution. The hull form optimization results are presented and discussed. The initial and optimized hulls are compared, allowing the optimization framework and its effectiveness to be verified.

2   Case study description

The optimization flow and the method of icebreaker design are introduced in section 2.1. The advantages of the various existing optimization algorithms are discussed through a literature review. The method of calculating the resistance is presented in section 2.2, including the resistance in open water and level ice. Fully parameterized modeling is used to model the stem of the icebreaker. The parameterization and deformation methods are introduced in section 2.3. To ensure that the optimized icebreaker can break level ice, the average pressure in the ice zone of the bow is used as a constraint.

2.1   Optimization platform

This article focuses on the optimization of China's icebreaker Xuelong 2, which meets the PC3 level of the IACS polar ice zone regulations. The waterline length L is 116 m, the ship width B is 22 m, and the draft D is 7.8 m. The buttock angle is 20° and the waterline incident angle is 40°. Under an ice thickness of 1.5 m, the speed of breaking level ice is 1.5 m/s. The present optimization problem for a trimaran hull form can be expressed as follows:

where f (X) and g (X) are objective functions relating to the navigation resistance of icebreakers and the average pressure in the ice zone, respectively. S ⊆ R^N is the set of feasible solutions, with the design constraints limiting the feasible design space. The N-dimensional vector of design variables X corresponds to the geometry reconstruction in the optimization process.

The computer-assisted design (CAD) optimization process is illustrated in Figure 1. The optimization loop consists of three components: CAD parametrization, a resistance solver, and an optimizer. Geometry reconstruction is conducted through a fully parameterized modeling approach based on CAESES, and the new geometry is introduced to the ice resistance solver and nonlinear maneuvering model. Finally, new sample points (i.e., new geometries) and corresponding design variables are generated by the DE algorithm and introduced to the CAD environment.

Figure  1  Flowchart of CAD-based optimization strategy

Download: Full-Size Img

For many problems, DE-based algorithms achieve superior performance over the corresponding genetic algorithms (Tušar and Filipič, 2007). An initialization method based on opposition learning has been proposed for generating the initial population of candidate solutions (Rahnamayan et al., 2007). DE has demonstrated its efficacy in resolving intricate challenges, encompassing multi-objective problem-solving, clustering weight determination, and neural network training (Plagianakos et al., 2008). In comparing the performance of genetic algorithms and DE for the optimal design of a radial flux permanent magnet generator, experimental results show that DE effectively improves the convergence speed and solution quality (Lilla et al., 2013). A DE algorithm has been developed to enhance search quality while mitigating premature convergence and stagnation issues (Ali et al., 2017), and a hybrid technique using genetic algorithms and DE has been proposed to prevent premature convergence by improving the seed initialization quality of the K-means algorithm (Mustafi and Sahoo, 2019). After selecting 30 benchmark functions to measure the performance of genetic algorithms, ant colony optimization, DE, artificial bee colony optimization, and other methods, DE was found to be the best or equal best algorithm for 24 of the functions. DE performs very well on multimodal functions, achieving the optimal results for 11 out of 12 such functions (Ab Wahab et al., 2015). DE has high practicality and versatility, providing a powerful tool for complex optimization problems without any assumptions. In many cases, DE can produce better and more stable solutions than genetic algorithms (Ahmad et al., 2022). DE only uses vector operations and a random number generator, allowing it match the speed of simpler but limited genetic algorithms. Unlike genetic algorithms, DE adopts a greedy strategy in selecting the next-generation population. The parent individuals compare their fitness values with those of their offspring, and individuals with higher fitness form the next generation. In terms of crossover, DE first generates mutation vectors based on the parent generation and then performs crossover operations between the parent individuals and their mutation vectors. DE also eliminates the encoding operation required in genetic algorithms. This gives DE the advantages of fast convergence speed, few controlled parameters, and easy implementation. The method of automatic optimization iteration connects the optimization algorithm with automatic modeling. This approach outputs and receives models, and simulates them comprehensively. In addition to the iterative process of the DE algorithm, it is necessary to write instruction files for parameterized modeling based on individual information, process the simulation results to obtain the objective function, and handle individuals that do not satisfy the constraints.

The process of parametric modeling in project files and defining parameters as design variables can be represented by modifiable values in instruction files. The waterline points and auxiliary waterline points of the icebreaker are output in the form of two model files, named line 1.sat and line 2. sat. Line 1.sat includes information such as the area and volume of the stem. The. sat model file is a text schema of the geometry kernel, which can store model information such as points, lines, faces, and volumes.

2.2   Navigating resistance of icebreakers

The hydrodynamic loads are estimated by the nonlinear maneuvering model, in which the hydrodynamic loads are functions of the longitudinal, lateral, and yaw velocities. The hydrodynamic coefficients are estimated by the block coefficient, draft, and beam (Zhou et al., 1983). Icebreakers sailing in level ice are subject to the icebreaking load generated by the compression and destruction of the ice by the hull, as well as the immersion resistance generated by the relative sliding of broken ice along the hull and immersion in water after it falls off the level ice.

The icebreaking loads are calculated according to the cyclic process of contact compression bending failure. We assume the pure bending of sea ice, and neglect the shear effect between the icebreaker and the sea ice boundary. The idealized continuous icebreaking process is as follows: first, the hull contacts the sea ice and exerts a squeezing effect. As the icebreaker advances, the force exerted by the hull on the sea ice gradually increases. When it reaches the bending limit, the sea ice bends and breaks, and wedge-shaped broken ice falls off from the intact and level ice. To ensure that the numerical value of the icebreaking force is close to the actual period of the ice load, the secondary fracture of sea ice must be considered. The icebreaking force is expressed as a piecewise function, and the residual ice force after the ultimate icebreaking force load at the initial fracture is calculated based on the contact surface inclination angle and ship speed at discrete points of the waterline. The ultimate icebreaking force during secondary fracture is calculated as the difference between the loading and unloading icebreaking forces in a ratio related to speed and angle (Gao et al., 2019).

The immersion resistance is calculated by the Lindqvist formula. The total immersion resistance R_s is

where δ_p is the difference between the density of pure water and that of sea water, g is the gravitational acceleration, h is the thickness of the level ice, μ is the frictional coefficient, θ is the buttock angle, α is the waterplane angle, and γ is the flare angle. The relationship between the speed and immersion resistance is

where v is the velocity component and V is the ship speed (Gao et al., 2019).

The proposed approach is verified through its application to the Swedish Tor Viking Ⅱ icebreaker. When the ship is sailing at a speed of 5 m/s in level ice of thickness 0.6 m, the average icebreaking force calculated by the proposed method is 0.43 MN. The result calculated by the Lindqvist method is 0.404 MN, a relative error of 6.4%. In the process of turning at a 45° rudder angle in level ice of thickness 0.5 m, the maximum turning diameter of the full-scale test is 495 m and the simulation result is 463 m, a relative error of 6.5%.

2.3   Fully parameterized modeling of ship stem

The fully parameterized modeling approach has been extended to enable the optimization of complex ship designs. The parametric design of a twin stern fin has been carried out using the full parametric modeling method, with the optimal design of the twin stern fin ship completed using CFD simulation technology combined with the NSGA-Ⅱ optimization algorithm (Zhang et al., 2021a). The key parameters affecting the hydrodynamic performance of a double tail fin have been explored, whereby the optimal parameters were determined to establish a fully parametric model of the tail (Xu et al., 2023). A fully parametric SWATH model has been established to enable automatic geometric transformations, with the Sobol algorithm used for space sampling and sensitivity analysis conducted on the SWATH resistance parameters, along with a correlation study between design variables and resistance (Guan et al., 2021). By integrating NSGA-Ⅱ, successful automatic optimization of SWATH was achieved, allowing the hull resistance to be minimized across various speeds (Guan et al., 2022). The windshield of a ship's bow bulwark should have a smooth, large curved surface designed by Bezier curves. This approach ensures precision and smoothness in the windshield design (Wang et al., 2024).

The parametric modeling software CAESES is used to model the stem of the icebreaker. The stem surface of the icebreaker is generated by stretching a cubic non-uniform rational B-spline (NURBS) curve in the depth direction. Using a cubic NURBS curve and four control points (P₀–P₃), the buttock angle parameter is converted to a direct parameter related to the tangential angle parameter of the curve endpoint and the position parameter of the control points. For the bow waterline in Figure 2, the coordinates of the control points are

Figure  2  Bow waterline and its control points

Download: Full-Size Img

where α and β are the tangential angles of the curve's ends relative to the longitudinal direction, D is the distance between P₀ and P₃, and ρ₁ and ρ₂ are position parameters of the control points.

The parameterized model of the stem below the ice contact height on the prototype is shown in Figure 3. Both line 1 and the scanning curve are expressed using the curve parameters in Figure 2. Line 1 is the mid-longitudinal section line of the bow and line 2 is the first horizontal boundary line. On any horizontal plane at any height, the intersection of the scanning curve and two lines is the endpoint of the NURBS curve, which must satisfy a tangent relationship with the parallel middle body.

Figure  3  Parameterized model of the stem below the ice contact height on the prototype

Download: Full-Size Img

The expression for a k-order NURBS curve P (τ) is

where ω_i is the weight factor corresponding to the n+1 control vertices d_i, N_{i, k} (k) is the basis function of k-order B-spline curves, n is the degree of the Bezier curve, and τ is the knot partition. The tangential angle of the first end of the scanning curve, the position parameters ρ₁ and ρ₂ of the two control points, and the variations in the weight factors ω₁ and ω₂ in the pattern depth direction are shown in Figure 4.

Figure  4  Changes in curve parameters in the depth direction

Download: Full-Size Img

The buttock angle φ at the upper ice waterline and waterline angle θ are chosen as the optimization parameters. The spline curve in Figure 2 changes significantly when the control points P₁ and P₂ are adjusted, even if the tangential angle between the end and endpoints remains unchanged. To consider the changes in bow profile details, the position parameters ρ₁ and ρ₂ and weight factors ω₁ and ω₂ of P₁ and P₂ on the scanning curve are taken as the bow optimization parameters. To reduce the network identification and translation of parameters at the bottom of the bow, the right end of each curve in Figure 4 changes accordingly, and the left side remains almost unchanged.

2.4   Average pressure in the ice zone of the bow

According to the IACS standard, the bow section of a ship is divided into four sub-regions of equal length based on the length of the waterline. For the bow of a PC3 class icebreaker, the shape factor a_i is

For PC4 and higher-grade icebreakers, x = 1.5 m. The force F_i on the stem is

The aspect ratio AR_i of the load surface is

the line load Q_i is

and the pressure P_i is

where i = 1, …, 4, α is the waterline angle, β′ is the normal frame angle at the upper ice waterline and the quarter width of the bow of the ship, CF_F is the bending failure level factor, CF_C is the squeezing failure level factor, CF_D is the load surface aspect ratio level factor, and D_UI is the displacement. The average pressure in the ice zone is

where the ice zone height is b = Q_i /P_i and the load plate width is w = F_i /Q_i. The four sub-regions of Xuelong 2 (from bow to rear) have average pressures in the ice zone of 8.615, 8.583, 8.521 and 8.337 MPa, respectively.

3   Sensitivity analysis

The optimization design of icebreakers must not only consider their resistance performance in ice areas, but also meet the design requirements of the specifications, which can be used as the objective function for subsequent optimization. However, the correlation and degree of correlation between each objective and the bow parameters will vary. Before the optimization, the influence and importance of the modeling parameters on the resistance and average pressure in the ice zone should be discussed to determine the optimization objectives and parameters.

3.1   Ice resistance

To determine the relationship between the objective function and optimization parameters, a correlation analysis was conducted to examine the resistance of longitudinal navigation icebreakers in level ice at various parameter values. This section focuses on the influence of the position parameters P₁ and P₂ and the weight factors on the icebreaking load. Table 1 presents the initial, maximum, and minimum values of the sampling parameters to ensure a reasonable and feasible bow model. Using the Latin hypercube method, 500 sample points were collected within the parameter range, and two ice thicknesses of 1.2 m and 1.5 m were considered. The dashed and solid lines in Figure 5 show the ranking of the ice resistance in longitudinal navigation with an ice thickness of 1.5 m and 1.2 m, respectively. The two ice thickness conditions produce a similar trend, and the calculation result for the 1.2-m ice thickness give a smaller range of fluctuations. This indicates that the trend of resistance performance under different ice thicknesses is the same. Therefore, when optimizing the icebreaking resistance in flat ice areas, the influence of the ice thickness can be ignored.

Figure  5  Resistance under different ice thicknesses

Download: Full-Size Img

Figure 6 shows the proxy model for a ship speed of 1.5 m/s with an ice thickness of 1.5 m. The polynomial proxy model has a complex phase relationship value of 0.953 01 and takes 1 s to calculate. After adding control point parameters, the overall objective function increases, which indicates that changing the control points leads to an increase in icebreaking force. This is because changing the control point parameters increases the complexity of the curve and produces local shapes that are unfavorable to ice resistance.

Figure  6  Proxy model for a ship speed of 1.5 m/s with ice thickness of 1.5 m

Download: Full-Size Img

The influence of the control point parameters on the icebreaking force is a complex function of various parameters, and an effective proxy model cannot be directly established. Thus, the influence of the control point parameters on the icebreaking resistance is expressed as a variable of the waterplane coefficient. Keeping the bow and waterline angles the same as for the original ship, 100 sample points covering the parameter ranges in Table 1 were sampled using the Latin superelevation method. The relationship between the waterplane coefficient of the stem and the icebreaking force is shown in Figure 7. The icebreaking force has a significant linear relationship with the stem waterplane coefficient, with a linear correlation of − 0.53. The linear relationship is more pronounced when the waterplane coefficient of the bow is close to the minimum or maximum region because the waterline tightens at the extremum and the shape becomes more stable. The shape of the waterline in the middle section of the waterplane coefficient varies greatly, resulting in the most unstable icebreaking force.

Figure  7  Relationship between the waterplane coefficient of the bow and the icebreaking force

Download: Full-Size Img

To identify which parameter has the greatest impact on resistance, the correlation between the buttock angle, waterline angle, and resistance in longitudinal navigation was investigated. A multi-layer full-factor sampling method was adopted because of the limited number of parameters, with a buttock angle range of 16° ‒ 30° and a waterline angle range of 36°‒50°. Each factor was divided into 15 layers, giving a total of 225 sample points. The correlation analysis considered i) a ship speed of 1.5 m/s with an ice thickness of 1.5 m and ii) a ship speed of 2.7 m/s with an ice thickness of 1.2 m. The ice resistance is the same under these two operating conditions. A fitting surrogate model was constructed using polynomial functions. As shown in Figure 8, when the ice thickness is 1.2 m, the icebreaking load decreases with increasing waterline angle and with decreasing buttock angle. The change in immersion resistance exhibits the opposite trend to that of the icebreaking load. The load trend when the ice thickness is 1.5 m is the same as for 1.2 m. The complex correlation coefficients, calculated using the leave-one-out method, are 0.959 47 (1.2 m) and 0.939 36 (1.5 m). The computation time for both ice thicknesses is approximately 1 s. The calculation time for the complex correlation coefficient of the polynomial proxy model is relatively short, indicating a strong correlation between the parameters and the objective function.

Figure  8  Response surfaces of different force component

Download: Full-Size Img

The linear correlations (Corr) between the parameters x_i and resistances y_i at n sampling points were calculated as follows using the average of the parameters x and resistances y:

We considered 15 samples near each of the sample points with a waterline angle of θ = 40° and a buttock angle of φ = 20°. Table 2 presents Corr for an ice thickness of 1.5 m and ship speed of 1.5 m/s, because the ice thickness does not affect the trend between the resistance and parameters. Both the waterline angle and buttock angle exhibit a strong correlation with the ice resistance (absolute values of Corr are greater than 0.6). The correlation with the incidence angle of the waterline is slightly weaker than that with the buttock angle, so changes in the mid-longitudinal section will have a greater impact on the icebreaking process than changes in the waterline.

3.2   Average pressure in the ice zone

The average pressure on the stem in the ice zone was calculated according to IACS specifications. The relationship between the maximum average pressure in the ice zone and the two main parameters of the stem was calculated using the multi-layer full-factor method (each factor is divided into 15 layers, giving a total of 225 sample points), as shown in Figure 9. The average pressure in the ice zone of the stem is approximately linearly related to the two angles. The complex correlation coefficient values obtained by using the leave-one-out method is approximately 1. With a buttock angle of 30° and a waterline angle of 36°, the minimum average pressure at the sample point is 8.482 MPa. With a buttock angle of 16° and a waterline angle of 50°, the maximum average pressure is 8.662 MPa. Fifteen sample points were selected for scenarios with a buttock angle of 20° and a waterline angle of 40°. The linear correlations between the average pressure of the four sub-regions (from bow to rear) and the buttock angle were calculated to be −0.998, −0.998, −0.979 and −0.744, respectively. The linear correlations with the waterline angle were 0.989, − 0.317, − 0.898 and − 0.712, respectively. Therefore, the first sub-region has the highest correlation.

Figure  9  Response surface between the maximum average pressure in the ice zone and the angle parameters

Download: Full-Size Img

The relationship between the maximum average pressure in the ice zone and the buttock angle, waterline angle, and control point parameters was calculated at 500 sample points selected using the Latin superelevation method. The results are shown in Figure 10; the range of parameter changes is listed in Table 1. The control point parameters have a limited influence on the average pressure in the ice zone. The maximum average pressure in the ice zone increases as the waterline angle becomes larger and the buttock angle becomes smaller. At 55% of the sample points, the values are greater than those obtained by only changing the main parameters of the bow, indicating that randomly selecting larger control point parameters may lead to a decrease in the ultimate performance of the bow.

Figure  10  Response surface between the maximum average pressure in the ice zone and the angle parameters, position, and weight factors

Download: Full-Size Img

4   Optimization

According to the results from the proxy model, the influence of various parameters on the performance of the icebreaker can be comprehensively considered. This section first analyzes the objective function and performs multi-objective optimization on the resistance and average pressure. For this discrete optimization problem, a modified differential evolution algorithm with fast convergence is used, and the weighted summation method is applied to decompose the multi-objective optimization problem.

The key to optimizing the performance of icebreakers navigating in level ice is the shape of the stem. The main basis for evaluating the quality of a ship's bow is its longitudinal resistance while navigating in level ice. The correlation between the strength performance of the bow structure and the main parameters of the bow exhibits the opposite trend to that of the resistance. Therefore, the strength and resistance cannot be combined into a single objective, and multi-objective optimization must be employed. We take the buttock angle φ, waterline angle θ, control point position parameters P₁ and P₂, and weight factors ω₁ and ω₂ as the optimization parameters. The constraint condition is that the surface area and volume change of the icebreaker cannot exceed 3%. The main parameters of the stem have a significant impact on the optimization objective function. Reducing the buttock angle and increasing the waterline angle reduces the ice resistance, but has a significant impact on the weight distribution, internal configuration, hydrodynamic performance, and longitudinal stability of the ship. An excessive increase in the buttock angle and a reduction in the waterline angle to reduce the load in the ice zone negatively affects the overall arrangement. The expected optimization result is an icebreaker with better objective function values without significant changes in the main parameters of the bow.

The optimal point on the Pareto solution set appears when there is more than one objective function and there is a contradictory relationship between the objective functions that cannot simultaneously satisfy the optimal conditions. At least one objective function value on the solution set is better than at all other points. The optimization objectives are the longitudinal navigation resistance of the icebreaker in a level ice zone and the average pressure in the ice zone in the ice zone of the stem. The DE algorithm is adopted for the simulation calculations. There is no need to obtain Pareto solutions for the entire parameter range because the desired optimized bow parameters are similar to those of the original ship. Therefore, the weighted sum method is used to transform the multi-objective optimization into a single-objective problem. A one-dimensional approximate relationship is established between the two objective functions and the main parameters of the bow based on the proxy model described in section 3, as shown in Figure 11. The sampling points (φ, θ) represent the bow and waterline angles. To ensure that the two objective functions are of equal magnitude, the resistance at the point based on the original ship parameters (20, 40) is adjusted to be equal to the average pressure.

Figure  11  Relationship between design variables and objective functions after balancing the order of magnitude of the objective functions

Download: Full-Size Img

When the weight coefficients of both objective functions are set to 0.5, the trend of the transformed single-objective function with respect to parameter changes is the same as the total load (see Figure 11). The sum of the two objective functions increases monotonically on (φ, θ), and the optimization result must have the form (min, max). However, a small buttock angle and a large waterline angle do not conform to the overall design of the stem. Hence, it is necessary to further improve the objective function to ensure that the design variables of the expected final optimization result are closer to those of the original ship. Therefore, the average pressure function f_pavg of the load plate is changed to

This equation causes all points to the left of (20, 40) to become a decreasing function, and the growth trend on both sides is the same. As a result, the point (20, 40) becomes close to the minimum value of the objective function. At this point, the Pareto solution forms a convex set. The final optimization result is somewhere around (20, 40), but produces a better icebreaker.

In the optimization process, the method of randomly generating parameters is first used to generate the initial population as the parent generation. After each generation of simulation calculation is complete, the objective function is obtained and offspring are generated for the next generation. The DE algorithm stores the parent generation with the objective function and the offspring without the objective function. After each generation is used to obtain the objective function, the result is compared with that for the parent generation, and the parent generation for the next iteration is determined. The mutation factor of the DE algorithm is 0.6 and the selected population size is 60; the CPU utilization rate of each simulation process is about 2%. The WaitForMultipleObjects function is used to create multiple threads, with 30 threads created for each iteration.

The optimization result of the objective function is shown in Figure 12. The maximum average pressure in the ice zone in the ice zone at the stem of the original ship is 8.615 MN. The resistance value after dimensionless processing is also 8.615. The weight coefficients of both objectives are 0.5, so the objective function of the original ship is 8.615. The points in the figure represent the solution with the smallest objective function value in each generation of the population. The minimum value of the initial population is 8.432. After 49 iterations, the minimum value reaches 8.337, a reduction of approximately 3.2% compared with the original ship. To increase the global search capability of the optimization algorithm, a larger population size can be used, meaning that better individuals are likely to be sampled in the initial population.

Figure  12  Objective function iteration results

Download: Full-Size Img

The multi-objective optimization results of the population, producing a significant decrease in resistance, are shown in Figure 13. The dashed line represents the Pareto front of multi-objective optimization. The final location of the aggregate population is the tangent point between the line perpendicular to the weight coefficient vector and the Pareto front. The average pressure at the stem in the ice zone is approximately 8.61 MPa in the final generation, which is 0.6% lower than for the original ship. The dimensionless ice resistance has decreased significantly, by 4.2% compared with the original ship. The reason is that the maximum average pressure in the ice zone occurs at the forefront of the bow, and is thus closely related to the waterline angle. The navigation resistance in a straight line is related to the shape of the entire waterline, and can be further optimized by adjusting the control points of the waterline spline. Equation (13) also affects the optimization results. To prevent the buttock angle from becoming too small, the main parameters of the ship stem are constrained to prevent an uncontrolled reduction in the average pressure in the ice zone.

Figure  13  Multi-objective optimization results and the Pareto front

Download: Full-Size Img

The changes in the main parameters of the bow during the iteration process are shown in Figure 14. The sample points do not include large angles because of the displacement constraint, although the waterline angle varies from 36° to 50° and the buttock angle ranges from 16° to 30°. The iterative process largely discards the initial population because it does not satisfy the constraint conditions. The buttock angle of the last generation is around 21° and the waterline angle is around 44°. There are unlikely to be any solutions with a buttock angle less than 20°, which means the main parameters of the bow are similar to those of the original ship. This demonstrates the effectiveness of the optimization design and methods used in this study. The clustering of sample points also indicates that the solution is globally optimal.

Figure  14  Optimization trend of main parameters of bow in five generations

Download: Full-Size Img

The optimized ship types are summarized in Table 3. The position parameters and weight factors are both different from those of the original ship. The displacement of the icebreaker has increased by about 1% because of the increase in the main parameters of the bow. However, a stem with a larger displacement is more conducive to depressing the ice layer. The weight factor ω₂ is consistently negative, corresponding to the better resistance performance of the bow waterline in Figure 15. The control point parameters are listed in Table 4. The waterplane coefficient in the figure is 0.646. Although the area enclosed by the two waterlines and coordinate axes is the same, waterline 1 has a uniform curvature, with a certain outward convex curvature in the first half and a smooth transition between the second half and the midship. The first half of waterline 2 tends towards a straight line, while the second half connects quickly with the midship section, resulting in a significant change in curvature. Waterline 1 has a better icebreaking effect at the bow of the ship because the first half of the ice zone has a smaller contact angle with the ice and the flat ice mainly bends and breaks. However, the second half of the ice zone has a larger contact angle with the flat ice and the ice mainly undergoes compression breakage. The convex shape is beneficial for increasing the outward inclination angle of the hull, which enhances the ship's ability to exert downward pressure on the ice layer. A smaller value of ω₂ produces a smoother connection between the rear half of the bow waterline and the parallel body.

Figure  15  Relationship between waterline shape and control points

Download: Full-Size Img

Taking Optimized Ship No. 4 as an example, the ship's straight-line resistance has decreased by 4.5% while ensuring that the average pressure values of the four load plates at the bow of the ship are lower than those of the original ship. The wet surface area of the ship is 3 523 m², an increase of 0.44% compared with the original ship. The horizontal profile lines of the bow on the optimized ship and the original ship are compared in Figure 16. The first half of the optimized waterline is full because of the waterline angle. The latter half ensures a smooth transition with the parallel body. The relationship between ice thickness and speed is shown in Figure 17 for a mooring column tension of 254.9 t, power of 15 000 kW, and open-water propeller speed of about 5.9 m/s. When the flat ice is thin, the ship moves faster and the water resistance increases significantly. However, as the ice thickness increases, the speed of the optimized ship is significantly improved compared with that of the original ship.

Figure  16  Stem lines plan of waterlines before and after optimization

Download: Full-Size Img

Figure  17  Relationship between ice thickness and speed before and after optimization

Download: Full-Size Img

5   Conclusion

This paper has described a multi-objective optimization process for enhancing the total resistance and average pressure in the ice zone of an icebreaker. A fully parameterized model-based optimization platform combining the cyclic process of contact compression bending failure and the Lindqvist formula was proposed. A sensitivity analysis was applied to obtain information on the trends regarding the significant parameters of the stem, the impact of ice thickness on the icebreaking load, and the most important region for the pressure. Based on numerous sample points, the weighting coefficients of each objective were determined, and the objective function of pressure was modified to make the optimized bow parameters similar to those of the original ship. Geometry parametrization focused on the waterline angle, buttock angle, two feature lines, and scanning curves for shape deformation. The position and weight factors of cubic NURBS were defined for parametric modeling, allowing a comprehensive design to be acquired. Accordingly, in the second phase of optimization, CAD was combined with the DE algorithm and a resistance solver in a platform environment, and multi-threaded optimization calculations were performed using the WaitForMultipleObjects function. Process automation made it possible to manage the design using an optimization algorithm, ultimately yielding the global optimum design without intervention.

The optimal values of the design variables were obtained by investigating the Pareto optimal solution set. The optimization result is a hull form that reduces the resistance by 4.2% and decreases the average pressure in the ice zone by 0.6%. The buttock angle and the first sub-region are the most important components in enhancing the icebreaker's performance. Resistance is positively correlated with the buttock angle and negatively correlated with the waterline angle, while the pressure in the ice zone exhibits the opposite trends. A single control point parameter has no direct relationship with the performance of the icebreaker, but affects the icebreaking resistance by influencing the waterplane coefficients and the shape of the waterline. In the optimized design, the buttock angle and waterline angle of the icebreaker have both increased by 6%–9%. A smaller weight factor makes the connection between the rear half of the bow waterline and the parallel hull smoother and has little effect on the immersion resistance caused by broken ice. The curvature of the waterline is greater, which makes it easier for level ice to bend and break. The optimized waterline curvature changes uniformly in the direction of the ship length, which pertains to lower ship hull resistance. A comparison between the initial and optimized hulls and the optimization process results demonstrates the validity of the proposed optimization design strategy.

This study can be further improved. The resistance of polar ships can be further optimized by considering both open-water and icebreaking performance. Using historic ice data and ship routes, the navigation state recognition model aids voyage assessments, enabling hull optimization for polar ships in terms of minimizing resistance in both open-water and ice-covered conditions, while maintaining the desired icebreaking capability. The energy efficiency design index can be used as the objective of the optimization function and the main dimensions of the ship can be designated as design variables.