An Approach to Determine Optimal Bow Configuration of Polar Ships under Combined Ice and Calm-Water Conditions

Selecting an optimal bow configuration is critical to the preliminary design of polar ships. This paper proposes an approach to determine the optimal bow of polar ships based on present numerical simulation and available published experimental studies. Unlike conventional methods, the present approach integrates both ice resistance and calm-water resistance with the navigating time. A numerical simulation method of an icebreaking vessel going straight ahead in level ice is developed using SPH (smoothed particle hydrodynamics) numerical technique of LS-DYNA. The present numerical results for the ice resistance in level ice are in satisfactory agreement with the available published experimental data. The bow configurations with superior icebreaking capability are obtained by analyzing the sensitivities due to the buttock angle γ, the frame angle β and the waterline angle α. The calm-water resistance is calculated using FVM (finite volume method). Finally, an overall resistance index devised from the ship resistance in ice/water weighted by their corresponding weighted navigation time is proposed. The present approach can be used for evaluating the integrated resistance performance of the polar ships operating in both a water route and ice route.


Introduction
The bow configuration of a polar ship is determined by a basic idea, i.e., the minimum power to move forward in the ice. From the design point of view, level ice resistance is also one of the most significant factors. There are many cases of polar ships with selficebreaking ability sailing in level ice, for example, PUGNAX (polar deck carrier) [1], MV Xue Long (polar scientific research vessel) and MV Xue Long II. At the current stage of polar exploration, this kind of ship sailing in level ice is also considered as an important condition. The primary requirement for these polar ships is a good performance in level ice. Good performance means low ice resistance, high propulsive efficiency and guaranteed continuity in icebreaking [2]. Some researches [3][4][5] suggested that ice resistance is the main contribution to the total resistance in level ice. Apart from the ice properties and the power of polar ships, the ice resistance is significantly influenced by the bow configurations. In the evolution of the bow configuration, the frame angle has been increasing, while the waterline angle has decreased, and the buttock angle has reduced from about 30 • to 20 • in the early years [6]. From the perspective of ship-ice interaction, a description about the changes of the bow angles is explained below. When the ship is operating in ice region, the ice is crushed at the stem post of the bow. The buttock angle has been reduced because a smaller buttock angle induces more bending moment and less horizontal force. The waterline angle mainly affects the crushed ice-pushing capability of the bow. The decrease et al. [24] combined an elastoplastic softening constitutive model with cohesive element method (CEM) to simulate the continuous icebreaking process in level ice. This method captured well the main features, including the local crushing and the bending failure, of the ship-ice interactions.
In the past few years, the smoothed particle hydrodynamics (SPH) method has been used to effectively model large deformation and failure behavior of solids, including ice. Das and Ehlers [25] carried out a numerical simulation of crushing and bending failure of ice using SPH. In their study, the numerical bending results for force, displacement and failure time were compared with earlier simulations of in situ four-point bending test results and finite element simulations. The crushing failure was compared with experiments conducted by Häusler [26], and it was found that the SPH approach was in good agreement with the experiment. Furthermore, the SPH method is extended so as to simulate the actual behavior of sea ice as ships progress through level ice. Zhang et al. [27] reported that the SPH results agreed well with the experimental data of three-point bending test. This implies that the present SPH model can produce accurate results for simulating the ice failure problem. Zhang et al. [27] also showed a good agreement between the results of the SPH and the experimental data for the ice-ship interaction. They indicated that the numerical accuracy and stability of SPH were satisfactory.
Usually, polar ships also navigate in water ways, but the calm-water performance is ignored in the researches above. Especially for polar ships converted from the merchant vessel, which, most of the time, are operated in both water ways and ice regions. It is necessary to conduct integrated analyses that account for the resistance both in ice and in water. However, the literature on the optimization of bow configuration of polar ships navigating in combined routes is limited. Polach et al. [28] proposed a method which includes the ship performance evaluation based on ship merit factor (SMF). This method combined SMF with a route-specific ship and allows one to compare the technical and economic performance of ships operating in open water and ice. The performance of different ship designs that operated along the route Rotterdam to Yokohama through the Suez Canal and the Northern Sea Route was studied.
This study aims at devising an approach to determine the optimal bow of polar ships operating in both water route and ice route. In the paper, a numerical simulation of an icebreaking vessel going straight ahead in level ice is performed using the SPH numerical technique. The present numerical results are compared with the experimental data by Zhou et al. [29]. A polar research vessel, MV Xue Long, is used as a reference ship to demonstrate the usage of the proposed approach. Comparative analyses are carried out to evaluate the sensitivities of the buttock angle γ, the frame angle β and the waterline angle α on the icebreaking capability. Integrated numerical simulations are conducted to study the icebreaking capability and calm-water performance of the polar ship with different bow configurations. An overall resistance index C r is proposed to determine the optimal bow configuration of polar ships navigating in a combined route.

Validation of Ice Material Model
According to the mechanical properties of ice, the constitutive model of ice is established using the SPH method. The numerical model is validated by comparing the results with the experiment by Kim et al. [30]. Furthermore, a ship-ice-water interaction model is proposed for the simulation of a ship moving in level ice. The present numerical results for the icebreaking patterns and the ice resistance are then compared with the experimental data by Zhou et al. [29].

