船舶力学

Definition	Symbol	Value	Definition	Symbol	Value
Water drag force coefficient		0.20	Thickness of sea ice	h	1.0 m
Drag moment coefficient		0.14	Normal bond strength		0.1, 0.5, 1.0, 1.5 MPa
Ice density	ρ_ice	920 kg/m³	Seawater density	ρ_w	1.025 kg/m³
Young's modulus	E	1.0×10³ MPa	Number of elements	N	4000
Poisson's ratio	ν	0.3	Current drag coefficient	c_w	0.0055

Definition	Symbol	Value	Definition	Symbol	Value
Water drag force coefficient		0.20	Thickness of sea ice	h	1.0 m
Drag moment coefficient		0.14	Normal bond strength		0.1, 0.5, 1.0, 1.5 MPa
Ice density	ρ_ice	920 kg/m³	Seawater density	ρ_w	1.025 kg/m³
Young's modulus	E	1.0×10³ MPa	Number of elements	N	4000
Poisson's ratio	ν	0.3	Current drag coefficient	c_w	0.0055

基于多面体离散元方法的非规则海冰动力特性分析

PDF下载

李继 , 王嗣强 , 刘璐 , 季顺迎

船舶力学 | 流体力学 2024,28(12): 1849-1863

收起

船舶力学 | 流体力学 2024, 28(12): 1849-1863

基于多面体离散元方法的非规则海冰动力特性分析

全屏

李继, 王嗣强, 刘璐, 季顺迎

作者信息

大连理工大学　工业装备结构分析优化及CAE软件全国重点实验室，辽宁　大连　116024

王嗣强（1993-），男，大连理工大学助理研究员

季顺迎（1972-），男，博士生导师，大连理工大学教授。

通讯作者:

通讯作者

Dynamic Characteristics of Irregular Ice Floes Based on Polyhedral Discrete Element Method

Ji LI, Si-qiang WANG, Lu LIU, Shun-ying JI

Affiliations

State Key Laboratory of Structural Analysis, Optimization and CAE Software for Industrial Equipment, Dalian University of Technology, Dalian 116024, China

LI Ji(1995-), male, Ph.D. candidate

LIU Lu(1990-), male, associate professor, corresponding author, E-mail: liulu@dlut.edu.cn.

出版时间: 2024-12-20 doi: 10.3969/j.issn.1007-7294.2024.12.004

文章导航

摘要

收起

极地海冰在多种空间尺度上呈现出离散分布特征，本文基于Gilbert-Johnson-Keerthi（GJK）和Expanding Polytope Algorithm（EPA）接触算法发展一种用于极地海冰的离散元方法。该方法允许模拟浮冰的漂移和接触碰撞过程，并基于简化的刚体有限元方法表征海冰的冻结破坏。通过在平整冰场上进行单轴压缩模拟证实程序的稳定性，并研究压缩强度对浮冰场变形的影响。在此基础上，对穿越狭窄通道的冰块集合进行离散元模拟，研究浮冰尺度、初始密集度和圆形度对其流动特性的影响规律。结果表明，在幂律分布的浮冰尺度分布函数中增大幂指数、增加圆形度以及降低初始密集度都会导致浮冰场的流动速度加快。

关键词

离散元方法 / GJK-EPA算法 / 海冰动力学 / 浮冰几何特征

Abstract

收起

In polar regions, floating ice exhibits distinct characteristics across a range of spatial scales. It is well recognized that the irregular geometry of these ice formations markedly influences their dynamic behavior. This study introduces a polyhedral Discrete Element Method (DEM) tailored for polar ice, incorporating the Gilbert-Johnson-Keerthi (GJK) and Expanding Polytope Algorithm (EPA) for contact detection. This approach facilitates the simulation of the drift and collision processes of floating ice, effectively capturing its freezing and fragmentation. Subsequently, the stability and reliability of this model are validated by uniaxial compression on level ice fields, focusing specifically on the influence of compression strength on deformation resistance. Additionally, clusters of ice floes navigating through narrow channels are simulated. These studies have qualitatively assessed the effects of Floe Size Distribution (FSD), initial concentration, and circularity on their flow dynamics. The higher power-law exponent values in the FSD, increased circularity, and decreased concentration are each associated with accelerated flow in ice floe fields. The simulation results distinctly demonstrate the considerable impact of sea ice geometry on the movement of clusters, offering valuable insights into the complexities of polar ice dynamics.

Key words

discrete element method / GJK-EPA algorithm / sea ice dynamics / floe geometry

引用本文

