1. Introducción
Ice loads are crucial control loads that impose safety limitations on offshore equipment structures in high-latitude marine regions. They constitute vital input parameters for designing marine structures with ice-strengthened features, such as lighthouses, conductor support platforms, and wind turbine foundations, which can be approximated as vertical cylindrical structures. The interaction between ice and cylindrical structures results in substantial ice loads and dynamic responses in the structures. Icebreakers are commonly employed to break level ice into broken pieces as a risk mitigation measure. Broken ice exerts notable effects on structures, such as collision, accumulation, rotation, and friction. Consequently, investigating the interaction loads between broken ice and cylindrical structures carries considerable practical importance.
Recently, there has been global attention toward the prediction and analysis techniques of ice loads. Timco [
1,
2,
3] summarized experimental data related to the impact of floating ice on offshore structures. Experimental data regarding the impact of floating ice on offshore structures were summarized. By taking into account the kinetic energy of floating ice, empirical formulas were derived based on ice mass and velocity to determine the load on structures caused by broken ice. The discrete element method (DEM), based on the discontinuous medium approach, has been utilized to capture the nonuniform, discontinuous, and large deformation characteristics of sea ice. Significant progress has been made in this method in recent years, owing to advancements in computer hardware. Ji et al. [
4,
5,
6] developed a FEM–DEM coupled method that allows for the simulation of the dynamic response of the interaction between structures and ice layers. They employed a GPU parallel algorithm to enhance computational efficiency and developed a rapid generation model for fragmented ice. By utilizing a fast search algorithm among polyhedra and a cohesive-fracture model, they improved the accuracy of analyzing structural ice loads and established a convenient computational tool. Jou et al. [
7] modeled sea ice using the local cohesive discrete element method. They investigated the collision between level ice and inclined icebreaking cones at various angles and compared the results with other DEM simulations.
The cohesive zone model (CZM) has been used to simulate the interaction between level ice and structures [
8,
9,
10,
11,
12,
13]. Theoretical methods were utilized to determine the initial stiffness, reducing the dependence on grid resolution in the cohesive element approach. In a previous study, a numerical model using the FEM–SPH adaptive method was developed to investigate the interaction between conical objects, icebreakers, and level ice. By comparing the model results with experimental data, the study demonstrated the potential of this method in accurately simulating ice fragmentation resistance [
14].
Research on icebreakers has made significant advancements due to the growing demand for Arctic scientific expeditions. Xue [
15] provided a comprehensive summary of the latest experimental techniques in studying ice loads, particularly emphasizing their application to ships, and conducted extensive foundational studies in this area. Dynamic methods have been employed by researchers to examine the impact of factors, including ice thickness and navigation speed, on average ice loads [
16,
17]. Wang developed a time-domain numerical calculation method to assess the icebreaking capabilities of ships by studying continuous icebreaking loads and ship motions [
18]. Furthermore, utilizing LS-DYNA simulations of collision-induced icebreaking contact loads, additional investigations were carried out to analyze the bending moments and structural loads on the ship’s cross-section [
19].
Several researchers have investigated ice loads in broken ice environments through experimental methods. Zong et al. conducted scaled model towing experiments employing a polypropylene ice model to examine the influence of broken ice concentration, ice block shape, and size on ship resistance. Additionally, they also confirmed the viability of employing the polypropylene ice model for studying ice resistance [
20,
21]. Islam et al. [
22,
23] employed model testing methods to investigate ship motion, loads, and maneuvering performance in ice-covered waters. Furthermore, they analyzed the loads on the ice–structure interface, thus affirming the potential of discrete element methods in engineering applications.
Sea ice is composed of discrete floating blocks that exhibit substantial variations in size. Roach et al. [
24] introduced a joint distribution model to characterize the thickness and block size of sea ice, and this model effectively captures the fluctuations in ice thickness, leading to enhanced simulations of ice–structure interaction loads. Li et al. employed measured data from Antarctic research vessel expeditions to derive probability distributions that yield the best fit for ice loads through statistical analysis. These probability distributions form a basis for long-term statistical analysis, thereby facilitating the investigation of ice-induced loads [
25,
26]. Huang [
27] conducted dynamic ice load tests at varying ice velocities to examine the vibration process of ice-induced flexible structures and investigated the mechanisms behind vibration suppression. Their findings suggest that the control mechanism for steady-state vibration in ice-induced flexible column structures can be attributed to the periodic impact of ice loads. Moreover, the interaction between the structure and fragmented ice, along with the periodic destruction of the ice layer, exacerbates the influence of structural vibration.
Considering all the abovementioned studies, the analysis of ice loads on floating structures operating in icy environments highlights the superior accuracy of broken ice loads in representing the actual load characteristics, owing to their inherent stochastic nature. However, predicting these loads poses significant challenges. To gain further insights into the mechanisms of broken ice loads and explore the methods for predicting the interaction between broken ice and structures, this study specifically investigates the interaction between a cylindrical structure and broken ice. By employing model similarity criteria, experiments were conducted in a wave flume at room temperature using an ice model made of polypropylene to assess the resistance of the cylindrical structure to broken ice. The experiments measured the structural loads under different ice velocities while maintaining a fixed density of broken ice. Leveraging high-performance discrete element technology, a simulation method was developed to estimate the resistance of the cylindrical structure to broken ice, and a sensitivity analysis was conducted to examine the impact of discrete element density on the resistance results. The structural resistance corresponding to the ice velocity employed in the experiments was predicted, and the findings were compared with the model experiments.
3. Numerical Method
3.1. DEM Theoretical Approach
- a
Ice Discrete Element Model
Sea ice is conceptualized as an assembly of discrete particles, wherein individual particles exhibit a relatively independent motion but engage in interactions upon contact. To simulate the collision process between particle units, a spring–damper–slider model [
5] was employed, as depicted in
Figure 3. This modeling approach allows for a comprehensive representation of the dynamic behavior of sea ice, capturing the effects of particle interactions and their subsequent motion.
The forces
generated by each unit are calculated by summing the normal contact force
and the tangential contact force
.
The normal contact force
can be expressed as the sum of the normal elastic force
and the normal viscous force
, that is,
where
is the normal contact stiffness coefficient;
is the normal overlap between particles;
is the normal damping coefficient; and
represents the relative normal velocity.
where
is the effective elastic modulus;
is the effective radius;
is the damping coefficient;
represents the effective mass; and
is the restitution coefficient of the unit.
The tangential contact force
can be expressed as the sum of the tangential elastic force
and the tangential viscous force
, that is,
In these equations,
represents the friction coefficient;
is the tangential overlap;
is the Poisson’s ratio;
is the relative tangential velocity;
represents the tangential damping coefficient; and
represents the maximum tangential overlap, which can be determined using the following equation:
- b
Coupled method of DEM and FEM
In this study, the coupled DEM–FEM algorithm proposed by Ji [
28] was applied to facilitate parameter transfer between the two computational domains in the coupled DEM–FEM model established for investigating the interaction between ice and a cylindrical structure. The collision detection between fragmented ice and the cylindrical structure was simplified by considering the contact between spherical and triangular elements, with the theoretical algorithm elaborated in [
29]. The conversion of particle interaction loads on shell elements to nodal loads was accomplished using the shape function method.
Figure 4 illustrates a schematic representation of the collision contact between a single particle and a specific triangular element on the outer surface of the cylindrical structure. The calculation of the equivalent nodal forces for triangular elements is presented as follows:
In the above equation, denotes the interpolation matrix for the element contacted by particle a, p represents the number of contact points exerted on the triangular element by particles, denotes the forces and moments exerted by particles on the triangular element, and represents the transformation matrix.
Figure 4.
Contact force–moment conversion diagram of triangle element.
Figure 4.
Contact force–moment conversion diagram of triangle element.
3.2. Numerical Model
- a
Structure and Broken Ice Model
In order to compare the obtained results with experimental data, a vertical cylindrical structure with a diameter of 0.55 m was selected as the subject of analysis in this study. The outer surface of the cylindrical structure was discretized using triangular elements to comply with the internal contact algorithm utilized in this study, as depicted in
Figure 5. The total element number of the cylindrical structure was 1400.
In this research, the SDEM computational program was employed to simulate the load interaction between ice and the structure. We disregarded the influence of waves and exclusively focused on the interaction between broken ice and the structure under steady flow velocity conditions. In order to mitigate the influence of boundary effects on the numerical simulation results, a computational domain measuring 100 m in length and 40 m in width was established. The density of the broken ice loads was set to 70%.
Figure 6 depicts the discrete element coupled model that elucidates the interaction between the ice and the cylindrical structure. In order to ensure the reliability of the computational results, a minimum increment step of 10
−6 s was used.
- b
Discrete Density Comparison
In discrete element algorithms, the term “discrete density” refers to the quantity or distribution density of particles within a specified region. Variations in discrete densities can significantly influence computational results and computation time. In order to assess the impact of discrete element density distribution on computational outcomes, a comparative analysis was conducted involving various particle distribution configurations, namely two layers, three layers, and four layers, with respect to a standardized ice thickness of 0.2 m. The number of discrete elements was 169,282, 533,264, and 1,220,996, respectively. The time-domain curves of the ice load are shown in
Figure 7.
To calculate the statistical values of random loads, three distinct statistical parameters were considered: the mean ice load, the significant ice load, and the maximum ice load. The mean ice load signifies the average value derived from all ice force measurements, providing an overall representation of the data. In contrast, the significant ice load specifically pertains to the ice force magnitude throughout the entire observation period. The maximum ice load refers to the highest recorded ice force value observed over the complete duration of the study.
The statistical method employed for significant wave height was utilized and adapted for the significant ice load. The step-by-step procedure was as follows:
- (1)
Time-history data pertaining to ice forces, obtained through either measurements or calculations, were acquired;
- (2)
The peak ice forces, corresponding to the maximum values within the dataset, were determined;
- (3)
The filtered ice forces were arranged in descending order and organized from the highest value to the lowest;
- (4)
The significant ice force was calculated by averaging the top one-third of the sorted ice forces.
The resultant ice force values for each computational scenario are succinctly presented in
Table 3.
Based on the findings presented in
Table 3, it is evident that the statistically obtained ice load values exhibit minimal discrepancies across the three density distributions. However, there is a notable contrast in computational time for simulations of the same duration, and the computational times for the three different discretization densities are 13 h, 60 h, and 144 h, respectively.
Nonetheless, this outcome is intricately linked to the scale of the analyzed sea ice. Enhancing the discrete density has the potential to yield more precise simulation outcomes. By augmenting the number or density of particles, the discrete element algorithm can aptly capture the interactions and behavior among particles within the system. Consequently, conducting a sensitivity analysis of the grid scale becomes imperative within the realm of engineering research.
In order to enhance computational accuracy and optimize computational resources, a two-layer discrete element approach was effectively employed in this specific case to simulate the resistance of fragmented ice for an ice thickness of 0.20 m.
4. Results and Discussion
4.1. The Results of the Ice–Structure Test
- a
Experimental phenomena
The experimental phase of the ice–cylindrical structure load is illustrated in
Figure 8. Initially, in the experiment, broken ice approached the cylindrical structure at a constant velocity influenced by the flow load. Throughout the experiment, the broken ice underwent rotation or experienced friction upon colliding with the structure, leading to a subsequent alteration in its original direction of movement and dispersion toward the sides or rearward. This interaction resulted in the accumulation and clustering of broken ice. The ice–structure interaction had a brief duration, displaying periodic impact characteristics. Downstream, particularly at the rear of the cylindrical structure, a relatively unobstructed water pathway developed, as depicted in
Figure 8c, exhibiting a notable contrast to the initial configuration. The interaction between the ice model and the structure exhibited evident random phenomena, consistent with relevant studies in the literature [
20].
At higher ice velocities, the frequency of collisions between the broken ice and the cylinder intensified, causing the broken ice to swiftly pass the cylindrical structure. Consequently, some of the broken ice accumulated behind the cylinder, as illustrated in
Figure 8d.
- b
Ice load characteristic analysis
Figure 9,
Figure 10 and
Figure 11 illustrate the time-domain curves of the ice load determined based on scaled experimental data under three distinct velocities: 0.05 m/s, 0.075 m/s, and 0.125 m/s. Taken together, these experimental findings reveal the prominent characteristics of random loading in the ice load behavior. In the case of cylindrical structures, the experimental results predominantly indicate collisions. Conversely, for upright structures lacking slopes, the occurrence of load accumulation or climbing phenomena was relatively infrequent. This can be attributed to the dispersion of broken ice around cylindrical structures under the influence of flow velocity, inhibiting the manifestation of accumulation. Notably, the ice load exhibited a distinct triangular loading–unloading pattern.
The distribution of load curves at varying ice velocities revealed that higher velocities led to increased load frequencies and a more densely distributed set of time-history curves, aligning with experimental observations. Despite the significant randomness of ice loads and the probabilistic nature of collisions between broken ice and the structure, the overall amplitude distribution of the loads remained consistent. The time history of ice loads demonstrated distinct pulse characteristics and periodicity, consistent with the research findings of numerous scholars.
Table 4 presents the statistical findings of ice loads under various experimental conditions. The mean, significant, and maximum values of the overall ice loads were positively correlated with the experimental ice velocity. This is a consequence of the increased kinetic energy of broken ice at higher velocities, resulting in a stronger impact on the structure. Specifically, in LC2, where the ice velocity was 1.5 times that of LC1, the corresponding statistical load values were 2.94 times, 2.64 times, and 2.45 times greater than those in LC1. Similarly, in LC3, with an ice velocity 2.5 times higher than LC1, the corresponding statistical load values were 8.06 times, 6.02 times, and 6.17 times greater than those in LC1. This comparison demonstrates the nonlinear growth of loads, tending to follow a square relationship with velocity. This relationship can be attributed to the power-law distribution of characteristic ice sizes, with the power exponent ranging between 1 and 2 [
30]. The fractal distribution of broken ice in space is a widely recognized characteristic.
To delve deeper into the characteristics of broken ice load amplitudes, a provisional relationship coefficient of 2 was assumed. The dimensionless coefficients of broken ice loads for different conditions were then calculated using Equation (13) to explore their inherent attributes. A comparison of the ice force coefficients revealed that while there were slight differences in the ice force coefficients for the same statistical sample under different ice velocities, the overall distribution range remained relatively constant. The ice force coefficient for average ice force was approximately 1.37, whereas the coefficient for significant ice force was approximately 2.68, and the coefficient for maximum ice force was around 4.67.
This highlights the feasibility of conveniently defining the coefficient of broken ice loads in simulations, enabling the rapid prediction of broken ice loads on multiscale ocean structures and holding significant engineering implications.
In this equation, indicates the ice load coefficient, is the total ice load, is the ice velocity, is the ice density, and represents the area of interaction.
To verify the experimental results with actual sea ice, we referred to the literature [
31] and relevant research findings to conduct experiments on cylindrical and flat ice. By calculating the ice force coefficient, it was found that, at an ice velocity of 0.125 m/s, the overall statistical ice force coefficient was 2.95. Compared to the LC3 condition in this study, which had a meaningful ice force coefficient of 2.19, the difference was 25%. Considering the fact that we used a 70% density of fragmented ice, it can be seen that the overall ice load is relatively close, which demonstrates the validity of the experimental results in this study.
4.2. The Results of Numerical Simulation
Figure 12 presents a comparative analysis of experimental and numerical simulation results under an ice velocity of 0.05 m/s. The graph reveals that the 5–7 s interval signifies the initial contact between fragmented ice and the structure, indicating a relatively good agreement between the experimental and numerical findings. However, as the contact progressed, a certain degree of deviation emerged between the numerical simulation and experimental results. This discrepancy can be attributed to the inherent randomness in the fragmented ice load, leading to variations in the distribution patterns observed in the experimental and numerical data. Nonetheless, the overall peak loads demonstrated a satisfactory level of correspondence.
Figure 13 depicts the fragmented ice region diagram for both numerical simulation and experimental phenomena. The numerical simulation results consistently aligned with the experimental observations, demonstrating the presence of a relatively open area behind the cylindrical structure with less ice formation. This alignment highlights the accuracy of the numerical simulation method.
Figure 14 and
Figure 15 illustrate the time-domain curves of fragmented ice resistance at ice velocities of 0.075 m/s and 0.125 m/s, respectively. To further validate the accuracy of the numerical simulation method,
Table 5 presents the statistical data on the average ice load, the significant ice load, and the maximum ice load for different ice velocities.
The numerical simulation results demonstrate that the temporal profiles of the overall broken ice loads exhibit significant pulsation characteristics, similar to the experimental results. Under the computational condition of 0.05 m/s, the average, significant, and maximum values of the overall loads were 0.17 kN, 0.38 kN, and 0.76 kN, respectively. These values differed from the corresponding experimental data (0.16 kN, 0.36 kN, and 0.64 kN) by 8.31%, 6.81%, and 18.13%, respectively. Although the maximum load displayed some randomness, the differences between the other stable loads and experimental data were within 10%, indicating a high level of agreement. Under the computational condition of 0.075 m/s, the average, significant, and maximum values of the overall loads were 0.41 kN, 0.83 kN, and 1.51 kN, respectively. These values differed from the corresponding experimental data (0.47 kN, 0.95 kN, and 1.57 kN) by 12.34%, 12.63%, and 3.82%, respectively. All load differences with the experimental data were within 15%, which may be attributed to the vibration on the structure’s surface as the flow velocity increased. Under the computational condition of 0.125 m/s, the average, significant, and maximum values of the overall loads were 1.06 kN, 2.39 kN, and 3.90 kN, respectively. These values differed from the corresponding experimental data (1.29 kN, 2.17 kN, and 3.95 kN) by 17.83%, 10.14%, and 1.27%, respectively. Overall, the load data demonstrated good agreement with the experimental data.
The comparison of the experimental phenomena with the aforementioned computational results reveals that the ice load prediction technique based on the discrete element method can accurately predict ice-breaking patterns and ice loads, and therefore, it has strong potential for practical applications, and the results highlight the effectiveness of the numerical simulation model.
To investigate the influence of the ice friction coefficient on the applied load, we conducted simulations under the conditions of 0.556 m/s ice velocity with friction coefficients of 0.2, 0.3, and 0.4.
Figure 16 illustrates the time-history curves of ice fragmentation for each simulation scenario. During the initial stage (0 s~1.5 s), there were minor differences in the ice load among the three cases. However, as time progressed, the movement and contact of fragmented ice varied to different degrees due to random collisions. Consequently, the peak magnitude and occurrence time of the ice load also differed. Overall, when the friction coefficient was 0.4, the applied load was slightly higher, but the difference was not significant. Therefore, further research on ice friction coefficients should consider factors such as ice fragmentation density and ice layer type.