Verification of Ice Material Model
During the process of deformation to failure of ice, the ductile-brittle transition occurs due to the effect of the loading rate. The behavior of ice at high strain rates (greater than 10 −3 s −1 ) is similar to the linear elastic material with brittle failure [31]. The material model takes into account high strain rates and is widely used to model the dynamic behavior of ice including elastic deformation and brittle fracture. Strain rate (s −1 ) is the change in strain of a material with respect to time. In this study, the strain rate effect was transformed to the strain effect included in the material model. An isotropic elastic-plastic material model with a von Mises yield surface was used in the numerical simulation.
The von Mises yield condition is given by: where the second stress invariant, J 2 , is defined in terms of the deviatoric stress components as and the yield stress, σ y , is a function of the effective plastic strain, ε p e f f , and the plastic hardening modulus, E p : Pressure is given by the expression where n is the number of time steps, K is the bulk modulus, and V is the relative volume defined as the ratio of the current volume over the initial volume. When an element reaches the plasticity phase under the action of compression, it follows the yield surface until failure. An element is removed by erosion if the pressure surpasses the failure pressure. Meanwhile, failure by element erosion is activated if the equivalent plastic strain reaches the plastic failure strain. The values including plastic failure strain and failure pressure for the cone-shaped ice were also used in the ice sheet simulations. When the ice breaks, compression failure acquirement is met, as shown in Table 1. In order to verify the accuracy of the above constitutive model of ice, a numerical simulation was carried out on the compressive cone-shaped ice experiments [30]. The models of cone-shaped ice and steel plate for the present numerical simulation are shown in Figure 1a. The SPH method is used to simulate the cone-shaped ice with a diameter of 100 mm and a cone angle of 120 degrees. In the bottom plane of the numerical cone-shaped ice, the movement of SPH particles is constrained to six degrees of freedom to simulate the ice fixed on the experimental device. A rigid body is used to simulate the steel plate. The cone-shaped ice is vertically compressed by the steel plate at the speed of 1 mm/s and 100 mm/s. The brittle fracture was observed in the numerical calculations, as shown in Figure 1b,c.   Figure 2a gives the convergence tests on compression force, in which N is the total particle number. Three particle numbers of 27,974; 37,196 and 51,192 are tested in the compressive cone-shaped ice. The relative error Erf is defined as the errors between experimental results [30] and SPH results, which are calculated by: where fi is the compression force of SPH with different particle numbers at t = ti; and f0 is the compression force of experimental data at t = ti. The results of Figure 2a indicate that the errors of the compression force decrease with the increasing particle number and demonstrate a roughly first-order convergence ratio. Therefore, 37,196 particles are deemed adequate for the simulation of ice cone. The force-displacement curves of the numerical simulation at the compression speed of 1 mm/s and 100 mm/s are presented in Figure 2b,c. Numerical results are compared with the experimental data [30]. It is found that, in both the high-speed loading and the low-speed loading, the ice force increases with the displacement. As shown in Figure 2b,c, the force pattern is captured by the numerical simulation with satisfactory accuracy. The numerical method agrees better with the experiment in the amplitude of ice force for lower loading speeds. It also can be seen that more fluctuations are present in the solution for higher loading speed. When the ice  Figure 2a gives the convergence tests on compression force, in which N is the total particle number. Three particle numbers of 27,974; 37,196 and 51,192 are tested in the compressive cone-shaped ice. The relative error Er f is defined as the errors between experimental results [30] and SPH results, which are calculated by: where f i is the compression force of SPH with different particle numbers at t = t i ; and f 0 is the compression force of experimental data at t = t i . The results of Figure 2a indicate that the errors of the compression force decrease with the increasing particle number and demonstrate a roughly first-order convergence ratio. Therefore, 37,196 particles are deemed adequate for the simulation of ice cone. The force-displacement curves of the numerical simulation at the compression speed of 1 mm/s and 100 mm/s are presented in Figure 2b,c. Numerical results are compared with the experimental data [30]. It is found that, in both the high-speed loading and the low-speed loading, the ice force increases with the displacement. As shown in Figure 2b,c, the force pattern is captured by the numerical simulation with satisfactory accuracy. The numerical method agrees better with the experiment in the amplitude of ice force for lower loading speeds. It also can be seen that more fluctuations are present in the solution for higher loading speed. When the ice is broken to a certain length, the ice force does not increase but decrease due to the unloading of the pressure. In Figure 2b,c, the strain rate is 10 −2 s −1 for v = 1 mm/s, and 1 s −1 for v = 100 mm/s. Both fall into the range of high strain rates (>10 −3 s −1 ), where the ice has the behavior as brittle materials. Moreover, the same cut-off stress value is used in the present numerical ice model. In fact, the cut-off stress decreases slightly with the increase in loading speed in the brittle range. This is why the peak contact forces are somewhat close (35 and 38 kN, respectively), and the present numerical method shows a little overestimation of the amplitude for a higher loading speed. The comparisons show that the present numerical ice model can give a similar trend of the compression process of ice with that of the experiment, although there are some differences in details.
is broken to a certain length, the ice force does not increase but decrease due to the un-loading of the pressure. In Figure 2b,c, the strain rate is 10 −2 s −1 for v = 1 mm/s, and 1 s −1 for v = 100 mm/s. Both fall into the range of high strain rates (＞10 −3 s −1 ), where the ice has the behavior as brittle materials. Moreover, the same cut-off stress value is used in the present numerical ice model. In fact, the cut-off stress decreases slightly with the increase in loading speed in the brittle range. This is why the peak contact forces are somewhat close (35 and 38 kN, respectively), and the present numerical method shows a little overestimation of the amplitude for a higher loading speed. The comparisons show that the present numerical ice model can give a similar trend of the compression process of ice with that of the experiment, although there are some differences in details.

