Parametric Study on Strength Characteristics of Two-Dimensional Ice Beam Using Discrete Element Method

: Strength characteristics of a two-dimensional ice beam were studied using a discrete element method (DEM). The DEM solver was implemented by the open-source discrete element method libraries. Three-point bending and uniaxial compressive tests of the ice beam were simulated. The ice beam consisted of an assembly of disk-shaped particles with a particular thickness. The connection of the ice particles was modelled using a cuboid element, which represents a bond. If the stress acting on the bond exceeded the bond strength criterion, the bond started to break, explaining the cracking of the ice beam. To ﬁnd out the effect of the local parameters of the contact and bond models on the ice fracture, we performed numerical simulations for various bond Young‘s modulus of the particles, the bond strength, and the relative particle size ratio. and the mechanical properties of the simulated ice was investigated. The effects of the bond Young’s modulus ( E b ), the bond strength ( σ b ), and the relative particle size ratio ( h / d ) on the mechanical properties of the simulated ice were studied. The simulated Young’s modulus ( E s ) increased in proportion to the bond Young’s modulus ( E b ) , while the ﬂexural strength was not affected by the bond Young’s modulus. This was because that the normal stiffness of the bond was proportional to the bond Young’s modulus. The simulated ﬂexural strength ( σ f ) and the compressive strength ( σ c ) increased linearly with respect to the bond strength ( σ b ). The normalized strength was constant for given simulation conditions. This was the key relationship for predicting the simulated ﬂexural strength. The inﬂuence of the relative particle size ratio ( h / d ) was investigated to obtain the converged value. The simulated Young’s modulus for the three-point bending tests tended to converge to the constant ratio ( E s / E b ) of 1.5 at a sufﬁciently high h / d . The simulated Young’s modulus for the uniaxial compressive tests increased as the bond Young’s modulus increased. As per the results, the normalized ﬂexural strength and compressive strength were 1.55 and 4.04, respectively. From the ice modelling in the Bohai Sea, the obtained parameters of the contact and bond models could be used to assess and predict the mechanical properties of the sea ice. In future research, ice breaking by a ship will be carried out.


Introduction
As the sea ice area coverage in the Arctic Ocean shrinks over the years due to climate change, the operation of ships in the Arctic Ocean have been issued. To operate a ship in the sea ice area, accurate prediction of ice breaking performance is needed. Studies for an ice breaking load estimation have been carried out using empirical, analytical, and numerical methods [1]. The empirical method used measured data in full-scale trials and model-scale experiments. Formulations based on full-scale data produce highly reliable methods for the ice load prediction, while there is a limitation with regard to obtaining definitive data on properties, e.g., thickness, strength, and friction [2,3]. The model tests in the ice model basin have the advantage of being able to measure the ice load under various operating conditions compared to the full-scale measurements, but it is very difficult to evaluate the ice performance with various design factors due to cost and time issues. Consequently, there is an increasing need for numerical models to predict the accurate ice load with regard to various sea ice conditions in the initial design stages of Arctic offshore structures [4].
Numerical approaches for ice modeling could be divided into a continuous method, such as the finite element method (FEM) and finite discrete method (FDM), and a discontinuous method, such as the discrete element method (DEM). The FEM is a dominant numerical tool in modelling continuous materials both in the linear and nonlinear range of deformation. It has some drawbacks when simulating macrocracks or fragmentation of the material [5,6]. On the other hand, the DEM easily generates realistic macrocrack patterns and material fragments given its discontinuous nature [5,7]. The DEM is known to better simulate the propagation of an ice crack and fracture behavior because the connection between the particles can be modelled [7]. The DEM is widely applied to ice modelling, ice breaking, and ice-structure interaction issues [8]. To produce physical deformabil-ity and strength of ice by the DEM, researchers require extensive careful calibration of parameters [5,9].
In the DEM, each individual particle that contains properties of the ice could be described as several shapes such as a disk, a sphere, and a polyhedron. The DEM could simulate various ice conditions, e.g., ice floes, level ices, and ice ridges by modeling the ice as individual particles or an assembly of particles [10][11][12][13]. For the ice floes, studies on an interaction between ships or offshore structures and the ice floes using disk-shaped particles that assumed unbreakable ice were carried out [14][15][16][17][18]. For the level ices, some studies used bonds between particles to simulate contacts and cracks in the level ice [19][20][21]. The bond between two particles was broken when the maximum force acting on the bond exceeded a criterion, which could explain the crack and fracture phenomenon [22].
The ice breaking load in the DEM was highly dependent on the mechanical properties of ice [23][24][25]. The bond Young's modulus, flexural strength, and compressive strength were related to parameters of contact and bond models. It is necessary to define the parameters of the models that affect the mechanical properties of ice and to figure out the relationship between the parameters and the mechanical properties of ice [23,26,27].
In the present study, the DEM was selected and applied to the three-point bending test and the uniaxial compressive test. For the simulations, the open-source discrete element method libraries were used [13,28,29]. The relationship between the mechanical properties of the simulated ice and the parameters associated with the contact and bond models was investigated.
The present paper is organized as follows. Section 2 describes numerical modeling including the governing equation for particle contact and bond models, and the parameters for ice modeling. Section 3 presents the results and discussion for the parametric study. Finally, in Section 4, concluding remarks are provided.

