Prediction of Ice-Resistance Distribution for R/V Xuelong Using Measured Sea-Ice Parameters

: Increasing human activity in polar areas is making ice-going ships more indispensable in multiple operation scenarios. An improvement in ice-resistance prediction, which cannot be performed without accurate ice parameters, will promote the development of hull design and operational safety in ice-infested waters. The Nataf transformation is applied to generate correlated pseudo-random numbers which represent ice parameters; then, as a numerical method, the circumferential crack method is introduced to calculate the ice resistance of R/V Xuelong in level ice. The main factors which may have a large inﬂuence on simulated ice load are studied. The simulation results show that the Burr distribution is the most suitable model to describe the distribution of ice resistance calculated and ice-force amplitude concentrated at a lower level. The statistical results are also discussed and compared with similar research through empirical formulas and Monte Carlo methods. The present simulation can obtain more detailed information during the icebreaking process compared to similar research: the ice force at each time step is achieved; the key ice-force amplitude can be collected, which can beneﬁt studies on hull structure; and potential stress generated by sea ice can be predicted. The present numerical tools and simulation results can provide a reference for ice-going ships sailing in level ice in most scenarios with regard to ice resistance and operational safety.


Introduction
With the increase in human activity in polar areas, ice-going ships play a significant role in commercial transportation, logistics supply, and personnel transportation. These ships cannot avoid colliding with sea ice when navigating in cold regions; therefore, such ships require the capability to break ice or be equipped with an ice-strengthened hull in order to sail in waters covered by ice and snow. An accurate prediction of ice resistance encountered by a ship in various situations is necessary at the design stage or during transportation planning, thereby promoting the development of structural design and operational safety in polar areas.
To date, a large number of studies on the ice-breaking process have been carried out. Lindqvist [1] and Riska et al. [2] proposed some semi-empirical and analytic forecast models, which provided effective methods for the early design of icebreakers. Su et al. [3] simulated ship maneuvers in ice areas involving the surge, sway, and yaw motions; however, only a single ice-failure pattern was considered. Lubbad and Løset [4] applied PhysX to simulate the motions of ice floes in the calculation domain. Zhou et al. [5,6] studied the relationship between the hull slope angle, ice friction coefficient, and sea-ice failure pattern, and analyzed the mathematical model of each component of ice load. In this paper, we apply a numerical simulation method based on a semi-empirical formula to simulate the straight operation of a ship in level ice, combining the advantages of the numerical simulation and empirical formula to make sure that accuracy and efficiency is guaranteed. This numerical tool can calculate ice resistance based on the interaction between the hull and sea ice. Geometric contact between the hull and ice was carried out at each time step in order to judge the ice-failure pattern and whether the sea ice was broken, and the ice force at each hull node was calculated and collected. Key ice force at certain time steps can be detected and extracted, and such results can be used as a reference resource for ship design and ship operation in polar areas.
Ice parameters such as physical and mechanical properties are required as calculation input in these methods, e.g., ice thickness, ice density, ice flexural strength, ice crushing strength, temperature, and elastic modulus [7]. These parameters show a high degree of randomness in nature [8], which have a great influence on ice-hull interaction, breaking the pattern of sea ice and ice resistance during interactions with different kinds of structures. Sodhi and Chin [9] used Gaussian, Rayleigh, and Weibull distributions for indentation tests in urea-model ice. The ice was then failed by crushing, and the distribution of peak ice forces followed the Weibull distribution. Takeuchi et al. [10] conducted research on the spatial distribution of the uniaxial compressive strength of ice, as well as the correlation with temperature and density. Uniaxial strength and other mechanical parameters fitted well with a normal distribution. Johnston et al. [11] reviewed the measured data of multiyear ice thickness in polar areas and proposed a lognormal distribution for the data. In 2012, Strub-Klein [12] analyzed the spatial and temporal distribution of uniaxial compressive strength, temperature, density, and salinity of first-year ice, then applied a methodology for statistical analysis of uniaxial compressive strength to full-scale ice data in 2017 [13]. Ni [14] applied a one-way computational fluid dynamics-discrete element coupling method to simulate the motion of a ship sailing in level ice. The ice resistance showed strong oscillation and was mainly affected by the ship's velocity and ice thickness. Wang [15] studied the trend of ice resistance under crushed ice. The ice resistance increased with the increase in speed, but when the ship's speed reached a certain value, the ice resistance no longer increased, and even reduced the trend. Zong [16] studied an ice-ship collision. Hydrodynamic influence turned out to be significant for elastic and elastoplastic modes of ice. Chai [17] studied the statistics of ice thickness and the ice strength of first-year sea ice; a scheme based on empirical formulas was proposed to investigate the correlation between different ice parameters. A novel probabilistic model was introduced to describe the distribution of ice flexural strength. However, statistical studies of long-term analysis and estimations of ice resistance are rare. Measured short-term ice-induced loads show that the Weibull probability distribution [18][19][20][21] or the exponential or logarithmic normal probability distribution [22,23] match well with measured data after later statistical research. Based on measured data collected by R/V Xuelong, Ji [24] found that ice-induced hull vibration amplitude increases with both ice thickness and ship speed. Wang et al. [25] studied the probability distribution of ice parameters based on measurements collected by R/V Xuelong in scientific voyages, fixed ice sites, or empirical equations. It was found that ice resistance follows the distribution of Gamma. Ice resistance was calculated based on Lindqvist's semi-empirical formulae [1], which provided only an approximate estimate of ice resistance. However, this method ignored geometric contact and weakened the interaction between the hull and the ice. The distribution of ice resistance needs to be analyzed by a more refined method.
On the basis of the probabilistic model proposed by Wang [25], a method which can generate pseudo-random numbers is adopted to produce the distribution and correlation of ice parameters. Ice data are then used as input for the empirical numerical simulation tool mentioned above in order to calculate the ice resistance of R/V Xuelong [26,27]. The ice resistance is obtained after the whole sailing motion is accomplished. Finally, the probability distribution of the ice resistance is analyzed and compared with related research.

