A combined experimental and numerical approach to predict ship resistance and power demand in broken ice

Despite its remoteness and hostile environmental conditions, the Arctic holds significant shipping lanes, such as the Northern Sea Route (NSR) and the Northwest Passage (NWP). Typically, merchant ships operate along these routes in summer only, when the dominating type of ice is broken ice. A challenge of operating in such ice conditions is that there is no cost-and time-efficient method for predicting the resulting ice resistance


Introduction
The Arctic holds many promising shipping lanes, including the Northern Sea Route (NSR) and the Northwest Passage (NWP) (Bergström et al., 2020;Li et al., 2020).Nevertheless, the volume of Arctic shipping has remained relatively low, in part due to various associated technical challenges (Bergström et al., 2022;Li et al., 2023).Typically, merchant ships operate in the Arctic in summer only, so that they stay within the so-called marginal ice zone, characterized by broken ice consisting of ice floes with diameters of a few meters (Williams et al., 2013).One of the challenges associated with operating in such ice conditions is the prediction of the resulting ice resistance, making route planning challenging, among others (Sazonov, 2011).
Available theoretical approaches for estimating hull ice resistance are generally based on the assumption that ice resistance is the sum of ice forces acting against the movement of the ship.Aboulazm (1989) proposed an analytical approach in two dimensions for a ship navigating through a broken ice field, consisting of a "micro model" and a "macro model".The "micro model", which is intended for application on ships sailing at medium speed in a broken ice field with a relatively low ice concentration, estimates the ice resistance considering the collision force between the hull and each individual ice floe.The "macro model", in turn, which is intended for application on ships sailing in a broken ice field with a high ice concentration at a relatively low speed, estimates the ice resistance in terms of the work done by the broken ice on the path of removal.It's worth noting that, because both of its models are two dimensions, the overturning motion of the broken ice is not considered.Idrissova et al. (2019) used a collision-energy-based method to predict ice loads on ships, which agrees well with full-scale measurements, but the implementation relies on the use of multiple assumed coefficients concerning the ship-ice interaction process.Using a theoretical approach Zong and Zhou (2019) considered part of fluid effects, such as the added mass of ice moving in a constant flow field, and verified the relationship between vessel speed, ice concentration, and ice resistance based on experimental results.They found that fluid effects play an important role in estimating hull ice resistance in a broken ice field.However, their approach is based on thin-ship and small ice floe hypotheses and it remains unclear if the approach is valid in case those assumptions are not met.Based on the above examples, we conclude available theoretical approaches are typically developed for assessing specific types of ship-ice interactions, and that they are often dependent on assumptions as well as multiple empirical coefficients, which in principle makes them semi-empirical.Combined, these features significantly limit their practical use.
Semi-empirical formulas provide a convenient and fast way to estimate ice resistance.In general, these are based on a combination of fullscale measurements, model tests, and theoretical studies (Nogid, 1959;Keinonen et al., 1998).Well-established semi-empirical formulas for predicting level ice resistance are presented by Riska et al. (1998), Enkvist (1972), and Lindquist (1989).In terms of ice resistance in broken ice fields, the formulas consider two main categories of influencing factors.The first category consists of parameters describing the hull, such as width and prismatic coefficient.The second category consists of parameters describing the operating condition, such as the Froude number, ice thickness, as well as floe size and ice concentration (Dubrovin, 1970;Aboulazm, 1989;Colbourne, 2000).In practice, unlike the brash ice channel or broken ice field, navigating in the newly created broken ice channel is becoming a common ice condition for ships following an icebreaker.A semi-empirical equation for estimating hull resistance in a broken ice channel was proposed (Sazonov and Dobrodeev, 2021).Their formula takes the effect of ice channel width into account, which is also an attempt to consider the factor of specialization.It has been shown that fluid effects play a key role in the prediction of broken ice resistance (Tsarau et al., 2014) and a semi-empirical formula to consider such effects is presented by Zong and Zhou (2019), Huang et al. (2021).A limitation of these formulas is that they rely on empirical coefficients derived from existing data and that they can only take into account a limited number of main influencing parameters.However, verification is still required for special applications, such as new ship types or special ice conditions.
A well-proven, albeit costly and time-consuming method for assessing ice resistance is to conduct ice-tank model tests, which allow for detailed consideration of a ship's hull parameters (Hellmann et al., 2005;Huang et al., 2018;Zhou et al., 2019).To enable less costly tests, attempts have been made to use synthetic non-refrigerated materials such as wax, polypropylene, and other mixed materials as substitutes for refrigerated ice in conventional towing tanks (Polojärvi et al., 2012;Zong et al., 2020).These approaches have been validated through comparisons with ice tank tests (Kim et al., 2013) and calculations based on empirical formulas (Van der Werff et al., 2015).The results of experiments utilizing synthetic ice have also served as references for the development of empirical formulas for ship resistance (Huang et al., 2021;Yang et al., 2021).
In recent years, numerical simulations have emerged as an alternative way of predicting ship ice resistance (Lau et al., 2011;Suominen et al., 2013;Ni et al., 2020a;Zheng et al., 2022).To predict ship resistance in broken ice fields, the Discrete Element Method (DEM) offers advantages by enabling the calculation of combined forces arising from simultaneous contacts between numerous elements, including ice-ice and hull-ice collisions (Hopkins and Tuhkuri, 1999;Hansen and Løset, 1999).Ji et al. ( 2013) applied a method based on three-dimensional disk elements to calculate the hull resistance of an icebreaker operating in a broken ice field.To simulate ice fields with a reduced number of cells, researchers have developed methods to study the different shapes of cells that constitute the ice field, such as the non-smooth discrete element method (NDEM) (Van den Berg et al., 2019) and the Dilated-polyhedron-based DEM (Liu and Ji, 2021).Recently, Polojärvi et al. ( 2021) compared full-scale measurements of ice loads induced by floe fields with corresponding loads calculated by DEM simulations.In general, DEM simulations simplify fluid effects.However, there are methods to consider fluid variations by employing the CFD-DEM (computational fluid dynamics and discrete element method) coupling method.Huang et al. (2020) demonstrated the significant effect of waves created by a ship's bow, which push ice floes away from the hull, on the resulting ship resistance.Luo et al. (2020) used the commercial software Star-CCM + to study ship resistance in a brash ice channel, modelling the accumulation of small ice pieces in a shipping channel.They investigated both one-way and two-way coupling schemes, which produced similar results, with the former requiring fewer computational resources.Furthermore, Ni et al. (2020b) examined ship resistance in level ice using CFD-DEM.
In the design of some new ships or exploration of special ice conditions, where existing empirical parameters may not meet the needs.Since numerical methods can be applied to a wide range of cases and require less time and resources than model testing, we believe that a well-validated numerical method could be used as a supplement to model testing to enable the consideration of multiple different ice conditions and operating scenarios without significant time and resource expenditures.However, as a validation, the cost of model test in ice-tank is consuming, and the use of a resource-efficient modeling test method can save time and costs.To explore this idea, in this study, we employ a combined numerical and experimental approach to predict ship resistance in broken ice.The aim is to find a way to predict a ship's resistance and power demand in broken with agreeable accuracy in a time-and cost-efficient manner, at the same time, details of individual differences can be taken into account as fully as possible.The proposed combined approach is composed of both experimental and numerical ones.The experimental tool makes use of a type of non-refrigerated synthetic model ice made of polypropylene, which makes it possible to conduct ice model tests in a conventional non-refrigerated towing tank rather than an ice tank.The numerical tool, in turn, is based on the state-of-the-art CFD-DEM method that accounts for fluid effect while the ship is moving ahead with affordable computation efforts.Besides, validation calculations against established empirical approaches are presented to demonstrate the accuracy of the proposed method.
The rest of this article is organized as follows: in Section 2, the experimental setup and the case study vessel are introduced.It also describes in detail the characteristics of the non-refrigerated model ice in use.Section 3 presents the theoretical frame of the coupled CFD-DEM method.Also, the setup of the numerical simulations is described.In Section 4, the empirical formulas for broken ice resistance are introduced and the three established formulas used in this study are described in detail.Section 5 presents the results of the model tests and the numerical simulations, which are compared with the validation calculations from the empirical formulas.The conclusion and discussion are presented in Section 6.

