A Data-Driven Model for Ice-Breaking Resistance of Structure Based on Non-Smooth Discrete Element Method and Artificial Neural Network Method

Abstract

In this paper, a data-driven model based on the Non-smooth Discrete Element Method (NDEM) and Artificial Neural Network Method (ANN) is proposed for the computation of the ice-breaking resistance of the structure. The idea of so-called “meta-modelling”, which means establishing an Artificial Neural Network (ANN) model based on a pre-computed ice failure database to avoid the time-consuming direct resolving of the ice fracture process, was integrated in the non-smooth discrete element method (NDEM). The developed model was validated by simulating the ice-breaking process of the cone structure, and the computational results match well with the experimental ones in the literature. After that, the effects of various parameters on the ice-breaking resistance were analyzed by the developed model. It was found that the factors that have great influence on the resistance of cone structure in level ice condition are the cone angle, navigation velocity and ice thickness.

1. Introduction

Ships and other types of marine structures in the Arctic area have long been a popular issue in the field of Naval Architecture and Ocean Engineering. The reliable prediction of the ice resistance of ships and marine structures in ice-covered areas are related to their safe and efficiency.

The interaction of ships and marine structures with ice has been studied by experimental and numerical methods. Field tests and model scale experiments are important tools to understand the underlining physics of the ice-structure interaction process, some representative work include: Riska et al. [1], Frederking [2], Suyuthi et al. [3], Chernov [4], Kujala et al. [5], and model tests: Huang et al. [6], Jones [7,8], Aksnes and Bonnemaire [9], Woolgar and Colbourne [10], and Jeong et al. [11]. Some experiments have been carried out in the Hamburg Ship Model Basin (HSVA) [12,13], including tests of cones breaking ice in layers of ice. Based on the analysis of field measurement data and test results, the ice load and ice—structure interaction process of the structure in the ice region can be obtained.

Alternatively, the rapid development of various numerical methods, including Finite Element Method (FEM) [14], Discrete Element Method (DEM) [15] and Peridynamics method (PD) [16], have achieved the direct simulation of physical process such as ice failure during ice-structure interaction. Additionally, Smoothed Particle Hydrodynamics method (SPH) models can also be applied in the simulation of ice-structure interaction, such as using DynaRICE [17]. Based on the above-mentioned numerical simulation methods, a more subdivided numerical simulation method is derived.

Discrete element method (DEM) can be divided into Smooth Discrete Element Method (SDEM) and Non-smooth Discrete Element Method (NDEM), which are distinguished by the methods of contact calculation and time integration [18]. Compared with SDEM, NDEM allows larger time step, due to its implicit time integration.

The programs that are based on NDEM is commonly called Physic Engines in the field of Computational Graphics and many of them are open-source libraries such as Bullet, ODE, Box2D, have been developed based on NDEM, and a comparison between them can be seen in [19]. Some of them have been used in ice-related studies. Lubbad and Løset [20] proposed a numerical model to simulate the ship–ice interaction by using the physics engine “PhysX”. Metrikin [21] proposed a further non-smooth discrete element method, including a computational framework for the Bullet physics engine.

However, physics engines are mainly used to simulate collision rigid body motion. It is difficult to use physics engines to simulate deformation and stress and strain. Considering the accuracy of calculation, the cohesion element method, such as Kuutti et al. [22], can be used to solve the problem of level ice failure, but it leads to a considerable increase of computational time. Alternatively, the data-driven models such as the Artificial Neural Network (ANN) method can establish a simple function between the input and output parameters of a very complex physical problem based on the data of huge amounts of precise simulations. The resultant simple function can be used to replace the time-consuming simulation for other complicated scenarios, which significantly improves the simulation efficiency. By using an Artificial-Neural-Network-based data-driven method, Li et al. [23,24] established a model for ice-breaking computation and promising results have been reported. This data-driven model is so-called “Meta-Modelling” [25].

Therefore, physics engines and “Meta-Modelling” can be combined by adding the ice-breaking calculation model established by the data-driven method into the physics engine to complete level ice fracture simulation. This method can be used to enrich the computational scenarios of the physics engine, such as simulating level ice fracture and the collision of shattered ice generated after fracture.