Description of Model Test
Zhou et al. [29] conducted a series of ice resistance tests where the model of an icebreaking vessel, MT Uikku, was towed in an ice basin. A scale model of 1:31.6 was built without any appendices such as rudders and propellers. The full-scale primary dimensions of MT Uikku are given in Table 2. The model was towed with the carriage to simulate the ice-hull interaction process. In the model tests, the ice drift angle and the hull heading angle were constant at 0°. The ice force was measured and digitally sampled on a computer and recorded as an analogue signal as a back-up. The measured data were

Description of Model Test
Zhou et al. [29] conducted a series of ice resistance tests where the model of an icebreaking vessel, MT Uikku, was towed in an ice basin. A scale model of 1:31.6 was built without any appendices such as rudders and propellers. The full-scale primary dimensions of MT Uikku are given in Table 2. The model was towed with the carriage to simulate the ice-hull interaction process. In the model tests, the ice drift angle and the hull heading angle were constant at 0 • . The ice force was measured and digitally sampled on a computer and recorded as an analogue signal as a back-up. The measured data were sampled at a frequency of 107 Hz. The test program and measured ice properties are summarized in Table 3, where h i is ice thickness; σ b and σ c are bending strength and compressive strength of ice; E i is the elastic modulus of ice; and V i is ship speed.  Table 3. Test matrix with the model speed and ice properties.

Numerical Simulation
A finite element ship model of the MT Uikku is created in the numerical simulation. The ship form parameters are presented in Table 2. The icebreaking pattern and the ice resistance are investigated regardless of the response of ship structures. Therefore, only the hull of the ship is established and simplified to a rigid body in the present numerical simulation. The model is going straight ahead in the x-direction at predefined constant speeds. Based on the SPH numerical technique of LS-DYNA hydrocode, the problem of a ship moving in level ice is simulated. The dimensions of domain are 600 m × 150 m × 75 m. The level ice floats on water, and its two long sides are constrained in the x-direction.
In the mathematical description of fluid, the equation of state describes the relationship between volume deformation and pressure. The null material is used to describe the water flow. The Gruneisen equation is adopted to describe the state of water. The parameters used in the numerical simulation are shown in Table 1. The force data are sampled every 0.01 s. Figure 3 shows the typical configuration of the numerical model applied to calculate the ice resistance. sampled at a frequency of 107 Hz. The test program and measured ice properties are summarized in Table 3, where hi is ice thickness; σb and σc are bending strength and compressive strength of ice; Ei is the elastic modulus of ice; and Vi is ship speed.

Numerical Simulation
A finite element ship model of the MT Uikku is created in the numerical simulation. The ship form parameters are presented in Table 2. The icebreaking pattern and the ice resistance are investigated regardless of the response of ship structures. Therefore, only the hull of the ship is established and simplified to a rigid body in the present numerical simulation. The model is going straight ahead in the x-direction at predefined constant speeds. Based on the SPH numerical technique of LS-DYNA hydrocode, the problem of a ship moving in level ice is simulated. The dimensions of domain are 600 m × 150 m × 75 m. The level ice floats on water, and its two long sides are constrained in the x-direction.
In the mathematical description of fluid, the equation of state describes the relationship between volume deformation and pressure. The null material is used to describe the water flow. The Gruneisen equation is adopted to describe the state of water. The parameters used in the numerical simulation are shown in Table 1. The force data are sampled every 0.01 s. Figure 3 shows the typical configuration of the numerical model applied to calculate the ice resistance.

Comparisons and Discussions of Results
Numerical simulations are performed to simulate the Test 1 and Test 2 in Table 3. Figure 4 shows the icebreaking stress and the icebreaking pattern in the numerical simulation of Test 2. As shown in Figure 4a, the stress is generated by the deformation of the ice, which mainly occurs in the contact area between ice and hull. A concentrated area of high stress is present near the bow, and larger ice pieces are observed in the bow and shoulder area. As the hull goes straight ahead in level ice, the open channel is generated in Figure 4b. As can be seen in Figure 4a,b, the width of the channel is approximately