李继, 王嗣强, 刘璐, 季顺迎. 基于多面体离散元方法的非规则海冰动力特性分析. 船舶力学, 2024 , 28 (12) : 1849 -1863 . DOI: 10.3969/j.issn.1007-7294.2024.12.004

Ji LI, Si-qiang WANG, Lu LIU, Shun-ying JI. Dynamic Characteristics of Irregular Ice Floes Based on Polyhedral Discrete Element Method[J]. Journal of Ship Mechanics, 2024 , 28 (12) : 1849 -1863 . DOI: 10.3969/j.issn.1007-7294.2024.12.004

正文

收起

0　Introduction

收起

Covering up to 10% of Earth's ocean surface, sea ice is a critical component of the Earth's climate system^[1]. The polar oceans exhibit multiple and diverse ice distribution patterns, necessitating a detailed nomenclature for various sea ice types^[2]. These characteristics are often described as discrete, irregular and rapidly changing. Traditionally, sea ice has been modeled on a geological scale as a continuum, incorporating the floe size distribution (FSD) and ice thickness distribution (ITD) into the governing dynamics and thermodynamics equations^[3-4]. This approach statistically represents floe clusters' properties, with properties within the mesh determined by solving these governing equations. In recent years, the quest for high-precision simulations has driven the development of finer grids and enhanced parameterization schemes for the physical processes of sea ice.

Hilber's first dynamic-thermodynamic coupling ice model^[5], modified from Coon's elastic-plastic model^[6], has been proven effective and is widely utilized. An explicit elastic-viscous-plastic (EVP) numerical scheme was later proposed to enhance computational efficiency and facilitate parallel computing^[7], yielding ice motion results similar to the original model^[8]. As demands for high-precision ice characterization increase, more detailed studies have emerged, including anisotropic constitutive^[9], ice ridge parameterization^[10] and air-ice drag coefficient^[11], further refining the constitutive model for sea ice. Initially, sea ice dynamics and thermal models had a spatial resolution of 125 km due to computational limitations. However, advances in computing have enabled simulations at the 5 km grid level. Yet, in the continuum model, sea ice material properties are uniformly distributed across the computational grid, making the continuity hypothesis challenging to establish at higher spatial resolutions, as randomly generated Linear Kinematic Features (LKFs) like crevasses and ridges cannot be considered uniformly distributed^[12].

The use of the same Euler grid solution idea as in atmospheric and ocean current models in sea ice modeling restricts the accuracy of simulations. Consequently, an alternative approach has been employed to describe the dynamic motion and thermal evolution of all ice floes in Lagrangian coordinates^[13]. First developed by Cundall^[14] and widely applied in fields like geotechnical mechanics^[15], chemical engineering^[16] and ice engineering^[17-18], the discrete element model (DEM) calculates the motion of each element, offering a more natural depiction of granular materials. While early discrete element descriptions focused on the material and engineering scales due to computational demands^[19], recent advancements in computer performance and parallelization technology have bolstered confidence in using DEM to describe full-scale polar sea ice dynamics^[20-22].

Initial models in sea ice research utilized the circular element discrete element developed by Shen^[17] and the polygonal discrete element developed by Hopkins^[23]. Shen deduced the collision rheology of idealized two-dimensional flows in crushing rinks, while Hopkins first applied DEM to simulate Arctic region dynamics and subsequently to explore ice ridge formation mechanisms, mesoscale floating ice collisions, and lead formations^[24]. Full-scale triangular element discrete element simulations of the Arctic region by Hopkins demonstrated the feasibility of DEM for geological-scale numerical simulations^[25]. Like the continuum model's evolution from plane to spherical coordinates, the SIKU model projects two-dimensional polygonal elements onto a spherical coordinate system to describe sea ice motion on the Earth's surface^[21]. Additionally, the variability of the Smoothed Particle Hydrodynamics (SPH) method has been considered, with plastic disk elements describing convergence and divergence^[26]. Two-dimensional disk elements based on LAMMPS describe the process of sea ice jams^[20], and any breakable polygon element is used to describe the rapid change in geometry during the collision of sea ice materials^[27].

DEM has evolved to represent various particle shapes to depict the real physical world more accurately. Discrete element particle morphology has achieved arbitrary uneven polyhedron shapes and irregular shapes with smooth surface functions^[28-30]. Furthermore, DEM can be combined with other numerical methods to develop models for large-deformation soft materials^[31], plastic materials^[32], breakable materials^[33], and more. Among these, the GJK-EPA method has emerged as one of the most widely-used polyhedron contact techniques in DEM. The dispersion of LKFs means that the frozen ice sheet is cut into broken ice floes with a logarithmic size distribution. The natural geometry of ice floes is typically polyhedral on basin scales, with the horizontal spatial scale much larger than the scale in the ice thickness direction. Compared with spherical elements, polyhedral elements achieve a balance between computational efficiency and scale.

