Numerical simulation method of sea ice motion response under wind wave current action based on CFD-DEM coupling

Abstract

The present application relates to the field of sea ice movement simulation, in particular to a numerical simulation method for sea ice movement response under the action of wind wave flow based on CFD-DEM coupling, which comprises: forming a multi-scale heterogeneous numerical base model of wind wave flow-sea ice; based on the multi-scale heterogeneous numerical base model, dividing and matching the mesoscale unit through the mesoscale bridging algorithm, completing the flow field force transmission of the mesoscale bridging unit by fusing the Liutex vorticity mechanism, and embedding the percolation mechanism to solve the DEM particle behavior characterization model, and obtaining the numerical simulation result after the solution is terminated. The present application realizes the matching of CFD fluid phase and DEM by means of the mesoscale bridging algorithm, at the same time, the Liutex vorticity mechanism is fused to improve the accuracy of the flow field force transmission to the sea ice particle phase, and the percolation mechanism is embedded in the mesoscale bridging unit to realize the bidirectional dynamic coupling of the sea ice particle phase and the fluid phase.

Description

Technical Field

[0001] This invention relates to the field of sea ice motion simulation, and more specifically to a numerical simulation method for sea ice motion response under the action of wind, waves and current based on CFD-DEM coupling. Background Technology

[0002] Current methods for simulating sea ice motion under the influence of wind, waves, and currents suffer from numerous technical deficiencies in the core aspect of handling the interaction between the fluid phase and the sea ice particle phase. In terms of macro- and micro-scale integration, the lack of effective meso-scale bridging methods leads to discontinuities in the transfer of macro- and micro-scale parameters and force chain interactions, preventing the accurate transmission of fluid phase forces to the sea ice particle phase. Furthermore, existing methods are significantly inadequate in identifying and quantifying vortex patterns at ice-water and wind-water interfaces. They fail to fully consider the actual impact of vortex forces on sea ice motion, and the force transmission chain in the flow field is incomplete. They simply calculate drag and lift, neglecting the fluid forces under vortex action, resulting in significant deviations in the calculation of the excitation load on sea ice particles.

[0003] Furthermore, existing simulation methods lack effective force chain and particle force analysis feedback mechanisms, making it difficult to achieve dynamic coupling feedback between the sea ice particle phase and the fluid phase. The coupled field model cannot respond in real-time to the mechanical changes of sea ice particles. Ultimately, this results in significant deviations between the simulated sea ice motion response under the influence of wind, waves, and currents and the actual sea ice motion state in the sea area. It fails to accurately reflect the actual behavior of sea ice, such as displacement, velocity, moment, and breakup evolution, and cannot provide reliable numerical references for polar engineering protection and sea ice disaster early warning. Summary of the Invention

[0004] This invention addresses the technical problems existing in the prior art by providing a numerical simulation method for sea ice motion response under the action of wind, waves and current based on CFD-DEM coupling.

[0005] The technical solution of this invention to solve the above-mentioned technical problems is as follows: a numerical simulation method for sea ice motion response under the action of wind, waves and current based on CFD-DEM coupling, the method comprising: S1. Based on the energy distribution data of wind, wave and current monitored in the target sea area, identify areas with intense energy exchange, and construct the CFD grid of the numerical water tank using the energy weighted grid densification method. Based on the sea ice morphology of broken ice and layered ice, and the sea ice particle size and layered ice thickness data obtained by a third-party engine, construct a DEM hierarchical particle system. Set the boundary dynamic loading in the numerical water tank, and fuse the constructed CFD grid and DEM hierarchical particle system to form a multi-scale heterogeneous numerical basis model of wind, wave and current-sea ice. S2. Based on the multi-scale heterogeneous numerical basis model, the mesoscale unit is divided and matched by the mesoscale bridging algorithm to obtain the mesoscale bridging unit. After integrating the Liutex vorticity mechanism to complete the flow field force transmission construction of the mesoscale bridging unit, the percolation mechanism is embedded to perform force chain and particle force analysis feedback in each mesoscale bridging unit to obtain the dynamic coupled field interaction model. S3. The obtained dynamic coupled field interaction model is coupled with particle mechanical properties and particle dynamic properties to transform the dynamic coupled field interaction model into a DEM particle behavior characterization model with multi-physical property coupling. Then, the fluid phase and force balance of the DEM particle behavior characterization model is solved and synchronized with the time step. After the solution is terminated, the numerical simulation results of sea ice motion response under the action of wind, waves and current are obtained.

[0006] In a preferred embodiment, S1 is based on the wind, wave and current energy distribution data of the target sea area monitored by a third-party engine. By statistically analyzing the energy density values ​​of each region in the wind, wave and current energy distribution data, regions with energy density values ​​exceeding the average energy density value are identified as regions of intense energy exchange. Based on the grid refinement criteria of multi-scale simulation, the grid scale of the region with intense energy exchange is set to one-tenth of the minimum vortex scale of the fluid. Regions with energy density below the average value are identified as far-field energy stable regions, and their grid scale is set to ten times that of the region with intense energy exchange. The grid scale is gradually transitioned from the grid scale of the region with intense energy exchange to the grid scale of the far-field energy stable region by linear interpolation, thus completing the CFD grid construction of the numerical water tank. After constructing the CFD mesh for the numerical water tank, S1 constructs the DEM-level particle system based on sea ice grain size and layer ice thickness data obtained from a third-party engine, including the following specific settings: The spherical particle classification model is used to classify the ice fragments. Based on the particle size distribution range of sea ice fragments, three continuous intervals are selected from the particle size distribution range of sea ice fragments according to the principle of equal interval division. These intervals are defined as the first particle size class, the second particle size class, and the third particle size class, respectively. The quantity distribution ratio of the first particle size class, the second particle size class, and the third particle size class is set according to the proportion of ice fragments in the three continuous intervals of sea ice fragments, thus completing the simulation of the natural randomness of ice fragment particle size. The ice layer was modeled using a hexagonal close-packed bonded particle model. Based on standard spherical DEM particles, the ice layer was bonded together using a parallel bonding model to form an ice layer plate. The ice thickness was set according to the common monitoring range of the target sea area, including three continuous ranges: thin ice layer, medium-thickness ice layer, and thick ice layer. The thin ice layer, medium-thickness ice layer, and thick ice layer were formed by stacking corresponding numbers of spherical particles along the thickness direction.