Comparisons and Discussions of Results
Numerical simulations are performed to simulate the Test 1 and Test 2 in Table 3. Figure 4 shows the icebreaking stress and the icebreaking pattern in the numerical simulation of Test 2. As shown in Figure 4a, the stress is generated by the deformation of the ice, which mainly occurs in the contact area between ice and hull. A concentrated area of high stress is present near the bow, and larger ice pieces are observed in the bow and shoulder area. As the hull goes straight ahead in level ice, the open channel is generated in Figure 4b. As can be seen in Figure 4a,b, the width of the channel is approximately equal to the beam of the model ship. The icebreaking patterns obtained in the present numerical simulation is similar to that in Zhou et al. [29]. equal to the beam of the model ship. The icebreaking patterns obtained in the present numerical simulation is similar to that in Zhou et al. [29].
(a) (b)  Figure 5a gives the convergence tests on ice resistance, in which M is the total particle number. Three particle numbers of 210,555; 262,709 and 307,601 are tested in the icebreaking simulations. The relative error Err is defined as the errors between experimental results [29] and SPH results, which are calculated by: where r is the ice resistance of SPH with different particle numbers from t = 0 s to t = 30 s, r0 is the ice resistance of experimental data from t = 0 s to t = 30 s. The results of Figure 5a indicate that the errors of the ice resistance decrease with the increasing particle number and demonstrate a roughly first-order convergence ratio. Hence, 262,709 particles are deemed adequate for the simulation of icebreaking. Time histories of ice forces in the longitudinal direction of Test 1 and Test 2 (Table 3) from the present simulation and the experiment [29] are presented in Figure 5b,c. Since the ice resistance commonly refers to the time average of all longitudinal forces due to ice acting on the ship, a comparison of results in a quantitative way is conducted. The mean values of the ice forces at the steady stage are also presented in Figure 5b,c. The differences between numerical results and model test data are identified to a certain extent. The ice force fluctuates during the icebreaking process while the force is never zero. This is because that there are always multiple ice contacts along the hull as the ship moves forward. These contacts are a series of events during the process, and the ice forces are the summation of all the contact forces. In addition, ice submergence resistance is a significant contribution. Some previous numerical simulation, e.g., [3], shows that submersion plays a significant role even at slow speed. The existence of submersion is also one of the reasons that the time history of ice forces rarely drops to zero. The mean values of the present simulation results for the ice forces, i.e., ice resistance, approach those of the experimental data from [29]. The fluctuation amplitudes of the simulated ice force are slightly higher than those of the experimental results. Figure 4b shows that the pieces of large ice cusps generated in the present numerical simulation are more than those observed in the experiment [29]. Data processing methods in Wang et al. [24] are referenced in Figure 5b,c. In order to eliminate the noise from the results, the numerical force curves are plotted by processing the original data through the filter used in the experiment. Higher fluctuation amplitudes of calculated ice resistance have been reduced. The results of the present numerical simulations are in better agreement with the experimental data. The ice resistances of numerical simulations and experiments are denoted by Ri num and Ri exp in Table 4 Figure 5a gives the convergence tests on ice resistance, in which M is the total particle number. Three particle numbers of 210,555; 262,709 and 307,601 are tested in the icebreaking simulations. The relative error Er r is defined as the errors between experimental results [29] and SPH results, which are calculated by: where r is the ice resistance of SPH with different particle numbers from t = 0 s to t = 30 s, r 0 is the ice resistance of experimental data from t = 0 s to t = 30 s. The results of Figure 5a indicate that the errors of the ice resistance decrease with the increasing particle number and demonstrate a roughly first-order convergence ratio. Hence, 262,709 particles are deemed adequate for the simulation of icebreaking. Time histories of ice forces in the longitudinal direction of Test 1 and Test 2 (Table 3) from the present simulation and the experiment [29] are presented in Figure 5b,c. Since the ice resistance commonly refers to the time average of all longitudinal forces due to ice acting on the ship, a comparison of results in a quantitative way is conducted. The mean values of the ice forces at the steady stage are also presented in Figure 5b,c. The differences between numerical results and model test data are identified to a certain extent. The ice force fluctuates during the icebreaking process while the force is never zero. This is because that there are always multiple ice contacts along the hull as the ship moves forward. These contacts are a series of events during the process, and the ice forces are the summation of all the contact forces. In addition, ice submergence resistance is a significant contribution. Some previous numerical simulation, e.g., [3], shows that submersion plays a significant role even at slow speed. The existence of submersion is also one of the reasons that the time history of ice forces rarely drops to zero. The mean values of the present simulation results for the ice forces, i.e., ice resistance, approach those of the experimental data from [29]. The fluctuation amplitudes of the simulated ice force are slightly higher than those of the experimental results. Figure 4b shows that the pieces of large ice cusps generated in the present numerical simulation are more than those observed in the experiment [29]. Data processing methods in Wang et al. [24] are referenced in Figure 5b,c. In order to eliminate the noise from the results, the numerical force curves are plotted by processing the original data through the filter used in the experiment. Higher fluctuation amplitudes of calculated ice resistance have been reduced. The results of the present numerical simulations are in better agreement with the experimental data. The ice resistances of numerical simulations and experiments are denoted by R i num and R i exp in Table 4 relative errors (Er) between Ri num and Ri exp are 9.04% for Test 1 and 10.56% for Test 2. Although this observation indicates that the present numerical simulation cannot capture all the details of the broken ice floes around the hull, the approach adopted here, as shown in the comparative analyses above, is appropriate for evaluating the ice resistance performance. Thus, the present approach can serve as a tool to study the bow configuration of the polar ship.