This paper presents a simulation method for sea ice dynamics based on polygonal discrete elements, designed to simulate the drift and collision of floating ice on mesoscale and larger scales. Compared with spherical elements, polyhedral elements are more efficient in calculating ice floes, because they can represent the complex geometric shape of a single ice floe with a smaller number of elements^[18]. GJK algorithm is a fast convex geometry contact detection algorithm with less computer memory^[34], and is widely used in the field of robot collision detection. The EPA algorithm developed on this basis can calculate the overlap^[35]. Utilizing the GJK-EPA contact algorithm, the method allows for the collision detection and contact force calculation of any irregular polygon. The subsequent chapter introduces the specific implementation process of this method, including the contact detection of floes, the DEM bond-fracture model, and the drifting process. To verify feasibility, uniaxial compression and flow tests are conducted to study the drift process under varying sea ice concentrations, FSDs, and circularities. The final chapter summarizes this paper and outlines future work.

1　Discrete element method for sea ice simulation

收起

The GJK algorithm is utilized to evaluate the overlap between floating ice elements, while the EPA quantifies this overlap. Subsequently, the Hertz contact force model is utilized to calculate the contact force based on the determined overlap. The external forces acting upon these elements are classified into three types: volume forces, surface forces and forces acting at surface contact points. Additionally, a breakable bonded floe ice model is integrated, which leverages multiple polyhedrons to accurately represent a level ice sheet.

1.1　Contact detection based on the GJK-EPA algorithm

In the GJK method, the Minkowski difference is employed to generate a new geometric construct, forming the basis for calculating the overlap between two elements. For clarity and ease of understanding, this concept is illustrated through the two-dimensional scenario of polygon intersection. As shown in Fig. 1, the Minkowski difference A ⊖ B between elements A and B can be expressed as

(1)

where x and y are the vertices of polygons A and B, respectively. The Minkowski difference A ⊖ B comprises all novel points derived by subtracting the geometric positions of each vertex in A from those in B. The criterion for determining the presence of an overlap between elements A and B is whether the convex hull of these new points encompasses the origin O. It is worth noting that the geometries of A and B must be convex for the criterion to be applicable, since only convex hulls can fulfill the condition expressed as

(2)

The surface overlap between the two geometric bodies is considered, specifically, scenarios where they are merely in contact expressed as O ∈ ∂(A ⊖ B). The criterion for calculating the contact force depends on the degree of overlap. Thus, if two elements are merely touching, the overlap is considered to be zero. As a result, no contact force will be generated between the two geometric bodies, making their state indistinguishable from a non-contact scenario.

The GJK algorithm determines the point with the maximum projection distance in a specified direction among all surface points of a convex body. This process is facilitated by constructing the support function, denoted as SF (A,v), where A represents a given convex body and v signifies any chosen direction on the convex body. At the commencement of the iteration process, the point is arbitrarily selected, with the initial query vector set as v₀ = (1, 1). The support function satisfies the following conditions:

(3a)

(3b)

EPA determines the depth of overlap by using the final simplex obtained from the GJK algorithm and extending it into polygons. After the GJK-EPA algorithm, the final normal overlap δ_n and tangential displacement δ_s can be obtained. The normal contact force F_n and tangential force F_s between two elements are determined as

(4)

(5)

where, k_n denotes the effective normal contact stiffness, defined as

. In this context, E^* represents the equivalent elastic modulus, and R^* signifies the equivalent radius. δ_n is the normal overlap,

is the vector of the normal relative velocity, and C_n denotes the normal damping coefficient. The constant κ is assigned a fixed value of 1.5. In Eq.(5), k_s represents the effective tangential stiffness between elements, defined as k_s = r_snk_n, where r_sn = 1/(2 + 2ν), and ν is the Poisson's ratio. δ_s is the tangential displacement, and μ is the coefficient of sliding friction.

A simple test between a rigid bar and a plate is employed to verify the accuracy and stability of the contact force calculation of our DEM model. As shown in Fig.2, a cube with an initial rotation angle α is placed above the plate and then moves downward at a constant velocity of v⁰. After the collision process, the cube's velocity changes to v' and the angular velocity around the Y axis is ω'. Friction and gravity are ignored during the simulation. The analytical forms of velocity v' and angular velocity ω' can be expressed as follows^[36]:

(6)