[0007] In a preferred embodiment, after constructing the CFD grid and DEM hierarchical particle system, acquisition equipment including wind speed sensors, wave monitors, and underwater current meters is deployed in the target sea area. The acquired wind speed values ​​are directly used as the calculation input for wind field shear stress, and the acquired wave height values ​​are converted into the velocity inlet parameters of the wave generator. The flow velocity values ​​are converted into the pressure inlet parameters of the flow field. Boundary dynamic loading is set at the inlet of the numerical water tank, including loading the waves using the velocity inlet method, the flow field using the pressure inlet method, and the wind field using the shear stress inlet method. The initial time step of loading is uniformly set. Preferably, it can be set to one percent of the fluid wave period. This step can ensure calculation accuracy and control the consumption of computing resources. After completing the boundary dynamic loading, the constructed CFD mesh and DEM hierarchical particle system are integrated to form a multi-scale heterogeneous numerical basis model of wind, wave and current-sea ice.

[0008] In a preferred embodiment, S2, based on a multi-scale heterogeneous numerical basis model, constructs a one-to-one mesoscale bridging unit by connecting the smallest micro-element of the CFD mesh fluid phase with the particle clusters of the DEM sea ice phase. Specifically: The number of particles in the ice fragment phase is determined according to the particle size class. The first particle size class particle cluster consists of three particles, the second particle size class particle cluster consists of four particles, and the third particle size class particle cluster consists of five particles. The spatial size of the particle cluster is consistent with the size of the smallest micro-element of the fluid. The particle clusters of the layered ice phase are a cohesive particle unit of the smallest structural unit of the corresponding layered ice. Within each mesoscale bridging unit, the smallest micro-element of the fluid phase and the particle cluster of the DEM sea ice phase are in the same spatial coordinate system, and the boundary of the smallest micro-element coincides with the outer boundary of the particle cluster, thereby achieving precise scale matching between the flow field micro-element and the sea ice particle cluster.

[0009] In a preferred embodiment, S2 uses the Liutex vortex mechanism to identify flow field vortices at the ice-water interface and wind-water interface within each mesoscale bridging unit. Specifically, it calculates the velocity components of each discrete point within the smallest infinitesimal element of the fluid phase, then calculates the velocity gradient tensor of each discrete point using the central difference method, decomposes the velocity gradient tensor into a symmetric part and an antisymmetric part, uses the antisymmetric part as the vortex tensor, and extracts the eigenvalues ​​of the vortex tensor. The magnitude of the eigenvalues ​​of the vortex tensor is used as the vortex rotation angular velocity of the corresponding discrete point, and the average value of the vortex rotation angular velocities of all discrete points is used as the average rotation angular velocity of the vortex within the corresponding mesoscale bridging unit. Twice the average rotation angular velocity is used as the rotation intensity of the vortex, thus completing the vortex identification and rotation intensity calculation of the mesoscale bridging unit. Next, the projected area of ​​the DEM sea ice phase particle clusters is calculated. Specifically, the three-dimensional model of the DEM sea ice phase particle clusters is projected vertically towards the direction of the flow field vortex to obtain a two-dimensional projected profile. The boundary of the projected profile is discretized into several vertices using the polygon area integration method. The vertices are arranged in a clockwise order. The area of ​​each small triangle is calculated using the coordinates of adjacent vertices. The areas of all small triangles are summed to obtain the projected area of ​​the sea ice phase particle clusters. After obtaining the projected area of ​​the sea ice particle cluster, the vortex rotation intensity is multiplied by the projected area of ​​the sea ice particle cluster, and then multiplied by the fluid density to complete the Liutex vortex force calculation. The fluid density is the median value of the standard density range of seawater. Furthermore, under general circumstances, when performing vortex force calculation in this application, the median value of the standard value range for similar vortex-particle interaction problems in the field of fluid mechanics can also be introduced as an empirical correction coefficient to further correct the vortex force. Finally, the direction of the vortex core is determined by the direction of the eigenvector of the vortex tensor, and the direction of the vortex force is perpendicular to the vortex rotation axis.

[0010] The drag force of the fluid on sea ice was calculated based on Archimedes' principle and the relative velocity of the fluid. The lift force of the fluid on sea ice was calculated based on the Magnus effect. The drag force, lift force, and Liutex vortex force were combined into an intermediate resultant force according to the vector triangle rule. The intermediate resultant force was then combined with the Liutex vortex force to obtain the total fluid force, thus establishing a flow field force transmission chain. The drag force is positively correlated with the fluid density, particle cluster volume, and the square of the fluid velocity relative to the particle cluster. The lift force is positively correlated with the fluid density, particle cluster volume, fluid relative velocity, and vortex rotation angular velocity. All relevant parameters were acquired in real time by the data acquisition equipment, and the calculation was based on Archimedes' principle and the Magnus effect.

[0011] In a preferred embodiment, S3 embeds a percolation mechanism within each mesoscale bridging unit, specifically as follows: The total number of sea ice particles in each mesoscale bridging unit is counted, and the number of possible contact pairs between all particles is calculated. The calculation logic for the number of possible contact pairs is to multiply the total number of particles by the difference between the total number of particles and the total number of particles minus one, and then divide by two. The actual number of contacting particle pairs is determined by the DEM particle contact detection algorithm, and the ratio of the actual number of contacting particle pairs to the number of possible contact pairs is used as the basic contact probability. For each actual contacting particle pair, the contact force and the maximum allowable contact force are calculated. The contact force and the maximum allowable contact force are determined by the compressive strength of the sea ice in the target sea area to obtain the contact strength coefficient. The basic contact probability is multiplied by the average of all contact strength coefficients to obtain the percolation probability of the force chain. A path search algorithm is used to find the connected force chains within the mesoscale bridging unit. The continuous chain formed by particle contact from one end to the other within the mesoscale bridging unit is taken as the connected force chain. The length of all connected force chains is counted. The length refers to the sum of the particle center-to-center distances contained in each connected force chain. The average length of all connected force chains is calculated, and the average length is divided by the side length of the mesoscale unit to obtain the connectivity coefficient of the force chain. The compressive strength of sea ice is obtained, and the maximum allowable contact force between particles is calculated based on the contact area of ​​the particles. The maximum allowable contact force is equal to the compressive strength multiplied by the contact area. The maximum allowable contact force is used as the force chain breakage threshold. When the actual contact force between particles in the mesoscale bridging unit exceeds the force chain breakage threshold, it is determined to be a force chain breakage. The parameters of the force chain are transferred to the mesoscale bridging end of the CFD fluid phase. When the percolation probability of the force chain of ice fragments exceeds the critical percolation probability of the fluid-particle coupling, the viscosity coefficient of the fluid phase in the corresponding mesoscale unit is adjusted according to the magnitude of the force chain connectivity coefficient. When the contact force of the force chain at the bonding surface of layered ice particles exceeds the fracture threshold and percolation failure occurs, the distribution of the fractured layered ice particles is transferred to the CFD fluid phase, and the volume fraction of layered ice particles in the unit is statistically analyzed. The multiscale heterogeneous numerical basis model is transformed into a dynamic coupled field interaction model.