Approach to Determine the Optimal Bow Configuration of Polar Ships
Based on the dimensions of MV Xue Long, a full-scale hull model is established. The full-scale model is taken as the master model. Nine hull models are created with different bow configurations so as to analyze the icebreaking capability. Then, the models are moved in the calm water (without ice) at a speed of 15 knots. The workflow of the approach is shown in Figure 6.

Approach to Determine the Optimal Bow Configuration of Polar Ships
Based on the dimensions of MV Xue Long, a full-scale hull model is established. The full-scale model is taken as the master model. Nine hull models are created with different bow configurations so as to analyze the icebreaking capability. Then, the models are moved in the calm water (without ice) at a speed of 15 knots. The workflow of the approach is shown in Figure 6.

Ship Model and Bow Configuration Parameters
An icebreaking scientific research vessel, MV Xue Long, is modeled for the present study. MV Xue Long is an icebreaking research vessel owned by the Polar Research Institute of China. It was built at the Kherson Shipyard in Ukraine in 1993. It started as an icebreaking cargo and supply ship designed for the Russian Arctic. It was then converted from an Arctic cargo ship to a polar research and resupplying vessel by the mid-1990s. The ice class of MV Xue Long, assigned by the China Classification Society (CCS), is B1 [32].
The primary dimensions of MV Xue Long are presented in Table 5. As the present study focuses on the bow configuration instead of the structural response, the ship is modeled as a rigid body without any hull structures. The finite element model of the hull is shown in Figure 7. The bow configuration can be characterized by the buttock angle γ, the frame angle β and the waterline angle α [33], as illustrated in Figure 8.   Figure 6. General workflow for the approach to determine the optimal bow configuration.

Ship Model and Bow Configuration Parameters
An icebreaking scientific research vessel, MV Xue Long, is modeled for the present study. MV Xue Long is an icebreaking research vessel owned by the Polar Research Institute of China. It was built at the Kherson Shipyard in Ukraine in 1993. It started as an icebreaking cargo and supply ship designed for the Russian Arctic. It was then converted from an Arctic cargo ship to a polar research and resupplying vessel by the mid-1990s. The ice class of MV Xue Long, assigned by the China Classification Society (CCS), is B1 [32].
The primary dimensions of MV Xue Long are presented in Table 5. As the present study focuses on the bow configuration instead of the structural response, the ship is modeled as a rigid body without any hull structures. The finite element model of the hull is shown in Figure 7. The bow configuration can be characterized by the buttock angle γ, the frame angle β and the waterline angle α [33], as illustrated in Figure 8.  Figure 6. General workflow for the approach to determine the optimal bow configuration.

Ship Model and Bow Configuration Parameters
An icebreaking scientific research vessel, MV Xue Long, is modeled for the present study. MV Xue Long is an icebreaking research vessel owned by the Polar Research Institute of China. It was built at the Kherson Shipyard in Ukraine in 1993. It started as an icebreaking cargo and supply ship designed for the Russian Arctic. It was then converted from an Arctic cargo ship to a polar research and resupplying vessel by the mid-1990s. The ice class of MV Xue Long, assigned by the China Classification Society (CCS), is B1 [32].
The primary dimensions of MV Xue Long are presented in Table 5. As the present study focuses on the bow configuration instead of the structural response, the ship is modeled as a rigid body without any hull structures. The finite element model of the hull is shown in Figure 7. The bow configuration can be characterized by the buttock angle γ, the frame angle β and the waterline angle α [33], as illustrated in Figure 8.

Schemes of Bow Configuration
An analysis is conducted to analyze the effect of γ on the icebreaking capability, which is done by varying γ with α and β being fixed. In the simulations, the ship model is going straight ahead in level ice (hice = 0.5 m). The speed the ship can attain in level ice is calculated based on the present numerical simulation. The ice resistance equals the extra thrust available at different power levels and speeds. The net thrust curves used in the analysis are based on the theoretical calculation given by Equation (7) [34] and are considered as a preliminary curve. The average longitudinal speed the ship achieves in the designated ice conditions is determined from the intersection of the net thrust curve and the ice resistance curve.
where TB is the bollard pull; u is the initial ship speed; and vow is the maximum calm-water speed. Figure 9a shows the effect of γ on the average longitudinal ship speed. It can be seen that the relationship between the longitudinal ship speed in the continuous icebreaking process and γ is approximately a downward curve. The main reason is that a small γ induces a large bending moment and a small horizontal force. The effects of α and β on the average longitudinal ship speed are also analyzed using in the same manner. Figure 9b shows the average longitudinal ship speed for α = 17°, 21°, 25°, 29° and 33° with γ fixed to 23° and β to 53°. It is observed that the longitudinal ship speed decreases slightly with α increases. The reason is that the longitudinal crushed ice force Fx increases with α and results in the decrease in the longitudinal ship speed, as shown in Figure 9d. Unlike γ and α, β has a non-monotonic effect on the longitudinal ship speed when γ is fixed to 23° and α to 25°, as shown in Figure 9c. A fluctuating pattern of the influence of β on the average longitudinal speed is observed. However, the difference between the maximum and the minimum ship speed is less than 4%. It suggests that β has little effect on the total ice resistance and the ship speed.