(7)

where, m is the mass of element, e is the coefficient of restitution, and I_yy is the moment of inertia around the Y axis. In this simulation, the side length of the cube l is 2 m, the value of e is 0.5, and the density of the cube is 920 kg/m³ of the sea ice. Then the mass m and moment of inertia I_yy can be calculated.

The simulation is carried out with angle α of varying initial values, and the dimensionless analytical solution is compared with the simulation results as shown in Fig.3. Although there is a small error between the simulation and the analytical solution, the results can validate the soundness of the contact force model and the contact detection process.

Intersecting LKFs, such as leads, are typically formed under external forces acting upon level ice. As shown in Fig.4, the incorporation of a breakable bonded floe ice model enables the scenarios in which a single piece of floating ice fractures into multiple segments. In this model, each individual ice floe is subdivided into smaller triangular prism elements, which are then rendered in a frozen state. More details about this model can be found in the literature of dilated polyhedron DEM^[19]. In scenarios where the number of intact bonding points on the bonding surface falls below two, the surface is considered to have undergone complete failure attributable to a loss of stability.

1.2　Drift force model

The dynamic governing equation for each unit area of ice is expressed as follows:

(8)

where, m denotes the mass of ice per unit area in the horizontal direction, u represents the velocity of the floe, τ_a and τ_w are the stresses induced by air and water velocities respectively,

accounts for the stress due to Coriolis acceleration, σ signifies the internal stress of the ice, mg∇H_o refers to the stress caused by differences in ocean surface height, and τ_b indicates the seabed stress. Considering the basic discrete element as an ice floe, the governing equation for the element can be simplified as

(9)

where, F_c denotes the contact force of the element, F_a and F_w represent the integral forms of ocean currents and atmospheric stress respectively, and F_k is the Coriolis force. In DEM calculations, τ_b and mg∇H_o are disregarded. As shown in Fig.5, the ocean force acting on an element is calculated as follows:

(10)

(11)

(12)

(13)

where, τ_w denotes the stress acting on the bottom surface of the ice; A is the bottom area of the ice; F_d represents the lateral force on the sea ice element; M_d is the moment exerted on the floe side; ρ_w is the density of seawater; A_sub is the lateral area of the ice element submerged below the water surface; r_e is the equivalent radius of the element, calculated as

; ω signifies the rotational speed of the surface of the underwater part of the sea ice; r_c is the RO vector as shown in Fig.5, with point R located at the geometric center of the ABCD plane; u and u_c are the ice velocity and ocean current velocity, respectively; and V and V_w are the velocity vectors of some floe surface. The summation symbol indicates the summation over all surfaces of the underwater part of the sea ice.

and

are the drag force coefficient and drag moment coefficient, respectively.

In comparison to hydrodynamic calculations for sea ice with the geometry of the element being considered, the atmospheric forcing is obtained as

(14)

where, c_a represents the drag force coefficient, ρ_a denotes the air density, and u_a is the wind velocity.

2　DEM simulations of floe dynamics

收起

Due to the significant disparity between the horizontal dimensions and thickness at larger scales, the ensuing numerical computations will exclusively focus on rotational dynamics in the horizontal plane of the floating ice. Uniaxial compression tests were conducted on sea ice samples of varying densities to evaluate the stability and modeling efficacy of our sea ice model. This was followed by numerical simulations of narrow channel flows, which considered variations in concentration, size distribution, and the statistical characteristics of sea ice morphology. The flow velocity of the ice under various conditions was comprehensively analyzed. Throughout the study, the sizes of sea ice elements adhere to a power-law distribution. Each element is derived through Voronoi tessellation, with their initial placement and orientations meticulously controlled to attain a predefined concentration of floating ice elements within a designated square region.

2.1　Uniaxial compression

The uniaxial compression test is widely used to evaluate the compressive behavior of sea ice. In this study, the floes are consolidated to examine the deformation of the floe field under ideal boundary conditions, aiming to elucidate the distribution of the internal stress field. As shown in Fig.6(a), the idealized floe is positioned on the sea surface and bonded together. The geometric configuration of the sea ice consists of polygons generated through Voronoi tessellation, with an average area of 100 m². The sea ice remains in a frozen state, measuring 2 km in both length and width, and is subjected to uniaxial compression at a velocity of 1 m/s. External driving forces, such as wind and currents, are not considered in this setup. However, during motion, the sea ice experiences drag forces from ocean currents acting upon its surface. Further details are provided in Tab.1.