[0012] In a preferred embodiment, S3 characterizes the mechanical properties of ice fragments and layered ice particles based on a dynamic coupled field interaction model, specifically as follows: The HertzMindlin non-adhesive contact model is adopted. During the relative sliding of ice fragments, the product of the normal contact force between the particles and the coefficient of friction is used as the sliding friction force. The velocity decay of the ice fragments after the collision is calculated based on the sliding friction force. The velocity decay calculated by the sliding friction force is obtained from Newton's laws of mechanics and force decomposition, so it will not be elaborated here. The HertzMindlin parallel bonding model was used to calibrate the bonding strength and bonding stiffness of the layered ice particles. The tensile stress was obtained by dividing the tensile force on the bonding surface of the layered ice by the area of ​​the bonding surface. The tensile stress was compared with the bonding strength. When the actual tensile stress exceeded the bonding strength, the bonding of the layered ice particles was triggered to break. The broken layered ice particles automatically switched to the HertzMindlin non-bonded contact model of the crushed ice particles. The friction coefficient and the coefficient of restitution were adopted from the parameters of the crushed ice particles. The calibrated bonding strength and bonding stiffness of the layered ice particles were derived based on the Young's modulus of sea ice materials.

[0013] In a preferred embodiment, after completing the mechanical property characterization, S3 uses the wind field shear stress and total fluid force collected in S2 as external excitation loads for the DEM particles and applies them to the centroid of the particles. The total fluid force is the sum of the wind field shear stress and the flow field force generated by the wave generator. The external excitation loads are distributed to the particle cluster composed of multiple particles according to their mass ratio. Furthermore, the distribution mechanism in this application assigns a load to each particle equal to the total load multiplied by the ratio of the particle's mass to the total mass of the particle cluster. The initial motion state of the particles is determined, that is, the initial displacement, initial velocity, and initial angular velocity are all set to the initial stationary state. All particles are traversed, and the external loads on each particle, including the total fluid force and wind shear force, and the interparticle interaction forces, including the normal contact force, sliding friction force, and adhesive force, are calculated. All forces are accumulated in the vector direction to obtain the resultant force of the particles. The resultant force is divided by the mass of the particle to obtain the translational acceleration of the particle, which is obtained by Newton's laws. The translational acceleration is integrated over time, with the integration interval being the current time step, to obtain the velocity change. The initial velocity plus the velocity change gives the current velocity of the particle. The current velocity is integrated over time to obtain the displacement change. The initial displacement plus the displacement change gives the current position of the particle. All the torques acting on the particle are calculated, where torque refers to the torque generated by the interaction forces between particles. The total torque is obtained by accumulating them, and finally, the particle dynamics characteristics including displacement, velocity, and torque are obtained. The particle dynamics and mechanical properties are assigned to the dynamic coupled field interaction model, which is then transformed into a DEM particle behavior characterization model with multi-physical property coupling, thus completing the particle phase solution.

[0014] In a preferred embodiment, S3 discretizes the CFD mesh constructed by S1 using the finite volume method, decomposing the fluid continuity equation and momentum conservation equation into mass conservation and force balance equations for each mesh element, and then performing mass conservation, force balance, and particle phase solutions. Specifically: The inflow and outflow masses of each grid cell within the current time step are calculated by multiplying the flow velocity of each face of the grid cell by the area and the time step. The inflow mass minus the outflow mass is the mass change within the grid cell, thus completing the mass conservation solution. Based on the pressure difference between adjacent grid cells, the pressure, viscous force, and inertial force on each grid cell are calculated. By solving the force balance equations, the velocity and pressure of the grid cells are obtained through the vector sum of all forces being zero. By iterating through all grid cells, the flow field is solved once, and the numerical simulation parameters of the flow field, including velocity, pressure, and vortex, are obtained, thus completing the force balance solution. Furthermore, the viscous force is calculated based on the RNGk-ε model, and the inertial force is calculated by multiplying the fluid density by the acceleration. The particle phase solution is the S3 step, which involves solving the DEM particle model. The detailed steps for solving the dynamic parameters are strictly followed in step three to calculate the real-time motion state of each sea ice particle and obtain parameters such as displacement, velocity, and angular velocity of the sea ice particles.

[0015] In a preferred embodiment, step S3: After solving at each time step, the residual coefficient of the CFD flow field is calculated in real time. Specifically, for each grid cell, the theoretical velocity of the grid cell is calculated based on the fluid dynamics theory formula and derived from the Navier-Stokes equation. The difference between the calculated velocity and the theoretical velocity is calculated to obtain the absolute difference. The absolute difference is divided by the theoretical velocity to obtain the residual ratio of the grid cell. All grid cells are traversed, and the largest residual ratio is taken as the residual coefficient of the flow field. The solution step size is adjusted according to the residual coefficient.

[0016] The beneficial effects of this invention are as follows: relying on the mesoscale bridging algorithm, the one-to-one precise matching of the smallest micro-element of the CFD fluid phase and the DEM sea ice phase particle cluster is realized, and an effective channel for the transfer of macro- and micro-scale parameters and the interaction of force chains is established. This breaks the discontinuity problem of macro- and micro-scale connection in the traditional coupling method. At the same time, the Liutex vortex mechanism is integrated to complete the accurate identification of ice-water and wind-water interface vortices and the quantitative calculation of vortex forces, which improves the accuracy of the flow field force transfer to the sea ice particle phase. Furthermore, the percolation mechanism is embedded in the mesoscale bridging unit to capture the contact state, force chain evolution and fragmentation behavior of sea ice particles, and dynamically transfer the force chain parameters to the CFD fluid phase, realizing the bidirectional dynamic coupling between the sea ice particle phase and the fluid phase. The finite volume method was used to discretize the CFD mesh, achieving accurate solutions for the mass conservation and force balance of the fluid phase. Simultaneously, the dynamic solutions of the particle phase and the fluid phase were synchronized in time step, and the solution step was dynamically adjusted by calculating the flow field residual coefficient in real time. This effectively controlled the numerical calculation error, avoided the numerical oscillation problem that easily occurs in traditional solution methods, improved the stability and overall accuracy of the simulation calculation, and ensured the consistency and reliability of the fluid-solid two-phase solution results. Attached Figure Description

[0017] Figure 1 This is a flowchart of the present invention; Figure 2 This is a schematic diagram of the mesh generation for the numerical water tank model of the present invention; Figure 3 This is a schematic diagram illustrating the simulation of broken ice using a polyhedral particle model in the embodiment. Figure 4 This is a schematic diagram of the spherical particle model used in this invention to simulate broken ice. Figure 5 This is a schematic diagram of ice fragments bonded by DEM particles according to the present invention; Figure 6 This is a distribution map of crushed ice constructed based on the DEM-based particle classification system of this invention. Detailed Implementation

