A Numerical Prediction of the Resistance of Bulk Carriers in Brash Ice Channels

Abstract

Ship resistance increases significantly when navigating a brash ice channel. In this study, the numerical method is applied to predict the full-scale ship resistance of bulk carriers in brash ice channels. The viscous flow computational fluid dynamics (CFD) solver was coupled with the discrete element method (DEM) to establish the brash ice model. The Euler multiphase flow’s volume of fluid (VOF) model was applied to simulate the interaction between the ship and water. The ship–brash ice interaction was simulated. Predictions of ships’ total resistance based on the numerical method and the Finnish Swedish ice class rules (FSICR) method were compared with the experimental results carried out in Hamburg Ship Model Basin (HSVA) ice tank. The numerical resistance shows a good agreement with the HSVA experiment reports and a better performance than the FSICR method. The present study shows that the numerical method could provide reasonable and practical ice resistance predictions for engineering applications.

1. Introduction

The distance between Yokohama and Hamburg through Suez is 18,350 km, while the Arctic route is only 11,100 km. Navigating through the Arctic route makes it possible to reduce the travel time from 22 days to 15 days [1,2]. Therefore, navigating the polar route could dramatically reduce transportation costs and improve efficiency. However, in such polar routes, the ships will face a more hostile and complex environment than in general due to the presence of sea ice. Due to a lack of autonomous ice-breaking capability, the commercial vessels usually sail through brash ice channels guided by icebreakers to ensure their safety [3]. The Finnish-Swedish ice class rules (FSICR) [4] are a accepted rules to ensure that ships can navigate in polar regions. According to FSICR, ships shall maintain a minimum speed of 5 knots to ensure smooth traffic progress in the brash ice channel. It is reported that ice-induced resistance can account for more than 80% of the total resistance when sailing in a brash ice channel at 5 knots [5]. Therefore, predicting the ice resistance of a ship in a brash ice channel is very important for the ship design process.

Several theoretical analyses and empirical formulas have been developed to predict ice-induced resistance. Lindqvist [6] categorized ice resistance into bending, crushing, and submersion components and presented an empirical formula. Based on a series of measured data in the Baltic Sea, Riska [7] presented a formula by improving the Lindqvist formulas [6]. The FSICR method, modified from the Riska formula, has evolved progressively over the years [4]. It has become one of the most common ways to predict the ice resistance of a ship. Dobrodeev and Sazonov [8] developed a method for calculating the ice resistance of ships sailing in a brash ice channel. The theoretical methods have remarkable advantages in fast prediction. However, the empirical methods usually provide too conservative predictions [9], which may lead to the delivered power’s failure to meet the EEDI Regulation. Moreover, both theoretical and empirical methods require a large amount of measured data. The performance of such methods varies for different hull shapes [10,11,12].

The model test in ice tanks is an efficient way to predict ship ice resistance due to the merit of providing reliable data [10,12]. However, model tests also have some limitations. For example, the ice models may be replaced by paraffin wax-made or polypropylene material limited by the test conditions [13,14]. In addition, model tests are only affordable for a few researchers [15].

The rapid development of computer technologies makes it possible to predict ice-induced resistance through numerical methods. The discrete element method (DEM) is one of the most used numerical methods to solve ice–ship interaction [16], and it was first developed by Cundall and Strack [17]. This method treats the discrete medium as a group of discrete elements. A two-dimensional DEM method was used to evaluate the impact of level ice on objects by Selvaderai and Sepehr [18]. The influence of current velocity, current direction, ice thickness, and ice floe size on ship resistance based on the discrete element method was investigated [19]. Following Mucha [20], the CFD-DEM coupling method was employed to analyze the ships moving in brash ice, and the sensitivity to the material properties and degree of coupling were discussed. The interaction of hull, ice, and water on a bulk carrier was studied by using the CFD-DEM coupling method by Luo et al. [5], and they analyzed the effects of contact models and brash ice particle shape on the ship’s ice resistance. The numerical results also vary with the numerical setups, such as the numerical mesh and parameters setup in contact models [21]. The object of the previous studies is only one vessel, and there is a lack of knowledge of the performance of the numerical setups on other vessels.