Schemes of Bow Configuration
An analysis is conducted to analyze the effect of γ on the icebreaking capability, which is done by varying γ with α and β being fixed. In the simulations, the ship model is going straight ahead in level ice (h ice = 0.5 m). The speed the ship can attain in level ice is calculated based on the present numerical simulation. The ice resistance equals the extra thrust available at different power levels and speeds. The net thrust curves used in the analysis are based on the theoretical calculation given by Equation (7) [34] and are considered as a preliminary curve. The average longitudinal speed the ship achieves in the designated ice conditions is determined from the intersection of the net thrust curve and the ice resistance curve.
where T B is the bollard pull; u is the initial ship speed; and v ow is the maximum calm-water speed. Figure 9a shows the effect of γ on the average longitudinal ship speed. It can be seen that the relationship between the longitudinal ship speed in the continuous icebreaking process and γ is approximately a downward curve. The main reason is that a small γ induces a large bending moment and a small horizontal force. The effects of α and β on the average longitudinal ship speed are also analyzed using in the same manner. Figure 9b shows the average longitudinal ship speed for α = 17 • , 21 • , 25 • , 29 • and 33 • with γ fixed to 23 • and β to 53 • . It is observed that the longitudinal ship speed decreases slightly with α increases. The reason is that the longitudinal crushed ice force F x increases with α and results in the decrease in the longitudinal ship speed, as shown in Figure 9d. Unlike γ and α, β has a non-monotonic effect on the longitudinal ship speed when γ is fixed to 23 • and α to 25 • , as shown in Figure 9c. A fluctuating pattern of the influence of β on the average longitudinal speed is observed. However, the difference between the maximum and the minimum ship speed is less than 4%. It suggests that β has little effect on the total ice resistance and the ship speed. To determine the optimal bow configuration, schemes of combinations of α, γ and β are obtained using the orthogonal design method [35]. A comprehensive research with three factors (α, γ and β) at three levels (i.e., three values of each factor) requires twentyseven different cases. By using the orthogonal design method, there are only nine nonrepetitive cases. As shown in Figure 10, the nine cases are evenly distributed in the entire research domain. The bow configuration schemes and parameters used in the orthogonal design are listed in Table 6.   To determine the optimal bow configuration, schemes of combinations of α, γ and β are obtained using the orthogonal design method [35]. A comprehensive research with three factors (α, γ and β) at three levels (i.e., three values of each factor) requires twenty-seven different cases. By using the orthogonal design method, there are only nine non-repetitive cases. As shown in Figure 10, the nine cases are evenly distributed in the entire research domain. The bow configuration schemes and parameters used in the orthogonal design are listed in Table 6. To determine the optimal bow configuration, schemes of combinations of α, γ and β are obtained using the orthogonal design method [35]. A comprehensive research with three factors (α, γ and β) at three levels (i.e., three values of each factor) requires twentyseven different cases. By using the orthogonal design method, there are only nine nonrepetitive cases. As shown in Figure 10, the nine cases are evenly distributed in the entire research domain. The bow configuration schemes and parameters used in the orthogonal design are listed in Table 6.

Sensitivity of the Icebreaking Capability to the Bow Configuration Parameters
Sensitivity studies are carried out to analyze the influence of bow configurations on the icebreaking capability. Analysis of variance is employed to identify the critical affecting factor in the significance tests. The variance of the test results originates from the test errors and the variations of different factors: where S T denotes the total variance of the test results. Q T is the sum of squares of the indicators, and P is the revised value.
where n is the total number of tests, and x k is the indicator for each test. The variance of each factor can be expressed as where x ij (13) where K i denotes the sum of test results for factor A at level i. a is the number of tests for each level, and x ij is the indicator of j-th test for factor A at level i. Taking the ice resistance as the indicator of the icebreaking capability, Table 6 presents the variance analysis results for the significance test, where T is the sum of indicators. The variance of factors at three levels are ranked in the order of S γ , S β , and S α .
The total variance of the test results S T = Q T − P = 110,800, and the variance of the test error S E = S T − S α − S γ − S β = 2400. Table 7 shows the significance test results of the effects of α, γ and β on the ice resistance, where d.f. denotes the degrees of freedom. The total d.f. of the test f T = n − 1 = 8, the d.f. of each factor f A = n a − 1 = 2, and the d.f. of the F is defined as the mean squared deviation (Mean Sq.) of each factor divided by the mean squared deviation of test error, which represents the magnitude of the sensitivity of each factor to the test results. Here, the F values of α, γ and β are 1.70, 40.5 and 3.01, respectively. Compared with the F-critical value, i.e., F 0.05 (2, 2) = 19.0, which is obtained by a joint hypotheses test, the F value of γ is greater, indicating that γ has a significant effect on the test results. In other words, γ is the most important factor to the ice resistance, marked with *. In contrast, the F values of α and β are much less than 19.0, suggesting that the effect of α and β on the ice resistance is insignificant. Therefore, the bow configurations from Case 1, Case 4 and Case 7 have an excellent icebreaking capability, which is achieved with the minimum value of γ (i.e., γ = 22 • ).