[0018] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0019] As attached Figure 1 As shown, this embodiment provides a numerical simulation method for sea ice motion response under the influence of wind, waves and currents based on CFD-DEM coupling, including the following steps: S1. Based on the energy distribution data of wind, wave and current monitored in the target sea area, identify areas with intense energy exchange, and construct the CFD grid of the numerical water tank using the energy weighted grid densification method. Based on the sea ice morphology of broken ice and layered ice, and the sea ice particle size and layered ice thickness data obtained by a third-party engine, construct a DEM hierarchical particle system. Set the boundary dynamic loading in the numerical water tank, and fuse the constructed CFD grid and DEM hierarchical particle system to form a multi-scale heterogeneous numerical basis model of wind, wave and current-sea ice. S1 is based on the wind, wave and current energy distribution data of the target sea area monitored by a third-party engine. By statistically analyzing the energy density values ​​of each region in the wind, wave and current energy distribution data, the region with an energy density value exceeding the average energy density value is identified as a region of intense energy exchange. Based on the grid refinement criteria of multi-scale simulation, the grid scale of the region with intense energy exchange is set to one-tenth of the minimum vortex scale of the fluid. Regions with energy density below the average value are identified as far-field energy stable regions, and their grid scale is set to ten times that of the region with intense energy exchange. The grid scale is gradually transitioned from the grid scale of the region with intense energy exchange to the grid scale of the far-field energy stable region by linear interpolation, thus completing the CFD grid construction of the numerical water tank. After constructing the CFD mesh for the numerical water tank, S1 constructs the DEM-based particle size distribution system based on sea ice grain size and layer ice thickness data obtained from a third-party engine, including the following specific settings: The spherical particle classification model is used to classify the ice fragments. Based on the particle size distribution range of sea ice fragments, three continuous intervals are selected from the particle size distribution range of sea ice fragments according to the principle of equal interval division. These intervals are defined as the first particle size class, the second particle size class, and the third particle size class, respectively. The quantity distribution ratio of the first particle size class, the second particle size class, and the third particle size class is set according to the proportion of ice fragments in the three continuous intervals of sea ice fragments, thus completing the simulation of the natural randomness of ice fragment particle size. The ice layer adopts a hexagonal close-packed bonded particle model, based on standard spherical DEM particles, and is bonded into ice layer plates through a parallel bonding model. The ice thickness of the ice layer is set according to the common intervals monitored in the target sea area, including three continuous intervals: thin ice layer, medium-thickness ice layer, and thick ice layer. The thin ice layer, medium-thickness ice layer, and thick ice layer are formed by stacking the corresponding number of spherical particles along the thickness direction. The use of standard spherical DEM particles as the basis refers to a particle diameter that is one-fifth of the smallest microstructure scale of the ice layer. This is based on the microstructure data of ice layer in the global average ice layer. Furthermore, this application also reserves microcrack units between ice layer particles. The diameter of the microcrack unit is one-tenth of the diameter of the standard spherical particle. The distribution density of the microcrack unit is set according to the average density of microcracks in the ice layer, simulating the microcrack characteristics that naturally exist in polar ice layer.

[0020] After constructing the CFD grid and DEM hierarchical particle system, data acquisition equipment including wind speed sensors, wave monitors, and underwater current meters were deployed in the target sea area. The collected wind speed values ​​were directly used as the calculation input for wind field shear stress, and the collected wave height values ​​were converted into the velocity inlet parameters of the wave generator. The flow velocity values ​​were converted into the pressure inlet parameters of the flow field. Boundary dynamic loading was set at the inlet of the numerical water tank, including loading the waves using the velocity inlet method, the flow field using the pressure inlet method, and the wind field using the shear stress inlet method. The initial time step of loading was uniformly set. Preferably, it can be set to one percent of the fluid wave period. This step can ensure calculation accuracy while controlling the consumption of computing resources. After completing the boundary dynamic loading, the constructed CFD mesh and DEM hierarchical particle system are integrated to form a multi-scale heterogeneous numerical basis model of wind, wave and current-sea ice.

[0021] In some other specific embodiments, this application also includes using the fluid volume method to process the gas-water two-phase interface, selecting the RNG k-ε turbulence model to characterize the fluid turbulence characteristics, setting a field coupling transition layer at the interface of wind field, flow field and wave field, the thickness of the transition layer being twenty times the grid scale of the region with intense energy exchange, and the grid scale within the transition layer being linearly transitioned to avoid numerical oscillations caused by the direct superposition of multiple fields. Furthermore, to meet the simulation requirements of physical model experiments for complex hydrodynamic environments, this application innovatively designs a bottom-supported structure in the hydrodynamic circulation system to achieve efficient water circulation. Therefore, regarding the above technical solution, this application proposes the following specific experimental operations for the construction of the numerical tank and CFD mesh: The experimental setup uses an adjustable false bottom system at a height of 0.6m at the bottom of the tank, and wave-damping banks are configured at the front and rear ends of the flume, with the distal wave-damping bank reaching a length of 2.4m, effectively suppressing wave reflection interference. In the numerical simulation stage, to balance computational accuracy and efficiency, the geometric parameters of the numerical flume are optimized and determined as follows: longitudinal 14m, transverse 2m, and vertical 0.6m.

[0022] A refined strategy was adopted for boundary condition setting: a velocity inlet boundary condition was used at the inlet to simulate wave generation; slip boundary conditions were set on the sidewalls to minimize the influence of wall viscous dissipation; a pressure boundary condition was used at the outlet, and a 3.33m long wave-damping zone was configured, combined with an active absorption algorithm to eliminate reflected wave effects. (See attached...) Figure 2 As shown, in the process of mesh generation, in addition to the above-mentioned technical solutions, local densification processing was also implemented for the free liquid surface region, and a dynamic adaptive mesh was introduced to ensure high-precision capture of the gas-liquid two-phase interface. In some other specific embodiments, taking ice crushing simulation as an example, this application further discloses the technical characteristics and effects of using a spherical particle model to simulate ice crushing, as follows: In numerical simulation studies of ice fragment dynamics, when using a polyhedral particle model to discretize ice fragments, such as... Figure 3 Numerical simulation results show that under dynamic loads such as waves, the edges of ice fragments exhibit significant verticalization, and the free motion behavior of ice fragments under dynamic loads cannot be effectively reproduced. Please refer to... Figure 4 In contrast, the discrete element method system for broken ice, constructed based on a spherical particle model, simulates the internal mechanical properties of broken ice through the bonding effect between particles, successfully eliminating the non-physical phenomenon of the edges of broken ice standing upright.

[0023] In summary, this application incorporates a DEM-based ice block model into its numerical water tank model. The ice block shapes are defined as square and elliptical; please refer to [reference needed]. Figure 5 All of them are constructed from DEM particles with a diameter of 10mm through bonding. In the ice fragment group, square and oval ice blocks each account for 50% of the total.