As the floe field undergoes compression and lacks lateral constraints, it fractures and expands bidirectionally, as shown in Fig.6(b). This process results in the formation of numerous channels between the emergent cracks. To accurately capture the stress distribution within the floe field, the stress for each freezing face of the individual elements is calculated with the tension strength set to 1.5 MPa. Fig.7(a) depicts the profile of tangential stress distribution across the floe field at a strain of 0.55%, confirming that compressive stress does not induce fracture in our model. Under uniaxial compression, stress concentration is particularly notable along the direction of compression and near the center of the ice field. The orientation of fractures serves as a crucial assessment tool in analyzing sea ice deformation. Under varying strain rate conditions, the probability distribution function of fracture angles exhibits an unimodal pattern(Fig. 7(b)), which is consistent with the findings of Feltham^[37] regarding uniaxial compression of sea ice driven by wind in scenarios without refreezing.

To investigate the resistance to deformation under identical compression conditions across different sea ice freezing strengths, a statistical analysis of the driving force values is conducted, as illustrated in Fig.8. Compressing sea ice along a single axis leads to a continual increase in compressive stress as the strain of the sea ice evolves. Higher freezing strength indicates an enhanced resistance of the floe field to deformation, necessitating a more substantial external driving force to achieve an equivalent strain within the same duration. Notably, the curve displays marked fluctuations within the initial 0.01% strain, a phenomenon attributed to the time required for stress propagation within the ice. Upon reaching a strain of 0.25%, sea ice detachment occurs due to the absence of compression in opposing directions on both sides, leading to considerable damage within the frozen ice layer atop the smooth ice field. Consequently, there is a rapid decrease in overall compressive stress. While subsequent stress peak reductions continue to occur, the magnitude of each successive decrease diminishes. This pattern is ascribed to the lack of lateral compression in the floating ice field, which gradually precipitates the disintegration of sea ice on both sides.

2.2　Sea ice flow through narrow channel

Within high-resolution Earth system models, irregular and narrow geographical channels are common. For instance, sea ice flows from the central Arctic ice sheet through the Fram Strait, serving as a primary source of Arctic sea ice and freshwater output. Transport processes in narrow straits often interact with complex geographical boundaries, making them a valuable tool for evaluating the simulation capability of sea ice discrete element models^{[20, 27, 38]}. Sea ice drift in the Nares Strait exhibits significant seasonal and annual variations, and its complex and unique conditions lead to the formation of ice jams in the south and north of the Nares Strait, creating an ice arch. The formation of ice arches is influenced by the internal strength of sea ice and external forcing from the atmosphere and ocean. It is widely believed that ice arches prevent sea ice from flowing outward from the Arctic, thereby impeding the reduction of the Arctic sea ice extent. Employing polyhedral ice floe elements, rather than morphological elements like spheres, offers advantages due to their superior ability to capture geometric variations. This approach enables a more accurate representation of floating ice blocks, resembling those observable in satellite imagery.

In this section, a square region of floating ice is constructed, measuring 1 km by 1 km, with a predetermined concentration and floe size distribution function, as shown in Fig.9. Subject to the forces of wind and ocean currents, the floating ice proceeded at a constant velocity through a narrow channel, with the right-side opening spanned 250 m.

In this study, the characteristics of the floating ice within the designated region are regulated through three key parameters: the concentration (A), floe circularity (C) and the exponent (α) of the power-law Floe Size Distribution (FSD). Several studies have substantiated that the exponent α is typically less than 2. The values of current velocity u_c is 4 m/s; wind velocity u_a is 20 m/s; air density is 1.2 kg/m³ and air drag coefficient c_a is 0.0012. The same parameters and values as in Section 2.1 are not repeated here.

The initial sea ice concentration is defined as the ratio of the surface area occupied by sea ice to the total surface area of the square region. During the concentration analysis, the average number of edges of ice floes between 5 and 10 is maintained, and the exponent of the power-law size distribution function is set to 2 to mitigate their influence. Following the generation of initial sea ice using the Voronoi method, random rotations and placements within the region were applied to prevent initial contact between individual ice floes. This methodology will be consistently applied for the initial arrangement of sea ice in all subsequent simulations. As shown in Fig.10, the square sea ice region is categorized into initial densities of 50%, 55% and 60%, corresponding to element counts of 1243, 1394 and 1516, respectively.