Numerical Method
The ice-load model was first introduced, which refined ice-failure patterns into common bending-failure patterns, and crushing-failure patterns occurred at a vertical frame angle. Dynamic global ice force was calculated based on hull-ice interaction and semiempirical formulae; then, we used Nataf transformation to process statistical model and generate corresponding ice parameters as input into the simulation [28].

Circumferential Crack Method
Vertical compression and horizontal tension of ships act continuously on sea ice during voyages in polar areas; then, ice fails at bending and crushing on most occasions. Circumferential cracks parallel to the contact surface will take the lead under compression and tension. The sea-ice-failure mechanism adopted is the occurrence of circumferential crack based on this physical phenomenon observed from model tests and ice trials, which assumes that ice blocks break from ice sheets in a vertical direction, and that their contact surface is flat and their contact area can be determined by its length and depth [4].
The geometric shape of the ice floes broken from level ice can be assumed to be wedgeshaped, as shown in Figure 1. The shape of the crack is represented by a red curve, and the original intact ice sheet is represented in blue. The angle of ice wedge is θ, and the ice wedge is shaped based on the icebreaking radius R, the expression of which is given in reference [29]: where C l and C v are empirical parameters, v rel n is the relative normal velocity between ice and hull, and l is the characteristic length of ice; where E i is the young's modulus of ice, h i is the thickness of ice, υ is the poisson ratio of ice, ρ w is the density of water, and g is the acceleration of gravity.

Mechanical Model of Interaction between Ice and Hull
The localized crushing force F cr generated on the contact surface of hull and sea ice can be expressed as: where A c is the contact area between ice and hull, σ c is the crushing strength of ice. Figure 2 shows velocities and ice forces during the icebreaking process. The parameters v rel n and v rel τ are the normal and tangential components of the relative velocity between hull and ice along the hull contact nodes. The parameters v rel n,1 and v rel n,2 are the components of v rel n , along and perpendicular to the contact surface in the vertical plane, respectively. f H and f V are the horizontal and vertical components of the frictional force, respectively. The combined crushing force and the component of the frictional force can be divided into the horizontal component F H and vertical component F V . Forces mentioned above can be expressed in Su et al.'s study, where, µ i is the frictional coefficient between hull and ice, ϕ is the slope angle of hull, α is the water entrance angle. According to Zhou et al. [26], the failure patterns of ice vary under the influence of different hull slope angles. Friction coefficient will affect the slope angle to distinguish the failure patterns. The difference in failure patterns will lead to different sizes of ice wedge and different ice forces, especially at the shoulder parts of the ship where crushing failure of ice usually occurs in turning or transverse towing operations, which should be distinguished in the simulation. Considering the ice-structure friction coefficient as 0.0394, a limit angle of 87 • is chosen to separate the different ice-failure patterns.