[0024] Furthermore, in specific experiments, the distribution pattern of the ice shards was as follows: Figure 6As shown, the ice shard area is rectangular, measuring 4500mm × 2000mm, with its leading edge 1m from the inlet of the numerical water tank. To accurately obtain wave propagation characteristics and the wave-damping effect of the ice shard area, this invention utilizes CFD-DEM bidirectional coupling, fully considering the bidirectional interaction between DEM particles and fluid. Based on this, wave height monitoring points 1 and 2 are respectively located 500mm upstream and 500mm downstream of the ice shard area to collect wave height data at both points in real time.

[0025] S2. Based on the multi-scale heterogeneous numerical basis model, the mesoscale unit is divided and matched by the mesoscale bridging algorithm to obtain the mesoscale bridging unit. After integrating the Liutex vorticity mechanism to complete the flow field force transmission construction of the mesoscale bridging unit, the percolation mechanism is embedded to perform force chain and particle force analysis feedback in each mesoscale bridging unit to obtain the dynamic coupled field interaction model. S2, based on a multi-scale heterogeneous numerical basis model, constructs a one-to-one mesoscale bridging unit by linking the smallest micro-elements of the CFD mesh fluid phase with the particle clusters of the DEM sea ice phase. Specifically: The number of particles in the ice fragment phase is determined according to the particle size class. The first particle size class particle cluster consists of three particles, the second particle size class particle cluster consists of four particles, and the third particle size class particle cluster consists of five particles. The spatial size of the particle cluster is consistent with the size of the smallest micro-element of the fluid. The particle clusters of the layered ice phase are a cohesive particle unit of the smallest structural unit of the corresponding layered ice. Within each mesoscale bridging unit, the smallest micro-element of the fluid phase and the particle cluster of the DEM sea ice phase are in the same spatial coordinate system, and the boundary of the smallest micro-element coincides with the outer boundary of the particle cluster, thereby achieving precise scale matching between the flow field micro-element and the sea ice particle cluster. Furthermore, this application connects macroscopic continuous media and microscopic discrete particles through a mesoscale bridging algorithm, thereby realizing macro- and micro-scale parameter transfer and force chain interaction. After the mesoscale bridging unit is divided and the parameter transfer and force chain interaction are established, it is necessary to build the flow field force transfer. For this purpose, this application provides a further solution. S2 uses the Liutex vortex mechanism to identify vortices in the flow field at the ice-water and wind-water interfaces within each mesoscale bridging unit. Specifically, it calculates the velocity components of each discrete point within the smallest infinitesimal element of the fluid phase, then calculates the velocity gradient tensor of each discrete point using the central difference method. The velocity gradient tensor is decomposed into a symmetric part and an antisymmetric part. The antisymmetric part is used as the vortex tensor, and the eigenvalues ​​of the vortex tensor are extracted. The magnitude of the eigenvalues ​​of the vortex tensor is used as the vortex rotation angular velocity of the corresponding discrete point. The average vortex rotation angular velocity of all discrete points is used as the average rotation angular velocity of the vortex within the corresponding mesoscale bridging unit. Twice the average rotation angular velocity is used as the rotation intensity of the vortex, thus completing the vortex identification and rotation intensity calculation of the mesoscale bridging unit. Next, the projected area of ​​the DEM sea ice phase particle clusters is calculated. Specifically, the three-dimensional model of the DEM sea ice phase particle clusters is projected vertically towards the direction of the flow field vortex to obtain a two-dimensional projected profile. The boundary of the projected profile is discretized into several vertices using the polygon area integration method. The vertices are arranged in a clockwise order. The area of ​​each small triangle is calculated using the coordinates of adjacent vertices. The areas of all small triangles are summed to obtain the projected area of ​​the sea ice phase particle clusters. After obtaining the projected area of ​​the sea ice particle cluster, the vortex rotation intensity is multiplied by the projected area of ​​the sea ice particle cluster, and then multiplied by the fluid density to complete the Liutex vortex force calculation. The fluid density is the median value of the standard density range of seawater. Furthermore, under general circumstances, when performing vortex force calculation in this application, the median value of the standard value range for similar vortex-particle interaction problems in the field of fluid mechanics can also be introduced as an empirical correction coefficient to further correct the vortex force. Finally, the direction of the vortex core is determined by the direction of the eigenvector of the vortex tensor, and the direction of the vortex force is perpendicular to the vortex rotation axis.

[0026] The drag force of the fluid on sea ice was calculated based on Archimedes' principle and the relative velocity of the fluid. The lift force of the fluid on sea ice was calculated based on the Magnus effect. The drag force, lift force, and Liutex vortex force were combined into an intermediate resultant force according to the vector triangle rule. The intermediate resultant force was then combined with the Liutex vortex force to obtain the total fluid force, thus establishing a flow field force transmission chain. The drag force is positively correlated with the fluid density, particle cluster volume, and the square of the fluid velocity relative to the particle cluster. The lift force is positively correlated with the fluid density, particle cluster volume, fluid relative velocity, and vortex rotation angular velocity. All relevant parameters were acquired in real time by the data acquisition equipment, and the calculation was based on Archimedes' principle and the Magnus effect.

[0027] S3. The obtained dynamic coupled field interaction model is coupled with particle mechanical properties and particle dynamic properties to transform the dynamic coupled field interaction model into a DEM particle behavior characterization model with multi-physical property coupling. Then, the fluid phase and force balance of the DEM particle behavior characterization model is solved and synchronized with the time step. After the solution is terminated, the numerical simulation results of sea ice motion response under the action of wind, waves and current are obtained.

[0028] S3 embeds a percolation mechanism within each mesoscale bridging unit, specifically as follows: The total number of sea ice particles in each mesoscale bridging unit is counted, and the number of possible contact pairs between all particles is calculated. The calculation logic for the number of possible contact pairs is to multiply the total number of particles by the difference between the total number of particles and the total number of particles minus one, and then divide by two. The actual number of contacting particle pairs is determined by the DEM particle contact detection algorithm, and the ratio of the actual number of contacting particle pairs to the number of possible contact pairs is used as the basic contact probability. For each actual contacting particle pair, the contact force and the maximum allowable contact force are calculated. The contact force and the maximum allowable contact force are determined by the compressive strength of the sea ice in the target sea area to obtain the contact strength coefficient. The basic contact probability is multiplied by the average of all contact strength coefficients to obtain the percolation probability of the force chain. A path search algorithm is used to find the connected force chains within the mesoscale bridging unit. The continuous chain formed by particle contact from one end to the other within the mesoscale bridging unit is taken as the connected force chain. The length of all connected force chains is counted. The length refers to the sum of the particle center-to-center distances contained in each connected force chain. The average length of all connected force chains is calculated, and the average length is divided by the side length of the mesoscale unit to obtain the connectivity coefficient of the force chain. The compressive strength of sea ice is obtained, and the maximum allowable contact force between particles is calculated based on the contact area of ​​the particles. The maximum allowable contact force is equal to the compressive strength multiplied by the contact area. The maximum allowable contact force is used as the force chain breakage threshold. When the actual contact force between particles in the mesoscale bridging unit exceeds the force chain breakage threshold, it is determined to be a force chain breakage. The parameters of the force chain are transferred to the mesoscale bridging end of the CFD fluid phase. When the percolation probability of the force chain of ice fragments exceeds the critical percolation probability of the fluid-particle coupling, the viscosity coefficient of the fluid phase in the corresponding mesoscale unit is adjusted according to the magnitude of the force chain connectivity coefficient. When the contact force of the force chain at the bonding surface of layered ice particles exceeds the fracture threshold and percolation failure occurs, the distribution of the fractured layered ice particles is transferred to the CFD fluid phase, and the volume fraction of layered ice particles in the unit is statistically analyzed. The multiscale heterogeneous numerical basis model is transformed into a dynamic coupled field interaction model.

