Discrete Element Analysis of High-Pressure Zones of Sea Ice on Vertical Structures

In cold regions, ice pressure poses a serious threat to the safe operation of ship hulls and fixed offshore platforms. In this study, a discrete element method (DEM) with bonded particles was adapted to simulate the generation and distribution of local ice pressures during the interaction between level ice and vertical structures. The strength and failure mode of simulated sea ice under uniaxial compression were consistent with the experimental results, which verifies the accuracy of the discrete element parameters. The crushing process of sea ice acting on the vertical structure simulated by the DEM was compared with the field test. The distribution of ice pressure on the contact surface was calculated, and it was found that the local ice pressure was much greater than the global ice pressure. The high-pressure zones in sea ice are mainly caused by its simultaneous destruction, and these zones are primarily distributed near the midline of the contact area of sea ice and the structure. The contact area and loading rate are the two main factors affecting the high-pressure zones. The maximum local and global ice pressures decrease with an increase in the contact area. The influence of the loading rate on the local ice pressure is caused by the change in the sea ice failure mode. When the loading rate is low, ductile failure of sea ice occurs, and the ice pressure increases with the increase in the loading rate. When the loading rate is high, brittle failure of sea ice occurs, and the ice pressure decreases with an increase in the loading rate. This DEM study of sea ice can reasonably predict the distribution of high-pressure zones on marine structures and provide a reference for the anti-ice performance design of marine structures.


Introduction
Ice pressure exists in the process of collision between sea ice and marine structure in ice-covered waters. Sea ice pressure poses great risk for navigation and leads to ship besetting and damages [1,2]. By studying the ice pressure, people can strengthen the ice resistance of the structure and predict the ice load of the structure. Ice pressure mainly includes local ice pressure and global ice pressure [3]. Global ice pressure is the pressure on the whole contact surface of the structure under ice load. The local ice pressure is used to characterize the ice pressure that acts on a small local area of the structure during the destruction of sea ice. The area of concentrated stress caused by the non-simultaneous failure of sea ice and the structure is called the high-pressure zone [4][5][6]. Many experimental studies have proved that the maximum local ice pressure on a marine structure caused by sea ice can exceed 60 MPa, which is much higher than the global ice pressure [7,8]. Therefore, the existence of local high-pressure zones may cause local deformation and damage to marine structures and even affect the overall stability of structures.
Field tests and observational methods are often adopted to study the interaction process between sea ice and a vertical structure and to analyze the causes and distribution of high-pressure zones on the structure [9][10][11]. Through the observation of the sea ice 2 of 14 breaking process after the action with a vertical structure, it can be found that the sea ice breaking mainly occurs at the upper and lower edges, while the sea ice in the middle layer is not easy to break. The broken sea ice is constantly squeezed out of the contact surface, but the sea ice in the middle layer always keeps in contact with the structure, forming a continuous ice pressure on the structure. This steady ice pressure forms the high-pressure zone of the structure. With increases in the load, the ice in the middle layer will break, the internal microcracks in the ice will gradually expand, and the ice on the upper and lower surfaces will peel off again. These results show that the high-pressure zone is continuous and does not change with the ice load, and the distribution of the high-pressure zone is closely related to the sea ice failure mode. Therefore, the high-pressure zone often occurs near the center line of the contact surface and has a linear distribution, called a "line-like load" [12,13]. However, the distribution of the high-pressure zone is also affected by other conditions, including sea ice thickness, loading rate, and structure shape.
Through an experimental study, it was found that the local ice pressure on the structure decreases with an increase in the contact area between the sea ice and the structure [14][15][16]. This relationship is described in the form of an index and has been widely recognized. However, it is also affected by other loading conditions, such as ratio of structural width to ice thickness, sea ice strength, and loading rate. When the ice is thick, owing to the change in sea ice failure mode, the distribution of the high-pressure zone is no longer concentrated in the central line, and multiple high-pressure zones may appear in the vertical direction of the contact surface [17,18]. Sodhi [19] also pointed out that the loading rate affects the failure mode of sea ice and the distribution of the high-pressure zone. When the loading rate is low, the ice pressure increases with the increase in the loading rate. Meanwhile, the ice pressure is evenly distributed on the contact surface. When the loading rate is high, the ice pressure decreases with the increase in the loading rate, which is distributed along the middle line of the contact surface. The test method can analyze the influence of various factors on the distribution of the high-pressure zone from the data results, but it is difficult to explain the change mechanism of the high-pressure zone from the failure mode of sea ice, especially the internal breakup of sea ice.
The numerical method can not only describe the interaction process between sea ice and marine structures but also deeply analyze the transition mechanism of the sea ice failure mode and the influencing factors of ice load. Based on the wide application of the finite element method in engineering, this method has also been applied to the study of the fracture characteristics of sea ice [20,21]. However, this method has the disadvantages of a huge computational resource consumption and low computational efficiency. Since the 1980s, discrete element method (DEM) has been used to analyze the dynamic evolution of sea ice and the ice load on marine structures. The DEM is not only suitable for the simulation of ice collisions at small scales [22,23] but also for the numerical calculation of ice ridge formation, ice crevasse generation, and sea ice evolution at mesoscales [24]. In recent years, the DEM based on GPU (graphics processing unit) parallel computing has allowed great progress to be made in the study of the interaction of sea ice with offshore structures and ship hulls [25,26]. The results show that the DEM can simulate the different failure modes of sea ice and calculate the dynamic ice load during the interaction between sea ice and a structure. This method can describe the shape, size, and accumulation of sea ice after breaking, and it has excellent computational performance and high computational efficiency in large-scale sea ice. In order to simulate the crushing failure mode of sea ice, the accuracy and quantity of particle elements are required. The large number of particle elements on the interface between sea ice and a structure is helpful for calculating the ice pressure distribution on the structure. Therefore, the DEM also has the potential to simulate the ice pressure distribution of structures.
In this study, the DEM with bond and break functions was adopted to simulate the interaction between sea ice and vertical structures. The simulated results of the local pressure distribution and ice load were compared with field experiments [27] to validate the rationality of the DEM. The causes of the distribution characteristics in high-pressure zones were studied, and the influences of contact areas and load speeds on local ice pressures were analyzed.