Bending-Failure Pattern
When slope angle ϕ of hull is less than the limit angle, and the vertical component F V of crushing force and friction force is large enough to exceed the ice-bearing capacity, bending failure occurs. Ice will break from the ice sheet as an intact ice wedge. In the simulation, the waterline of ship is represented with a bunch of discrete nodes which form a closed polygon as the hull; level-ice boundary is also distinguished into discrete nodes. The interaction between ship and ice is transformed into the interaction between a closed polygon and discrete nodes. Whether the ice boundary consisting of discrete nodes enters the closed polygon at each time step in the simulation is detected; contact area and other parameters are collected through the program.
The vertical component F V increases with the increase in contact area. When F V meets the ice-bearing capacity, initial ice boundary around the contact area will be updated to a new ice wedge. Due to the distribution of low slope angle, bending failure always occurs around the bow area. Ice force at this time step is then calculated and stored in the program for calculation in the next step. According to the circumferential crack method mentioned above, when bending failure happens to level ice under the action of icebreaker hull, the opening angle and the radius of an ice wedge formed due to bending failure can use parameters θ and R in expression (1) to present. The top part of the ice wedge is under the load in vertical direction. The bearing capacity of ice P f can then be introduced based on reference [30]: where σ f is the flexural strength of ice, h i is the ice thickness, C f is the empirical parameter. According to the previous work conducted by Su et al. [31], the value of C f is 2.9.

Crushing-Failure Pattern
When the slope angle ϕ exceeds limit angle and the icebreaking condition is achieved, according to ISO/FDIS 19906-2019 [32] and without the limit mechanism, the localized crushing force acting on the hull can be expressed as: where P G is the global average ice pressure; L is the projected width of the hull; h is the ice thickness; h 1 is the reference thickness of 1 m; m and n are empirical parameters, the value of m is −0.16 and n is equal to −0.5 + h/5 for h < 1.0 m and to −0.3 for h > 1.0 m; and C R is the ice-strength coefficient, which changes under different regions. As shown in Figure 3, the length of the contact area projected on the hull structure indicated as a red line is calculated by program at each time step, and increases with the advance of icebreaker; when the crushing length L reaches the predefined limitation of crushing length, ice will break from the sheet. In the program, the initial ice boundary will be replaced by the updated ice boundary represented as red curve; ice force at this time step will be calculated and stored. According to the model test, the predefined limitation of crushing length we choose as the crushing failure limit condition is ch, where c is a constant. The limitation of crushing length increases with the increase in ice thickness, which indicates that ships need to reach wider crushing lengths in order to break thicker ice.

Excitation Forces and Moments
The excitation forces and moments involved in simulation mainly include the ice force F ice (t), generated due to the interaction between hull and level ice. The detailed calculation formula of each part of ice force is explained below.
Ice force includes the continuous ice-breaking forces F brk (t) and ice-submerging forces. The total continuous icebreaking forces are calculated based on the local icebreaking forces, which change under the different parts of the hull at each time step. Local icebreaking forces are calculated and collected with ice-failure mechanism summarized above. The ice-submerging forces are calculated according to ice-force formula [1]; then the ice force in general can be expressed as: where v rel is the relative velocity between ice and hull; v rel 1 is the forward component of v rel ; v rel 2 is the transverse components of v rel ; L is the waterline length of the ship; and R S is the submersion component of ice force [1], which can be written as: where ρ w is the sea density; ρ i is the density of the level ice; B is the beam; T is the draught; µ is the frictional coefficient between ice and hull; ϕ is the stem angle of the ship; and α is the water entrance angle, ψ = arctan tanϕ sinα . After obtaining each component of ice force, the contact area and position between hull and ice is processed by computer at each time step in the simulation. Ice force and updated level-ice boundary can be used as input for the next time step; then ice channels created by ship and ice resistance in level ice can be achieved, as well as the time series curves of ice force.