[0029] Furthermore, the critical value of the critical percolation probability of fluid-particle coupling is the standard critical value used in percolation theory to describe the connectivity of dense particle systems. When adjusting the viscosity coefficient of the fluid phase in the mesoscale unit according to the magnitude of the force chain connectivity coefficient, it should be noted that the larger the connectivity coefficient, the higher the proportion of viscosity coefficient increase, thus simulating the obstructive effect of sea ice particle groups on the fluid flow field.

[0030] Furthermore, the volume fraction of ice particles within a unit is determined by the ratio of the total particle volume to the unit volume. The higher the volume fraction, the greater the reduction in the local velocity of the fluid phase, thus realizing the dynamic influence of sea ice breakup on the fluid flow field.

[0031] S3 characterizes the mechanical properties of ice fragments and layered ice particles based on a dynamic coupled field interaction model, specifically as follows: The HertzMindlin non-adhesive contact model is adopted. During the relative sliding of ice fragments, the product of the normal contact force between the particles and the coefficient of friction is used as the sliding friction force. The velocity decay of the ice fragments after the collision is calculated based on the sliding friction force. The velocity decay calculated by the sliding friction force is obtained from Newton's laws of mechanics and force decomposition, so it will not be elaborated here. The HertzMindlin parallel bonding model was used to calibrate the bonding strength and bonding stiffness of the layered ice particles. The tensile stress was obtained by dividing the tensile force on the bonding surface of the layered ice by the area of ​​the bonding surface. The tensile stress was compared with the bonding strength. When the actual tensile stress exceeded the bonding strength, the bonding of the layered ice particles was triggered to break. The broken layered ice particles automatically switched to the HertzMindlin non-bonded contact model of the crushed ice particles. The friction coefficient and the coefficient of restitution were adopted from the parameters of the crushed ice particles. The calibrated bonding strength and bonding stiffness of the layered ice particles were derived based on the Young's modulus of sea ice materials.

[0032] After completing the mechanical property characterization, S3 uses the wind field shear stress and total fluid force collected by S2 as the external excitation load of the DEM particles and applies it to the centroid of the particles. The total fluid force is the sum of the wind field shear stress and the flow field force generated by the wave generator. The external excitation load is distributed to the particle cluster composed of multiple particles according to the mass ratio. Furthermore, the distribution mechanism in this application assigns a load to each particle equal to the total load multiplied by the ratio of the mass of the particle to the total mass of the particle cluster. The initial motion state of the particles is determined, that is, the initial displacement, initial velocity, and initial angular velocity are all set to the initial stationary state. All particles are traversed, and the external loads on each particle, including the total fluid force and wind shear force, and the interparticle interaction forces, including the normal contact force, sliding friction force, and adhesive force, are calculated. All forces are accumulated in the vector direction to obtain the resultant force of the particles. The resultant force is divided by the mass of the particle to obtain the translational acceleration of the particle, which is obtained by Newton's laws. The translational acceleration is integrated over time, with the integration interval being the current time step, to obtain the velocity change. The initial velocity plus the velocity change gives the current velocity of the particle. The current velocity is integrated over time to obtain the displacement change. The initial displacement plus the displacement change gives the current position of the particle. All the torques acting on the particle are calculated, where torque refers to the torque generated by the interaction forces between particles. The total torque is obtained by accumulating them, and finally, the particle dynamics characteristics including displacement, velocity, and torque are obtained. The particle dynamics and mechanical properties are assigned to the dynamic coupled field interaction model, which is then transformed into a DEM particle behavior characterization model with multi-physical property coupling, thus completing the particle phase solution.

[0033] S3 uses the finite volume method to discretize the CFD mesh constructed by S1, decomposing the fluid continuity equation and momentum conservation equation into mass conservation and force balance equations for each mesh element, and then performing mass conservation, force balance, and particle phase solutions. Specifically: The inflow and outflow masses of each grid cell within the current time step are calculated by multiplying the flow velocity of each face of the grid cell by the area and the time step. The inflow mass minus the outflow mass is the mass change within the grid cell, thus completing the mass conservation solution. Based on the pressure difference between adjacent grid cells, the pressure, viscous force, and inertial force on each grid cell are calculated. By solving the force balance equations, the velocity and pressure of the grid cells are obtained through the vector sum of all forces being zero. By iterating through all grid cells, the flow field is solved once, and the numerical simulation parameters of the flow field, including velocity, pressure, and vortex, are obtained, thus completing the force balance solution. Furthermore, the viscous force is calculated based on the RNGk-ε model, and the inertial force is calculated by multiplying the fluid density by the acceleration. The particle phase solution is step S3, which involves solving the DEM particle model. The detailed steps of solving the dynamic parameters are strictly followed in step three to calculate the real-time motion state of each sea ice particle and obtain parameters such as displacement, velocity, and angular velocity of the sea ice particles. S3: After solving at each time step, the residual coefficients of the CFD flow field are calculated in real time. Specifically, for each grid cell, the theoretical velocity of the grid cell is calculated based on the Navier-Stokes equations derived from the fluid dynamics theory formulas. The difference between the calculated velocity and the theoretical velocity is calculated to obtain the absolute difference. The absolute difference is divided by the theoretical velocity to obtain the residual ratio of the grid cell. All grid cells are traversed, and the largest residual ratio is taken as the residual coefficient of the flow field. The solution step size is adjusted according to the residual coefficients. In some other specific implementations, when the residual coefficients are all below the standard accuracy threshold of fluid numerical calculation, the current time step remains unchanged; when the residual coefficients are higher than the standard accuracy threshold but lower than ten times the standard accuracy threshold, the current time step is reduced to half of the original; when the residual coefficients are higher than ten times the standard accuracy threshold, the current time step is reduced to one-quarter of the original. After adjustment, the solution for the time step is repeated until the residual coefficients all drop below the standard accuracy threshold.

Claims