Integrated Evaluation Based on the Resistance Performance
For polar ships sailing in both water and ice regions, it is economically beneficial to consider both water resistance and ice resistance by using a typical route for conceptual design. In the present study, the calm-water resistance of Case 1, Case 4 and Case 7 is calculated using the FVM-based commercial code of STAR-CCM+. The fluid domain and mesh of the present numerical model are shown in Figure 11a. The ship speed is 15 knots for the estimation of calm water resistance, which represents the design speed during operation. A grid convergence study and time-step study are performed in order to verify the present numerical model. Table 8 indicates the grid information and the resulting average calm-water resistance. Three grid numbers of 1,357,071; 3,088,921, and 4,401,879 are tested, respectively. Figure 11b shows the plot of average calm-water resistance with varying grid numbers. As the grid number increases, the average calm-water resistances approach an asymptotic infinite-grid number value. Three simulations (coarse, medium and fine) are completed with a constant refinement ratio r = 2. The order of convergence, p, is calculated using: This value is also plotted on Figure 11b. The grid convergence index for the fine grid solution can now be computed. A factor of safety of F S = 1.25 is used since three grids were used to estimate the average calm-water resistance. The grid convergence index (GCI) for the medium and fine refinement levels is: which is approximately one and indicates that the solution is well within the asymptotic range of convergence. Therefore, 3,088,921 grids are sufficient to produce accurate results, which are used to calculate the calm-water resistance.
For the time-step convergence study, the time step is selected based on the numerical simulations in which a variety of Courant numbers, C = 0.075, 0.1, 0.25, and 0.5 are carried out. The Courant number describes the relationship between the time step and the space step. Intuitively speaking, the Courant number is the number of grids that a fluid particle can pass through in a time step. As shown in Figure 11c 0.091826% 2 0.002316% 0.99938 = which is approximately one and indicates that the solution is well within the asymptotic range of convergence. Therefore, 3,088,921 grids are sufficient to produce accurate results, which are used to calculate the calm-water resistance. For the time-step convergence study, the time step is selected based on the numerical simulations in which a variety of Courant numbers, C = 0.075, 0.1, 0.25, and 0.5 are carried out. The Courant number describes the relationship between the time step and the space step. Intuitively speaking, the Courant number is the number of grids that a fluid particle can pass through in a time step. As shown in Figure 11c, when the Courant numbers of 0.075 and 0.1, the calm-water resistance obtained with Courant number being 0.075 and 0.1 is insignificant. Therefore, the Courant number C of 0.1 is used in this present study.    Figure 12 shows the time histories of the calm-water resistance for Case 1, Case 4 and Case 7. The calm-water resistance amplitude in the stable range, from 80 to 200 s, is used to calculate the mean value, i.e., the average calm-water resistance. Comparison of the numerical results shows that Case 7 has the maximum calm-water resistance. As can be seen in Table 6, Case 7 has the fullest bow, corresponding to the largest α and β, resulting in the maximum calm-water resistance. Results for the average calm-water resistance (R water ) and the ice resistance (R ice ) are presented in Table 9. It was found that the average calm-water resistance for Case 4 is about 3.87% larger than Case 1, and the average calm-water resistance for Case 7 is 4.79% larger than Case 4. The ice resistance for Case 4 is 3.35% smaller than for Case 1, and the ice resistance for Case 7 is 8.09% smaller than for Case 4. R water and R ice are significantly influenced by the bow configuration. The present results indicate that full bows are favorable for the performance in the ice but not beneficial to the performance in the calm-water. In addition, Table 9 shows that total resistance in the ice region is more than four times higher than in water region, because it is assumed that the ship speed is fixed to 15 knots in the water route and the invariant ice thickness is 0.5 m in the ice region. In this case, the ice resistance is approximately four times higher than the water resistance. However, the speed changes when the ship navigates along the ice routes according to the ice condition, such as the variety of ice thickness. In addition, the ship probably does not sail at the same speed in the water region. Therefore, the ratio of the ice resistance to the water resistance is not constant.  Figure 12 shows the time histories of the calm-water resistance for Case 1, Case 4 and Case 7. The calm-water resistance amplitude in the stable range, from 80 to 200 s, is used to calculate the mean value, i.e., the average calm-water resistance. Comparison of the numerical results shows that Case 7 has the maximum calm-water resistance. As can be seen in Table 6, Case 7 has the fullest bow, corresponding to the largest α and β, resulting in the maximum calm-water resistance. Results for the average calm-water resistance (Rwater) and the ice resistance (Rice) are presented in Table 9. It was found that the average calmwater resistance for Case 4 is about 3.87% larger than Case 1, and the average calm-water resistance for Case 7 is 4.79% larger than Case 4. The ice resistance for Case 4 is 3.35% smaller than for Case 1, and the ice resistance for Case 7 is 8.09% smaller than for Case 4. Rwater and Rice are significantly influenced by the bow configuration. The present results indicate that full bows are favorable for the performance in the ice but not beneficial to the performance in the calm-water. In addition, Table 9 shows that total resistance in the ice region is more than four times higher than in water region, because it is assumed that the ship speed is fixed to 15 knots in the water route and the invariant ice thickness is 0.5 m in the ice region. In this case, the ice resistance is approximately four times higher than the water resistance. However, the speed changes when the ship navigates along the ice routes according to the ice condition, such as the variety of ice thickness. In addition, the ship probably does not sail at the same speed in the water region. Therefore, the ratio of the ice resistance to the water resistance is not constant.  In most ice regions, the ship speed changes when a ship navigates along the ice routes according to the ice condition, such as the variety of ice thickness. If the ship captain wants  In most ice regions, the ship speed changes when a ship navigates along the ice routes according to the ice condition, such as the variety of ice thickness. If the ship captain wants to achieve a certain speed in the ice region, the ship needs to overcome resistance in ice. For that purpose, propeller must provide sufficient thrust, i.e., the ship engine is selected based on this. Ship probably does not sail at the same speed in ice region and in region without ice due to significantly different total resistances. In this case, the entire route of a polar ship should be divided by different segments based on ship speed from the harbor to the ice region. Thus, an overall resistance index C r comprising the ship resistance and the corresponding weighted navigation time is proposed as follows, where k water i and k ice j are the weight factors of the i-th route segment in the water region and the j-th route segment in the ice region, respectively; R water i and R ice j are the resistance of the i-th route segment in the water region and the j-th route segment in the ice region, correspondingly.
In the present study, C r is calculated by only one uniform water segment and one uniform ice segment. It is assumed that the speeds of a polar ship in calm water and in level ice are both constant. For MV Xue Long, the weighted factor k water is assumed to be 0.8, and k ice is 0.2. C r of Case 1, Case 4 and Case 7 is shown in Table 9, where C r of Case 7 is the smallest. Therefore, Case 7 is selected as the optimal bow configuration for overall performance in both water and ice. It should be noted that k water and k ice are consistent with the actual navigation time of the polar ship of concern. For example, in August 2017, the Russian icebreaking LNG carrier Christophe de Margerie [36] passed through the Northern Sea Route (NSR) in six days and completed the entire 19-day journey from Hammerfest, Norway, to Boryeong, South Korea. In this case, k water = 0.68 and k ice = 0.32 can be used according to the article posted by Rachael [36]. The present method can be used to improve and optimize the bow configuration by assessing the resistance performance. The optimal route of a polar ship from the harbor to the ice region can also be studied using Equation (14). For the operation of vessels in the different routes, the mileage is considered in k water and k ice .