Therefore, the understanding of the prediction of ice-induced ship resistance through the numerical method still needs to be improved. The present work systematically studies the CFD-DEM method’s performance in predicting bulk carriers’ ship ice resistance in a brash ice channel. The DEM method establishes ship–ice interaction under the Lagrangian framework, while the interaction between ship and water is realized under the Eulerian framework. The ship–ice interaction process and mechanical behavior of five bulk carriers in the brash channel are simulated, considering different loading conditions and ice classes. The numerical results are compared with the experiments carried out in the Hamburg Ship Model Basin (HSVA) ice tank [22,23,24,25,26]. The calculated results show a good agreement with the experimental results. Additionally, the numerical resistance predictions provide better performance than the FSICR method.

2. Numerical Method

2.1. CFD Method

The present numerical simulations were all calculated by use of the Star-CCM+ software. The interaction between ship and water is realized under the Eulerian framework by the CFD method. This study uses the shear stress transport ($SST$) $k-\omega$ turbulence model with the two-layer all y+ method to do the numerical simulations. The VOF method with a high-resolution interface capture (HRIC) scheme [27] is used to capture the free surface. A wave-damping area is established at the inlet, the outlet, and the side boundary to avoid an unphysical reflection.

2.2. DEM Method

The DEM method establishes the interaction between ship and ice under the Lagrangian framework in the present numerical simulations. The ice particles in a channel can be considered as many discrete distribution units in numerical simulations [20,21]. The critical point of the DEM method is calculating the collision forces between the elements. The commonly used calculation models are the elastic contact model, elastic–plastic contact model, viscous-elastic contact model, etc. [28]. The elastic contact of discrete elements is assumed in this study, and the contact model between the cubic particle units is shown in Figure 1.

The contact force between the two cubic particles (particle i and particle j) can be described as follows:

$$F_{contact}=F_{nij}+F_{tij}$$

(1)

In Equation (1), $F_{nij}$ is the normal force and $F_{tij}$ is the tangential force. The normal force can be expressed as

$$F_{nij}=-K_n{d}_n-N_n{V}_n$$

(2)

where $K_n$ and $N_n$ are normal spring stiffness and normal damping, respectively. $V_n$ represents the velocity component in the normal direction at the contact point.

The tangential force is defined as