Nataf Transformation
According to ice-property measurements along R/V Xuelong's scientific sailing and ice sites in polar areas, it can be found that ice thickness follows distribution of Burr XII, while ice density follows Beta distribution. Gamma distribution and GEV distribution can be applied to describe the distribution of ice flexural strength and ice effective modulus, respectively [25]. Based on the distribution of ice parameters and correlation coefficient matrix, we can use Nataf transformation to generate pseudo-random numbers as corresponding ice parameters. These ice parameters can be used in the numerical simulation to calculate the ice resistance exposed to R/V Xuelong in polar areas.
The Nataf transformation is a tool to rebuild joint distribution once the distribution of input parameters is obtained. We can obtain ice properties by following specific distribution through inverse Nataf transformation. In Nataf transformation, X = [x 1 , x 2 , . . . , x n ] is the n-dimensional random variable. The probability density function and cumulative distribution function of random variable x i are f i (x i ) and F i (x i ), respectively. According to literatures [33,34], the relationship between random variable X = [x 1 , x 2 , . . . , x n ] and standard normal random variable Y = [y 1 , y 2 , . . . , y n ] can be described as: where φ(·) is the cumulative distribution function of standard normal distribution, assuming that correlation coefficient matrix of X and Y are ρ and ρ 0 , respectively. According to Nataf transformation, elements ρ ij and ρ 0ij (i = j) at corresponding positions in matrix follow the below relationship: where, µ i and µ j are the expectation of x i and x j , respectively; σ i and σ j are the standard deviation of x i and x j , respectively; and φ 2 y i , y j , ρ 0ij is the joint distribution function of two-dimensional standard normal random variables y i and y j with correlation coefficient of ρ 0ij . Correlation coefficient ρ ij of ice properties is shown in Table 1. Correlation coefficient ρ 0ij of standard normal random variable is calculated by solving the root of nonlinear function in program; the results are shown in Table 2.  The correlation coefficient matrix ρ 0 can be expressed using Cholesky factorization: where L 0 is lower triangular matrix of ρ 0 after Cholesky factorization, and T represents matrix transposition. Correlation standard normal random variables Y can be transformed into independent standard normal random variables U through L 0 : Nataf transformation is then established. Inverse Nataf transformation of independent standard normal random variables U transformed into correlation standard normal random variable Y can be expressed as: N-dimensional random variable X = [x 1 , x 2 , . . . , x n ] that presents different ice parameters can be generated and obtained; 1000 samples are generated for each ice parameter.
The distribution of each generated ice parameter follows the same distribution as revealed in [25], with a small margin of error. Ice thickness follows the distribution of Burr Due to sea-ice-failure mechanism considering the crushing failure of sea ice in ice model, the crushing strength of sea ice is required; however, many previous studies revealed only the distribution of crushing strength. Quantifying the relationship of crushing strength with other ice parameters is a difficult task and has not been clearly studied. We used other methods to generate this key parameter. According to Timco and Weeks's research [7], the ratio of crushing strength and flexural strength of sea ice is 0.5-5, while model scale measurements show a ratio of 1.5-2.7 in model test by Zhou et al. [5], which narrows this range. The ratio of crushing strength and flexural strength we choose is 2.12, which is the average ratio of 10 model test cases. Crushing strength follows the distribution of Gamma, with the parameter of a = 19.8, b = 27.8.

Numerical Simulation Setup
The hull used in the simulation is the icebreaker R/V Xuelong, which undertook the scientific expedition and the field measurements. This vessel was built in 1993 and has the capability to break level ice with a thickness of 1.2 m and snow of 0.2 m at 0.8 m/s. The main parameters of the ship and ice parameters are shown in Table 3. The ship was set to sail in level ice with a fixed forward velocity of 0.8 m/s. The ship sailed for 200 s in total, and the time interval was 0.1 s in the simulation. The ice thickness in one case was constant. At the initial stage, the ship kept a certain distance away from the ice and was going to collide at the boundary of level ice. The random ice boundary settled around the hull; however, there was no contact, as shown in Figure 4. The ice boundary is represented by the blue dotted line. It distinguishes the open water from the level ice, and between the two parallel blue curves is the ice channel created for the icebreaker. The waterline of the ship is represented by the black dotted line. With ship parameters and generated pseudo-random numbers, a simulation can be performed to predict ice-resistant ship encounters.