Contact Model
The translational and rotational motions of a particle with a constant mass (m) could be expressed by Newton's second law as follows: where the subscripts c and b represent the contact and the bond, respectively. The mass (m = ρh d πr 2 ) is calculated as the disk-shaped particle that has a certain thickness (h d ). u is translational velocity of the particle. F c and F b are the contact force and the bond force, respectively. mg is the gravitational force. I is the inertia moment of the particle, and ω is the angular velocity of the particle. The particle rotation is calculated by the sum of the torque (M c ) caused by collisions and the moment (M b ) by the bond that transferred to the particle. The contact force (F c ) by the collisions between particles could be simplified using a spring-dashpot model. The contact force consisted of the normal (F c,n ) and tangential (F c,t ) forces, which could be expressed as where the subscripts b and t represent normal and tangential components, respectively. k and γ are the spring stiffness and the damping coefficient, respectively. δ and . δ represent the overlap and the relatively velocity between two contact particles, respectively. The normal overlap (δ n ) was a function of the particle radius and the distance between the particle centroids at the current time step. The tangential overlap (δ t ) was calculated by the incremental distance in the tangential direction after an incidence of contacting. The maximum tangential force was limited by the Coulomb friction law with the friction coefficient (µ) and normal contact force (F c,n ). The spring stiffness and damping coefficient were dependent on the particle properties and selected contact law. The contact law could be divided into a linear contact law (Hook contact model) and a non-linear contact law (Hertz contact model). In this study, Hertz contact model was used [30].

Bond Model
The bonding method between particles was based on the parallel bond method [22]. Figure 1 shows two disk-shaped particles bonded by a cuboid element in the two-dimensional model [13]. Each bond was characterized by the thickness (h b ), width (2R b ), length (l b ), radius factor (λ R ), bond Young's modulus (E b ), and the normal stiffness (k b,n ) and shear stiffness (k b,t ) per unit area. Here, the bond thickness (h b ) and the disk-shaped particle thickness (h d ) exist in the z-direction, and both values are identical. The width and length of the bond could be expressed as where the radius factor (λ R ) ranged from 0 to 1. k b,n and k b,t depend on the bond Young's modulus (E b ), bond length (l b ), and the ratio of the normal to shear stiffness (λ ns ), k b,n and k b,t could be expressed as The forces (F b ) and torques (M b ) acting on the particles connected with a bond caused a relative displacement and rotation of those particles.
Before the start of the simulation, the forces (F b ) and moments (M b ) in the bond were initialized to zero. The relative displacement and rotation of the particles were governed by the linear elastic material and calculated as follows: where S b = 2R b h b is the cross-sectional area of the bond, and ∆u n and ∆u t represent the parallel and perpendicular components to the bond axis, respectively. The subscripts tw and bn indicate a twisting and a bending, respectively. ∆w n and ∆w t represent the normal and tangential components, respectively.
, and J b = I b,y + I b,z represent the moments of inertia. J denotes the polar moment of inertia. The torque could be expressed as follows: In the 2-D model, F b.t contributed to the particle's rotation around the z-axis. The twisting moment (M b,tw ) was zero, and only the z-component of the bending moment (M b,bn ) was not zero.
To simulate the ice fracture, we used the elastic beam theory. The normal bond strength (σ b ) and the shear bond strength (τ b ) could be expressed as a maximum value as follows [22]: where n is a unit normal vector, H is a 3 × 3 projection matrix from the 3-D space to the x-y plane, H 11 = H 22 = 1, and all remaining elements are zero. The negative sign in Equation (16) indicated that only tensile normal stress could cause bond failure. The normal stress acting on the bond element includes the bonding moment due to the relative rotation of the particle and the normal force. The bond began to break if at least one of the maximum stress components acting on that bond exceeded the bond strength as follows:

Parameters for Ice Modeling
The parameters for the contact model were the particle radius (r), density (ρ), Poisson ratio (ν), friction coefficient (µ), and restitution coefficient (e). Some parameters, namely, Poisson ratio (ν), friction coefficient (µ), and restitution coefficient (e) for the contact model, were less crucial to approximate the strength and stiffness of the continuous material with the inter-particle bonds [23,24]. The parameters of the bond model considered in this study were the bond Young's modulus (E b ), width (2R b ), length (l b ), thickness (h b ), ratio of the normal to shear stiffness (λ ns ), normal strength (σ b ), and shear strength (τ b ). The normal and shear bond strengths were applied to be identical [31][32][33]. In the bending problem, the effect of the contact point of the load was not significant [25]. Within a particle assembly, all particle properties were assumed to be distributed uniformly. The distance between particles was entirely covered by the bond. Thus, the bond length (l b ) and the bond width (2R b ) were regarded as the particle diameter (d = 2r). The radius factor (λ R ) was 1. For both the three-point bending tests and the uniaxial compressive tests, the particle diameters were set at d = 2.5, 5, 10 and 20 mm. The particle size effect was studied using the relative particle size ratio (h/d). Here, h indicates the cross-sectional height of the ice specimen in the three-point bending test and uniaxial compressive test as shown in Figure 2. h/d is proportional to the number of particle layers of the ice specimen. The parameters related to the contact and the bond of ice modeling are shown in Table 1.

Three-Point Bending Tests and Uniaxial Compressive Tests of Ice Beam
To investigate the relationship between the parameters for ice modeling and the simulated mechanical properties of ice, we simulated the three-point bending test and uniaxial compressive test of an ice beam. Figure 2 shows the dimensions and test setup of the ice specimen. In the three-point bending simulation, the distance (l) between the two fixed supporting points was 500 mm and the length of the ice specimen (L) was 700 mm. The width (b) and height (h) were set to be the same at 70 mm. A constant downward vertical load with a constant rate of 0.002 m/s was applied at the middle point on the top side of the ice beam. The supporting and loading points also were modelled by a disk-shaped particle.
In the uniaxial compressive tests, the distance (l) between top and bottom plates was 250 mm. The width (b) and height (h) were set to be the same at 100 mm. The bottom plate was fixed, and the constant downward load of 0.002 m/s was applied to the top plate. The bottom and top plates were modelled by a disk-shaped particle.
Sea ice is quasi-brittle heterogenous and anisotropic. In the present study, for simplicity, the sea ice was assumed to be homogeneous, anisotropic, and elastic brittle [24,25,32]. The ice beam was represented by the particle assembly with a regular arrangement such as the Hexagonal Close Packing (HCP) [24,25,32]. This arrangement leads to anisotropy but yields a less realistic crack pattern as compared to the randomized packing [27]. Despite the limitations of the regular arrangement, it could lead to a consistent and predictable mechanical behavior, which was beneficial for establishing the relationship between the parameters for ice modeling and the simulated mechanical properties of ice [20,[24][25][26]32].
In the modeling regarding the level ice for the ice-structure interaction issues, the critical mechanical properties were the bond Young's modulus, flexural strength, and compressive strength [34]. The three-point bending and uniaxial compressive tests were carried out to obtain the simulated Young's modulus (E s ), as well as the flexural strength (σ f ) and the compressive strength (σ c ) of the ice beam. The total contact force acting on the loading particle indicated the load applied to the ice beam, while the deformation of the ice beam was expressed by the displacement of the loading particle. The flexural strength and the compressive strength of the ice beam could be calculated as where P max is the maximum load when the ice beam is broken. The simulated Young's modulus (E s ) can be derived from the stress-deflection curve as where the subscripts A and B denote the two arbitrary selected points in the stressdeflection curve.
In the three-point bending and uniaxial compressive tests, the bond Young's modulus (E b ), the bond strength (σ b ), and the relative particle size ratio (h/d) were studied as the main parameters of the contact and bond models. Figure 3 shows the failure process of the three-point bending test. The compressive stress was increased at the upper part and the tensile stress was increased at the lower part of the ice beam until the crack appeared at t = 0.4792 s. It could be observed that the crack occurred near the lower part at t = 0.4794 s. As the compressive stress concentrated near the upper part at t = 0.4796 s, the ice beam broke at t = 0.4800 s. The fracture of the ice beam occurred at the middle point with a gradual increase of the normal stress at the top and bottom sides [33,35].    Figure 5 shows the obtained stress-deformation curve. During the short time after the deformation was started, the stress reached to the maximum stress, which denoted a brittle material behavior. The normal stress of the ice beam was increased until the fracture occurred. The simulated Young's modulus (E s ) could be obtained by the stress and deformation at the two arbitrary points A and B. The maximum stress at the point C indicated the flexural strength (σ f ) and the compressive strength (σ c ). In Figure 5a, the simulated Young's modulus (E s ) and the flexural strength (σ f ) were 1.47 GPa and 1.15 MPa, respectively. As shown in Figure 5b, the simulated Young's modulus (E s ) and the compressive stress (σ c ) were 1.258 GPa and 3.03 MPa, respectively. The ratio of the compressive stress to the flexural stress (σ c /σ f ) was 2.63, which was acceptable for that measured for the Bohai Sea [35]. To validate the bond model, we compared the result for the flexural strength with the DEM simulations [24,26,33] and experimental data of sea ice [35] as listed in Table 2. The DEM simulations [24,26,33] used the parallel bond [22] and the Mohr-Coulomb law for the breaking criteria. For the same bond Young's modulus, the bond strengths were similar. The real sea ice data represented the measured mean flexural strength of sea ice in the Bohai Sea [35]. The simulated flexural strength in the present study was slightly higher than the measured mean flexural strength of sea ice [35] and lower than the DEM results [24,26,33]. All simulated results (E s , σ c , σ f , σ c /σ f ) were within the range of the mechanical properties of the actual sea ice [35]. It could be seen that the simulated ice characteristic represented the mechanical properties of real sea ice to a reasonable level.  (19) and (20), respectively.

Effect of Bond Young's Moduls
The effect of the bond Young's modulus (E b ) was investigated. The relationship between the bond Young's modulus and simulated Young's modulus (E s ) is examined in Figure 6. The bond strength (σ b ) was 1.0 MPa and the relative particle size (h/d) was 7.06. The simulations were performed at three bond Young's modulus. Figure 6a presents the stress-deflection curve. The slope of the stress-deflection curve increased as the bond Young's modulus increased, while the maximum stress was almost consistent. Figure 6b shows the stress-strain curve. The slope of the stress-strain curve increased as the bond Young's modulus increased, while the maximum stress was almost consistent. The material stiffness was proportional to the bond Young's modulus, but the flexural strength (σ f ) and the compressive strength (σ c ) were not related to the bond Young's modulus. Figure 6c shows the ratio of the simulated Young's modulus to the bond Young's modulus (E s /E b ). The normal bond stiffness (k b,n ) was proportional to the bond Young's when the bond length (l b ) was constant. Therefore, a uniform ratio (E s /E b ) was obtained under the given simulation conditions. The ratio (E s /E b ) greater than 1 meant that the ice fracture occurred at a lower deflection than the targeted ice condition. However, the uniaxial compressive tests showed a lower result of the ratio (E s /E b ) than that of the three-point bending tests. This was due to the assumption of equal contact and bond modulus. When the contact and bond modulus were equal, the contact stiffness at a small overlap was much lower than the bond stiffness. Therefore, the simulated Young's modulus was lower in the uniaxial compressive tests where contact was a dominant part of the load-resisting mechanism [25].  Figure 7 shows the stress-deformation curve for various bond strengths (σ b ). The bond Young's modulus (E b ) and the relative particle size (h/d) were 1.0 GPa and 7.06, respectively. The flexural strength (σ f ) and the compressive strength (σ c ) increased proportionally with the bond strength. This was due to the breaking condition where the bond was breaking when the maximum stress exceeded the bond strength. However, the deformation of the ice beam according to the constant vertical load showed a constant behavior until the ice was fractured. This indicated that the simulated material stiffness was not affected by the bond strength.    Figure 9 presents the effect of the relative particle size (h/d) to the ratio of simulated Young's modulus to bond Young's modulus (E s /E b ). The bond Young's modulus and the bond strength were 1.0 GPa and 1.0 MPa, respectively. As h/d increased, E s /E tended to increase slightly in the uniaxial compressive tests, while converging in the three-point bending tests. In Figure 9a, E s /E b converged to 1.5 with h/d > 7.06, while it was low with h/d < 7.06. Following the results shown in Figure 9b, the simulated Young's modulus was proposed as a function of E s /E b and h/d and expressed as The relative particle size ratio (h/d) was dependent on the porosity and the number of particles. In other words, smaller particle sizes reduced a void between the particles and increased the number of particles per unit volume, resulting in a stiffer assembly. Therefore, if the relative particle size was large enough, the simulated Young's modulus converged [23,36]. Figure 10 shows the effect of the relative particle size (h/d) to the simulated ice strength. The normalized strength converged to be 1.55 for the flexural strength and 4.04 for the compressive strength at h/d > 7.06, indicating that the consistent behavior of ice could be predicted when the relative particle size ratio was higher than 7.06.

Ice Modeling in the Bohai Sea
The simulated flexural strength (σ f ) and compressive strength (σ c ) were compared to the actual strength of the Bohai Sea ice by the three-point bending test and the uniaxial compressive test [35]. The measured mean and the maximum flexural strengths were 1.08 MPa and 2.42 MPa, respectively. The measured mean and maximum compressive strengths were 2.56 MPa and 5.98 MPa, respectively. The Young's modulus in the Bohai Sea was assumed to be 1.0 GPa. The two particle sizes, 5 mm and 10 mm, were considered.  Table 3 shows the ice properties used in the simulation. Figure 11 presents the stress-deflection curves. Table 4 shows the simulated flexural strengths and the actual flexural strengths of ice. The simulated strength was evaluated to be close to the actual strength. From the results in Figure 11 and Table 4, we can confirm that the parameters determined by the relationships obtained in the parametric study (the E s /E ratio and the normalized strength) were able to achieve the desired flexural strength and compressive strength. This feature is crucial for the simulation of ice-structure interactions.

Conclusions
In the present paper, the three-point bending simulation was performed using the DEM. The open-source DEM libraries were used for the simulation [13,28,29]. The focus was to figure out the relationship between the model parameters and the mechanical properties of ice because there is a unique relationship in the regular particle arrangement [27]. This relationship is important for predicting an ice breaking load caused by the interaction with ships or offshore structures. To simulate the ice fracture, we used the bonds to connect the particles. The ice fracture occurred when the maximum stress acting on the bond exceeded the bond strength. The relationship between the parameters for ice modeling and the mechanical properties of the simulated ice was investigated. The effects of the bond Young's modulus (E b ), the bond strength (σ b ), and the relative particle size ratio (h/d) on the mechanical properties of the simulated ice were studied. The simulated Young's modulus (E s ) increased in proportion to the bond Young's modulus (E b ), while the flexural strength was not affected by the bond Young's modulus. This was because that the normal stiffness of the bond was proportional to the bond Young's modulus. The simulated flexural strength (σ f ) and the compressive strength (σ c ) increased linearly with respect to the bond strength (σ b ). The normalized strength was constant for given simulation conditions. This was the key relationship for predicting the simulated flexural strength. The influence of the relative particle size ratio (h/d) was investigated to obtain the converged value. The simulated Young's modulus for the three-point bending tests tended to converge to the constant ratio (E s /E b ) of 1.5 at a sufficiently high h/d. The simulated Young's modulus for the uniaxial compressive tests increased as the bond Young's modulus increased. As per the results, the normalized flexural strength and compressive strength were 1.55 and 4.04, respectively. From the ice modelling in the Bohai Sea, the obtained parameters of the contact and bond models could be used to assess and predict the mechanical properties of the sea ice. In future research, ice breaking by a ship will be carried out.
Author Contributions: Conceptualization, S.S. and S.P.; simulation, S.S. and W.J.; formal analysis, S.S. and W.J.; writing-original draft preparation, S.S. and S.P.; writing-review and editing, S.S. and S.P.; visualization, S.S. and W.J.; supervision, S.P. All authors have read and agreed to the published version of the manuscript.