Due to variations in concentration that alter the total mass of sea ice within the entire region, quantitatively comparing the resistance of fragmented ice to boundary conditions poses a challenge. Our findings also suggest that a higher initial ice concentration within the region leads to increased ice accumulation forces at the boundaries. Consequently, the velocity field of sea ice in a narrow channel is used to characterize floe motion. Fig.11 illustrates the evolution of the sea ice velocity field after 170 seconds of drifting and upon entering the narrow channel at different initial concentrations. It is apparent that sea ice with a higher density is more susceptible to congestion, leading to decreased drift speed and affecting the subsequent replenishment velocity within the channel. The floating ice's velocity significantly diminishes upon encountering the boundary. At lower densities, the speed reduction near the boundary exerts minimal effect on the ice velocity in the central channel. However, at higher densities, a pronounced reduction in channel speed occurs, influencing the velocity of sea ice as it exits the channel.

The size distribution of ice floes, known as Floe Size Distributions (FSD), is significantly influenced by dynamic processes within sea ice. These processes, which involve leads, ridges, and the size of sea ice, play pivotal roles in the dynamic-thermodynamic evolution of sea ice. As such, FSD at the sub-grid scale has garnered considerable scholarly attention. Numerous researchers, through the analysis of satellite or aerial images of sea ice, have found that FSD can be characterized by a power-law formula in a non-cumulative distribution format represented as f(x) as follows:

(15)

where, x denotes the size of the ice floes and α represents the power exponent. Practically, due to the existence of an upper limit for ice floe sizes within a particular scenario, the non-cumulative distribution coefficient c can be expressed as

(16)

where, a and b signify the lower and upper limits of the ice floe size range, respectively. In this study, a is set at 5 m and b is set at 100 m. In this section, the power-law Floe Size Distribution (FSD) functions with exponent values of α = 1.2, α = 1.6, and α = 2.0 are investigated. The concentration has been consistently set to 50% across all cases, resulting in ice floe field element counts of 512, 744 and 1243, respectively. The construction of the initial ice floe field is shown in Fig.12.

Under conditions of a lower power exponent, the ice floe field tends to produce larger-sized ice floes, resulting in a correspondingly reduced total number of floes when concentration is maintained constant. Owing to the maximum size limit imposed on ice floes, a higher power exponent, when larger ice floe sizes are held constant, correlates with a more fragmented ice floe field. After 300 seconds, observations revealed that large ice floes, generated under a lower exponent condition (α = 1.2), led to ice jams in narrow channels. Conversely, in the other two scenarios (α = 1.6 and α = 2.0), the ice floes navigated through the narrow channels, as shown in Fig.13.

The simulation results further demonstrate a correlation between the velocity of ice floes through narrow passages and the Floe Size Distribution function. Notably, larger exponents correspond to faster ice floe velocities. At the 200-second mark, the ice velocities for the three scenarios, in conditions where ice jams have not occurred, are shown in Fig.14.

Sphericity is an essential metric for evaluating the flow characteristics of particle motion. Considering the predominantly plane-like, two-dimensional motion characteristics of drifting sea ice, the circularity of ice floe elements can be expressed as

(17)

where, D_p denotes the perimeter of a given floe element, and D_s is defined as the perimeter of an area-equivalent circle. The ice floe fields are constructed by regulating the generation of corners in three scenarios: 5-10 corners, 10-15 corners, and 15-20 corners. Subsequently, the average circularity is calculated using Eq.(17). The sea ice concentration is uniformly maintained at 60%, and the power-law distribution exponent is consistently set at 2.0. As shown in Fig.15, the initially generated elements exhibit three distinct numbers of vertices, with corresponding average circularities of 0.83, 0.90 and 0.92, respectively.

At varying levels of circularity, significant variations in flow velocities are observed among the ice floes. Clusters of ice floes with higher circularity demonstrate markedly faster flow velocities compared to those with more irregular shapes, as shown in Fig.16.

3　Concluding remarks

收起

The present study introduces a discrete element method based on the GJK-EPA contact detection approach for simulating the drift and collisions of ice at the floe size scale. The inclusion of a bond-failure model between elements makes it more conducive for representing the phenomenon of freezing and breaking in level ice regions. The stability and effectiveness of the model were validated through uniaxial compression tests on level ice surfaces under ideal conditions. A simple sensitivity analysis of the bonding strength parameters was conducted. Increased bond strength results in ice layers that are more resistant to breaking.

While there has been considerable recognition and further development in dynamic simulation studies related to the scale of floe ice, research on the impact of floating ice geometry on its flow characteristics remains limited. Numerous discrete element methods based on geometric shapes such as circles and spheres have been employed for simulating ice floe scales. This study qualitatively explores the influence of ice floe roundness, concentration, and Floe Size Distribution on the drifting process of ice through a narrow channel. The results indicate that higher roundness, increased density, and lower power-law distribution exponent values result in difficulty for ice to flow even in the jam phenomenon.