Based on the simulation method mentioned above, this paper calculates the ice-breaking parameters by using the empirical formula fitted by artificial neural network, and simulated interaction between the cone and level ice by Bullet physics engine. The remaining parts of this paper are organized as follows: Section 2 introduces the calculation principle of physics engine and the setting of ice-breaking model. In Section 3, the simulation of cone–ice interaction and resistance was compared with the test to verify the accuracy of the model, and then the cone’s ice-breaking resistance was parameterized, including cone angle, cone motion velocity, cone waterline width and ice thickness. Finally, a discussion and direction for further improvement of this method were concluded in Section 4.

2. Methodology

In this study, one of the mainstream physics engines, i.e., a Bullet open-source library [26] is used for the simulation. The ice breaking was computed by a data-driven model based on a pre-calculated database and ANN. The details of both the Bullet and ice failure model are presented in this section.

2.1. Numerical Model

This section describes the technique details of the adopted physics engine, i.e., the Bullet open-source library and additional physics models used in simulation. The main physics model is from a Bullet physics engine [26] which includes collision detection, constraint computation and time integration. The fluid model is not included by Bullet physics engine, so the simulation adopted method of Yang et al. [27], using Morison’s equation [28] to compute the drag of rigid bodies.

2.1.1. Collision Detection

Bullet uses the same two-level search scheme for collision detection as most physics engines. It searches for objects that are far away from the calculated target through broad phase detection, and then uses the narrow phase detection to find the collision pairs and their corresponding contact information. Sweep and Prune (SAP) algorithm [29] is used in the broad phase detection, and Gilbert–Johnson–Keerthi (GJK) algorithm [30] is used in the narrow phase detection. The contact points and normal vectors between convex bodies of arbitrary shape are solved by this algorithm. Thin envelopes are added around the collision shapes by Minkowski and to allow small mutual penetration, and the Expanding Polytope Algorithm (EPA) [31] is applied to avoid large overlaps.

2.1.2. Constraint Computation

The momentum equation of body i in the multi-body collision system can be written as:

$$m_i\mathsf{\Delta }{u}_i=F_{ext,i}\mathsf{\Delta }t+{\displaystyle {\sum }_k{P}_{cont,k}}$$

(1)

$$I_i\mathsf{\Delta }{\mathsf{\Omega }}_i={\tau }_{ext,i}\mathsf{\Delta }t+{\displaystyle {\sum }_k{r}_{k,i}\times P_{cont,k}}$$

(2)

where $m_i$ and $I_i$ are mass and moment of inertia, respectively. $\mathsf{\Delta }{u}_i$ and $\mathsf{\Delta }{\mathsf{\Omega }}_i$ are the translational and rotational velocity change over a given time step $\mathsf{\Delta }t$. $F_{ext,i}$ and ${\tau }_{ext,i}$ are the generalized force vector of body i except the contact force. $P_{cont,k}$ is the contact impulse at the contact point k and $r_{k,i}$ represent the vector from the center of mass to the contact point k. For a contact pair, as shown in Figure 1, the contact impulse $P_{cont,k}$ at the contact point k can be expressed as:

$$P_{cont,k}^n=\frac{\left(\mathsf{\Delta }{u}_{rel}\cdot n\right)\cdot n}{1/m_A+1/m_B+\left(r_{k,A}\times n/I_A\right)\times r_{k,A}\cdot n+\left(r_{k,B}\times n/I_B\right)\times r_{k,B}\cdot n}$$

(3)

$$P_{cont,k}^t=\frac{\left(\mathsf{\Delta }{u}_{rel}\cdot t\right)\cdot t}{1/m_A+1/m_B+\left(r_{k,A}\times n/I_A\right)\times r_{k,A}\cdot t+\left(r_{k,B}\times n/I_B\right)\times r_{k,B}\cdot t}$$

(4)