$$F_{tij}=\left\{\begin{array}{cc}-K_t{d}_t-N_t{V}_t\hfill & \left(|K_t{d}_t| < |K_n{d}_n|C_{fs}\right)\hfill \\ |K_n{d}_n|C_{fs}{d}_s/\left|d_s\right|\hfill & \left(|K_t{d}_t| > |K_n{d}_n|C_{fs}\right)\hfill \end{array}\right.$$

(3)

where $K_t$ is tangential spring stiffness, and $N_t$ is tangential damping; $C_{fs}$ represents the friction coefficient; $d_t$ is the overlap of the normal and tangent directions at the collision point. $V_t$ represents the velocity component in the tangential direction at the contact point.

Different contact models could describe the coefficients mentioned above [29,30]. Zhang et al. [21] investigated the effect of three contact models on contact force: the linear spring contact model, the Hertz–Mindlin contact model, and the Walton–Braun contact model. The results showed a minor difference in the simulation for ship resistance in brash ice of the three contact models. The linear spring contact model is applied in this study, and details of this model can be found in the guidelines of STAR-CCM+ [31].

2.3. Governing Equation of CFD Coupling DEM

Fluid–ice interaction is mainly related to the mechanical characteristics of the interaction between ship and ice. The motion of an incompressible Newton fluid satisfies continuity equation and conservation of momentum equations [32,33]:

$$\frac{\partial \left({\rho }_f{ϵ}_f\right)}{\partial t}+\nabla ·\left({\rho }_f{ϵ}_f\overrightarrow{u}\right)=0$$

(4)

$$\frac{\partial \left({\rho }_f{ϵ}_f\overrightarrow{u}\right)}{\partial t}+\nabla ·({\rho }_f{ϵ}_f\overrightarrow{u}·\overrightarrow{u})=-{ϵ}_f\nabla p-\nabla ·\left({ϵ}_f\overrightarrow{{\tau }_f}\right)+{\rho }_f{ϵ}_f\overrightarrow{g}-\overrightarrow{F}$$

(5)

where ${\rho }_f$ and ${ϵ}_f$ are the fluid’s density and volume fraction. The relationship is that ${ϵ}_f=1-{ϵ}_p$, ${ϵ}_p={\sum }_{n_p}^{i-1}{V}_{pi}/\Delta V$, where ${ϵ}_p$ is the discrete ice term control volume’s integral number; $V_{pi}$ and $\Delta V$ are the discrete ice particle i’s volume and all control volumes within the region, respectively. $\overrightarrow{u}$ is the fluid’s average velocity; p is the pressure’s mean value; $\overrightarrow{F}$ is the particle resistance based on the volume average to the nearby fluid of the discrete ice unit in the control volume, which includes drag, pressure gradient, shear stress, and other interactions between fluid and discrete terms. ${\tau }_f$ is the fluid’s stress tensor. ${\tau }_f$ and $\overrightarrow{F}$ could be stated as:

$$\overrightarrow{{\tau }_f}=-u_f(\nabla \overrightarrow{u}+{(\nabla \overrightarrow{u})}^{{}^{\prime }})+\frac{2}{3}{\mu }_f(\nabla ·\overrightarrow{u})\overrightarrow{\delta }$$

(6)

$$\overrightarrow{F}=\frac{1}{V_{cell}}\sum _{i=1}^{n_p}\overrightarrow{F_i}=\frac{1}{V_{cell}}\sum _{i=1}^{n_p}(\overrightarrow{F_{di}}+\overrightarrow{F_{bi}}+\overrightarrow{F_{li}}+\overrightarrow{F_{ai}}+\overrightarrow{F_{pi}})$$

(7)

where $n_p$ is the discrete particle unit number in the control volume of each fluid; ${\mu }_f$ is the fluid’s dynamic viscous coefficient; $\overrightarrow{\delta }$ is the tensor of each unit; $\overrightarrow{F_{ai}}$,$\overrightarrow{F_{bi}}$, $\overrightarrow{F_{di}}$, $\overrightarrow{F_{li}}$, and $\overrightarrow{F_{pi}}$ are additional mass force, buoyancy force, drag force, lift force, and pressure-gradient force acting on the fluid term of ice particle i in the discrete term, respectively. More details can be found in the STAR-CCM+ user guide [31].

2.4. Modeling of Brash Ice Particles

In HSVA ice tests, the brash ice particles are made manually by breaking up the parent-level ice. Figure 2 shows the picture of brash ice particles in the ice tank of HSVA [23], which are in different shapes. Therefore, it is challenging to reproduce each ice piece’s shape, location, and distribution characteristics in numerical simulations [34]. The shape of ice particles has a remarkable influence on the numerical results. Following recent research [5,21], an irregular polyhedral particle is used in the present study; see Figure 2.

Despite the geometric shape of the brash ice particle, it is also necessary to determine its mechanical characteristics. In the present simulations, the ice density is 916.72 kg/m${}^3$, Young’s modulus is 9 GPa, and Poisson’s ratio is 0.3 [35,36,37]. The friction coefficient between hull and ice is 0.1, the same as the HSVA experiments [26]. The brash ice particle number ranges from 16,000 to 33,000 in the present simulations.

2.5. Computational Domain and Ship Model

The geometric model of full-scale ships is established. According to the experiment test conditions and previous research [5,21,26], the present numerical domain is set based on the ship’s length. It ranges from $-2L_{pp}$ to $6L_{pp}$ in axial direction. The range of spanwise direction is $-L_{pp}\le y\le L_{pp}$, and the range of vertical direction is $-2L_{pp}\le z\le L_{pp}$. The width of the brash ice channel is twice that of the ship breadth following the FSICR (Trafi, 2011). The simulation of the level ice zone on each side of the channel is realized by establishing wall areas, consistent with the brash ice model test in HSVA. The relative coordinate system is applied in this study. The ship is maintained still, while the water and brash ice particles move toward the ship at 5 knots. The computational domain is demonstrated in Figure 3.

This study selects five bulk carriers, namely 38K BC, 68K BC, 82K BC, 95K BC, and 108K BC.

Table 1. Characteristics of flow and sediment.

Ship ID	${\mathit{L}}_{\mathit{pp}}$ (m)	B (m)	${\mathit{TF}}_{\mathit{UIWL}}$ (m)	${\mathit{TA}}_{\mathit{UIWL}}$ (m)	${\mathit{TF}}_{\mathit{LIWL}}$ (m)	${\mathit{TA}}_{\mathit{LIWL}}$ (m)
38K BC	177.0	32.0	10.7	10.7	5.3	7.4
68K BC	223.4	32.3	13.3	13.3	5.1	8.3
82K BC	225.5	32.3	15.5	15.3	4.9	7.4
95K BC	225.5	38.0	15.0	15.0	5.1	8.3
108K BC	241.8	43.0	14.5	14.5	6.6	9.3

lists these ships’ main particulars. In Table 1, $L_{pp}$ is the length between perpendiculars, B is the breadth, $TF$ is the fore draught, and $TA$ is the aft draught. The upper ice waterline (UIWL) is the load draught condition, and lower ice waterline (LIWL) is the ballast draught condition.

2.6. Numerical Mesh Setup

In the present simulations, the numerical setup is all kept the same, including the strategy of the numerical mesh setup. The difference among each simulation of different vessels is the base size of the numerical mesh, which depends on the length of each vessel.

Firstly, the grid sensitivity analysis was performed to avoid the influence of the grid. The simulations were carried out on three grids, each with a different density, while the other settings are all fixed. The numerical mesh of 108K BC based on the medium size is shown in Figure 4.

It should be noted that even though the International Towing Tank Conference (ITTC) proposed the procedures and guidelines for ice testing, there has yet to be an individual proposal for the corresponding CFD simulations [21]. The ship-ice interaction simulations in this study were accomplished based on the general CFD Verification and Validation Procedures and Guidelines [38]. The numerical results of 108K BC based on different density grids are shown in

Table 2. Grid convergence results.

Grid Density	Grid Number	Results—Num	Convergence Ratio
Corase	1.41 million	1531.40 kN	0.29
Medium	2.05 million	1458.89 kN
Fine	2.97 million	1437.91 kN

.

The convergence ratio could be calculated by

$$\mathrm{Convergence} \mathrm{Ratio}=\frac{{\mathrm{Results}}_{\mathrm{Fine}}-{\mathrm{Results}}_{\mathrm{Medium}}}{{\mathrm{Results}}_{\mathrm{Medium}}-{\mathrm{Results}}_{\mathrm{Corase}}}$$

(8)

The estimated convergence ratio is 0.29, which indicates that the total resistance obtained by CFD simulation shows a monotonic convergence in this study. The numerical results in this study are all based on medium-density grids. The present study performed all the simulations based on a computer with 48 cores and a maximum frequency of 2.4 GHz, and the simulation time ranged from 35 to 40 h for each case.

3. Results and Discussions

3.1. Ship–Ice Interaction

Figure 5 shows the brash ice distribution of 108K BC under UIWL conditions with the 1A ice class at different times. The time–history curve of total resistance, combining the ship–water and ship–ice contact resistance, is shown in Figure 6. In the beginning, the brash ice particles move forward and begin to contact the ship bow at 16 s. The interaction between the ship fore part and the ice is the main factor contributing to the ice resistance before 28 s; see Figure 5a,b. In Figure 5c,d, the brash ice particles are constantly being pushed over the ship’s shoulders and sides. The whole hull is gradually immersed in the brash ice particles, including the fore part, middle body, and ship stern. At this stage (28 s–180 s), the number of ice particles colliding with the ship, as well as the surface area of the hull in contact with the ice, is increasing. Consequently, the total resistance shows a fluctuating upward trend, as shown by the blue line in Figure 6.

With the brash ice particles passing through the stern and moving downstream to the outlet, it turns to the second stage. The total resistance value gradually reaches equilibrium (red line in Figure 6). Unlike the calm water resistance, which can be calculated to a very stable value, the resistance in ice varies in a larger border range. Because of the limited time history of the total resistance data from the experiments, the numerical results are verified through the averaged total resistance. The averaged total numerical resistance (black dashed line in Figure 6) is the mean value between 180 s–350 s.

Figure 7 and Figure 8 show the 108K BC and the brash ice interaction from the bow and stern views under UIWL and LIWL conditions, respectively. Under UIWL conditions, the ice particles are pushed by the ship to the sides of the channel when the ship is navigating the channel. The numerical results are in good agreement with the model test results. Unlike the results under UIWL condition, the brash ice also passes through underneath the ship under LIWL condition. It should be pointed out that the other selected vessels also show similar behavior. This phenomenon is well represented by the numerical method in Figure 8. As mentioned, the numerical method performs well in simulating ship–ice interaction while sailing in the brash ice channel.

Figure 9 shows the experimental and numerical results of brash ice distribution at the end of the test for 108K BC, respectively. As mentioned above, the hull pushed the brash ice particles toward the sides of the channel during the ship’s navigation process. However, from Figure 8a, it is clear that the brash ice particles are refilled in the whole channel after the test run. In Figure 8b, ice particles can be found regrouped away from the ship. The numerical results are in good agreement with the experimental results.

3.2. Ice Resistance

There are three common ways to predict the ship’s ice resistance: the model test method, theoretical method (empirical formulas), and numerical method. The predictions of resistance based on each method are listed in the following section.

3.2.1. Experimental Predictions

In HSVA ice experiment reports, the predictions of full-scale ice resistance are derived from the measured model-scale results. This approach should consider several factors. For example, the measured model ice thickness usually differs from the target value. Due to the difference, the model experiences a higher or lower resistance. As a result, it is necessary to modify the actual measured results to obtain the target values. The scale effect should also be taken into account. Some correction methods have been adopted in the HSVA report [22,23,24,25,26] for such purposes. Even though the correction methods bring some uncertainties to the full-scale predictions, the model test method could provide the most reliable predictions among the three methods. The HSVA experimental results of full-scale resistance predictions are shown in

Table 3. Full-scale resistance predictions reported by HSVA experiments.

Ship ID	Scale	Loading Condition	Ice Class	Brash Ice Thickness (m)	Results—HSVA (kN)
38K BC	26.602	UIWL	1A	1.448	772.48
		UIWL	1B	1.248	648.79
		LIWL	1A	1.448	768.94
		LIWL	1B	1.248	652.30
68K BC	29.508	UIWL	1A	1.452	1131.60
		UIWL	1B	1.251	955.40
		LIWL	1A	1.452	1071.00
		LIWL	1B	1.251	896.70
82K BC	29.333	UIWL	1A	1.452	1008.70
		UIWL	1B	1.252	857.99
		LIWL	1A	1.452	1058.24
		LIWL	1B	1.252	896.27
95K BC	33.206	UIWL	1A	1.532	1682.70
		UIWL	1B	1.332	1442.90
		LIWL	1A	1.532	1804.30
		LIWL	1B	1.332	1525.40
108K BC	33.634	UIWL	1A	1.602	1562.79
		UIWL	1B	1.402	1357.51
		LIWL	1A	1.602	1471.09
		LIWL	1B	1.402	1261.44

. The target brash ice thickness of each case is also included in Table 3.

3.2.2. Theoretical Predictions

The most common way to predict ship ice-induced resistance in practical engineering is with theoretical methods. As the FSICR method [4] is one of the most widely used theoretical methods, the total resistance predicted by FSICR methods is presented in

Table 4. Full-scale resistance predictions based on HSVA experiments and FSICR.

Ship ID	Loading Condition	Ice Class	Results—HSVA (kN)	Results—FSICR (kN)	APE (%)
38K BC	UIWL	1A	772.48	1180.70	52.8
	UIWL	1B	648.79	925.20	42.6
	LIWL	1A	768.94	1132.60	47.3
	LIWL	1B	652.30	891.20	36.6
68K BC	UIWL	1A	1131.60	1357.10	19.9
	UIWL	1B	955.40	1081.10	13.2
	LIWL	1A	1071.00	1275.20	19.1
	LIWL	1B	896.70	997.20	11.2
82K BC	UIWL	1A	1008.70	1490.00	47.7
	UIWL	1B	867.99	1204.20	40.4
	LIWL	1A	1058.24	1288.40	21.7
	LIWL	1B	896.27	1043.50	16.4
95K BC	UIWL	1A	1682.70	1897.00	12.7
	UIWL	1B	1442.90	1520.80	5.4
	LIWL	1A	1804.30	1617.10	6.4
	LIWL	1B	1525.40	1274.60	16.4
108K BC	UIWL	1A	1562.79	2901.80	85.7
	UIWL	1B	1357.51	1903.80	40.2
	LIWL	1A	1471.09	2166.60	47.3
	LIWL	1B	1261.44	1707.60	35.4

. The experimental predictions and the absolute percentage error (APE) between the two results are also shown in this table. The APE is calculated by the following expression:

$$\mathrm{APE}=|\frac{{\mathrm{Results}}_{\mathrm{FSICR}}-{\mathrm{Results}}_{\mathrm{HSVA}}}{{\mathrm{Results}}_{\mathrm{HSVA}}}| \times 100\%$$

(9)

Several statistical indices are introduced to evaluate the performance of the FSICR method, including mean absolute percentage error (MAPE), root-mean-square error (RMSE), and discrepancy ratio (DR). The following expression calculates these statistical indices:

$$\mathrm{MAPE}=\frac{100\%}{\mathrm{n}}\sum _{\mathrm{i}=1}^{\mathrm{n}}\mid \frac{{\mathrm{Results}}_{\mathrm{FSICR}}-{\mathrm{Results}}_{\mathrm{HSVA}}}{{\mathrm{Results}}_{\mathrm{HSVA}}}\mid$$

(10)

$$\mathrm{RMSE}=\sqrt{\frac{1}{\mathrm{n}}\sum _{\mathrm{i}=1}^{\mathrm{n}}{({\mathrm{Results}}_{\mathrm{FSICR}}-{\mathrm{Results}}_{\mathrm{HSVA}})}^2}$$

(11)

$$\mathrm{DR}=\frac{1}{\mathrm{n}}\sum _{\mathrm{i}=1}^{\mathrm{n}}\frac{{\mathrm{Results}}_{\mathrm{FSICR}}}{{\mathrm{Results}}_{\mathrm{HSVA}}}$$

(12)

where n is the test number.

For MAPE and RMSE, zero shows perfect accuracy, and a significant positive value indicates poor prediction. A DR of 1 means the model is unbiased, a DR < 1 means that the predicted value is smaller than the actual value, and a DR > 1 means that the predicted value is larger than the actual value. As the predictions with an APE of less than 20% are regarded as reasonable and practical values for engineering, comparisons are also made by listing the number of samples within the ±20% error lines. The values of such statistical indices are summarized in

Table 5. Statistical indices based on results of FSICR and HSVA experimental predictions.

Ship ID	MAPE	RMSE	DR	within ±20%
All	31.1%	441.6 kN	1.28	45%
38K BC	44.8%	328.8 kN	1.45	0%
68K BC	15.8%	172.1 kN	1.16	100%
82K BC	31.6%	326.4 kN	1.32	25%
95K BC	11.2%	193.6 kN	0.98	100%
108K BC	52.1%	832.8 kN	1.52	0%

. Scatter plots of full-scale total resistance between experimental and FSICR results are shown in Figure 10. The ±20% error lines are also illustrated in this figure.

Results in Table 5 and Figure 10 show that the FSICR method provides conservative results, with DR = 1.28 and 90% of the predictions larger than the experimental data. This behavior is favorable for engineering applications due to its high safety margins. However, the differences between the FSICR predictions and the experimental reports are significant. The maximum APE reaches 85.7%. The MAPE of all the results is 31.1%, which is much more significant than 20%. Regarding the MAPE for each report, the results of only two reports (68K BC and 95K BC) are less than ±20%. From Figure 10, all the predictions of 38K BC and 108K BC are beyond the $\pm 20\%$ error lines. Furthermore, over 50% of FSICR predictions exceed ±20% error lines. Even though the FSICR method has advantages in efficiency and safety margins, it could not be satisfactory in precision.

3.2.3. Numerical Predictions

Both the numerical and experimental full-scale total resistance of each vessel at different loading conditions and ice classes are summarized in

Table 6. Full-scale resistance predictions based on HSVA experiments and numerical method.

Ship ID	Loading Condition	Ice Class	Results-HSVA (kN)	Results-Num (kN)	APE $(\%)$
38K BC	UIWL	1A	772.48	990.96	28.3
	UIWL	1B	648.79	663.38	2.2
	LIWL	1A	768.94	1051.99	36.8
	LIWL	1B	652.30	727.12	11.5
68K BC	UIWL	1A	1131.60	1164.31	2.9
	UIWL	1B	955.40	835.88	12.5
	LIWL	1A	1071.00	1131.92	5.7
	LIWL	1B	896.70	831.19	7.3
82K BC	UIWL	1A	1008.70	1137.21	12.7
	UIWL	1B	867.99	824.01	4.0
	LIWL	1A	1058.24	1344.26	27.0
	LIWL	1B	896.27	966.06	7.8
95K BC	UIWL	1A	1682.70	1598.02	5.0
	UIWL	1B	1442.90	1023.84	29.0
	LIWL	1A	1804.30	1517.18	15.9
	LIWL	1B	1525.40	1085.45	28.8
108K BC	UIWL	1A	1562.79	1458.89	6.6
	UIWL	1B	1357.51	1166.46	14.1
	LIWL	1A	1471.09	1647.38	12.0
	LIWL	1B	1261.44	1349.44	7.0

and Figure 11. Table 6 also includes the APE between numerical and experimental results. The values of statistical indices between numerical and experimental results are illustrated in

Table 7. Statistical indices based on results of numerical and HSVA experimental predictions.

Ship ID	MAPE	RMSE	DR	within ±20%
All	13.9%	200.8 kN	1.02	75%
38K BC	19.7%	182.8 kN	1.20	50%
68K BC	7.1%	76.4 kN	0.97	100%
82K BC	12.9%	161.5 kN	1.11	75%
95K BC	19.7%	338.7 kN	0.80	50%
108K BC	9.9%	146.7 kN	1.00	100%

.

The computed results agree well with the experimental data. Most APEs (75%) between numerical and experimental results are less than 20%. The minimum APE is only 2.2%, while the maximum APE is 36.8%. The MAPE of all the results is 13.9%, smaller than 20%, and much smaller than the FSICR method. Additionally, the MAPE of each report is less than 20%, with the maximum value being 19.7%. The RMSE of the numerical results is 200.8 kN, which is also much less than that of the FSICR method. In addition, regarding both MAPE and RMSE, the numerical method shows a much better performance on the resistance prediction of each vessel except 95K BC.

Similar to the FSICR method, the numerical simulations also provide conservative results with DR = 1.02. More than 50% of the predictions are larger than the experimental results. However, it should be noted that the highest underprediction of numerical results reaches 29.0%, more significant than the FSICR method.

The performance of the numerical method varies on individual vessels. The numerical predictions of 68K BC show the best agreement with the HSVA report, with the MAPE being 6.9% and RMSE being 80.7kN. The numerical results of 95K BC have the most significant deviation from the experimental results, with the MAPE being 19.7% and RMSE 338.7 kN. In addition, the numerical method shows an overestimation behavior on 38K BC, while all the numerical predictions of 95K BC are smaller than the value of experimental results.

4. Conclusions

This study analyzes the CFD-DEM method performance on the full-scale resistance prediction of bulk carriers in brash ice channels. The CFD solvers realized ship–water interaction, and numerical brash ice particles were built using the DEM method. The numerical results were compared to those derived from the HSVA experiments and FSICR method. The present study can draw the following conclusions:

(1) The interaction between ship and brash ice particles was numerically simulated. The calculated results are in good agreement with the observations of the HSVA model test. The simulated phenomena show that the contact force caused by the interaction between the parallel middle body and brash ice is the primary source of the total resistance.

(2) The FSICR was applied to calculate the ice resistance in this work. The FSICR method provides conservative results. However, the FSICR predictions derive dramatically from the HSVA experimental results, with the MAPE being 31.1% and the maximum APE being 85.7%. Additionally, 55% of predictions are beyond the $\pm 20\%$ error lines.

(3) The computed results agree well with the experimental results, and the numerical method performs much better than the FSICR method. The MAPE of total resistance between numerical and experimental results is 13.9%, while the maximum APE is 36.8%. 75% of numerical predictions are within ±20% error lines.

(4) The numerical method shows various performances in predicting the ice resistance of different ships. The numerical method provides conservative results on predicting resistance of 38K BC, while it shows an underestimated behavior on 95K BC.