1. A numerical simulation method for sea ice motion response under wind, wave and current conditions based on CFD-DEM coupling, characterized in that, The method includes: S1. Based on the energy distribution data of wind, wave and current monitored in the target sea area, identify areas with intense energy exchange, and construct the CFD grid of the numerical water tank using the energy weighted grid densification method. Based on the sea ice morphology of broken ice and layered ice, and the sea ice particle size and layered ice thickness data obtained by a third-party engine, construct a DEM hierarchical particle system. Set the boundary dynamic loading in the numerical water tank, and fuse the constructed CFD grid and DEM hierarchical particle system to form a multi-scale heterogeneous numerical basis model of wind, wave and current-sea ice. S2. Based on the multi-scale heterogeneous numerical basis model, the mesoscale unit is divided and matched by the mesoscale bridging algorithm to obtain the mesoscale bridging unit. After integrating the Liutex vorticity mechanism to complete the flow field force transmission construction of the mesoscale bridging unit, the percolation mechanism is embedded to perform force chain and particle force analysis feedback in each mesoscale bridging unit to obtain the dynamic coupled field interaction model. S3. The obtained dynamic coupled field interaction model is coupled with particle mechanical properties and particle dynamic properties to transform the dynamic coupled field interaction model into a DEM particle behavior characterization model with multi-physical property coupling. Then, the fluid phase and force balance of the DEM particle behavior characterization model is solved and synchronized with the time step. After the solution is terminated, the numerical simulation results of sea ice motion response under the action of wind, waves and current are obtained.

2. The numerical simulation method for sea ice motion response under wind, wave and current based on CFD-DEM coupling as described in claim 1, characterized in that, S1 is based on the wind, wave and current energy distribution data of the target sea area imported by a third-party engine. By statistically analyzing the energy density values ​​of each region in the wind, wave and current energy distribution data, the region with an energy density value exceeding the average energy density value is identified as a region of intense energy exchange. The grid scale of the region with intense energy exchange is set to one-tenth of the minimum vortex scale of the fluid. Regions with energy densities below the average value are identified as far-field energy stable regions, and their grid scale is set to ten times that of the region with intense energy exchange. The grid scale is then gradually transitioned from the region with intense energy exchange to the region with far-field energy stable region using linear interpolation to complete the CFD grid construction of the numerical water tank. After constructing the CFD mesh for the numerical water tank, S1 constructs the DEM-level particle system based on sea ice grain size and layer ice thickness data obtained from a third-party engine, including the following specific settings: The spherical particle classification model is used to classify the ice fragments. Based on the particle size distribution range of sea ice fragments, three continuous intervals are selected from the particle size distribution range of sea ice fragments according to the principle of equal interval division. These intervals are defined as the first particle size class, the second particle size class, and the third particle size class, respectively. The quantity distribution ratio of the first particle size class, the second particle size class, and the third particle size class is set according to the proportion of ice fragments in the three continuous intervals of sea ice fragments, thus completing the simulation of the natural randomness of ice fragment particle size. The ice layer was modeled using a hexagonal close-packed bonded particle model. Based on standard spherical DEM particles, the ice layer was bonded together using a parallel bonding model to form an ice layer plate. The ice thickness was set according to the common monitoring range of the target sea area, including three continuous ranges: thin ice layer, medium-thickness ice layer, and thick ice layer. The thin ice layer, medium-thickness ice layer, and thick ice layer were formed by stacking corresponding numbers of spherical particles along the thickness direction.

3. The numerical simulation method for sea ice motion response under wind, wave and current based on CFD-DEM coupling as described in claim 2, characterized in that, After constructing the CFD grid and DEM hierarchical particle system, data acquisition equipment including wind speed sensors, wave monitors, and underwater current meters were deployed in the target sea area. The collected wind speed values ​​were directly used as the calculation input for wind field shear stress, and the collected wave height values ​​were converted into the velocity inlet parameters of the wave generator. The flow velocity values ​​were converted into the pressure inlet parameters of the flow field. Boundary dynamic loading was set at the inlet of the numerical water tank, including loading the waves using the velocity inlet method, the flow field using the pressure inlet method, and the wind field using the shear stress inlet method. The initial time step of the loading was uniformly set. After completing the boundary dynamic loading, the constructed CFD mesh and DEM hierarchical particle system are integrated to form a multi-scale heterogeneous numerical basis model of wind, wave and current-sea ice.

4. The numerical simulation method for sea ice motion response under wind, wave and current based on CFD-DEM coupling as described in claim 1, characterized in that, The S2 model, based on a multi-scale heterogeneous numerical basis model, constructs a one-to-one mesoscale bridging unit by linking the smallest micro-element of the CFD mesh fluid phase with the particle clusters of the DEM sea ice phase. Specifically: The number of particles in the ice fragment phase is determined according to the particle size class. The first particle size class particle cluster consists of three particles, the second particle size class particle cluster consists of four particles, and the third particle size class particle cluster consists of five particles. The spatial size of the particle cluster is consistent with the size of the smallest micro-element of the fluid. The particle clusters of the layered ice phase are a cohesive particle unit of the smallest structural unit of the corresponding layered ice. Within each mesoscale bridging unit, the smallest infinitesimal element of the fluid phase and the particle cluster of the DEM sea ice phase are in the same spatial coordinate system, and the boundary of the smallest infinitesimal element coincides with the outer boundary of the particle cluster.

5. The numerical simulation method for sea ice motion response under wind, wave and current based on CFD-DEM coupling as described in claim 4, characterized in that, The S2 method employs the Liutex vortex mechanism within each mesoscale bridging unit to identify vortices in the flow field at the ice-water and wind-water interfaces. Specifically, it calculates the velocity components of each discrete point within the smallest infinitesimal element of the fluid phase, then calculates the velocity gradient tensor of each discrete point using the central difference method. The velocity gradient tensor is decomposed into a symmetric part and an antisymmetric part. The antisymmetric part is used as the vortex tensor, and the eigenvalues ​​of the vortex tensor are extracted. The magnitude of the eigenvalues ​​of the vortex tensor is used as the vortex rotation angular velocity of the corresponding discrete point. The average vortex rotation angular velocity of all discrete points is used as the average rotation angular velocity of the vortex within the corresponding mesoscale bridging unit, and twice the average rotation angular velocity is used as the vortex rotation intensity. This completes the vortex identification and rotation intensity calculation for the mesoscale bridging unit. Next, the projected area of ​​the DEM sea ice phase particle clusters is calculated. Specifically, the three-dimensional model of the DEM sea ice phase particle clusters is projected vertically towards the direction of the flow field vortex to obtain a two-dimensional projected profile. The boundary of the projected profile is discretized into several vertices using the polygon area integration method. The vertices are arranged in a clockwise order. The area of ​​each small triangle is calculated using the coordinates of adjacent vertices. The areas of all small triangles are summed to obtain the projected area of ​​the sea ice phase particle clusters. After obtaining the projected area of ​​the sea ice phase particle cluster, the vortex rotation intensity is multiplied by the projected area of ​​the sea ice particle cluster, and then multiplied by the fluid density to complete the Liutex vortex force calculation. Based on Archimedes' law and the relative velocity of the fluid, the drag force of the fluid on the sea ice is calculated, and the lift force of the fluid on the sea ice is calculated based on the Magnus effect. The drag force, lift force and Liutex vortex force are combined according to the vector triangle law to obtain the intermediate resultant force. The intermediate resultant force is then combined with the Liutex vortex force to obtain the total fluid force, thus establishing the force transmission chain of the flow field.