where $P_{cont,k}^n$ and $P_{cont,k}^t$ are the normal and tangential component of the contact impulse. Likewise, n and t = (t₁, t₂) are the unit vectors in the normal and tangential direction. $\mathsf{\Delta }{u}_{rel}=\left( u_{rel}^{+}-u_{rel}^{-}\right)$ is the change of the relative velocity at the contact point. $u_{rel}^{+}$ is the separate velocity while $u_{rel}^{-}$ is the approaching velocity. Considering the energy dissipation during the collision, $u_{rel}^{+}$ and $u_{rel}^{-}$ obey the Newton impact law:

$$u_{rel}^{+}=-e{u}_{rel}^{-}$$

(5)

where e is the restitution coefficient.

The normal velocity and the normal contact impulse constitute the complementarity condition that if one of them is positive and the other one must be zero. Meanwhile, the tangential impulse is limited by the Coulomb’s friction law:

$$-\mu \left|\left|P_{cont.k}^n\right|\right|\le P_{cont,k}^t\le \mu \left|\left|P_{cont,k}^n\right|\right|$$

(6)

where μ is the coefficient of friction. The dynamic and static coefficients of friction are assumed to be equal in the present simulation.

2.1.3. Time Integration

The velocities and positions of the bodies are updated by the simplistic Euler method:

$$u_i^{t+\mathsf{\Delta }t}=u_i^t+\mathsf{\Delta }{u}_i$$

(7)

$${\mathsf{\Omega }}_i^{t+\mathsf{\Delta }t}={\mathsf{\Omega }}_i^t+\mathsf{\Delta }{\mathsf{\Omega }}_i$$

(8)

$$x_i^{t+\mathsf{\Delta }t}=x_i^t+u_i^{t+\mathsf{\Delta }t}\mathsf{\Delta }t$$

(9)

$${\psi }_i^{t+\mathsf{\Delta }t}={\psi }_i^t+{\mathsf{\Omega }}_i^{t+\mathsf{\Delta }t}\mathsf{\Delta }t$$

(10)

where x_i and ${\psi }_i$ represent the position and rotation of body i. It should be noted that the contact force could not be calculated directly in Bullet, but it could be retrieved from the contact impulse. Dividing the impulse by the numerical time step yields the contact force.

As for the selection of the time step, the contact force conversion model of Yulmetov is adopted [32], and the collision time is calculated as the following:

$$T_{col}\approx \frac{\pi }{2}\sqrt{\frac{M^{*}}{K}}$$

(11)

here M* is reduced mass, K is contact stiffness.

In cone and ice interaction, the ice floes’ mass is much smaller than cone’s mass, so reduced mass is equal to ice floes’ mass. Another important influence is the added mass of the ice floe, according to Zong [33] and Isaacson [34]. The ice floes’ added mass will become very large when the ice floe is close to the cone, so we consider to use added mass coefficient to estimate the reduced mass:

$$M^{*}\approx \left(1+\alpha \right)M_{ice}$$

(12)

here α is additional mass coefficient and $M_{ice}$ is mass of ice floe in collision. Referring to the study of Yang [27], the additional mass coefficient of the triangular ice here is about 0.5, and the contact stiffness is 1.0 × 10⁶ N/m.

The collision time calculated by Equation (11) is about 0.012–0.061 s, and Lubbad [20] chose the 0.01s as time step, so the simulation adopts 0.01 s as time step.

2.1.4. Fluid Force Model

The water-induced force exerting on the ice floe was decomposed into hydrostatic and hydrodynamic. For the hydrostatic force, namely buoyancy force, Catto’s algorithm [35] was adopted in the present numerical model. In Catto’s algorithm, the free surface was assumed to be a plane. The underwater volume was calculated by cutting the immersion part of the floating body. In terms of the hydrodynamic force, only the drag force was considered. The hydrodynamic calculation was a very important part of the simulation, and its calculation results were related to the impact effect. However, accurate calculation of hydrodynamic forces required relatively high computational resources. For example, Computational Fluid Dynamics (CFD) [36] was used to calculate hydrodynamic forces during collision. Our model focused on the overall drag effect, and only considering water drag was a simplification for simulation efficiency. The Morison’s equation was introduced in the model, which had been widely applied in DEM or NDEM numerical simulation [21]:

$$f_{drag}=-0.5C_d{\rho }_w{A}_{proj}{U}_{rel}\left|U_{rel}\right|$$

(13)

$${\tau }_{drag}=-0.5C_a{\rho }_w{l}_{avg}^3{A}_{sub}\mathsf{\Omega }\left|\mathsf{\Omega }\right|$$

(14)

where $f_{drag}$ and ${\tau }_{drag}$ represented the drag force and moment, respectively. Likewise, $C_d$ and $C_a$ were the drag coefficients of force and moment. ${\rho }_w$ was the water density. Ω was the rotational velocity of the ice floe. $U_{rel}$ was the relative velocity between the ice floe’s velocity and water velocity at the buoyancy center. $A_{sub}$ was the submerged area of the ice floe and $A_{proj}$ was the projection of the submerged area $A_{sub}$ along the $U_{rel}$’s direction. The negative projection area was not included in the sum of $A_{proj}$, which meant that the drag force is not applied to the shaded surface. $l_{avg}$ was the average length which was set as the average of the length, breadth and thickness of the ice floe.

2.2. Ice-Breaking Modelling—“Meta Modelling”

The above-mentioned Bullet code can handle the multibody interaction process, but the ice-breaking process cannot be handled by the original version and the specific numerical models for the ice breaking have to be integrated. In this study, a data-driven model called “Meta-Modelling” [25] is adopted and the technique details are discussed in this section.

“Meta-Modelling” is essentially a data-driven model. More specifically, the whole structure-ice interaction process is firstly decomposed into several simple yet representative scenarios. The ice-breaking process for these simple scenarios can be simulated by state-of-the-art methods such as FEM (Finite Element Method). After that, the ice parameters (e.g., contact force, ice geometry, ice properties) and the ice breaking parameters (e.g., the stress in the field, location of the bending failure) are collected and fitted by an Artificial Neural Network. As a result, the complicated ice-breaking process under various conditions can be easily interpolated based on the results of typical scenarios, without the need of conducting the time-consuming simulation. Meta-modeling can maintain the accuracy of the local scale simulation while being computationally efficient on the global problem.

2.2.1. Crushing Force Consideration

According to the study of Lubbad et al. [20], the impact force on a two-dimensional plane can be calculated as follows:

$$F_n={\sigma }_c\cdot A\cdot R$$

(15)

Among them, ${\sigma }_c$ is the uniaxial compressive strength of ice; A is the contact area, which is taken as two-dimensional projection intersection area; R is the contact coefficient reflecting the non-uniformity of full-size contact area. In the calculation of physics engine, time step is introduced into the calculation of contact force. Therefore, when the time step is small, there is abnormal force, which does not conform to the physical law. Here, by limiting the upper limit impulse to limit the contact force calculation of abnormal phenomenon, the upper limit value $I_{limit}$ is set as [21]:

$$I_{limit}=\mathsf{\Delta }t\cdot {\sigma }_c\cdot A\cdot R$$

(16)

In addition, when the broken ice debris is sandwiched during the collision between the structure and the level ice, because the mass of the structure and the level ice is relatively large, the crushed ice debris has very large contact force. Therefore, the ice block was examined in the simulation program and its failure condition was defined as whether its stress is greater than the uniaxial compressive strength of the ice ($\sigma >{\sigma }_c$). Once ice blocks are in failure condition, they would be deleted.

2.2.2. Ice-Breaking Model

The bending breaking pattern of level ice are considered in this section. According to Li [25], based on the observation of the structure-level ice interaction process, the irregular ice edge profile after the falling off of fractured ice can always be characterized as a common shape, as shown in Figure 2. So, the approximate trapezoid can be described in 3 parameters, i.e., the wedge depth b, the contact length c and the wedge angle θ. Point B is the collision point, and point A is the max principal stress point.

• Empirical formula usage

The ice-breaking process is then calculated by the following formula yielded by an Artificial Neural Network based on a pre-calculated database [25]:

$$\mathrm{y}=W_2\mathrm{S}\left(W_{1}X+b_1\right)+b_2$$

(17)