Experimental setup
A series of towing tank model tests were carried out using a type of non-refrigerated synthetic model ice made of polypropylene.The experiments focus on the following two factors affecting the motion characteristics of a ship moving in a broken ice field: (a) the motion of the broken ice in the vicinity of the hull and (b) the total resistance of the ship.
The model tests were conducted in the outdoor towing tank at Harbin Engineering University (HEU), the principal dimensions of which are 20 m × 2 m × 1.2 m.As shown in Fig. 1, the tank is divided into three sections: (a) a 5 m long acceleration section with open water; (b) a 10 m long experimental section with synthetic broken ice, and (c) a deceleration section with open water.The ship model is towed by a moving crane whose speed can be adjusted between 0 and 2 m/s.
The ship model is a scaled model of an ice breaker with an overall length of 122.5 m.The scale λ is 60.0, meaning that the overall length of the model ship is 2.04 m.The main parameters of the ship model are presented in Table 1 and a photo of the model is presented in Fig. 2.
The speed of the ship model is determined by the Froude similarity criterion: where V and L WL represent the ship speed and length of the waterline, respectively, and g = 9.81 m/s 2 is the gravitational acceleration.The specific speeds used in the experiments were 0.5, 0.6, 0.7, and 0.8 m/s.In the model tests, the secondary break-up of the fragmented ice was not considered.It is assumed that the ice particles are not subject to bending failure.The inertial forces of the ice pieces and the surface friction in the ship-ice interaction are therefore the key factors to the ice resistance.Polypropylene is selected as the material for the synthetic ice, following Zong et al. (2020).The density of the synthetic ice is approximately 917.0 kg/m 3 .The friction coefficient between the wet hull and the synthetic ice was measured to be 0.138.It is noticeable that the friction coefficient between two synthetic ice surfaces is 0.29.Both of these values are close to the properties of actual sea ice.The kinetic friction coefficient between a ship's hull and sea ice is 0.14 to 0.04 because of the hull's smoothness and there is also a variation between 25 and 30% regardless of whether the hull is smooth like it is printed or only corroded steel (Timco and Weeks, 2010).Therefore, this friction coefficient from the model test is also applied to numerical simulation.
The synthetic ice pieces used in this study are square-shaped with a side length (l i ) of 0.067 m and thickness (H) of 0.01497 m, which corresponds to 0.9 m in thickness and 4.0m in side length in the full scale.Squared ice pieces of similar size and thickness were reported from Liu (2018), andCho et al. (2013), respectively.It's noted that this is a simplification of real conditions of broken ice fields which are composed of ice pieces of various sizes and random shapes following many model tests as a verification (Van den Berg et al., 2018).However, the ice shape can affect the mean as well as the standard deviation of ice resistance without considering the fracture of ice pieces, because the parallel opposite edges of model ice can lead to the "force chain" earlier than who has a smaller "roundness" (Van den Berg et al., 2019).In order to keep the test comparable, the shape of the broken ice in this paper is consistent with the test in the numerical simulations.The dimensions of the ice pieces should be equal to or below a specific critical length l defined by Lu et al. (2016) as shown in Eq. ( 2) to make sure that the ice pieces do not break during the test: where k is the elastic coefficient and k = ρ water g, where ρ water is the density of water.W is the flexural rigidity of the sea ice given by where E and υ are Young's modulus and Poisson ratio of sea ice, respectively, E = 1.0 GPa and υ = 0.3, and H is the thickness of the broken ice.Accordingly, considering Froude and Cauchy similarity criteria (Huang et al., 2018), l is defined as 0.1497m as per Eq.(2).Another important parameter is the ice concentration, which is defined as the percentage of the total water area covered by ice.In the real world, the ice concentration varies between seasons and regions.To assess the influence of such variations, we consider three different ice concentration (η) values, namely 60%, 70%, and 80%.In the model