Parallel Bond Model of Sea Ice
A regular arrangement of spherical particles is adopted to simulate sea ice using the DEM. The spherical particles are glued together with the parallel bond model, as shown in Figure 1 [28]. The forces and moments are transferred between two adjacent particles with an elastic bond disk, which can represent the elastic spring with constant stiffnesses in normal and shear directions. Hence, the bond force can be calculated on the bond disk and centered at the contact point. rationality of the DEM. The causes of the distribution characteristics in high-pressure zones were studied, and the influences of contact areas and load speeds on local ice pres sures were analyzed.

Parallel Bond Model of Sea Ice
A regular arrangement of spherical particles is adopted to simulate sea ice using the DEM. The spherical particles are glued together with the parallel bond model, as shown in Figure 1 [28]. The forces and moments are transferred between two adjacent particles with an elastic bond disk, which can represent the elastic spring with constant stiffnesses in normal and shear directions. Hence, the bond force can be calculated on the bond disk and centered at the contact point. , (3 where and are the normal and shear bonding strengths between particles, re spectively.

Verification of Sea Ice Discrete Element Method
To describe more accurately the mechanical properties of sea ice materials using the DEM, in this section, a uniaxial compression test of sea ice based on the mechanical tes data of the Bohai Sea ice is simulated [29]. The sea ice specimen in the uniaxial compres sion experiment is constructed using spherical elements of the same size and regular ar rangement. The spherical elements are arranged and bonded in a hexagonal close-packed (HCP) structure to form a cuboid-shaped sea ice specimen. The cohesive failure of the particle element is used to simulate the formation of sea ice cracks, and the crushing fail ure process of sea ice is then calculated. The calculated compression strength and crushing failure process of sea ice are compared with the test data to verify the accuracy of the discrete element parameters of sea ice. The force and moment in the normal and shear directions associated with the parallel bond are denoted by F n , F s , M n , and M s . The tensile and shear stresses acting on the bond disk are calculated based on the beam theory as where the variables A, I, and J are given by If the tensile stress exceeds the normal bond strength, or the shear stress exceeds the shear bond strength, then the parallel bond breaks.
where σ b max and τ b max are the normal and shear bonding strengths between particles, respectively.

Verification of Sea Ice Discrete Element Method
To describe more accurately the mechanical properties of sea ice materials using the DEM, in this section, a uniaxial compression test of sea ice based on the mechanical test data of the Bohai Sea ice is simulated [29]. The sea ice specimen in the uniaxial compression experiment is constructed using spherical elements of the same size and regular arrangement. The spherical elements are arranged and bonded in a hexagonal close-packed (HCP) structure to form a cuboid-shaped sea ice specimen. The cohesive failure of the particle element is used to simulate the formation of sea ice cracks, and the crushing failure process of sea ice is then calculated. The calculated compression strength and crushing failure process of sea ice are compared with the test data to verify the accuracy of the discrete element parameters of sea ice.
According to the sea ice specimen size in the tests, the uniaxial compressive specimen size (a × b × h) in this DEM simulation is 0.15 m × 0.15 m × 0.375 m, as shown in Figure 2. The bottom plate is fixed, and the top indenter loads the specimen at a speed of 0.01 m/s in the uniaxial compressive test. The main computational discrete element parameters of the sea ice specimen are listed in Table 1. In the uniaxial compression test, the pressure sensor on the top indenter measures the force on the loading surface. If the maximum load measured is P max , the uniaxial compressive strength of sea ice can be expressed as: According to the sea ice specimen size in the tests, the uniaxial compressive specimen size ( ℎ) in this DEM simulation is 0.15 m × 0.15 m × 0.375 m, as shown in Figure 2. The bottom plate is fixed, and the top indenter loads the specimen at a speed of 0.01 m/s in the uniaxial compressive test. The main computational discrete element parameters of the sea ice specimen are listed in Table 1. In the uniaxial compression test, the pressure sensor on the top indenter measures the force on the loading surface. If the maximum load measured is , the uniaxial compressive strength of sea ice can be expressed as:  In the DEM simulation of sea ice, the key to the accuracy of numerical simulation is to select the DEM parameters according to the known mechanical properties of sea ice. The results of discrete element simulation show that the compressive strength of sea ice is closely related to particle element size, bonding strength, and friction coefficient [30]. When the particle bonding strength is 0.7 MPa, the compressive strength of sea ice increases with the decrease of particle size, and the relationship between them is nonlinear, as shown in Equation (5). Meanwhile, the compressive strength of sea ice increases linearly with the increase of particle bonding strength. Therefore, the relationship between particle bonding strength and compressive strength of sea ice can be expressed as Equation (6), which considers the particle size effect.   In the DEM simulation of sea ice, the key to the accuracy of numerical simulation is to select the DEM parameters according to the known mechanical properties of sea ice. The results of discrete element simulation show that the compressive strength of sea ice is closely related to particle element size, bonding strength, and friction coefficient [30]. When the particle bonding strength is 0.7 MPa, the compressive strength of sea ice increases with the decrease of particle size, and the relationship between them is nonlinear, as shown in Equation (5). Meanwhile, the compressive strength of sea ice increases linearly with the increase of particle bonding strength. Therefore, the relationship between particle bonding strength and compressive strength of sea ice can be expressed as Equation (6), which considers the particle size effect. (6) When simulating the interaction between level ice and a marine structure, L can be regarded as sea ice thickness according to the direction of action between the sea ice and the structure, and L/D can be regarded as the number of particle layers that constitute the level ice. When sea ice thickness, particle size, and compressive strength of sea ice are determined in the simulation, the bonding strength between particles also can be determined. This method can effectively reduce the influence of particle size on the strength of sea ice in DEM simulations and provides a basis for the selection of discrete element parameters.
During the continuous compression of sea ice, the sea ice specimen is broken into many fragments, and the cracks spread from top to bottom throughout the specimen. The breaking process of sea ice shows a brittle failure at a fast loading rate, as shown in Figure 3. After a dominant crack is formed in the specimen, the sea ice splits rapidly along the crack. Figure 4 shows the results of the DEM simulation for the failure process of sea ice under uniaxial compression, and the color of the particles indicates its velocity. In the sea ice breaking process, both the specimens in the test and DEM simulation show a splitting failure mode. Although the number and location of the cracks in the specimens in the simulation results are different from those in the test results, the form and direction of the cracks are relatively close. = 0.7 / 3.77 3.41 / . (6) When simulating the interaction between level ice and a marine structure, L can be regarded as sea ice thickness according to the direction of action between the sea ice and the structure, and / can be regarded as the number of particle layers that constitute the level ice. When sea ice thickness, particle size, and compressive strength of sea ice are determined in the simulation, the bonding strength between particles also can be determined. This method can effectively reduce the influence of particle size on the strength of sea ice in DEM simulations and provides a basis for the selection of discrete element parameters.
During the continuous compression of sea ice, the sea ice specimen is broken into many fragments, and the cracks spread from top to bottom throughout the specimen. The breaking process of sea ice shows a brittle failure at a fast loading rate, as shown in Figure  3. After a dominant crack is formed in the specimen, the sea ice splits rapidly along the crack. Figure 4 shows the results of the DEM simulation for the failure process of sea ice under uniaxial compression, and the color of the particles indicates its velocity. In the sea ice breaking process, both the specimens in the test and DEM simulation show a splitting failure mode. Although the number and location of the cracks in the specimens in the simulation results are different from those in the test results, the form and direction of the cracks are relatively close.
When simulating the interaction between level ice and a marine structure, L can be regarded as sea ice thickness according to the direction of action between the sea ice and the structure, and / can be regarded as the number of particle layers that constitute the level ice. When sea ice thickness, particle size, and compressive strength of sea ice are determined in the simulation, the bonding strength between particles also can be determined. This method can effectively reduce the influence of particle size on the strength of sea ice in DEM simulations and provides a basis for the selection of discrete element parameters.
During the continuous compression of sea ice, the sea ice specimen is broken into many fragments, and the cracks spread from top to bottom throughout the specimen. The breaking process of sea ice shows a brittle failure at a fast loading rate, as shown in Figure  3. After a dominant crack is formed in the specimen, the sea ice splits rapidly along the crack. Figure 4 shows the results of the DEM simulation for the failure process of sea ice under uniaxial compression, and the color of the particles indicates its velocity. In the sea ice breaking process, both the specimens in the test and DEM simulation show a splitting failure mode. Although the number and location of the cracks in the specimens in the simulation results are different from those in the test results, the form and direction of the cracks are relatively close.   The stress-strain curve of sea ice can be obtained according to the load variation with the time measured by the top indenter. A comparison between the test results and the DEM simulation results is shown in Figure 5. With the constant loading of the indenter, the stress and strain of the sea ice specimen increase gradually, and there is no obvious crack in the sea ice specimen. When the indenter load reaches its peak value, many cracks appear in the specimen; then, the specimen breaks, causing the indenter to separate from the specimen, and the load drops rapidly to zero. The compressive strength of sea ice is 4.2 MPa in the test and 3.8 MPa in the DEM simulation. The slope of the straight-line OA in Figure 5 is the elastic modulus of sea ice. The elastic modulus of sea ice in the DEM simulation is 1.42 GPa, which is close to the measured elastic modulus of sea ice in the Bohai Sea. Although the test and DEM simulation results are not the same in value, the variation trends in the sea ice stress-strain curve and failure mode are similar. This shows that the DEM can be used to simulate the compression failure process of sea ice, and the parameters in the DEM simulation can reasonably describe the mechanical properties of sea ice. The stress-strain curve of sea ice can be obtained according to the load variation with the time measured by the top indenter. A comparison between the test results and the DEM simulation results is shown in Figure 5. With the constant loading of the indenter, the stress and strain of the sea ice specimen increase gradually, and there is no obvious crack in the sea ice specimen. When the indenter load reaches its peak value, many cracks appear in the specimen; then, the specimen breaks, causing the indenter to separate from the specimen, and the load drops rapidly to zero. The compressive strength of sea ice is 4.2 MPa in the test and 3.8 MPa in the DEM simulation. The slope of the straight-line OA in Figure 5 is the elastic modulus of sea ice. The elastic modulus of sea ice in the DEM simulation is 1.42 GPa, which is close to the measured elastic modulus of sea ice in the Bohai Sea. Although the test and DEM simulation results are not the same in value, the variation trends in the sea ice stress-strain curve and failure mode are similar. This shows that the DEM can be used to simulate the compression failure process of sea ice, and the parameters in the DEM simulation can reasonably describe the mechanical properties of sea ice.

Comparison with a Field Test
The interaction process between the sea ice and a vertical structure was simulated by the DEM to analyze the formation mechanism of high-pressure zones in the sea ice crushing failure process. The simulated results are compared with the field test data collected by the Japan Ocean Industries Association (JOIA) in February 1999 [27]. In this experiment, an indenter with a height of 0.5 m and width of 1.5 m was used to act vertically on the stationary ice sheet at a horizontal velocity of 0.003 m/s, as shown in Figure 6(a). The pressure sensors distributed uniformly on the indenter monitored the change in ice pressure at all times during the loading process and provided reliable verification data for the numerical simulation of the ice pressure distribution when sea ice interacts with a vertical structure. In addition, the thickness of sea ice was 0.168 m, and its strength was 2.46 MPa. The relevant discrete element parameters were determined according to the sea ice and structural parameters in the test, and the established calculation model is shown in Figure  6(b). The sea ice in the simulation was composed of spherical elements with the same mass and size and regular arrangement, and the structure was composed of triangular elements. To finely simulate the crushing of sea ice, ten element layers were arranged in the direction of the ice thickness, and the three sides not in contact with the structure were fixed by linear spring boundaries. In the simulation, the indenter was regarded as a rigid

Comparison with a Field Test
The interaction process between the sea ice and a vertical structure was simulated by the DEM to analyze the formation mechanism of high-pressure zones in the sea ice crushing failure process. The simulated results are compared with the field test data collected by the Japan Ocean Industries Association (JOIA) in February 1999 [27]. In this experiment, an indenter with a height of 0.5 m and width of 1.5 m was used to act vertically on the stationary ice sheet at a horizontal velocity of 0.003 m/s, as shown in Figure 6a. The pressure sensors distributed uniformly on the indenter monitored the change in ice pressure at all times during the loading process and provided reliable verification data for the numerical simulation of the ice pressure distribution when sea ice interacts with a vertical structure. In addition, the thickness of sea ice was 0.168 m, and its strength was 2.46 MPa. The relevant discrete element parameters were determined according to the sea ice and structural parameters in the test, and the established calculation model is shown in Figure 6b. The sea ice in the simulation was composed of spherical elements with the same mass and size and regular arrangement, and the structure was composed of triangular elements. To finely simulate the crushing of sea ice, ten element layers were arranged in the direction of the ice thickness, and the three sides not in contact with the structure were fixed by linear spring boundaries. In the simulation, the indenter was regarded as a rigid body with a fixed velocity, and the structural vibration, deformation, and other motion modes were not considered. The loading indenter also acted on the sea ice at a loading rate of 0.003 m/s. The bonding strength of particles can be calculated by Equation (6). The main discrete element parameters in the simulation are shown in Table 2, and the other parameters are shown in Table 1. body with a fixed velocity, and the structural vibration, deformation, and other motion modes were not considered. The loading indenter also acted on the sea ice at a loading rate of 0.003 m/s. The bonding strength of particles can be calculated by Equation 6. The main discrete element parameters in the simulation are shown in Table 2, and the other parameters are shown in Table 1.
(a) (b) Figure 6. Interaction between sea ice and indenter: (a) Field test model [27]; (b) DEM simulation. The interaction process between the sea ice and the indenter was calculated by the DEM, as shown in Figure 7. When the interaction between the sea ice and structure began, small cracks gradually appeared in the sea ice. Then, the sea ice on the upper and lower surfaces of the ice sheet was broken and gradually squeezed out of the contact surface. The broken sea ice above the surface of the water could not be removed and continued to accumulate. With the continuous loading of the structure, the height of sea ice accumulation increased gradually.  The interaction process between the sea ice and the indenter was calculated by the DEM, as shown in Figure 7. When the interaction between the sea ice and structure began, small cracks gradually appeared in the sea ice. Then, the sea ice on the upper and lower surfaces of the ice sheet was broken and gradually squeezed out of the contact surface. The broken sea ice above the surface of the water could not be removed and continued to accumulate. With the continuous loading of the structure, the height of sea ice accumulation increased gradually.
In Figure 8, the ice load time history curve calculated by the DEM is compared with the field test data. The changing trend of the ice load can be roughly divided into three stages. In the first stage, within 10 s, the ice load reaches its maximum and then decreases rapidly. When the ice load increases, the sea ice is in the elastic deformation stage. At this time, the sea ice does not show obvious breaks; however, internal cracks have formed, and the sea ice shows ductile failure. When the ice load reaches its maximum, large cracks appear in the ice and large-scale fractures occur, which leads to the rapid reduction of the ice load. In the second stage (10-70 s), the ice load peak increases gradually. The broken sea ice at the front of the ice sheet extrudes continuously out of the contact surface under the action of the subsequent sea ice, and the subsequent sea ice contacts the structure again, which gradually increases the peak ice load. The third stage is as follows: after 70 s, the ice load is in a stable period; the crushing failure mode of the sea ice after interacting with the structure is stable, which reduces the amplitude and frequency of the ice load. However, the peak value of the ice load in the increasing and stable periods is far less than the maximum value of the ice load at the initial contact. This is because the sea ice contact surface area is the largest when it first contacts the structure, and as there are no cracks inside the sea ice, its strength is at its highest. Once the sea ice is destroyed, there are always new cracks in the front of the ice, which will reduce its strength. In Figure 8, the ice load time history curve calculated by the DEM is compared with the field test data. The changing trend of the ice load can be roughly divided into three stages. In the first stage, within 10 s, the ice load reaches its maximum and then decreases rapidly. When the ice load increases, the sea ice is in the elastic deformation stage. At this time, the sea ice does not show obvious breaks; however, internal cracks have formed, and the sea ice shows ductile failure. When the ice load reaches its maximum, large cracks appear in the ice and large-scale fractures occur, which leads to the rapid reduction of the ice load. In the second stage (10-70 s), the ice load peak increases gradually. The broken sea ice at the front of the ice sheet extrudes continuously out of the contact surface under the action of the subsequent sea ice, and the subsequent sea ice contacts the structure again, which gradually increases the peak ice load. The third stage is as follows: after 70 s, the ice load is in a stable period; the crushing failure mode of the sea ice after interacting with the structure is stable, which reduces the amplitude and frequency of the ice load. However, the peak value of the ice load in the increasing and stable periods is far less than the maximum value of the ice load at the initial contact. This is because the sea ice contact surface area is the largest when it first contacts the structure, and as there are no cracks inside the sea ice, its strength is at its highest. Once the sea ice is destroyed, there are always new cracks in the front of the ice, which will reduce its strength. From the simulation results, the peak value and variation form of the ice load time history curve are close to the field test data. Therefore, the sea ice crushing failure simulated by the DEM is similar to the real situation, which proves that the selection of the discrete element parameter is reasonable for sea ice.

Formation of High-Pressure Zones
The formation mechanism of the high-pressure zone under the action of sea ice can be explained by analyzing the sea ice failure process with the DEM. When simulating the interaction between sea ice and a vertical structure, the surface of the indenter is evenly divided into 0.02 mm × 0.02 mm grids. The contact force between the ice and structure in each grid on the contact surface is calculated at intervals. The maximum value of the local ice pressure in all grids at each time is calculated to obtain its variation with time, as  From the simulation results, the peak value and variation form of the ice load time history curve are close to the field test data. Therefore, the sea ice crushing failure simulated by the DEM is similar to the real situation, which proves that the selection of the discrete element parameter is reasonable for sea ice.

Formation of High-Pressure Zones
The formation mechanism of the high-pressure zone under the action of sea ice can be explained by analyzing the sea ice failure process with the DEM. When simulating the interaction between sea ice and a vertical structure, the surface of the indenter is evenly divided into 0.02 mm × 0.02 mm grids. The contact force between the ice and structure in each grid on the contact surface is calculated at intervals. The maximum value of the local ice pressure in all grids at each time is calculated to obtain its variation with time, as shown in Figure 9. The maximum local ice pressure is 118.9 MPa when t = 266 s, and the average ice pressure is 60.2 MPa. It can be seen from Figure 8 that when t = 266 s, the global ice pressure on the structure is only 0.6 MPa, which is far less than the local ice pressure. The local ice pressure is not related to the global ice pressure. It shows that the local ice pressure does not change with the ice load, but is closely related to the failure mode of sea ice. history curve are close to the field test data. Therefore, the sea ice crushing failure simu lated by the DEM is similar to the real situation, which proves that the selection of th discrete element parameter is reasonable for sea ice.

Formation of High-Pressure Zones
The formation mechanism of the high-pressure zone under the action of sea ice can be explained by analyzing the sea ice failure process with the DEM. When simulating the interaction between sea ice and a vertical structure, the surface of the indenter is evenly divided into 0.02 mm × 0.02 mm grids. The contact force between the ice and structure in each grid on the contact surface is calculated at intervals. The maximum value of the loca ice pressure in all grids at each time is calculated to obtain its variation with time, a shown in Figure 9. The maximum local ice pressure is 118.9 MPa when t = 266 s, and the average ice pressure is 60.2 MPa. It can be seen from Figure 8 that when t = 266 s, the globa ice pressure on the structure is only 0.6 MPa, which is far less than the local ice pressure The local ice pressure is not related to the global ice pressure. It shows that the local ic pressure does not change with the ice load, but is closely related to the failure mode of sea ice. The damage and movement of the sea ice in the contact area can be observed more clearly by vertically dividing the contact surface of the sea ice and indenter. In Figure 10, the breaking state of the sea ice at different times on the vertical plane is shown, and the color indicates the stress condition of the particle element. With the continuous breaking of sea ice, the contact area and crushing degree between the sea ice and structure changes with time. The sea ice is in the elastic deformation stage when it comes into contact with the indenter. With the increase in the displacement of the indenter, the number of particle elements contacting the indenter increases gradually, and the number of internal bonding failure elements increases. However, because there is no obvious damage to the sea ice, the ice pressure is distributed over the entire contact surface and increases gradually. When the broken sea ice is continuously squeezed out of the contact surface, the ice sheet gradually forms an obvious middle bulge in the vertical direction. During the loading process, the middle bulge remains on the vertical section of the ice sheet. The broken sea ice accumulates on the ice sheet, and the ice load increases with the accumulation. When the ice breaking and accumulation process is stable, the ice pressure is mainly concentrated in the middle line of the ice sheet.
gradually forms an obvious middle bulge in the vertical direction. During the loading process, the middle bulge remains on the vertical section of the ice sheet. The broken sea ice accumulates on the ice sheet, and the ice load increases with the accumulation. When the ice breaking and accumulation process is stable, the ice pressure is mainly concentrated in the middle line of the ice sheet.  Figure 11 shows the distribution of the high-pressure zone on the contact surface and the change in ice pressure in the vertical direction obtained by the DEM simulation. When t = 59 s, the distribution of the high-pressure zone is mainly concentrated at the middle line of the contact surface. The width of the contact surface is the width of the indenter, the height of the contact surface is the ice thickness, and the color changes on the contact surface indicate the ice pressure. There is also a small amount of pressure in areas other than in the middle line of the ice sheet. The sea ice in the front of the ice sheet is not easy to break because of the compression caused by the accumulation of broken ice. Due to the slow loading rate, the sea ice fragments are not easily cleaned up, and the front of the ice sheet keeps in contact with the structure, forming a stable ice load. When t = 200 s, owing to the slow interaction between sea ice and structure, the broken sea ice is not removed in time, which hinders the generation of new cracks in the ice sheet. The high-pressure zone on the contact surface is not only distributed near the midline but also may appear in other places, forming multiple simultaneous high-pressure zones.  Figure 11 shows the distribution of the high-pressure zone on the contact surface and the change in ice pressure in the vertical direction obtained by the DEM simulation. When t = 59 s, the distribution of the high-pressure zone is mainly concentrated at the middle line of the contact surface. The width of the contact surface is the width of the indenter, the height of the contact surface is the ice thickness, and the color changes on the contact surface indicate the ice pressure. There is also a small amount of pressure in areas other than in the middle line of the ice sheet. The sea ice in the front of the ice sheet is not easy to break because of the compression caused by the accumulation of broken ice. Due to the slow loading rate, the sea ice fragments are not easily cleaned up, and the front of the ice sheet keeps in contact with the structure, forming a stable ice load. When t = 200 s, owing to the slow interaction between sea ice and structure, the broken sea ice is not removed in time, which hinders the generation of new cracks in the ice sheet. The high-pressure zone on the contact surface is not only distributed near the midline but also may appear in other places, forming multiple simultaneous high-pressure zones.
gradually forms an obvious middle bulge in the vertical direction. During the loading process, the middle bulge remains on the vertical section of the ice sheet. The broken sea ice accumulates on the ice sheet, and the ice load increases with the accumulation. When the ice breaking and accumulation process is stable, the ice pressure is mainly concentrated in the middle line of the ice sheet.  Figure 11 shows the distribution of the high-pressure zone on the contact surface and the change in ice pressure in the vertical direction obtained by the DEM simulation. When t = 59 s, the distribution of the high-pressure zone is mainly concentrated at the middle line of the contact surface. The width of the contact surface is the width of the indenter, the height of the contact surface is the ice thickness, and the color changes on the contact surface indicate the ice pressure. There is also a small amount of pressure in areas other than in the middle line of the ice sheet. The sea ice in the front of the ice sheet is not easy to break because of the compression caused by the accumulation of broken ice. Due to the slow loading rate, the sea ice fragments are not easily cleaned up, and the front of the ice sheet keeps in contact with the structure, forming a stable ice load. When t = 200 s, owing to the slow interaction between sea ice and structure, the broken sea ice is not removed in time, which hinders the generation of new cracks in the ice sheet. The high-pressure zone on the contact surface is not only distributed near the midline but also may appear in other places, forming multiple simultaneous high-pressure zones. In conclusion, the formation of high-pressure zones is due to the non-simultaneity of sea ice crushing. The size of the high-pressure zone does not change with the global ice pressure but is closely related to the failure mode of the sea ice. The high-pressure zone is concentrated near the middle line of the ice sheet because the ice in the middle layer always remains in contact with the structure after the ice sheet is broken by the structure. However, the midline distribution of the high-pressure zone will also be affected by the sea ice failure mode and loading rate. The results show that the DEM can effectively illustrate the occurrence and propagation of cracks in sea ice and explain the distribution of high-pressure zones in combination with the sea ice crushing process.

Influencing Factors of Ice Pressure
Many factors affect the ice pressure, including the contact area, shape of the contact In conclusion, the formation of high-pressure zones is due to the non-simultaneity of sea ice crushing. The size of the high-pressure zone does not change with the global ice pressure but is closely related to the failure mode of the sea ice. The high-pressure zone is concentrated near the middle line of the ice sheet because the ice in the middle layer always remains in contact with the structure after the ice sheet is broken by the structure. However, the midline distribution of the high-pressure zone will also be affected by the sea ice failure mode and loading rate. The results show that the DEM can effectively illustrate the occurrence and propagation of cracks in sea ice and explain the distribution of high-pressure zones in combination with the sea ice crushing process.

Influencing Factors of Ice Pressure
Many factors affect the ice pressure, including the contact area, shape of the contact surface, ice strength, and loading rate. In this study, the contact area and loading rate between the sea ice and structure were studied. The DEM was used to analyze the reason for their influence on ice pressure. The influence of the contact area on the global ice pressure (p G ) and local ice pressure (p L ) was calculated, and the accuracy of the results was verified through a comparison with the formula in the relevant ISO 19906 [31]. The formula for the global ice pressure of the structure is given in the ISO and is limited to the calculation of global ice pressure of fixed structures with the ratio of structure width to ice thickness greater than 2.0: where w is the structure width; h is the ice thickness; h 1 is the reference thickness, 1 m; m is the empirical coefficient, −0.16; n is the empirical coefficient, when h < 1.0 m, n = −0.50 + h/5, when h ≥1.0 m, n = −0.30; and C R is the strength coefficient of sea ice (MPa), which is affected by environmental factors and directly related to the sea ice strength. The reference value of C R is given in the ISO, where C R = 2.8 in the Beaufort Sea and C R = 1.8 in the Baltic Sea.
In ISO 19906, the formula for the local ice pressure is also given as where the range of contact area (a) is 0.1-10.0 m 2 . When a > 10.0 m 2 , p L = 1.48 MPa.
In the DEM simulation, the width of the structure was set at 0.5-1.5 m, while the ice thickness and loading rate remained unchanged. A comparison between the simulation results and the ISO formula is shown in Figure 12. The results show that the global and local ice pressures decrease with an increase in the contact area, and the variation law of the simulation is consistent with that of the ISO formula. The global ice pressure calculated by the DEM simulation was close to the ISO standard. However, the local pressures in the simulations are approximately twice that of the standard. This is because the formula in the ISO standard is not applicable to all ice conditions, especially for thin ice. Another important factor affecting the ice pressure is the loading rate. The influence of the loading rate on the ice load and pressure is mainly caused by the change in the sea ice failure mode [32,33]. When the loading rate is low, flow deformation (creep) of sea ice occurs, and a large number of cracks appear in the sea ice without obvious fractures. At this time, the sea ice is broken by ductile extrusion. When the loading rate is high, a large amount of sea ice breakage occurs after a small deformation, considered as the brittle crushing of sea ice. When the loading rate changes slowly, the two sea ice failure modes can change into each other in the ductile brittle transition zone. At this time, the sea ice Another important factor affecting the ice pressure is the loading rate. The influence of the loading rate on the ice load and pressure is mainly caused by the change in the sea ice failure mode [32,33]. When the loading rate is low, flow deformation (creep) of sea ice occurs, and a large number of cracks appear in the sea ice without obvious fractures. At this time, the sea ice is broken by ductile extrusion. When the loading rate is high, a large amount of sea ice breakage occurs after a small deformation, considered as the brittle crushing of sea ice. When the loading rate changes slowly, the two sea ice failure modes can change into each other in the ductile brittle transition zone. At this time, the sea ice undergoes a plastic deformation with many microcracks, which may cause shear failure, extrusion failure, and tensile failure at the same time.
In the discrete element simulation, the phenomenon in which the sea ice failure mode changes with the loading rate is also observed, as shown in Figure 13. With different loading rates and the same pressure head displacement, the degree of sea ice breakage is very different. When the loading rate is low, the crack formation speed is slow, and there is less broken sea ice. The contact surface between the sea ice and structure is relatively intact, which leads to a low local ice pressure. When the loading rate is high, the sea ice is destroyed quickly after contact with the structure, and the sea ice on the contact surface will be destroyed completely. Therefore, it is difficult for a sustained high-pressure zone to be formed, and the local ice pressure is reduced. Thus, when the loading rate increases, the local ice pressure first increases and then decreases, as shown in Figure 14. When the local ice pressure increases with the loading rate, sea ice ductile failure occurs; however, when the loading rate increases, sea ice brittle fracture occurs. The influence of the loading rate on the failure mode of sea ice can be verified by discrete element simulation, which explains why the local ice pressure also changes.

Conclusions
In this study, the interaction between sea ice and vertical structures was simulated with a spherical DEM model considering inter-particle bonding and breaking. A comparison with a field test was provided to validate the rationality of this numerical model. Both the ice load and local ice pressure distribution characteristics from the DEM results were in agreement with the experimental data. Moreover, a line-like distribution of high-pressure zones was also obtained in this simulation. Through the discrete element analysis, it can be seen that the generation and distribution of the high-pressure zones is closely related to the failure mode of sea ice. Two factors influencing the local pressures were investigated, including the contact area and load speed. The local pressure decreases with increasing area, which is consistent with the relationship in the ISO standard. On the other hand, the load rate has a positive influence on the local ice pressure, which could lead to a change in the failure modes.

Conclusions
In this study, the interaction between sea ice and vertical structures was simulated with a spherical DEM model considering inter-particle bonding and breaking. A comparison with a field test was provided to validate the rationality of this numerical model. Both the ice load and local ice pressure distribution characteristics from the DEM results were in agreement with the experimental data. Moreover, a line-like distribution of high-pressure zones was also obtained in this simulation. Through the discrete element analysis, it can be seen that the generation and distribution of the high-pressure zones is closely related to the failure mode of sea ice. Two factors influencing the local pressures were investigated, including the contact area and load speed. The local pressure decreases with increasing area, which is consistent with the relationship in the ISO standard. On the other hand, the load rate has a positive influence on the local ice pressure, which could lead to a change in the failure modes.
This paper shows that the DEM simulation can calculate the ice load of marine structures and analyze the breaking process of sea ice by observing the internal changes of sea ice. This is an advantageous way to reasonably explain the generation and variation of the high-pressure zones on a structure. Therefore, DEM simulations can be applied in the study of ice pressure to reasonably reinforce structures. However, other factors, such as the ratio of structure width to ice thickness and strength, are also necessary for a detailed analysis of ice pressure in the future.