where $W_1$ and $W_2$ are weights, matrices of 10 × 5, 10 × 1, respectively; $b_1$ and $b_2$ are biases; X is the normalized input vector; y is the normalized output; $\mathrm{S}$ is a sigmoid function representing a nonlinear relationship.

According to the model of Li [25], the input vector X is associated with the non-dimensional parameters $\mathrm{b}/l_c$, $\mathrm{c}/l_c$, θ, ${\mu }_1$, ${\mu }_2$ while the outputs y is described by the non-dimensional parameters ${\sigma }_{Maxprinc}/p$ and $l_{break}/l_c$. As shown in Figure 2: $\mathrm{b}/l_c$ is the non-dimensional wedge depth, $\mathrm{c}/l_c$ is the non-dimensional contact length, θ is the wedge angle,${\mu }_1$ is the ratio between the out-of-plane force and the in-plane normal force, and ${\mu }_2$ is the friction coefficient for the in-plane friction force.

The characteristic length $l_c$ is calculated as:

$$l_c=\sqrt[4]{\frac{E{h}^3}{12\rho g(1-{\nu }^2)}}$$

(18)

where h is ice thickness, ρ and ν are ice density and Poisson’s ratio, respectively.

Substituting $\mathrm{b}/l_c$, $\mathrm{c}/l_c$, θ, ${\mu }_1$, ${\mu }_2$ into Equation (17), we can obtain the maximum value of the maximum principal stress ${\sigma }_{Maxprinc}/p$ and the non-dimensional fracture length $l_{break}/l_c$. The stress ${\sigma }_{Maxprinc}$ calculated at the maximum stress point in Equation (17) is compared with the ultimate bending stress of sea ice ${\sigma }_f$. If the maximum principal stress is greater than the ultimate bending stress ${\sigma }_f$, the circumferential crack forms at the maximum stress point. The breaking length $l_{break}$ represents the distance between the midpoint of the fracture edge and the circumferential crack, as shown in Figure 2. The breaking pattern is set as a circle with a radius of $l_b$.

According to Nevel [37], an intact level ice can be divided into multiple ice wedges, then calculate the stress of each wedge to obtain ice wedge fragment amount when ice breaking happened. To simplify breaking ice wedge fragment calculation, the ice wedge amount is randomly taken to be 3–5 from Nevel [38].

2.

Obtaining input parameter

The extraction of input parameters can be divided into three parts. The first step is to recognize a trapezoid such as a polygon from the collision point. Second, gain b, c, and θ from the polygon. At last, get ${\mu }_1$ and ${\mu }_2$. After obtaining all input parameters, they are substituted into Equation (17).

This step needs projection on the waterplane of the level ice, it named as 2d ice. With 2d ice, we could obtain two-dimensional intersection pattern; the recognition algorithm is shown in Figure 3.

Generally, the recognition of most ice edges is feasible, but there are still a few edges which are too complex to identify the desired trapezoid patterns accurately. Therefore, some distorted or non-conforming shapes are re-identified, and these shapes are replaced by two-dimensional intersecting shapes. Intersecting graphics extraction algorithm is shown in Figure 4.

Then, ${\mu }_1$ and ${\mu }_2$ are the ratios of out-of-plane forces to in-plane normal forces and friction coefficients of in-plane frictional forces, which can be obtained by extracting the values of contact forces by ${\mu }_1=f/N$, ${\mu }_2=P/N$. The tangent line between trapezoid and cone is taken as the y direction, and the force in this direction is f. N is the horizontal force in vertical direction, and P is and the normal upward force.

However, the values of $\mathrm{b}/l_c$, $\mathrm{c}/l_c$, θ, ${\mu }_1$, ${\mu }_2$ obtained by the method above may exceeds the range of database of Equation (17). Because the extrapolation ability of the artificial neural network is not ideal, the parameter is set to the minimum input value of the database when the value is too small. As for the condition of exceeding the maximum input value, a shorter time step can avoid it from happening.

3.

Parameter adjustment

Because the setting of this work is different from that of Li [25], the adjustment calculation of the stress parameters should be added to the calculation of the maximum principal stress according to the scale ratio [24]:

$$\sigma \propto \frac{p*E^{0.25}}{h^{0.25}}$$