Table 1
The main parameters of the ship model.tests, broken ice fields of different concentrations were generated by randomly placing the corresponding number of synthetic ice blocks in the towing tank, as shown in Fig. 3.
During the experiments, several cameras were adopted to record the ship-ice interaction process from different angles.The connection arrangement between the model ship and the carriage is shown in Fig. 4. Fig. 4 (a) shows the overall overview of connections, there is a tow bar between the carriage and the model ship, which is fixed to constrain in the rest of the directions except for the heave motions, which are achieved by two sets of pulleys at the front and rear of the tow bar.Under the tow bar, as shown in Fig. 4(b), a force sensor is connected to the hull through a cross-axis.The hull can achieve pitch and roll motions through the cross-axis.Force sensor was used to measure the ship resistance for various ship speeds and ice concentrations.The capability of the force sensor is 1 kN and its combined error margin is ±0.023%.To account for the impact of random factors, three repeated tests were carried out for each of the three considered ice concentrations.

Numerical methods
The proposed CFD-DEM coupling method is implemented using the commercial software Star-CCM+ (Siemens PLM Software, 2019).The CFD-DEM coupling method is a kind of Euler-Lagrangian coupling approach in which the fluid phase is considered a continuum phase modeled using the Euler method, and ice particles are modeled as unbreakable discrete elements following the Lagrangian method.For solving the fluid phase, the finite volume method (FVM) is used.A detailed description of the numerical modeling process is presented in the following.

The governing equations of the fluid domain
In the CFD-DEM coupling scheme, the continuity and momentum conservation equations are satisfied in an incompressible fluid domain (Norouzi et al., 2016): where ρ f and ε f are the density and volume fraction of fluid, respectively, ρ water = 1025kg/m 3 and ρ air = 1.184kg/m 3 ; t is time, u f is the velocity of the fluid, p is pressure, is the fluid viscous stress tensor, where μ is the dynamic viscosity of the fluid, δ ij is the Kronecker delta symbol, and F OM represents the average volume force of fluid acting on the solid particle in a fluid cell.The standard K-Epsilon turbulent model, which has the advantage of computing the course grid condition (Siemens PLM Software, 2019), is also applied.
The volume of fraction (VOF) method is adopted to treat the free surface (Wang and Wan, 2020;Ghamari et al., 2022).A volume fraction ε f,i is defined to represent the volume ratio occupied by the fluid in a cell, where i = 1 and 2 denote water and air, respectively, so that ε f,1 + ε f,2 = 1 and 0 ≤ ε f,1 ≤ 1.The water-air interface is captured by solving the following equation (Hirt and Nichols, 1981):

The motion equations of ice particles
The ice particles floating in the water are modeled using DEM, whose motion is governed by the equilibrium equations of force and angular momentum (Siemens PLM Software, 2019): where m and I represent the mass and inertia moment of the ice particle, respectively; v p and ω p represent the velocity and angular velocity of the ice particle, respectively; F g represents the gravity force, F f and M f represent the fluid force and torque on the particle, respectively; which are discussed in Section 3.2.1.F c and M c represent the contact force and torque from the other particles or the hull, respectively, which are discussed in section 3.2.2.

Ice-fluid interaction force and torque
Considering the ice-fluid interaction, the fluid force acting on each ice particle is F f , which here includes drag force F d , pressure gradient force F p , and lift force F l .The fluid forces acting on the ice particles are calculated based on the fluid grid.Each ice particle is projected onto the fluid grid, while the surrounding fluid grid contains fluid forces that are interpolated onto the ice particles to contribute to the total force and torque (Huang et al., 2020).
The drag force F d and drag torque M d are expressed as: where A p represents the projected area of the particle, D represents the particle diameter, and when the DEM particle shape is a polyhedron, the size of that polyhedron is specified by equating the volume of that polyhedron to a sphere, and thus by changing the diameter of the equivalent sphere.v r and ω r represent the relative velocity and angular velocity between the ice particle and fluid, respectively.C D and C R represent the coefficients of drag force and drag torque, detailed expressions of which can be found in Haider and Levenspiel (1989).
The pressure gradient force is determined according to Siemens PLM Software (2019): where V p represents the volume of ice particles and ∇p s represents the gradient of the static pressure in the fluid.The lift force F l is divided into shear lift force F ls and spin lift force F lr as per the following: where ω represents the angular velocity of the fluid, and the coefficients C lr and C ls are determined as per Saffman (1965) and Sommerfeld (2000).