Sensitivity Analysis
In order to deeply reveal the relationship of each ice parameter (ice thickness, flexural strength, crushing strength, effective modulus, and ice density) with ice resistance, sensitivity analysis on parameters was carried out. Based on the initial case, we changed the value of one parameter at a time; meanwhile, other parameters remained the same. It should be noted that the ice parameters selected in the initial case are summarized by the ship's trial in the Arctic in order to achieve appropriate representativeness, while such a goal is hard to reach via pseudo-random numbers. A total of eight cases were conducted for each parameter. The parameters of the initial case and the corresponding changing intervals are listed in Table 4.

Ice-Thickness Sensitivity Analysis
We conducted sensitivity analysis on ice thickness first. Ice thickness is one of the most intuitive parameters affecting the icebreaking process. The changing range of ice thickness is from 0.6 m to 1.3 m, with intervals of 0.1 m. The trend of ice resistance under rising ice thickness is shown in Figure 5. As shown in Figure 5, the trend of ice resistance roughly increases with the increase in ice thickness; the trend fluctuation is small. Ice resistance increases rapidly when ice thickness exceeds 0.9 m. Ice thickness has a great influence on ice resistance; ice resistance ranges from 0.34 MN to 1.01 MN-a great span of 0.67 MN-while ice thickness ranges from 0.6 m to 1.3 m.
The reason why a strict upward trend was not presented may be that the change in ice thickness will lead to a change in failure size and bearing capacity of sea ice; this will make interactions between the hull and the ice more complex, so certain unclear trends in simulation are acceptable.

Flexural Strength Sensitivity Analysis
The sensitivity analysis of flexural strength was then carried out, the parameters of which are involved in the bearing capacity of ice in the ice model. The changing range of flexural strength is from 300 kPa to 650 kPa, with intervals of 50 kPa. The trend of ice resistance under rising flexural strength is shown in Figure 6. Ice resistance increases strictly with the rise in flexural strength, as shown in Figure 6. Ice resistance ranges from 0.30 MN to 0.70 MN. The span of ice resistance is 0.40 MN. This parameter is related to the bearing capacity of ice in the ice model; the increase in flexural strength will make sea ice harder to break, which matches the total trend of resistance.

Crushing Strength Sensitivity Analysis
The changing range of crushing strength is from 1500 kPa to 2900 kPa, with intervals of 200 kPa. The trend of ice resistance under rising crushing strength is presented, shown in Figure 7. Ice resistance decreases as crushing strength increases. Crushing strength has little impact on ice resistance compared to other parameters; the presented curve is flat. The changing range of ice resistance is from 0.632 MN to 0.473 MN; the span of ice resistance is 0.159 MN.

Effective Modulus Sensitivity Analysis
The changing range of effective modulus is from 3.8 GPa to 6.6 GPa, with intervals of 0.4 GPa. The trend of ice resistance under rising effective modulus is shown in Figure 8. The trend of ice resistance under changing effective modulus is hard to describe and an overall fluctuating trend is presented, as shown in Figure 8. The influence of a different effective modulus is quite small. Ice resistance oscillates at around 0.52 MN, ranging from 0.506 MN to 0.543 MN. The span of ice resistance is 0.037 MN.

Ice-Density Sensitivity Analysis
Finally, the sensitivity analysis of ice density was studied. The changing range of ice density is from 0.80 g/cm 3 to 0.94 g/cm 3 , with intervals of 0.02 g/cm 3 . The trend of ice resistance under rising ice density is shown in Figure 9. Ice resistance decreases as ice density increases, as shown in Figure 9. The changing rate of this parameter is stable. Ice resistance ranges from 0. 582 MN to 0.468 MN, with a span of 0.114 MN; the impact of ice density on ice resistance is stable, and this parameter only affects the ice-submerging force; such a trend can be explained in simulation.