(19)

where p is the pressure at the point of collision, E is the elastic modulus of ice, and h is the thickness of ice.

3. Results and Discussion

3.1. Validation

3.1.1. Simulation Setup

The developed model is validated by the cone structure-level ice interaction problem reported in Lu et al. [12]. The level ice is set 50 m (Length) × 50 m (Width) in size, and the initial distance from the cone center to the ice edge is 7 m. The calculation duration is 20 s. The rest parameters are as given in

Table 1. Simulation setup.

Parameter	Value
Ice thickness (m)	0.33
Ice field length (m)	50
Ice field width (m)	50
Cone waterline width (m)	13.6
Cone angle (deg)	61.4
Cone height (m)	9.08
Cone height under water (m)	5.68
Interaction velocity (m/s)	1.202
Cone density (kg/m3)	7830
Ice flexural strength (kPa)	539
Young’s modulus of ice (GPa)	0.35

:

The present model has difficulty in solving the interaction between structure and intact level ice, so we set an initial radial and circumferential crack at the initial stage of the simulation. The initial failure pattern is shown in Figure 5.

3.1.2. Parameter Choice

The following parameters are adopted in all the simulations in this paper.

• Uniaxial compressive strength of ice

According to Luo et al. [39], the uniaxial compressive strength of ice is 0.36–3.67 Mpa at −0.7–7.5 °C. So, in Equation (15), uniaxial compressive strength of ice is taken as 1 Mpa in simulation. R is generally taken as 0.5 according to Valanto [40].

2.

Friction coefficient

Friction coefficient is same as Lu et al. [41]. Ice structure friction coefficient is 0.027, obtained by standard test procedures in HSVA [42]. Ice–ice friction coefficient is 0.3, taken from Paavilainen et al. [43] and Sukhorukov et al. [44].

3.

Time step

To validate the choice of time step used in simulation, the impact loading curves under different time steps are analyzed. The interaction between cone and ice under the condition presented in Section 3.1.1 (ice not breaking) is simulated. Through the comparison of impact forces with different time steps (see Figure 6), we can see that ice loads behave similarly with different time steps, which means the impact process is not sensitive with the different choice of time step. The choice of 0.01 s time step can achieve the balance between accuracy and efficiency of this problem.

3.1.3. Validation Simulation

It can be seen from Figure 7 that the oscillation pattern of ice resistance curve is similar compared with the original experiment result reported in Lu [12]. The mean resistance obtained by the numerical model in Lu [12] is 110 kN and the experiment result is 111.8 kN. The mean resistance obtained by the proposed model is 112.93 kN, the difference is about 1.01%. In terms of the frequency and local peaks of the resistance time history, the numerical and experimental results indeed present a certain level of difference. However, this is not only the case for our results; actually, a similar difference between Lu [12]’s simulation and this experimental results can also be observed in ref. [12]. This could be caused by the randomness of the ice properties, as well as the complicated interaction between crushing and bending breaking patterns. Further refinement of the numerical model is planned in the near future. The crack path behaves similarly, as shown in Figure 8.

3.2. Parametric Analysis

The study of ice force on conical structure is useful in offshore structure design. The parametric study of cone ice-breaking can also be applied to ship ice-breaking process. It is significant that making the ship have better ice-breaking performance and safety by adjusting the parameters. Therefore, the research on the cone and ice thickness with different parameters and velocities can benefit ship design and navigation.

The research of cones with different parameters and ice thickness and velocity can bring some reference value for ship design and navigation. Based on Section 3.1, several parameters of cone numerical calculation were analyzed, including: cone angle, cone motion velocity, cone waterline width and ice thickness. Except the parameters for parametric analysis, the parameters in 3.1 are used for other parameters.

The examples for parametric analysis are as

Table 2. Parametric analysis simulation setup.

Parameter	Setting Groups
Cone angle (deg)	50, 55, 65, 70
Cone motion velocity (m/s)	0.5, 1.0, 1.5, 2.0
Cone waterline width (m)	10.88, 12.24, 14.96, 16.32
Ice thickness (m)	0.3, 0.4, 0.5

:

The cone angle here refers to the intersection angle between the cone surface and water, which is shown in Figure 9:

3.2.1. Influence of Cone Angle

The results of mean resistance varying with cone angles are shown in Figure 10:

The comparison of resistance calculation shows that the ice-breaking resistance increases gradually with the increase of cone angle and larger angle generates higher increase rate of resistance, which is consistent with the results in many experiment [45,46,47]. Figure 11 shows the time history of resistance under different cone angles. It can be seen that the results with higher angles tend to show more high frequency oscillation components.

3.2.2. Influence of Cone Motion Velocity

The results of the mean resistance varying with cone velocity are shown in Figure 12:

The results show that the ice-breaking resistance increases with the increase of velocity. Figure 13 shows the distribution of broken ice pieces around cone structure under different velocities. At low velocity, the broken ice pieces tend to accumulate around the waterline area, while the ice pieces will stretch to larger areas below and behind the cone structures.

Similarly, to the cases with larger cone angle, time history with larger velocity show more high frequency components as well (see Figure 14).

In the case of certain ice thickness and waterline width (see Figure 15), the influence of different velocities on different cone angles is also different:

Several groups of calculations are carried out for different velocities and cone angles under the same ice thickness. It can be seen from the results in Figure 15 that, at the condition of 50°, the velocity of cone movement has little influence. With the increase of cone angle, the influence of velocity on resistance becomes more significant. More specifically, when the angle reaches 60°, the influence of velocity at the same cone angle begins to increase, and the positive correlation between high velocity and high resistance is strengthened. As for the cone moving at the same velocity, cone angle has little influence at low velocity. The cone angle has little influence and can only slightly reduce the ice-breaking resistance at 0.5 m/s. However, at velocities up to 2.0 m/s, the resistance of 50° cone is 49.1% of 70° cone.

3.2.3. Influence Analysis of Cone Waterline Width

The results of the cone waterline width group are shown in Figure 16:

The results show that the cone radius does not have great influence on the ice-breaking resistance (horizontal). From the ice fragment formation scene (see Figure 17), as the radius of the waterline increases, the accumulation of ice fragments at the bottom of the cone due to breakage increases, resulting in the force on the cone being mainly concentrated in the Z-direction (vertical) rather than in the X-direction (horizontal). According to Figure 16, the Z-direction(vertical) mean ice resistance increases with the increase in the waterline width, and the overall resistance also presents an upward trend.

3.2.4. Influence of Ice Thickness

The results of the calculation of the ice thickness group are shown in Figure 18:

Through the comparison of the resistance (see Figure 18), it can be seen that the breaking resistance increases with the increase in the ice thickness. The resistance increases linearly with the thickness from 0.2 m to 0.4 m. When the thickness reaches 0.5 m, the breaking resistance increases sharply. In the time history curve (see Figure 19), the resistance fluctuates little when the ice sheet is thin. With the thickness increasing to 0.4 m, the resistance fluctuations begin to intensify, and the peak value becomes obvious. It is obvious that the ice thickness has a great influence on the ice-breaking resistance during the ice-breaking process.

4. Conclusions

In this paper, a data-driven model based on the Non-smoothed Discrete Element Method (NDEM) and Artificial Neural Network (ANN) method was used for the simulation of level ice and structure interaction problem. Using the NDEM method, the collision can be accurately simulated by large time steps. At the same time, the ANN method can be used to approximate the solution of level ice fracture without direct simulation.

The application of the NDEM-ANN method can draw several conclusions:

(1)

The method of calculating level ice failure by physics engine was realized. This method got rid of the slow and complex calculation of equation solution and made the calculation results close to the high precision numerical simulation. Moreover, it is in good agreement with the experimental results.

(2)

By using the proposed method, the influence of various parameters such as cone angle, navigation velocity and ice thickness on ice resistance were investigated. It was found that the influence of velocity on ice resistance ice resistance increased greater in large cone angle than small cone angle.

The simulation method in this paper considered simulating efficiency for a simple interaction scenario of level ice and structure, and the improvement of simulation in further works are aiming to simulate more complicated ice–structure interaction scenario such as ice floe interaction.