Ice-ice and ship-ice contact force and torque
The contact forces F c of the ice-ice or ship-ice interactions are calculated using a linear spring contact model, which is a linear simplification of the Hertz contact model that is widely used in DEMbased models to solve contact forces (Renzo and Alberto, 2004;Luo et al., 2020).The linear spring contact model is essentially a kind of spring-dashpot model, in which the spring provides an elastic force and the dashpot provides viscous damping (Johnson, 1987).As a result, the contact force F c can be divided into two components along the normal and tangential directions of the touching surface of the two spheres, as shown in Fig. 5.In other words, the normal direction is along the connecting line of the centers of the two spheres, and the tangential direction is perpendicular to it.The normal component F n and the tangential component F t can be presented as: where the subscripts n and τ represent the normal and tangential components, respectively, and n and τ are the unit vectors of the normal and tangential directions, respectively.The particle-particle contact model and the contact force are shown in Fig. 5, in which the direction of incidence indicates the relative velocity direction between these two particles.When the ship-ice contact force is calculated, the ship particle can be seen as a particle with an infinite diameter.Under this circumstance, the normal and tangential directions are normal and tangential to the ship's hull, respectively.
Similarly, the contact torque can be expressed as where M n and M t are normal and tangential torque components, respectively.For a detailed background of Eqs. ( 15) and ( 16), see Ni et al. (2020b).

Coupling scheme and iteration process
There are three components of the two-way CFD-DEM coupling scheme: (a) CFD-model, (b) DEM-model, and (c) coupling model.Following the initialization, the calculation iteration of these three components is conducted at each time step.First, the porosity in each fluid cell is calculated before calculating the fluid force F f acting on each particle.The volume average force acting on each fluid cell F OM can also be obtained.The forces acting on both the fluid cell and particles are calculated in the coupling part.Second, in the DEM part, the contact force F c acting on each particle is calculated.Combining the fluid force acting on each particle, F f is calculated in the coupling part following the motion equations ( 7) and ( 8), in which the motion and position of each particle are updated.The iteration of the DEM component is repeated m times with the existing fluid information.Third, in the CFD part, equations ( 4)-( 6) of the fluid phase are solved.Subsequently, the velocity field and the pressure field are obtained.This process is repeated until the termination condition is met (Norouzi et al., 2016).
In the two-way coupling scheme, the ice-fluid interaction force mutually affects the continuous phase and discrete elements, making the required calculations highly resource-and time-consuming.On the other hand, in the one-way coupling scheme, ice particles are affected by the fluid, but the related feedback is not considered.For a low-speed simulation (V = 0.464 m/s), Luo et al. (2020) demonstrated that one-way coupling is recommended for the less influential CFD-DEM case, resulting in the same calculation effect as two-way coupling at a     2020) explored the accuracy of the one-way coupling approach and showed that the one-way coupling approach was sufficient to provide accurate resistance predictions for V ≤ 1.18 m/s.Therefore, because the one-way coupling scheme can simulate ship-ice interactions with sufficient accuracy and at a low calculation cost, we choose to adopt and apply it in the present study.

Computational domain set-up
In the computational domain setup, the ship model is defined to correspond to that used in the model test.Accordingly, the flow domain (including both water and air) is set as 8.65 L in length, 2.40 L in width, and 1.00 L in height.Out of the 1.0 L height of the calculation domain, the heights of the water and air domains are 0.543 L and 0.457 L, respectively.The length of the broken ice field is set as 4.0 L. Following Huang et al. (2020), the relative velocity method is adopted, which means that the ship model is kept stationary whereas the fluid together with the ice particles exerts a specific velocity astern.
A sketch map of the computational domain with boundary conditions is shown in Fig. 7, where the inlet boundary is set as the 'velocity inlet' with a given velocity setting, and the outlet is set as the 'pressure outlet' with the hydrostatic pressure distribution.In addition, the top, bottom, and sidewalls are also set as 'velocity inlets', but they are set as a field function that can obtain the velocity distribution on the boundary from its surrounding velocity field instead of from any specific velocity (Siemens PLM Software, 2019).
The applied mesh layout in the computational domain is presented in Fig. 8.The total number of cells is taken as 0.862 million following the cell independence verification in Fig. 9  element setup in DEM, the size and shape of each ice element is consistent with the model test, with each broken ice containing only a single square element.The relevant parameters were kept consistent with the model tests in the DEM setup, with ice particle density ρ i = 917.0kg/m 3 , the friction coefficient is 0.138 for ship-ice.For Young's modulus of ice particles, there is λ times difference between real ice in the full-scale and model scale (Luo et al., 2020), and other parameters in the DEM setup are shown in Table 2.A plane part injector is used to produce the ice particles in the simulation.Specifically, a plane injector with a specific size is set to cover the target area before the ice particles are injected into the fluid domain.The distribution of broken ice particles is random.All the broken ice particles are inserted into the fluid domain at the first-time step in the computing process.The initial state of the simulation is shown in Fig. 6.