Conclusions
An approach to determine the optimal bow of polar ships based on present numerical simulation and available published experimental studies was proposed. In order to validate the constitutive model of ice, numerical simulations were performed using SPH, and the results were compared with the experiment by Kim et al. [21]. Then, a ship-ice-water interaction model was proposed for the simulation of a ship moving in level ice. The present numerical results for the ice resistance in level ice were in satisfactory agreement with the experimental data by Zhou et al. [20]. This method was then used to analyze the ice resistance of a polar ships with various bow configuration parameters, including the buttock angle γ, the frame angle β and the waterline angle α. Sensitivity analyses of the ice resistance to these parameters were evaluated by performing the analysis of variance. It was found that the effect of γ on the ice resistance is much more significant than that of α and β; thus, small γ values are desirable.
To assess the overall resistance performance of a ship that is operated in both ice and waters regions, the resistance in water was calculated using the FVM-based method. Then, an overall resistance index C r devised from the ship resistance in ice/water weighted by their corresponding weighted navigation time was proposed. The formula to calculate the index takes a simple form but is particularly useful in practice as the resistances and weights can be easily obtained. Since the bow geometry is the most important factor to the overall resistance performance, the proposed formula provides a convenient approach to determine the optimal bow configuration of polar ships. A polar research vessel, MV Xue Long, was used as a reference ship to demonstrate the usage of the proposed formula. In most cases, the entire route of a polar ship should be divided by different segments based on ship speed from the harbor to the ice region. In the present study, C r was calculated by one uniform water segment and one uniform ice segment. It was assumed that the speeds of a polar ship in calm water and in level ice are both constant. More calculations by Equation (14) are needed when the route is divided by different segments based on ship speed. In addition, the present approach is applied to the evaluation based on the resistance performance of polar ships moving in level ice, which is a typical condition for polar ships with self-icebreaking ability. For merchant polar ships that may be operated in pack ice, improvement and validation of the approach shall be made for the desired accuracy.