Influence of Different Ice Parameters on Resistance
The influence of different ice parameters on resistance is studied based on the initial case. The dimensionless number is introduced to study the extent of influence, shown in Figure 10. As shown in Figure 10, different ice parameters are represented by different colors and shapes; their relationships with ice resistance are revealed by dimensionless numbers. The extent of influence of ice parameters can be reflected directly by the changing rate of curves. Ice thickness, represented by the black line, has the greatest positive influence on resistance among all parameters. Ice density, represented by the green line, has a great effect on resistance, second only to ice thickness; however, such effect is negative. Flexural strength, represented by orange triangles, shows the third-greatest impact on ice resistance. According to the figure, effective modulus, represented by the red line, has the least effect on the result; the curve is almost flat. Finally, crushing strength, represented by the blue dotted line, has the second-smallest impact on the results.

Analysis of Numerical Simulation Results
Considering the calculation time required for simulation and efficiency, 400 sets of cases using generated ice parameters were carried out in order to study the distribution pattern of ice resistance. After completing all 400 cases, the ice parameters and corresponding ice resistance of five random cases were collected and are listed in Table 5 to study individual cases. According to the sea-ice-failure mechanism, a difference in ice properties will not only affect ice resistance directly, but also lead to different ice-failure shapes created by the ship or the breaking frequency of sea ice. The detailed icebreaking process is studied to reveal such a phenomenon. The ice channel created during the icebreaking process and the corresponding ice force encountered by the ship in the surge direction are shown in Figures 11 and 12.  As shown in Figures 11 and 12, two typical cases are studied. The simulation results of two cases are significantly different, with a span of 1.85 MN. The radius of the ice wedge induced by ice-hull interaction in case 2 is larger than that of case 3, according to Figures 11a and 12a. The number of ice wedges created during the contact in case 3 is more than that of case 2. One period in simulation is identified by the corresponding maximum value of ice force; then, ice-force time history curves show that ice force duration in the period of case 3 is similar to that of case 2, which is about 8.6 s in case 2 and 9 s in case 3. However, one period of duration in case 3 can be separated into several load events due to the increasing number of ice wedges and the weakening ice-bearing capacity. Both cases indicate that the peak value of ice force will flatten over time, and the icebreaking process will follow a regular pattern.
The ice-resistance distribution of 400 sets of cases is shown in Figure 13. Ice resistance ranged from 0.07 to 3.42 MN, and the average ice resistance was 0.74 ± 0.40 MN with a median value of 0.70 MN. The Burr distribution turned out to be the most suitable model to describe the distribution of ice resistance obtained, which was also accepted by K-S test. The parameters of Burr distribution are estimated as a = 0.90, c = 2.92, and k = 1.91, with standard errors of 0.09, 0.18, and 0.38, respectively. In addition, based on the rain flow counting method [35], we processed all of the recorded ice force in the surge direction and obtained the amplitude of ice force. A total of 66,599 ice force amplitudes were extracted from 400 cases. These ice-load amplitudes are more convincing than ice resistance in regard to hull safety, since ice resistance only reflects the average value of ice force in the surge direction during one sail, whereas a large amount of ice force exceeds that value in the icebreaking process. If we focused only on the average value and ignored the high value of ice force, the ice force could be underestimated and cause damage to the ship during sailing in polar areas. The ice-force amplitudes and probability density curves of all ice-force amplitudes are shown in Figures 14 and 15.
As shown in the figure, 91.3% of ice-force amplitudes are concentrated at a lower level of 0-1 MN. The number of ice-force amplitudes decreases rapidly when the amplitude exceeds 1 MN, and the number of ice-force amplitudes exceeding 2 MN is quite rare, close to 2.3%. The maximum ice-force amplitude is 8.845 MN in the study; extreme ice force amplitude may be caused by extreme ice parameters. Such situations cannot be ignored, because these situations still can happen in real life, albeit with a low probability. It was concluded that the ice force encountered by ships seldom exceeds 2 MN in most navigation scenarios; however, some extreme cases are still worth a deep study in the future.   Figure 16 shows a comparison of the probability density curve and cumulative frequency curve, presented by Wang et al. [25] and simulated with the present method. The ice resistance calculated with the present numerical tools, represented by the red line, ranges from 0.07 to 3.42 MN, and the average ice resistance is 0.74 ± 0.40 MN with a median value of 0.70 MN, while ice resistance ranges from 0.06 to 2.69 MN, and the average ice resistance is 0.73 ± 0.34 MN with a median value of 0.68 MN by the Monte Carlo method, represented as black line. The differences in average value and median value between two results are quite small, and two probability density curves show the same positively skewed bell shape. Both distributions reveal that ice resistance tends to be concentrated in the area of low resistance, around 0.4-0.8 MN. The probability decreases significantly when the resistance exceeds about 1.5 MN. In addition, the upper limit of ice resistance is higher than that of the Monte Carlo method, which is also reflected in the figures; the value of the upper limit is 3.42 MN in the present numerical simulation. Ice resistance obtained through numerical simulation has high probability in the range of 0.7-0.8 MN, while ice resistance calculated by Wang has high probability in the range of 0.5-0.6 MN.