Empirical formulas
The total resistance (R) of a ship navigating in a broken ice field can be broadly divided into open water resistance (R OW ) and broken ice resistance (R i ).Most empirical formulas focus on the latter, which can be obtained directly from the formula.Based on the foregoing introduction, in general, the ice resistance formula would consist of three main components: the ship's hull parameters (e.g., ship's width, bow characteristics, etc.), operational parameters (e.g., ship speed, ice thickness, ice concentration, etc.), and empirical coefficients.The derivation of the formulas is based on assumptions about the influencing parameters, and then empirical coefficients are given based on available results from experimental or numerical simulation results, etc.In the following, three empirical formulas concerning navigation in the broken ice field are introduced.
Colbourne (2000) presents a formula given by Equations ( 17)-( 19).A series of modeling test results support it to obtain them.Two different model ship were selected and 40%-100% concentration broken ice fields were prepared in which the broken ice was cut into similar rectangular blocks.In Coulbourne's formula, the coefficient of ice resistance (C i ) is introduced, which relies on the ice Froude number determined by the ship's speed, ice thickness (H) and concentration of broken ice (η), as shown in Equations ( 18) and ( 19).Besides, two main impact parameters of ice condition, the η and H are also directly considered when calculating R i .
where n = 2 for ice-breaker (Molyneux and Kim, 2007) 2021) for ships navigating through floating ice floes.For the ice floe region, the shape of ice floe is simplified to a disk, with different sizes but the same thickness.In terms of the effect on the hull parameters, three different ships were selected.A series of hull parameters were taken into account, and specifically, in order to further consider the effect of ship-water-ice coupling, two characteristic parameters at the bow of the ship were analyzed in detail for the effect of ship-water-ice interaction in vertical and transverse direction: buttock angle (γ), and waterline angle (α).They analyzed the relationship between the relevant parameters based on the ice resistance of three model ships, verified by the experimental data and CFD-DEM approach, and proposed the formula to estimate the R i .
where α is the waterline angle at 1/4B, γ is the buttock angle, d i is a parameter introduced to indicate the size of broken ice pieces, modeled as the equivalent diameter of the upper surface of the ice and can be taken 10 times H in practice; a more detailed description of the equivalent diameter of the ice piece is referred to Huang et al. (2021).
A further formula to estimate R i was proposed by Sazonov and Dobrodeev (2021), which is based on an earlier semi-empirical method for the calculation of ice resistance in the newly created broken ice channel.In the derivation of their formulas, they mainly relied on their model test data (Dobrodeev, 2018).Similarly, five different ships were modeled for the tests, and the effects of a range of hull parameters were considered, as detailed in Eq. ( 21).It is noteworthy that they took the effects of the width of the broken ice channel into account as well.However, since the data source for the formula focuses more on the effect of narrow channels on resistance (B C /B WL up to 2.0), its application to wider areas of broken ice remains to be explored.The formula is expressed as follows: where B C is the channel width, here is the same setting as in the experimental and numerical simulation is 2.0 m, Froude number for a given ice thickness, β = arctan(tan γ /sin α 0 ).

Ship-broken ice interaction
In this section, the ship-ice interaction processes observed during the towing tank tests and numerical simulations are analyzed against each other.In the numerical simulation, all the parameter values are defined to correspond to those of the towing tank tests.Fig. 10 shows the corresponding experimental and numerical simulation results side-by-side at three different time instances in which t = 0 s represents the moment when the ship bow enters the broken ice field.

Table 3
The total resistance (in Newtons) as well as the uncertainty from the experiments, the numerical simulations, and the empirical formulas, the differences of the resistances from the simulations and the formulas against the experimental results.
Case Nr.
V (m/s) η Experiment Numerical simulation Formula-1 Formula-2 Formula-3 Fig. 16.Comparison of open water resistance and total resistance from the experiment, the numerical simulations, and the formulas.
Fig. 10 (a) shows the situation when the bow contacts the ice pieces for the first time, and some ice pieces are pushed towards the bow, as highlighted by the red dotted line.In the following, the ice pieces that are touching the hull surface are referred to as being in the 'contact mode', whereas the ice pieces that are not in direct contact with the hull surface are referred to as being in the 'noncontact mode'.As shown in Fig. 10 (b), as the ship continues to move ahead, the ice pieces being in 'contact mode' slide from the bow to the ship side, as observed in a moving coordinate system.While the side of the mid-ship section of the hull is vertical, the ice pieces being in 'contact mode' is also upright in the water.Once the entire hull has entered the broken ice field, as shown in Fig. 10 (c), the ice pieces being in 'contact mode' is pushed aside by the ship stern and turned on its side.Because this ice piece cannot move quickly towards the central line of the ship, an open water channel with a width corresponding to the beam of the ship is formed, as denoted by the red dotted line.We find that the outcomes of the numerical simulation generally agree well with those of the experiments, including the contact, turnover, and sliding of ice pieces in 'contact mode'.However, some ice pieces being in the 'noncontact mode' demonstrate an unrealistic turning phenomenon at a distance from the hull.This is because, in the CFD-DEM coupling simulation, the fluid forces are exerted on the center of the DEM particles instead of on the boundaries.This saves calculation resources and enables the consideration of a large number of particles (Norouzi et al., 2016).However, because the boundaries of DEM particles are not captured in this numerical algorithm, ice pieces are easily pushed upright in the water in a collision and may overturn.As shown in Fig. 10 (c), even once the ice pieces have been separated from the hull, they may still have difficulty recovering to their initial floating state.Similar phenomena were observed in previous numerical simulations, see Huang et al. (2020).Future research is needed to address this issue.
Although unrealistic, the turning phenomena indicate that some ice pieces being in the 'noncontact mode' are significantly affected by the ship's motion.To check this influence in the experiment, Fig. 11 presents three instances of the experiment, in which the yellow dotted line denotes the extension cord of the intersecting line of the hull side and water surface.For clarity, in Fig. 11 we select and track five individual ice pieces.In the first frame (Fig. 11 (a)), the extension cord cuts through the numbered model ice pieces.When the ship's shoulder passes through the model ice pieces (Fig. 11 (b)-(c)), they are pushed to the same side of the extension line due to their interaction with the model ice pieces being in the 'contact mode'.The movement of these model broken ice in the y-direction can be observed from these frames.It should be noted that the model ice pieces being in the 'contact mode' are in turn also affected by those being in the 'noncontact' mode, as these affect the level of ice concentration around the ship side.
The velocity of ice pieces is difficult to measure in the experiments, but it can be easily obtained from numerical simulations.The relative velocities between the ice pieces and the ship in the directions of the yaxis (V ry ) and x-axis (V rx ) are shown in Fig. 12.As the ship moves forward, V ry at both sides of the ship is increased, indicating that the ice pieces are pushed outwards (Fig. 12 (a)).V rx is negative, as its direction is opposite to the x-axis.It can be observed that V rx decreases sharply in the vicinity of the bow area, mainly because of the interaction with the ship bow, as shown by the 'velocity variation zone' in Fig. 12 (b).