In future work, exploring the discrete element method applied to arbitrary concave and convex polygonal shapes is among the proposed research directions. Additionally, it is essential to incorporate considerations for changes in floe thickness and size resulting from thermodynamic to enable simulations spanning seasonal time scales. Although the study does not provide a numerical investigation specific to geographical boundary conditions, the results are sufficient to underscore the importance of considering the geometric features of floating ice in the study of ice drift at the geophysical scale. Further research is warranted to investigate the influence of additional geometric features of floating ice on its flow characteristics. For instance, although the scale distribution of sea ice conforms to a power-law distribution in different seasons, the maximum size of ice floes undergoes significant variations during the processes of freezing, fragmentation and complete disappearance. This variability is particularly noteworthy in satellite observations and poses a research question, especially concerning the volume of ice passing through narrow channels. Furthermore, the LKFs of floating ice often lead to the formation of elongated ice floes with a high aspect ratio, which should also be considered in future studies.

基金

收起

国家自然科学基金资助项目(52192693; 52192690; 42176241; 12102083)
工业和信息化部专项(2021-342)

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

[1]

Leppäranta

. The drift of sea ice[M]. Springer Science & Business Media, 2011.

[2]

WMO sea-ice nomenclature: Terminology, codes and illustrated glossary[R]. Secretariat of the World Meteorological Organization, 1970.

[3]

Rothrock

D A

, Thorndike

A S

. Measuring the sea ice floe size distribution[J]. Journal of Geophysical Research, 1984, 89(C4): 6477.

[4]

Thorndike

A S

, Rothrock

D A

, Maykut

G A

, et al. The thickness distribution of sea ice[J]. Journal of Geophysical Research, 1975, 80: 4501-4513.

[5]

Hibler

W D

. A dynamic thermodynamic sea ice model[J]. Journal of Physical Oceanography, 1979, 9(4): 815-846.

[6]

Coon

M D

, Maykut

G A

, Pritchard

R S

, et al. Modeling the pack ice as an elastic-plastic material[J]. AIDJEX Bulletin, 1974, 24: 1-105.

[7]

Hunke

E C

, Dukowicz

J K

. An elastic-viscous-plastic model for sea ice dynamics[C]//International Symposium on Representation of the Cryosphere in Climate & Hydrological Models, British Columbia (Canada), Los Alamos National Laboratory, 1996.

[8]

Hunke

E C

, Zhang

. A comparison of sea ice dynamics models at high resolution[J]. Monthly Weather Review, 1999, 127(3): 396-408.

[9]

Wilchinsky

A V

. Rheology of discrete failure regimes of anisotropic sea ice[J]. Journal of Physical Oceanography, 2012, 42.7 (2012): 3042.

[10]

Hopkins

M A

. Four stages of pressure ridging[J]. Journal of Geophysical Research, 1998, 103(C10): 21883-21891.

[11]

Fujisaki

, Yamaguchi

, Toyota

, et al. Measurements of air-ice drag coefficient over the ice-covered Sea of Okhotsk[J]. Journal of Oceanography, 2009, 65(4): 487-498.

[12]

Mohammadi

A M

, Goessling

H F

, Losch

, et al. Predictability of Arctic sea ice on weather time scales[J]. Scientific Reports, 2018, 8(1): 6514.

[13]

Hopkins

M A

. A discrete element Lagrangian sea ice model[J]. Engineering Computations, 2004, 21(2/3/4): 409-421.

[14]

Cundall

P A

, Strack

O D L

. Discrete numerical model for granular assemblies[J]. Geotechnique, 1979, 29: 47-65.

[15]

Liu

Guangyu

, Xu

Wenjie

, Nicolin

, et al. Simulation of rock fracture process based on GPU-accelerated discrete element method[J]. Powder Technology, 2021, 377: 640-656.

[16]

Wang

, Fan

, Ji

. Interaction between super-quadric particles and triangular elements and its application to hopper discharge[J]. Powder Technology, 2018, 339: 534-549.

[17]

Shen

H H

, Hibler III

W D

, Leppäranta

. The role of floe collisions in sea ice rheology[J]. Journal of Geophysical Research:Oceans, 1987, 92(C7): 7085-7096.

[18]

Liu

, Ji

.Comparison of sphere-based and dilated-polyhedron-based discrete element methods for the analysis of ship-ice interactions in level ice[J]. Ocean Engineering, 2022, 244: 110364.

[19]

Liu

, Ji

. Bond and fracture model in dilated polyhedral DEM and its application to simulate breakage of brittle materials[J]. Granular Matter, 2019, 21(3): 41.1-41.16.