Comparison and Discussion
Although the probability density curves of the two results are quite similar, results follow different distributions due to the diversity in results. Gamma distribution is the most fitting model for the Monte Carlo method's results, with parameters of a = 4.6 and b = 0.2, and corresponding standard errors of 0.06 and 0.002, respectively. The Burr distribution turned out to be the most suitable model to describe the distribution in the present research. Reasons for such differences can be summarized as follows: crushing failure occurrences on the steep hull surface and geometric contact between the ice and hull were considered in the numerical simulation, while the Monte Carlo method was based on semi-empirical formulae taken before. In addition, 400 cases were performed due to the computational efficiency of the program, which may have led to errors caused by a small number of samples. Although there are no full-scale measured results to compare with the simulated results, the ice resistance of M/V Tianen based on measurements was estimated to verify the results [25]. The probability density curve shows a positively skewed bell shape, which matches with the probability density curves of both results. With the ice parameters of one area, the ice resistance of the area may be obtained and the ship safety operation could be achieved.
Sensitivity analysis of ice parameters conducted indicates that ice resistance shows a rising trend with the increase in ice thickness and flexural strength in general. Ice resistance decreases with the decrease in crushing strength and ice density, while ice resistance shows an unclear and low-changing range trend under the influence of effective modulus. Changing in ice thickness and flexural strength will affect ice resistance significantly.
Crushing strength and effective modulus shows a small impact on ice resistance. Moreover, the rain flow counting method is adopted to obtain the amplitude of ice force, which provides more detailed information to study the actions of ice force on structures than the Monte Carlo method.

Conclusions
Ice parameters following certain distributions are generated. A numerical simulation method based on circumferential crack assumption is then applied to calculate the ice load. The ice-failure pattern is refined into bending and crushing effects; then, the trigger conditions of failure patterns are introduced. Simulated results show more detailed information than the Monte Carlo method; ice force at each time step is obtained, key ice force amplitude is extracted, and significant differences in sea-ice break cracks are also observed.
The main findings based on numerical simulations and comparisons with another researcher are summarized as follows: (1) The change in ice thickness and flexural strength will affect ice resistance significantly, while effective modulus has the least effect on ice resistance among all the parameters; (2) Ice resistance obtained by pseudo-random numbers tends to be concentrated in the area of low resistance, around 0.4-0.8 MN. The probability decreases significantly when the resistance exceeds about 1.5 MN; (3) Ice-force amplitude is concentrated at 0-1 MN. Only 2.3% of ice-force amplitudes exceed 2 MN. The maximum ice-force amplitude is 8.845 MN. The results of ice amplitude could be a useful reference for ship design and structural analysis, and potential stress generated by sea ice could be predicted at an early stage; (4) The distribution of ice resistance follows the Burr distribution in the simulation. The probability density curve shows the same positively skewed bell shape with results obtained by the Monte Carlo method, a trend matched well with measurements of another ship.
In the future, the combination of this study in ice resistance and real-time image processing may provide more detailed ice data for simulation in field operation, continuous accurate ice-force predictions could be achieved, and serious damage to ships may be avoided.
It should be noted that some assumptions made in the empirical formula, or errors occurred in the generation of ice parameters, will affect the simulation results. An insufficient number of samples also may lead to certain deviations.