Processing and analysis of results
For the ship resistance (R) from the experiment, Fig. 13 shows the recorded results of the force sensor signal in a whole experiment process, when V = 0.5 m/s, η = 60%.Here, time zero corresponds to the time when the hull enters the broken ice field.The course of the experiment was divided into four stages, which were roughly separated and marked with dotted lines in the figure.In stage 1, the towing system accelerates together with the ship model to the required speed.Stage 2 represents the period when the ship model has reached its required speed but has not yet made any contact with the ice pieces.Stage 3 represents the period when the ship model passes through the broken ice field, which starts with the bow touching the broken ice field and ends with it leaving the broken ice field.Stage 4 represents the deceleration period.It can be seen from the figure that in stages 1 and 4 there are positive force signals that accelerate the ship model and negative force signals that slow down the ship model, respectively.Stage 2 represents the towing resistance of the ship model in the open water area, and its resistance value is slightly above zero at most times.On this basis, there are a few peak values.In stage 3, because the ship model has entered the broken ice field, the number and magnitude of the peak values are significantly larger than that in stage 2. Fig. 14 shows the time-history curves of ship resistance as obtained from the numerical simulation and stage 3 of the experiment.The dotted lines in the figure show the average value of resistance at 0.2s intervals.It is clear that there are many periodic spikes in experimental results, which may be due to the "surfing phenomenon" of the ship's hull, which was mentioned by Guo et al. (2018) in a previous model test, which described how the bow was supported by the upward force of the broken ice, resulting in increased pitch motion of the ship.In this experiment, the icebreaker bow makes it easier for the broken ice to move to the bottom of the bow, and a large amount of broken ice provides support to the bow (shown in Fig. 15), which further increases the reciprocal motion of pitching.As a result, the peak value is affected by this reciprocal pitch.A similar phenomenon was also reported by Ettema et al. (1987) and von Bock und Polach (2010).They both found that the peak and mean value of resistance is greater in free mode than in constrained mode and the difference in mean value can be decreased for medium and high ship speeds.A similar performance can be seen in Fig. 14, where the ship model is constrained mode in the numerical simulation and generally has a smaller peak than the experimental results.However, the significance of peak values has not been evaluated.This paper is more concerned with the assessment of the mean value of ship resistance at medium and high speeds, and the significance of these peaks needs to be explored by conducting tests that are more focused on it.
The results from numerical simulation consist of two components, one of which is the ice resistance, which is shown in many fluctuating peaks.These peaks are the recorded results of the x-axis collision force between the hull and the ice particles at each timestep.The other is the relatively stable water resistance, the sum of which with the ice resistance represents the total ship resistance (R) in the numerical simulation.However, since the pitch of the ship was not modeled in the numerical simulation, the "surfing phenomenon" is not obvious, and the peak value is also smaller than the experimental result.Nevertheless, the impact of increased pitching on the peak in broken ice regions needs more detailed studies.
To account for variations in the initial conditions of broken ice field arrangement in model test, three repeated tests were carried out for each of the considered test conditions.Subsequently, the analysis process recommended by ITTC ( 2005) is applied to estimate the experimental results as well as their uncertainty.The methodology of this analysis process is separating the single run test into several segments without considering the unstable stage, the mean resistance of each segment is a data sample to calculate the random uncertainty (U).After the segments are ready, the Chauvenet's criterion is applied to identify the valid samples, where the (Chauv#) mean is the Chauvenet number, F T Mean is the mean resistance of each segment, Mean F T mean and STD F T mean are the mean and standard deviation of F T Mean .The Chauvenet number can be found in ITTC ( 2008) and if the (Chauv#) mean of segment is larger than the Chauvenet number, it should be disregarded.The valid segments are selected, and the random uncertainty (U) can be calculated as where t can be found in ITTC (2008) and N is the number of valid segments.The length of each segment is recommended as 1.5-2.5 times the length of model ship (Jeong et al., 2021), and in this experiment, the segment is assumed as 2.0 times L. Take that the whole ship has already entered in broken ice field and before it exits as the stable stage, one run of test can be separated into 2 segments, and there is a total of 6 segments in total after three repeated runs.The mean value of F T mean valid segments is taken as the ship resistance in model test, which as well as the random uncertainty are shown in Table 3.
Based on this data processing process, ship resistance from the experiment and numerical simulation can be obtained.Table 3 and Fig. 16 show the results of ship resistance (R) in different cases from the experiment, the numerical simulations, and the empirical formulas, which include four ship speeds (V) and three ice concentrations (η).The uncertainty of model tests as well as the resistance differences from the simulations and the formulas against the experimental results are also listed in Table 3. From the uncertainties in the test results, most of the uncertainties due to the different initial ice fields are below 20%, but in very few cases they exceed 20%, which is related to the distribution of the broken ice uniformly at the initial ice arrangement.Most of the uncertainties at 80% concentration are smaller than others.Compared to 70% and 60%, the ice arrangement at η = 80% makes it easier to achieve a more uniform distribution.Nevertheless, the results from the numerical simulation are overall in good agreement with the experimental results.However, one thing that is quite obvious is that the numerical results overwhelmingly give lower estimates than the model tests, and this may be due to the above-mentioned unrealistic turning phenomenon of broken ice in simulation, which indicates that the outer broken ice is unable to provide a restraining force on the ice in "contact mode" for as long a sustained period of time as they do in the model test.
In the estimation by the empirical formulas, the open water resistance (R OW ) corresponds to that used in the numerical simulation.For brevity, the sum of estimations of R i from Coulbourne's formula and R ow is named formula-1, the sum of R i from Huang et al. (2021) and R ow is named formula-2, and the sum of R i from Sazonov and Dobrodeev (2021) and R ow is named formula-3.The values of the water resistance for each condition are also shown in Fig. 16, and their approximate components of the total resistance are likewise evident.It can be clearly seen that the open water resistance as a percentage of the total resistance increases significantly as the velocity increases.As can be seen from Fig. 16, in most cases the results from each approach agree well with each other.However, looking at them separately, in most cases, formula-1 provides results similar to the model test, with most of the differences being under 10%; formula-2 generally gives a lower estimate, while formula-3 generally gives the largest estimate, with the other prediction methods fluctuating between them.
From the background comparison of the empirical formulas, it is clear that differences between the data on which the formulas are based as well as the target scenarios can have some impact on the results.Formula-1 has the most similarity to the numerical simulation and experiments in the ice arrangement, and also achieves similar results; for formula-2, the estimate is smaller because it simplifies the ice floe to a disk, which is less likely to produce the "force chain" phenomenon common in rectangular ice; for formula-3, it takes into account the effect of narrow ice channels, but the forecast is slightly larger for the broken ice field, which lacks channel boundary constraints.

Effects of ship speed and ice concentration
In order to analyze the effect of ship speed and ice concentration in the broken ice field, we carried out simulations at six different speeds between 0.3 and 0.9 m/s, at 0.1 intervals, and for three different ice concentrations, between 60% and 80%, at 10% intervals.Model tests were carried out at four different speeds between 0.5 and 0.8 m/s.Fig. 17 presents V rx and V ry of ice pieces at different ship speeds for η = 60%.The region in the red ellipse is the 'velocity variation zone' and it grows along with the ship speed, especially in the y-axis direction.This is because the collisions become heavier at larger V, and ice pieces obtain more energy from the ship.However, the impenetrable condition on the ship significantly limits the motion of the x-direction around the bow, resulting in submerged ice pieces interacting with the bottom of the hull.Fig. 18 presents a comparison of the relationship between ship resistance and speed for different ice concentrations (η = 60%, 70%, and 80%).The endpoints of the vertical line in the experiments are the results of uncertainties.The results of the experiments and numerical simulations agree well in all three cases of ice concentration, both of them tend to increase significantly with ship speed.This also proves the applicability of numerical simulation and modeling tests in the cases of this paper.The collision forces become stronger at larger V, and more energy is transferred from the ship to the ice pieces, resulting in higher ice resistance.On the other hand, R ow in the ship resistance also increases rapidly with increasing V.
For the results from the empirical formulas, which show a velocitydependent multiplying power trend and have a similar growth pattern increases with V.At each velocity point, the results from formula-1 are closest to the numerical results overall and maintain a good agreement for all three ice concentration cases.Formula-2 and formula-3 both show good agreement with the numerical results, but the difference increases as the speed increases, the formula-2 has some downward deviation of the estimates compared to numerical results and the formula-3 is opposite.As mentioned earlier, the differences brought about by different target ice fields widen as the velocity increases, and the predictions of Formula-3, in particular, suggest that the effect of fairways is more pronounced at high speeds.
For the analysis of ice concentration, Fig. 19 shows a bottom-view comparison of ships at different ice concentrations for V = 0.5 m/s.
Here the initial distributions of the broken ice field of η = 70% and η = 80% are different from that of η = 60%.For η = 60%, the random distribution is adopted.However, for η = 70% and η = 80%, a random distribution cannot be used as there would be a significant risk of overlapping broken ice.In the numerical simulation, if the overlap occurs, it is necessary to terminate the calculations and select another distribution pattern.Through an iterative process, we find that the regular distribution of broken ice can satisfy the high ice concentration and that the ice resistance obtained is close to that of a random distribution.Based on this finding, we select a regular distribution for ice concentrations exceeding 60%.
As per Fig. 20, for each considered speed, ship resistance increases with η.At constant speed, the impact forces between the hull and the ice pieces and the water resistance can be assumed independent of the ice concentration.Higher ice concentration results in more ice pieces being submerged by the hull (shown in Fig. 18), an increase in the ship-ice contact area and increased friction force as well as more frequent collisions on the hull increasing the overall resistance of the ship.On the other hand, in terms of the comparison of the results from the empirical formulas and numerical simulation, they all seem to be linearly related in proportion to the η.The results of formula-1 provide a good comparison at all ice concentrations.However, for formula-2, there is generally a good contrast at η = 60%, but from η = 70% onwards, its predicted values start to be generally lower than the numerical results, most significantly at η = 80%, while this gap is clearly greatest at V = 0.8 m/s.Because the shape of the ice floe is simplified to a disk in formula-2, it is difficult to have "force chains" between ice-ice, which is very common at high ice concentrations.This also suggests that Formula-2 lacks sensitivity to η. Formula-3 always provides a larger estimate, which also shows an almost linear increase with increasing η.
The slopes of their predictions are significantly smaller than the experimental and numerical results.However, the predictions are significantly larger from η = 60% onwards than the experimental and numerical results, because the channel limitation effect is taken into account.The difference between them gradually decreases as η increases.However, for the cases studied in this paper, the effect of channel restrictions is not obvious and may lead to overestimation.

Estimation of propulsion power demand
In the design stage of polar ships, the estimation of ship resistance is mainly intended to be applied to the estimation of the engine output under various ice conditions.The required propulsion power of the ship for different speeds and ice concentrations is estimated by the previous ship resistance from numerical results and the empirical formulas in this section.
To obtain the required propulsion power in full-scale, the ship resistance in model-scale needs to be converted to full-scale.Ice resistance (R i ) can be scaled to full-scale to comply with ITTC-recommended guidelines, R i−full = λ 3 R i , where R i−full is the ice resistance in full-scale.The water resistance obtained from the numerical results is extrapolated to full-scale using the method proposed by Hughes (1954).The full-scale ice resistance is added to the corresponding water resistance to obtain the predicted ship resistance in full-scale (R full ).
The propulsion power demand for the design and construction stage of a new ship is estimated based on the in Finnish-Swedish ice class rules (TRAFI, 2021).The estimation process is based on three factors, the type of propulsion system, the diameter of the propeller and the ship resistance.Among them, the ship resistance is estimated based on the brash ice condition.In order to satisfy the propulsion power demand for the broken ice condition, it is estimated based on the resistance estimation values of several methods used in this paper, and the specific formula is as follows: where P is its required propulsion power, K e is taken as 1.44 according to its type of power equipment and the number of propellers, D p is the diameter of the propeller of 4.2 m.However, this estimation approach is only an attempt to broaden the scope of its original application, and more practice is needed to verify.Table 4 and Fig. 21 show the results of the estimated propulsion power demand, which tends to increase similarly with V at all η.As can be seen from Fig. 21, the output power predicted by these methods has a similar increasing trend with V, with good agreement at low speeds.However, the difference between them increases significantly as the speed increases.The numerical results are most similar to those predicted by formula-1, formula-2 predicts the slowest increase of required power with speed, in contrast to formula-3, which increases most rapidly, especially when V = 0.9 m/s.It is noted that we estimate the required propulsion power demand roughly, without considering detailed propeller characteristics or mechanical losses.

Conclusions and discussion
This study presents and analyses two complementary approaches to predict ship resistance in broken ice, of which one is experimental and the other numerical.The experimental approach makes use of a type of non-refrigerated synthetic model ice made of polypropylene, which makes it possible to test how a ship behaves in broken ice using a conventional non-refrigerated towing tank rather than an ice tank.Naturally, this makes the tests significantly less time and resource-intensive.On the other hand, it can well-validate the numerical methods, and make the results more reliable.The numerical approach is based on the CFD-DEM method and can be used to consider fluid effects, such as the changes in fluid velocity and ship waves, while the ship is moving ahead.
It is easy to conduct and can capture results that are difficult to measure in model tests, e.g., velocity variation in the broken ice field.To some extent, it can be used as a complement to model tests.Some conclusions can be presented as follows.
(1) The results obtained by different methods: modeling tests, numerical simulations, and empirical formulae predictions were compared, and in most cases they are in good agreement with each other and have similar trends, proving the feasibility of using this time-and cost-efficient manner in this paper for predicting the resistance.
(2) The differences between the empirical formula predictions, and numerical as well as experimental results are compared.
Combining the sources and backgrounds of the empirical formulas, it is found that when an empirical formula is applied, the predictions differ noticeably depending on the settings of the broken ice fields from which the data are derived.When the existing empirical formulas do not satisfy the application requirements, the combined method in this paper might be applied as a low-cost alternative with acceptable accuracy.(3) Combined approaches in this paper may, for instance, be used to estimate a ship's propulsion power demand for different speeds and ice concentrations.This approach can take into account the operational conditions as well as the various broken ice conditions to be faced, such as ice size, concentration, etc.On this basis, more detailed power estimates for various cases can be given.Thus, we believe that the proposed approaches may support the design of polar ships, and voyage planning in Arctic waters, among others.However, this study is only an attempt to apply the approach, and more practice is still needed to verify it before it can be used in real-life engineering applications.
Limitations of the work include the following.First, the applied nonrefrigerated model ice differs from natural ice in that it is non-breakable, meaning that it cannot be used to replicate ship-ice interactions resulting in the breaking of ice floes.Therefore, the utilization of the approach is limited to broken ice fields composed of relatively small ice floes, in which the ice resistance is mainly caused by ship-ice collisions and friction rather than by the breaking of ice floes.Second, for the numerical simulation, the applied one-way coupled ship-water-ice interaction model neglects the forces between the ice floes to the flow field.Nevertheless, this simplification can be considered a trade-off between computation accuracy and costs, as considering the flow field action would require significantly higher computational efforts.Thirdly, the present study considered a single type of ship.Finally, we applied an established methodology to scale up the predicted resistance value from the model to full-scale, but the accuracy of the scaled-up value is not evaluated.Future research is needed to evaluate the accuracy of model test predictions.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Three parts of the towing tank in the experiment.

Fig. 6 .
Fig. 6.The initial arrangement of the computational domain.

Fig. 8 .Fig. 9 .
Fig. 8.The mesh layout of the computational domain and the ship.
Fig. 10.Experimental and numerical results at three different time instances at V = 0.5 m/s, η = 60%.

Fig. 12 .Fig. 13 .
Fig. 12.The contours of the relative velocities of broken ice to ship in the y-axis and x-axis V ry and V rx in the numerical simulation at V = 0.5 m/s, η = 60%.

Fig. 17 .
Fig. 17.Velocity of broken ice in the x-axis (the top row) and y-axis (the bottom row) directions at different speeds.

Fig. 18 .
Fig. 18.Relationship between ship speed and resistance in various ice concentrations.

Fig. 19 .
Fig. 19.Speed of ice pieces in the direction of the x-axis (top row) and y-axis (bottom row) for different ice concentrations.

Fig. 20 .
Fig. 20.Relationship between ship resistance and ice concentration for various speeds.

Fig. 21 .
Fig. 21.The results of estimation of required propulsion power.

Table 2
Main parameters of the DEM simulation (in model scale).
Y. Xue et al. lower calculation cost.Furthermore, Huang et al. (

Table 4
The required propulsion power (P) from numerical simulation and formulas.