6. The numerical simulation method for sea ice motion response under wind, wave and current based on CFD-DEM coupling according to claim 5, characterized in that, The S3 embeds a percolation mechanism within each mesoscale bridging unit, specifically as follows: The total number of sea ice particles in each mesoscale bridging unit is counted, the number of possible contact pairs between all particles is calculated, and the number of actual contacting particle pairs is determined by the DEM particle contact detection algorithm. The ratio of the actual contacting particle pairs to the number of possible contacting particle pairs is used as the basic contact probability. For each actual contacting particle pair, the contact force and the maximum allowable contact force are calculated to obtain the contact strength coefficient. The basic contact probability is multiplied by the average of all contact strength coefficients to obtain the percolation probability of the force chain. A path search algorithm is used to find the connected force chains within the mesoscale bridging unit. The continuous chain formed by particle contact from one end to the other within the mesoscale bridging unit is taken as the connected force chain. The length of all connected force chains is counted, the average length of all connected force chains is calculated, and the average length is divided by the side length of the mesoscale unit to obtain the connectivity coefficient of the force chain. The compressive strength of sea ice is obtained, and the maximum allowable contact force between particles is calculated based on the contact area of ​​the particles. The maximum allowable contact force is used as the force chain breakage threshold. When the actual contact force between particles in the mesoscale bridging unit exceeds the force chain breakage threshold, it is determined that the force chain has broken. The parameters of the force chain are transferred to the mesoscale bridging end of the CFD fluid phase. When the percolation probability of the force chain of ice fragments exceeds the critical percolation probability of the fluid-particle coupling, the viscosity coefficient of the fluid phase in the corresponding mesoscale unit is adjusted according to the magnitude of the force chain connectivity coefficient. When the contact force of the force chain at the bonding surface of layered ice particles exceeds the fracture threshold and percolation failure occurs, the distribution of the fractured layered ice particles is transferred to the CFD fluid phase, and the volume fraction of layered ice particles in the unit is statistically analyzed. The multiscale heterogeneous numerical basis model is transformed into a dynamic coupled field interaction model.

7. The numerical simulation method for sea ice motion response under wind, wave and current based on CFD-DEM coupling as described in claim 1, characterized in that, The S3, based on a dynamic coupled field interaction model, characterizes the mechanical properties of ice fragments and layered ice particles, specifically as follows: The HertzMindlin non-adhesive contact model is adopted. During the relative sliding of ice fragments, the product of the normal contact force between the particles and the coefficient of friction is used as the sliding friction force. The velocity decay of the ice fragments after the collision is calculated based on the sliding friction force. The HertzMindlin parallel bonding model was used to calibrate the bonding strength and bonding stiffness of the layered ice particles. The tensile stress was obtained by dividing the tensile force on the bonding surface of the layered ice by the area of ​​the bonding surface. The tensile stress was compared with the bonding strength. When the actual tensile stress exceeded the bonding strength, the bonding of the layered ice particles was triggered to break. The broken layered ice particles were automatically switched to the HertzMindlin non-bonded contact model of the crushed ice particles. The friction coefficient and the coefficient of restitution were adopted from the parameters of the crushed ice particles.

8. The numerical simulation method for sea ice motion response under wind, wave and current based on CFD-DEM coupling as described in claim 7, characterized in that, After completing the mechanical property characterization, S3 uses the wind field shear stress and total fluid force collected by S2 as the external excitation load of the DEM particles, loads them onto the centroid of the particles, and distributes the external excitation load to the particle cluster composed of multiple particles according to the mass ratio. The initial motion state of the particles is determined, all particles are traversed, and the external loads on each particle, including the total fluid force and wind shear force, as well as the inter-particle interaction forces, including the normal contact force, sliding friction force, and adhesive force, are calculated. All forces are accumulated in the vector direction to obtain the resultant force of the particles. The resultant force is divided by the mass of the particle to obtain the translational acceleration of the particle. The translational acceleration is integrated over time, with the integration interval being the current time step, to obtain the velocity change. The initial velocity plus the velocity change gives the current velocity of the particle. The current velocity is integrated over time to obtain the displacement change. The initial displacement plus the displacement change gives the current position of the particle. All the torques acting on the particle are calculated and summed to obtain the total torque. Finally, the particle dynamics characteristics including displacement, velocity, and torque are obtained. The particle dynamics and mechanical properties are assigned to the dynamic coupled field interaction model, which is then transformed into a DEM particle behavior characterization model with multi-physical property coupling, thus completing the particle phase solution.

9. The numerical simulation method for sea ice motion response under wind, wave and current based on CFD-DEM coupling as described in claim 1, characterized in that, S3 uses the finite volume method to discretize the CFD mesh constructed by S1, decomposing the fluid continuity equation and momentum conservation equation into mass conservation and force balance equations for each mesh element, and then performing mass conservation and force balance solutions, specifically: The inflow and outflow masses of each grid cell within the current time step are calculated by multiplying the flow velocity of each face of the grid cell by the area and the time step. The inflow mass minus the outflow mass is the mass change within the grid cell, thus completing the mass conservation solution. Based on the pressure difference between adjacent grid cells, the pressure, viscous force, and inertial force on each grid cell are calculated. By solving the force balance equations when the vector sum of all forces is zero, the velocity and pressure of the grid cells are obtained. By iterating through all grid cells, the flow field is solved once, and the numerical simulation parameters of the flow field, including velocity, pressure, and vortex, are obtained, thus completing the force balance solution.

10. The numerical simulation method for sea ice motion response under wind, wave and current based on CFD-DEM coupling according to claim 9, characterized in that, S3: After the solution is completed at each time step, the residual coefficient of the CFD flow field is calculated in real time. Specifically, for each grid cell, the theoretical velocity of the grid cell is calculated according to the fluid dynamics theory formula. The difference between the calculated velocity and the theoretical velocity is calculated to obtain the absolute difference. The absolute difference is divided by the theoretical velocity to obtain the residual ratio of the grid cell. All grid cells are traversed, and the largest residual ratio is taken as the residual coefficient of the flow field. The solution step size is adjusted according to the residual coefficient.

Images