[20]

Damsgaard

, Adcroft

, Sergienko

. Application of discrete element methods to approximate sea ice dynamics[J]. Journal of Advances in Modeling Earth Systems, 2018, 10(9): 2228-2244.

[21]

Kulchitsky

, Hutchings

, Johnson

J B

. Siku sea ice discrete element method model[R]. Coastal Marine Institute, College of Fisheries and Ocean Sciences, University of Alaska Fairbanks, 2017.

[22]

Åström

J A

, Vallot

, Schäfer

, et al. Termini of calving glaciers as self-organized critical systems[J]. Nature Geoscience, 2014, 7(12): 874-878.

[23]

Hopkins

M A

, Hibler III

W D

. On the shear strength of geophysical scale ice rubble[J]. Cold Regions Science and Technology, 1991, 19(2): 201-212.

[24]

Hopkins

M A

. The effects of individual ridging events on the ice thickness distribution in the Arctic ice pack[J]. Cold Regions Science and Technology, 1996, 24(1): 75-82.

[25]

Hopkins

M A

. On the mesoscale interaction of lead ice and floes[J]. Journal of Geophysical Research: Oceans, 1996, 101(C8): 18315-18326.

[26]

, Li

, Liu

, et al. A modified discrete element model for sea ice dynamics[J]. Acta Oceanologica Sinica, 2014, 33(1):56-63.

[27]

Manucharyan

G E

, Montemuro

B P

. SubZero: A sea ice model with an explicit representation of the floe life cycle[J]. Journal of Advances in Modeling Earth Systems, 2022, 14(12): e2022MS003247.

[28]

Delaney

G W

, Cleary

P W

. The packing properties of superellipsoids[J]. Europhysics Letters, 2010, 89(3): 34002.

[29]

Liu

, Ji

. A new contact detection method for arbitrary dilated polyhedra with potential function in discrete element method[J]. International Journal for Numerical Methods in Engineering, 2020, 121(24): 5742-5765.

[30]

Wang

, Xu

, Ji

. A novel Minkowski sum contact algorithm for arbitrarily shaped particles constructed by multiple dilated DEM models[J]. International Journal of Solids and Structures, 2023, 280: 112409.

[31]

André

, Iordanoff

, Charles

, et al. Discrete element method to simulate continuous material by using the cohesive beam model[J]. Computer Methods in Applied Mechanics and Engineering, 2012, 213: 113-125.

[32]

Nezamabadi

, Ghadiri

, Delenne

J Y

, et al. Modelling the compaction of plastic particle packings[J]. Computational Particle Mechanics, 2021: 1-8.

[33]

Harmon

J M

, Arthur

, Andrade

J E

. Level set splitting in DEM for modeling breakage mechanics[J]. Computer Methods in Applied Mechanics and Engineering, 2020, 365: 112961.

[34]

Feng

Y T

, Tan

. On Minkowski difference-based contact detection in discrete/discontinuous modelling of convex polygons/polyhedra[J]. Engineering Computations, 2019, 37(1): 54-72.

[35]

Van Den Bergen

. Proximity queries and penetration depth computation on 3d game objects[C]//Game Developers Conference, 2001, 170: 209.

[36]

Park

. Modeling the dynamics of fabric in a rotating horizontal drum[D]. Purdue University, 2003.

[37]

Wilchinsky

A V

, Feltham

D L

, Hopkins

M A

. Modelling the reorientation of sea-ice faults as the wind changes direction[J]. Annals of Glaciology, 2011, 52(57): 83-90.

[38]

Rabatel

, Labbé

, Weiss

. Dynamics of an assembly of rigid ice floes[J]. Journal of Geophysical Research: Oceans, 2015, 120(9): 5887-5909.

2024年第28卷第12期

PDF下载

113

引用本文

BibTeX

文章信息

doi: 10.3969/j.issn.1007-7294.2024.12.004

接收时间：2024-06-13
首发时间：2026-03-26
出版时间：2024-12-20

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2024-06-13

基金

National Natural Science Foundation of China(52192693; 52192690; 42176241; 12102083)

国家自然科学基金资助项目(52192693; 52192690; 42176241; 12102083)

Special Project of Ministry of Industry and Information Technology of China(2021-342)

工业和信息化部专项(2021-342)

作者信息

大连理工大学　工业装备结构分析优化及CAE软件全国重点实验室，辽宁　大连　116024

通讯作者:

通讯作者

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/cblx/CN/10.3969/j.issn.1007-7294.2024.12.004

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT