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

Abstract

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.

1. 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 self-icebreaking 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 in waterline angle leads to the decrease in longitudinal crushed ice force. However, the bow needs more contact area and extrusion force to bend the ice with the decrease in frame angle. Therefore, the purpose of these changes of the bow angles is to achieve a superior icebreaking capability of a polar ship. It seems that the bow configurations of modern polar ships are quite similar. The ice resistance may nevertheless show notable differences. Consequently, all the changes of the bow configuration aim to obtain a superior icebreaking capability of polar ships navigating in ice region.

Several methods for ice resistance prediction have been proposed, which take into account the influence of geometric variations of the hull. In practical ship design, empirical or semi-empirical prediction methods may be applied to estimate the resistance in the preliminary design phase. Erceg and Ehlers [7] reviewed six semi-empirical level ice resistance prediction methods (Vance [8]; Lewis et al. [9]; Zahn and Phillips [10]; Keinonen [11,12]; Lindqvist [13]; Riska et al. [14]). Ice resistances of four ships with different sizes were calculated and compared based on the six semi-empirical prediction methods. The methods show significant discrepancies for larger ships, with differences in the predicted ice resistance as large as 100–200%. In the study of Lindström [15], most of the relevant ice parameters and hull configurations including the bow angles were regarded in the ice resistance prediction. Based on the study of Lindqvist [13], Riska et al. [14] proposed a formula for calculating the level ice resistance. This formula takes into account the particulars of the ship and the ice thickness. The empirical coefficients are derived from full-scale tests performed in the Baltic Sea. Jeong [16] decomposed the ice resistance into ice buoyancy resistance, ice-clearing resistance and icebreaking resistance. Using the component method, Jeong [16] conducted an ice model test of MOERI standard icebreaker model, and non-dimensional coefficients in ice resistance prediction formula are obtained based on the data of model test. Lubbad and Løset [17] established new analytical closed-form solutions to represent the icebreaking process of ship–ice interaction in real-time. These works have offered good insights into the method of estimating the ship performance accounting for ship–ice interaction.

Because the experimental research can reveal real physical effects of ship–ice interaction, some detailed analyses have been carried out to study the bow configurations based on the full-scale experiment as well as the ice model tests. Varges [18] investigated the bow shape of an icebreaker based on the full-scale experiment. The vessel of concern was the Soviet polar icebreaker MUDYUG, which was built in 1972 and given a new Thyssen/Waas forebody in 1986. The conventional angled bow was replaced with a form more resembling a sled. The sleigh-shaped bow has better performance and less vibration in the icebreaking process. Based on the ice model test results of HSVA’s database from 1996 to 2014, Myland and Ehlers [19] focused in particular on the icebreaking patterns in level ice with regard to relevant hull shape parameters. The ship models with similar scales were chosen for analysis, so as to have model ice conditions, model speed values and ship model dimensions in the same range.

More recently, numerical simulations, e.g., discrete element method (DEM) and finite element method (FEM), have been widely used to analyze the influence of ship hull forms on the ice performance. Lu et al. [20] applied FEM to calculate the ice load acting on conical structures colliding with level ice. The mesh size dependency of the bending failure of level ice was eliminated by taking into account the random spatial distribution of material properties. Wang and Zou [21] used FEM to investigate the effects of the wedge stem and the shovel stem on the damage of the ice layer and the resistance during the icebreaking process. The result of wedge-shaped bow and scoop bow in the process of icebreaking were compared, and some suggestions for the design of icebreaker hull lines were given. Li et al. [22] estimated the ice load acting on the shoulder and midship of a polar ship by applying the extended finite element method (XFEM). Long et al. [23] adopted DEM with the bond and failure model to simulate the breaking process of sea ice colliding with a conical structure and vertical pile. The results showed that the bending failure occurs as the ice cover moves upward or downward, depending on the water level. Wang 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.

2. 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].

2.1. 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⁻³ s⁻¹) 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⁻¹) 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:

$$\varphi =J_2-\frac{{\sigma }_y^2}{3}$$

(1)

where the second stress invariant, J₂, is defined in terms of the deviatoric stress components as

$$J_2=\frac{1}{2}{s}_{ij}{s}_{ij}$$

(2)

and the yield stress, ${\sigma }_y$, is a function of the effective plastic strain, ${\epsilon }_{eff}^p$, and the plastic hardening modulus, E_p:

$${\sigma }_y={\sigma }_0+E_p{\epsilon }_{eff}^p$$

(3)

Pressure is given by the expression

$$p^{n+1}=K\left(\frac{1}{V^{n+1}}-1\right)$$

(4)

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. The parameters used in the numerical simulation.

Materials	Items	Values
Cone-shaped ice	Mass density	920 kg/m3
	Yield stress	2.12 MPa
	Shear modulus	2.20 GPa
	Plastic failure strain	0.35
	Failure pressure	−4 MPa
Ice sheet	Mass density	906 kg/m3
	Thickness	0.77 m/0.76 m
	Bending strength	724 kPa/844 kPa
	Compressive strength	1748 kPa/2192 kPa
	Elastic modulus	929 MPa/984 MPa
Water	Mass density	998 kg/m3
	Acoustic velocity	1650 m/s
	Cut-off pressure	10 N/m2
	Dynamic viscosity coefficient	1.0 × 10−3 N·s/m2
	Volume correction	1.39
	Gruneisen gamma	0.49
	Internal energy	0.25 Kg/m2s2

.

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 Er_f is defined as the errors between experimental results [30] and SPH results, which are calculated by:

$$E{r}_f=\sqrt{{\displaystyle \sum _{i=1}^n{(f_i-f_0{}_i)}^2}}/\sqrt{{\displaystyle \sum _{i=1}^n{f}_0{}_i^2}}$$

(5)

where f_i is the compression force of SPH with different particle numbers at t = t_i; and f₀ 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⁻² s⁻¹ for v = 1 mm/s, and 1 s⁻¹ for v = 100 mm/s. Both fall into the range of high strain rates (>10⁻³ s⁻¹), 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.

2.2. 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. Full-scale primary dimensions of MT Uikku.

Items	Values
Length	150 m
Length of parallel	65 m
Breadth	21 m
Tested draft	9.50 m

. 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. Test matrix with the model speed and ice properties.

Test No.	hi (m)	σb (kPa)	σc (kPa)	Ei (MPa)	Vi (m/s)
1	0.77	724	1748	929	0.2
2	0.76	844	2192	984	0.5

, 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.

2.3. 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.

2.4. 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].

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:

$$E{r}_r=\left|r_0-r\right|/r_0$$

(6)

where r is the ice resistance of SPH with different particle numbers from t = 0 s to t = 30 s, r₀ 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. Ice resistances of numerical simulations and experiments.

Test No.	Ri num (kN)	StDev of Ri num (kN)	CV of Ri num	Ri exp (kN)	Er
1	503	156	0.31	553	9.04%
2	618	99	0.16	691	10.56%

, respectively. The standard deviations (StDev) of R_{i num} are 156 kN for Test 1 and 99 kN for Test 2. The standard deviations of ice resistance versus the mean value (CV) for R_{i num} are 0.31 for Test 1 and 0.16 for Test 2. The relative errors (E_r) between R_{i num} and R_{i 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.

3. 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.

3.1. 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. Primary dimensions of MV Xue Long.

Items	Values
Length overall	167.0 m
Length between perpendiculars	147.2 m
Breadth	22.6 m
Draught	9.0 m
Displacement	21,025 tons
Deadweight	8759 tons
Propulsion power	13.2 MW
Bollard pull	1650 kN
Maximum calm-water speed.	18 knots

. 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.

3.2. 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.

$$T_{net}=T_B\left[1-\frac{1}{3}\cdot \frac{u}{v_{ow}}-\frac{2}{3}\cdot {\left(\frac{u}{v_{ow}}\right)}^2\right]$$

(7)

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 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. Analysis of variance for the ice resistance.

Case	α	γ	β	Indicator (kN)
1	24°	22°	52°	895
2	24°	24°	54°	945
3	24°	26°	56°	1065
4	26°	22°	54°	865
5	26°	24°	56°	1025
6	26°	26°	52°	1155
7	28°	22°	56°	795
8	28°	24°	52°	1025
9	28°	26°	54°	1095
K1	2905	2555	3075	T = 8865
K2	3045	2995	2905	T2 = 78,588,225
K3	2915	3315	2885	P = 8,732,025
K12	8,439,025	6,528,025	9,455,625	QT = 8,842,825
K22	9,272,025	8,970,025	8,439,025	—
K32	8,497,225	10,989,225	8,323,225	—
QA	8,736,091.67	882,9091.67	8,739,291.67	—
SA	4066.67	97,066.67	7266.67	—

.

3.3. 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:

$$S_T=Q_T-P$$

(8)

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.

$$Q_T={\displaystyle \sum _{k=1}^n{x}_k^2}$$

(9)

$$P=\frac{1}{n}{\left({\displaystyle \sum _{k=1}^n{x}_k}\right)}^2$$

(10)

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

$$S_A=Q_A-P$$

(11)

where

$$Q_A=\frac{1}{a}{\displaystyle \sum _{i=1}^{n_a}{K}_i^2}$$

(12)

$$K_i={\displaystyle \sum _{j=1}^a{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. Significance test of the bow configuration parameters to the ice resistance.

Source	α	γ	β	Error
Sum Sq.	4066.67	97,066.67	7266.67	2400
d.f.	2	2	2	2
Mean Sq.	2033.33	48,533.33	3633.33	1200
F Value	1.70	40.50 *	3.01	—
F0.05(2,2)	19.00	19.00	19.00	—
F0.01(2,2)	99.01	99.01	99.01	—

Note: * is the significant factor affecting the indicator.

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 test error f_E = f_T $-$ f_α $-$ f_γ $-$ f_β = 2.

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°).

3.4. 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. Grid information and the resulting average calm-water resistance.

Grid	Grid Number	Average Calm-Water Resistance (kN)
Coarse (3)	1,357,071	193.875
Medium (2)	3,088,921	188.520
Fine (1)	4,401,879	188.385

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:

$$p=\ln \left[\left(193.875-188.520\right)/188.520-188.385\right]/\ln \left(2\right)=5.31$$

Richardson extrapolation is performed to predict an estimate of the value of the average calm-water resistance at infinite grid number,

$$f_{n\to \infty }=188.385+\left(188.385-188.520\right)/\left(2^{5.31}-1\right)={188.382}_{}(\mathrm{kN})$$

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:

$$GC{I}_{23}=1.25\left|\left(188.520-188.385\right)/188.385\right|/\left(2^{5.31}-1\right)\cdot 100\%=0.002316\%$$

The grid convergence index (GCI) for the coarse and medium refinement levels is:

$$GC{I}_{12}=1.25\left|\left(193.875-188.520\right)/188.520\right|/\left(2^{5.31}-1\right)\cdot 100\%=0.091826\%$$

Grids are ensured in the asymptotic range of convergence by checking:

$$0.091826\%/\left(2^{5.31}0.002316\%\right)=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. C_r of the polar ship navigating in a combined route.

Case	Rice (N)	Rwater (N)	kice	kwater	Cr (N)
1	8.95 × 105	1.81 × 105	0.2	0.8	3.24 × 105
4	8.65 × 105	1.88 × 105	0.2	0.8	3.23 × 105
7	7.95 × 105	1.97 × 105	0.2	0.8	3.17 × 105

. 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.

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,

$$C_r={\displaystyle \sum _{i=1}^m{k}_{water}{}_i{R}_{water}{}_i}+{\displaystyle \sum _{j=1}^n{k}_{ice}{}_j{R}_{ice}{}_j}$$

(14)

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.

4. 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.