Method and device for calculating ice load under action of ship and ice ridge based on discrete element

Abstract

The invention provides an ice load calculation method and device under the action of a ship and an ice ridge based on discrete elements, and relates to the technical field of ice load calculation methods, and the method comprises the steps: carrying out the fluid particle swarm division of a seawater monitoring region through employing a cubic array method, calculating the acceleration of each fluid particle, and updating the speed of each fluid particle; calculating the total dragging force of all adjacent fluid particles on the ice particles in the preset search area; and calculating the hot spot stress of each grid unit, and calculating the ship body fatigue damage coefficient of the grid unit to be detected. Through real-time speed exchange of the fluid particles and the ice particles, the dragging force of the fluid to the ice is accurately calculated, and the motion state of the ice particles is reversely updated; a discrete element ice particle model generated based on CT scanning and a fractal algorithm can truly reflect a heterogeneous microstructure of an ice ridge and mechanical properties related to temperature and salinity.

Description

Technical Field

[0001] This invention relates to the field of ice load calculation methods, specifically to a method and apparatus for calculating ice loads under the action of ships and ice ridges based on discrete element method. Background Technology

[0002] In the analysis of the interaction between polar vessels and sea ice ridges, traditional methods often treat the ice as a homogeneous continuous medium or use simplified mechanical models, making it difficult to accurately characterize the complex, heterogeneous, porous, and discretized structure and mechanical properties inside the ice ridge. This results in insufficient accuracy in predicting the spatiotemporal distribution of ice loads on the ship structure, and in particular, it is impossible to effectively simulate the motion of ice particles under hydrodynamic forces and the evolution of local high-frequency impacts and fatigue damage to the hull.

[0003] In the prior art, document CN117556530A proposes a method for rapid modeling of ice ridges and calculation of ice loads under dynamic interaction between the ship hull and the ice ridge using discrete element method (DEM) on the CUDA platform. However, this method lacks a discrete element ice particle model generated by CT scanning and fractal algorithms, failing to accurately reflect the heterogeneous microstructure and temperature- and salinity-related mechanical properties of the ice ridge, thus hindering accurate calculation of the contact force between ice particles and the ship hull. Furthermore, it does not utilize bidirectional coupling of SPH-DEM to dynamically capture the drag effect of water flow on ice particles and the feedback of ice particle motion on the flow field, thus failing to accurately simulate the time-varying nature and spatial distribution of ice loads. Finally, it does not incorporate updated local loads for hotspot-based fatigue damage assessment and early warning, significantly improving the accuracy of predicting the local strength and fatigue life of ship structures in ice-covered areas. Therefore, there is an urgent need for a method and apparatus for calculating ice loads under the action of ships and ice ridges based on discrete element method.

[0004] The information disclosed in the background section is only intended to enhance the understanding of the background of this disclosure, and therefore may include information that does not constitute prior art known to those skilled in the art. Summary of the Invention

[0005] The purpose of this invention is to provide a method and apparatus for calculating ice loads under the action of ships and ice ridges based on discrete element method, so as to solve the problems mentioned in the background art.

[0006] To achieve the above objectives, the present invention provides the following technical solution: The method for calculating ice loads under the interaction of ships and ice ridges based on the discrete element method includes the following steps: S1: Construct a 3D model of the ship to be tested and establish a global coordinate system. Delineate the seawater monitoring area under the global coordinate system. Divide the ship to be tested into a finite element mesh and mark the mesh in contact with the ice ridge. Obtain the CT scan results of the ice ridge in contact with the ship and use a fractal algorithm to generate discrete element ice particles consistent with the geometric boundary. S2: Collect temperature, salinity, and porosity parameters for each ice particle according to the grid position and calculate the elastic modulus parameter of each ice particle; establish a dynamic stiffness model based on the elastic modulus parameter and calculate the normal force of each ice particle using the dynamic stiffness model; then calculate the tangential force of each ice particle based on the normal force of each ice particle using the Coulomb friction method. S3: Count the number of contact points between the hull mesh cells to be tested and the ice particles, and calculate the normal load and tangential friction force of the mesh cells to be tested respectively; divide the seawater monitoring area into fluid particle groups, calculate the acceleration of each fluid particle, and update the velocity of each fluid particle in the seawater monitoring area; based on the velocity of each fluid particle at the corresponding ice particle and the velocity of the ice particle, calculate the total drag force of all neighboring fluid particles on the ice particle in the preset search area. S4: Based on the drag force of fluid particles on ice particles, calculate the acceleration of ice particles and update the velocity of ice particles. Based on the updated velocity of ice particles, recalculate the normal force and tangential friction force of ice particles. Simultaneously update the normal load and tangential friction force of the mesh element under test, and complete the update of the hull and ice ridge load status at the current time step. S5: For each grid cell to be tested, the stress components are solved based on the updated normal load and tangential friction force. The hot spot stress is calculated and corrected to obtain the equivalent alternating stress amplitude. The fatigue damage coefficient is calculated based on the equivalent alternating stress amplitude and the number of contacts between the grid cell to be tested and the ice particles. The fatigue damage coefficient is classified according to the preset threshold.

[0007] The ice particles that generate discrete elements are as follows: A three-dimensional model of the ship to be tested is constructed and a global coordinate system is established. The positive direction of the Z-axis of the global coordinate system is defined as vertically upward. The X-axis is set to point towards the bow and stern of the ship, and the Y-axis points towards the port and starboard sides of the ship. The seawater monitoring area is delineated under the global coordinate system. The monitoring area extends forward from the foremost point of the bow by 1.0 to 1.5 times the ship length, to the stern by 0.5 to 1.0 times the ship length, from the bottom of the ship downward by 2.5 to 3.0 times the ship's draft, to the deck above the deck by 0.5 to 1.0 times the ship's width, and includes the total width of the port and starboard sides of the ship, which is 2.0 to 3.0 times the ship's width. The hull under test was meshed using finite element methods, and the meshes in contact with the ice ridge were marked. Specifically, the outer surface geometric data were extracted based on the hull CAD model, and meshes were generated for the bow, weld seams, deck, and stern. A denser mesh was used at the bow and weld seams, while a sparser mesh was used at the deck and stern. All meshes were then filtered and marked to identify those meshes in contact with the ice ridge. The CT scan results of the ice ridges in contact with the ship's hull are obtained, and discrete element ice particles are generated using a fractal algorithm. Specifically, the CT scan results are denoised, the porosity distribution characteristic parameters of the ice ridge CT scan results are extracted, and multiple fractal surfaces are generated in three-dimensional space using the random midpoint displacement method. These surfaces divide the space into ice matrix regions and porosity regions. The Voronoi cutting algorithm is used to generate polyhedral discrete element ice particles for the ice matrix regions.

[0008] Furthermore, the normal and tangential forces of the ice particles are calculated separately. The specific steps are as follows: The temperature, salinity, and porosity parameters of ice particles at each grid point are obtained, and the elastic modulus parameter of each ice particle is calculated. Specifically, based on 10.2, an exponential function value with the natural constant e as the base is multiplied. The exponent of the exponential function value is calculated by adding 273 to the temperature value of the ice particle, dividing by 120, and taking the negative value to obtain the first calculation term. The first calculation term is then multiplied by a salinity correction factor obtained by subtracting 0.15 from 1 and multiplying by the square root of salinity to obtain the second calculation term. The first and second calculation terms are multiplied together, and the product is then multiplied by a porosity correction factor obtained by subtracting the cube of the porosity from 1, thereby obtaining the elastic modulus of the ice particles in that grid cell. Based on the elastic modulus parameter, a dynamic stiffness model is constructed to calculate the normal force of each ice particle in the grid cell under test. Specifically, the value of the normal force is obtained by multiplying a calibration coefficient by the elastic modulus of the ice particle, then by the radius of the ice particle, and finally by the modulus of the normal velocity component of the ice particle. The normal velocity component is obtained by performing a dot product operation between the ice particle velocity vector and the contact point normal vector. Furthermore, based on the normal force of each ice particle, the tangential force of the ice particles in the grid cell under test is calculated using the Coulomb friction method. Specifically, the tangential velocity vector is obtained by subtracting the normal velocity component from the velocity vector of the ice particle; the tangential velocity vector is then divided by the sum of its magnitude and a minimal near-zero constant to obtain the unit tangential vector; the tangential force is finally determined by multiplying the product of the friction coefficient and the normal force by the unit tangential vector, and then by the hyperbolic tangent function value of the ratio of the tangential velocity magnitude to the ideal reference velocity.

[0009] Furthermore, the normal load and tangential friction force of the mesh element under test are calculated using the force region average distribution method. The specific steps are as follows: The number of contact points between the tested grid cells and ice particles on the hull was counted. The normal load and tangential friction force of the tested grid cells were calculated using the force region average distribution method. Specifically, the normal load was calculated by summing the normal forces of the ice particles corresponding to all contact points on the grid cell and then dividing the sum by the surface area of ​​the grid cell. The tangential friction force was calculated by summing the tangential forces of the ice particles corresponding to all contact points on the grid cell and then dividing the sum by the surface area of ​​the grid cell.

[0010] Furthermore, the velocity of each fluid particle is updated using the Leapfrog integral method, specifically through the following steps: The cubic array method is used to divide the seawater monitoring area into fluid particle swarms. Specifically, a three-dimensional monitoring cubic domain is defined according to the spatial range of the seawater monitoring area. It is regarded as a regular hexahedral computation space. Within the cubic domain, the entire water body is equally cut into a series of regularly arranged tiny cubic grids according to a preset uniform spatial interval. The vertex position of each cubic grid is used as the spatial coordinate of an initial fluid particle. The SPH method is used to calculate the acceleration of each fluid particle within a preset neighboring fluid particle region. Specifically, it involves: first, identifying all neighboring fluid particles within the preset neighboring region except for itself; for each neighboring particle, multiplying its mass by the sum of the quotient of the current particle's pressure and the square of its own density, and the quotient of the neighboring particle's pressure and the square of its own density, then multiplying by the gradient value of the kernel function between these two particles, and summing the calculation results for all neighboring particles within the region. This negative summation, plus the gravitational acceleration, is the acceleration of the fluid particle at the current time step. The pressure of the current particle is calculated using a state equation: the state equation is based on the ratio of particle density to reference density, where the particle density is obtained by dividing the sum of the weighted contributions of neighboring particles to its kernel function by the reference density, and then raising the calculated ratio to a specified power and subtracting one, finally multiplying by the stiffness coefficient of seawater. The velocity of each fluid particle is updated using the Leapfrog integral method. Specifically, the velocity value of each fluid particle in the previous half-time step is added to the product of the acceleration value of the fluid particle in the current integer time step and the preset integration time step, so as to obtain the velocity value of the fluid particle in the new half-time step.

[0011] Furthermore, the drag force of the fluid particles on the ice particles is calculated, and the specific steps are as follows: Based on the velocity of each fluid particle at its corresponding ice particle and the velocity of the ice particle itself, the SPH-DEM bidirectional coupling method is used to calculate the total drag force exerted by all neighboring fluid particles on the ice particle within a preset search area. Specifically, for each neighboring fluid particle within the preset search area, the vector difference between the fluid particle's velocity and the ice particle's velocity, as well as its modulus, are first calculated. Then, half of this difference is multiplied by the drag coefficient, then by the seawater density, then by the projected area of ​​the ice particle, then by the modulus, and finally by the vector difference to obtain the drag force contribution of that single fluid particle on the ice particle. Finally, the drag force contributions of all neighboring fluid particles within the search area are vector-summed to obtain the total drag force exerted by all fluid particles on the ice particle.

[0012] Furthermore, the normal load and tangential friction force of the measured mesh element are updated, specifically through the following steps: Based on the drag force of fluid particles on ice particles, the acceleration of ice particles is calculated using Newton's second law. Specifically, the total drag force on the ice particle is vector-added with the gravity obtained by multiplying the mass of the ice particle by the acceleration due to gravity. The total force on the ice particle is then divided by the mass of the ice particle, and the result is the acceleration of the ice particle. At the same time, the velocity of the ice particles is updated by vector addition of the current velocity of the ice particles and the calculated acceleration of the ice particles multiplied by the preset integration time step. The sum is the updated velocity of the ice particles. Based on the updated ice particle velocities, the normal force and tangential friction force of the ice particles are recalculated. Specifically, the updated ice particle velocity is subtracted from the ship's velocity, and the resulting relative velocity vector is projected onto the normal direction of the contact point between the ice particle and the ship. The magnitude and direction of the normal relative velocity are determined by calculating the cosine of the angle between the projection and the normal vector. The normal force of the ice particle consists of two parts: the first part is calculated by multiplying the calibration coefficient by the elastic modulus of the ice particle, then by the radius of the ice particle, and finally by the modulus of the normal relative velocity; the second part is calculated by multiplying the normal damping coefficient by the normal relative velocity. The updated normal force of the ice particle is obtained by adding these two parts together. The relative velocity obtained by subtracting the ship's speed from the updated ice particle's velocity is further subtracted by its projection onto the normal direction at the contact point between the ice particle and the ship, yielding the tangential velocity component. The tangential force is calculated by multiplying four parts consecutively: the first part is the product of the friction coefficient and the updated normal force; the second part is the quotient obtained by dividing the tangential velocity component by its modulus and the sum of a minimal near-zero constant; the third part is the hyperbolic tangent function value of the ratio of the tangential velocity modulus to an ideal reference velocity; the product of these three parts is then multiplied by a unit tangential vector used to determine the direction, finally yielding the updated tangential frictional force of the ice particle. Based on the normal force and tangential friction of the ice particles, the updated normal load and tangential friction of the measured mesh element are further calculated. Specifically, for each mesh element, the updated normal forces of the ice particles corresponding to all its contact points are accumulated, and the sum is divided by the surface area of ​​the mesh element to obtain the updated normal load; at the same time, the updated tangential friction of the ice particles corresponding to all its contact points are accumulated, and the sum is divided by the surface area of ​​the mesh element to obtain the updated tangential friction.

[0013] Furthermore, the equivalent alternating stress amplitude is calculated, and the specific steps are as follows: Based on the updated normal load and tangential friction of the grid cells to be tested, the hot spot stress of each grid cell is calculated using the Von Mises criterion method. Specifically, the normal load is squared, and then three times the square of the tangential friction is added. The sum is then squared, and the result is the hot spot stress value of the grid cell. Furthermore, Goodman correction is applied to the hot spot stress to obtain the equivalent alternating stress amplitude. Specifically, the hot spot stress is divided by the quotient of the normal load and the ultimate tensile strength of the marine steel, and the result is the equivalent alternating stress amplitude of the mesh element.

[0014] Furthermore, based on the threshold results, corresponding alarms are issued. The specific steps are as follows: The number of contacts between the grid cell under test and the ice particles is counted. Based on the equivalent alternating stress amplitude, the hull fatigue damage coefficient of the grid cell under test is calculated. Specifically, the total number of contacts between the grid cell under test and the ice particles is counted. Based on the equivalent alternating stress amplitude generated by each contact, the hull fatigue damage coefficient of the grid cell is calculated. In the calculation, the equivalent alternating stress amplitude corresponding to each contact of the grid cell is taken as the specified power of the slope of the SN curve. Then, the power of the results of all contact numbers is accumulated. Finally, the accumulated sum is divided by the SN curve constant of the hull material. The quotient is the fatigue damage coefficient of the grid cell. Threshold classification is applied to the hull fatigue loss coefficient, and corresponding alarms are issued based on the classification results. Threshold classification is also applied to the hull fatigue damage coefficient, and corresponding alarms are issued based on the classification results. When the calculated fatigue damage coefficient is greater than or equal to 0.3, it indicates that the damage accumulation of the grid cell has entered the accelerated stage, and the system will issue a reminder to prompt increased monitoring. When the fatigue damage coefficient is greater than or equal to 0.7, it indicates that the remaining life of the grid cell has been significantly reduced, and an immediate overhaul is required. When the fatigue damage coefficient is greater than or equal to 1.0, it indicates that the grid cell has entered a state of fatigue failure risk, and it must be immediately stopped and overhauled.

[0015] The present invention also provides a discrete element method-based ice load calculation device for ships interacting with ice ridges, the device being used to perform the above-described calculation method, comprising: The ice particle construction module is used to construct a 3D model of the ship under test and establish a global coordinate system. Under the global coordinate system, the seawater monitoring area is delineated, the ship under test is meshed using finite element methods, and the meshes in contact with the ice ridges are marked. The CT scan results of the ice ridges in contact with the ship are obtained, and the discrete element ice particles consistent with the geometric boundaries are generated using a fractal algorithm. The force calculation module is used to collect temperature, salinity, and porosity parameters for each ice particle according to the grid position and calculate the elastic modulus parameter of each ice particle; a dynamic stiffness model is established based on the elastic modulus parameter and the normal force of each ice particle is calculated using the dynamic stiffness model; then, the tangential force of each ice particle is calculated using the Coulomb friction method based on the normal force of each ice particle. The fluid force calculation module is used to count the number of contact points between the hull mesh cells to be tested and the ice particles, and to calculate the normal load and tangential friction force of the mesh cells to be tested; to divide the seawater monitoring area into fluid particle swarms, calculate the acceleration of each fluid particle, and update the velocity of each fluid particle within the seawater monitoring area; based on the velocity of each fluid particle at the corresponding ice particle and the velocity of the ice particle, to calculate the total drag force of all neighboring fluid particles on the ice particle within the preset search area. The parameter update module is used to calculate the acceleration of ice particles and update the velocity of ice particles based on the drag force of fluid particles on ice particles. Based on the updated velocity of ice particles, the normal force and tangential friction force of ice particles are recalculated, and the normal load and tangential friction force of the mesh element under test are updated synchronously to complete the update of the hull and ice ridge load status at the current time step. The classification alarm module is used to solve the stress components for each grid cell under test based on the updated normal load and tangential friction force, calculate the hot spot stress and correct it to obtain the equivalent alternating stress amplitude, calculate the fatigue damage coefficient based on the equivalent alternating stress amplitude and the number of contactes between the grid cell under test and the ice particles, and classify the fatigue damage coefficient according to the preset threshold.

[0016] Compared with the prior art, the beneficial effects of the present invention are: By integrating the discrete element method, smoothed particle hydrodynamics, and finite element analysis, high-precision simulation of the multiphysics coupling effect between the ship, ice ridge, and fluid was achieved. The discrete element ice particle model generated based on CT scanning and fractal algorithms can realistically reflect the heterogeneous microstructure of the ice ridge and its temperature and salinity-related mechanical properties, thereby accurately calculating the contact force between the ice particles and the ship hull. The drag effect of water flow on ice particles and the feedback of ice particle motion on the flow field were dynamically captured, accurately simulating the time-varying nature and spatial distribution of ice loads. Combined with updated local loads, hotspot-based fatigue damage assessment and early warning were performed, significantly improving the accuracy of predicting the local strength and fatigue life of ship structures in ice-covered areas, providing strong technical support for the safe design and real-time risk management of ships. Attached Figure Description

[0017] Figure 1 This is a schematic diagram of the overall method flow of the present invention; Figure 2 This is a graph showing the relationship between the normal load on the mesh element and the corresponding hot spot stress. Figure 3 This is a schematic diagram of the overall device of the present invention. Detailed Implementation

[0018] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments.

[0019] It should be noted that, unless otherwise defined, the technical or scientific terms used in this invention should have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in this invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" mean that the element or object preceding the word encompasses the elements or objects listed following the word and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0020] Example: Please see Figures 1-2 The present invention provides a technical solution: The method for calculating ice loads under the interaction of ships and ice ridges based on the discrete element method includes the following steps: S1: Construct a 3D model of the ship to be tested and establish a global coordinate system. Delineate the seawater monitoring area under the global coordinate system. Divide the ship to be tested into a finite element mesh and mark the mesh in contact with the ice ridge. Obtain the CT scan results of the ice ridge in contact with the ship and use a fractal algorithm to generate discrete element ice particles consistent with the geometric boundary. The ice particles that generate discrete elements are as follows: A three-dimensional model of the ship to be tested is constructed and a global coordinate system is established. The positive direction of the Z-axis of the global coordinate system is defined as vertically upward. The X-axis is set to point towards the bow and stern of the ship, and the Y-axis points towards the port and starboard sides of the ship. The seawater monitoring area is delineated under the global coordinate system. The monitoring area extends forward from the foremost point of the bow by 1.0 to 1.5 times the ship length, to the stern by 0.5 to 1.0 times the ship length, from the bottom of the ship downward by 2.5 to 3.0 times the ship's draft, to the deck above the deck by 0.5 to 1.0 times the ship's width, and includes the total width of the port and starboard sides of the ship, which is 2.0 to 3.0 times the ship's width. In the above process, the bow needs to be extended by 1.0-1.5 times the ship's length. This is because when a ship is sailing, the bow will have a pre-compression effect on the water in front of it. 1.0 times the ship's length is the basic distance at which pressure disturbances decay to a negligible level. If the ship speed is high or the bow is full, such as oil tankers or bulk carriers, the pressure propagation distance will be even greater. Therefore, it needs to be extended to 1.5 times the ship's length to ensure that the high-pressure area at the bow is completely included in the calculation domain and to avoid the bow pressure being underestimated due to truncation. The wake field needs to extend 0.5-1.0 times the ship's length behind the stern. This is to fully capture the development of the wake field behind the stern. When the water flows past the hull and converges at the stern, it will form a complex vortex structure and a velocity loss area, i.e., the wake field. The length of the wake field is related to the ship type. The wake field of slender ships such as container ships is relatively short, while the wake field of large ships can extend to a considerable distance behind the stern. 0.5 times the ship's length is the basic length of the wake core area, while 1.0 times the ship's length is to capture the complete attenuation of the wake. The required draft of 2.5-3.0 times the hull depth is derived from the theoretical analysis of the shallow water blockage effect. When the hull depth is too close to the seabed, a narrow channel forms between the hull and the seabed, accelerating the water flow and causing additional sinking force and drag on the hull—the shallow water effect. Studies show that the shallow water effect becomes significant when the ratio of water depth to draft is less than 3.0; when the ratio is less than 2.5, the impact on the hull is already non-negligible. Therefore, to simulate deep-water conditions—open water unaffected by the seabed—it is essential to ensure that the seabed boundary is sufficiently far from the hull depth to allow the water flow to develop freely below the hull without compression. A draft of 2.5-3.0 times the hull depth is precisely the safe distance to ensure that the seabed boundary has no significant impact on the flow field around the hull. The deck needs to be 0.5-1.0 times the ship's beam upwards to accommodate the upwelling height of ship-borne waves. When a ship is sailing, water accumulates at the bow, forming wave crests higher than the calm water surface. The height of these wave crests depends on the ship's speed and hull type, and can typically reach several meters. Furthermore, in wave environments, wave height must also be considered. 0.5 times the ship's beam is the upper limit for typical ship-borne wave heights, while 1.0 times the ship's beam considers extreme wave conditions. If this height is set too low, the wave crest will impact the top boundary of the computational domain, generating non-physical pressure reflections that interfere with the actual flow field at the bow.

[0021] The hull under test was meshed using finite element methods, and the meshes in contact with the ice ridge were marked. Specifically, the outer surface geometric data were extracted based on the hull CAD model, and meshes were generated for the bow, weld seams, deck, and stern. A denser mesh was used at the bow and weld seams, while a sparser mesh was used at the deck and stern. All meshes were then filtered and marked to identify those meshes in contact with the ice ridge. The CT scan results of the ice ridges in contact with the ship's hull are obtained, and discrete element ice particles are generated using a fractal algorithm. Specifically, the CT scan results are denoised, the porosity distribution characteristic parameters of the ice ridge CT scan results are extracted, and multiple fractal surfaces are generated in three-dimensional space using the random midpoint displacement method. These surfaces divide the space into ice matrix regions and porosity regions. The Voronoi cutting algorithm is used to generate polyhedral discrete element ice particles for the ice matrix regions.

[0022] In the above process, filtering techniques are applied to the original CT images to reduce noise. Non-local mean filtering or Gaussian filtering is usually used to remove noise generated during the scanning process, while preserving the edge details of the pores as much as possible. The pore distribution feature parameters of the ice ridge CT scan results are extracted as follows: the proportion of pores in the entire ice ridge volume is statistically analyzed, then each connected pore is marked, their equivalent diameter and volume are calculated, and the frequency distribution of pore size is statistically analyzed. Multiple fractal surfaces are generated in 3D space using the random midpoint displacement method. The aim is to reconstruct a 3D porous structure of ice ridges with realistic statistical characteristics based on extracted distribution feature parameters. First, a cubic mesh containing multiple grid points is established in 3D space, much like drawing many equally spaced grid lines within a cube. Then, the random midpoint displacement method is used to generate fractal Brownian motion surfaces. Specifically, initial height values ​​are assigned to the eight corner points of the cube. Then, the heights of the midpoints of each edge and the center points of each face are calculated. The height of the new point is equal to the average height of the surrounding points, plus a random offset. The mesh size is then halved, and this process is repeated until the desired mesh resolution is achieved. The magnitude of the random offset at each step is controlled by the previously calculated fractal dimension; a higher fractal dimension results in a rougher surface, and vice versa. By repeating this process, multiple intersecting and nested fractal surfaces can be generated. These surfaces collectively divide the 3D space into ice matrix regions and porous regions, where the region below the fractal surface corresponds to the porous region, and the region above corresponds to the ice matrix. The Voronoi slicing algorithm is used to generate polyhedral discrete element ice particles in the ice matrix region. The purpose is to discretize the continuous ice matrix region into independent ice particles represented by polyhedra. Within the ice matrix region, a series of seed points are randomly arranged according to the distribution characteristics of ice crystal size. The density of the seed points determines the size distribution of the final generated ice particles. Then, using these seed points as the core, Voronoi slicing is performed on the entire ice matrix region. This involves drawing a dividing surface at the midpoint of adjacent seed points to form a series of closely spaced polyhedral elements, each corresponding to a potential ice particle. Next, the fractal surface generated in the previous step is used as the cutting boundary to trim the Voronoi elements. The portion above the fractal surface is retained as ice particles, while the portion below the fractal surface is removed to form pores. This introduces pores into the ice ridge model. For each generated ice particle, the local salinity value is extracted from the original CT image based on its centroid position, and the particle temperature parameter is assigned in conjunction with the ambient temperature field. The porosity is determined by the proportion of the pore volume inside the particle.

[0023] S2: Collect temperature, salinity, and porosity parameters for each ice particle according to the grid position and calculate the elastic modulus parameter of each ice particle; establish a dynamic stiffness model based on the elastic modulus parameter and calculate the normal force of each ice particle using the dynamic stiffness model; then calculate the tangential force of each ice particle based on the normal force of each ice particle using the Coulomb friction method. The specific steps for calculating the normal and tangential forces on the ice particles are as follows: The temperature, salinity, and porosity parameters of ice particles at each grid point are obtained, and the elastic modulus parameter of each ice particle is calculated. Specifically, based on 10.2, an exponential function value with the natural constant e as the base is multiplied. The exponent of the exponential function value is calculated by adding 273 to the temperature value of the ice particle, dividing by 120, and taking the negative value to obtain the first calculation term. The first calculation term is then multiplied by a salinity correction factor obtained by subtracting 0.15 from 1 and multiplying by the square root of salinity to obtain the second calculation term. The first and second calculation terms are multiplied together, and the product is then multiplied by a porosity correction factor obtained by subtracting the cube of the porosity from 1, thereby obtaining the elastic modulus of the ice particles in that grid cell. The formula used in the above process is: in, Indicates the first The first grid cell of the first grid cell Elastic modulus parameter of an ice particle; Indicates the first The first grid cell of the first grid cell The temperature of each ice particle; Indicates the first The first grid cell of the first grid cell The salinity of each ice particle; Indicates the first The first grid cell of the first grid cell Porosity of each ice particle; Indicates the grid cell number index; Indicates the index number of ice particles; In the above process, the dependent variable of the formula is the elastic modulus parameter of the ice particles. Specifically, it reflects the stiffness characteristics of ice materials under a given microscopic physical state. The technical advantage lies in its ability to quantify the heterogeneous physical properties of ice into mechanical parameters that can participate in discrete element method (DEM) calculations, providing accurate and physically meaningful material constitutive inputs for subsequent dynamic calculations of ice-ship contact forces. The independent variable of this formula is temperature. ,salinity With porosity These factors collectively determine the microstructure and composition of natural sea ice, and are key physical variables affecting its macroscopic mechanical properties: temperature directly influences the elastic modulus by affecting the thermal vibration and phase state of the ice lattice; salinity determines the brine content and distribution within the ice, affecting its internal defects and strength; porosity characterizes the proportion of void volume within the ice, directly related to its effective load-bearing area; and the elastic modulus exhibits an exponential decay with increasing temperature. The increase in salinity is reflected in a linear correction term. And the increase in porosity is reflected in the cubic attenuation term. The decrease is consistent with the actual physical law that the mechanical strength of ice materials decreases when they are warmer, saltier, and more porous. The constant 10.2 is the reference value for the elastic modulus of ice under ideal low-temperature, salt-free, and non-porous conditions, providing the basis for dimensional matching and numerical normalization; the value 273 is used to convert Celsius temperature to absolute temperature to ensure the physical continuity of the exponential function; the value 120 is a characteristic parameter of temperature sensitivity, reflecting the rate at which the elastic modulus of ice decays exponentially with increasing temperature, and its value is based on the fitting of experimental data on the thermal vibration energy levels of the ice lattice; the coefficient 0.15 characterizes the degree of attenuation of the elastic modulus by salinity, which is correlated with the volume effect of brine inclusions within the ice through the square root of salinity, and its magnitude is determined by empirical statistics on the salt precipitation pattern during seawater freezing; the exponent 3 reflects the significant nonlinear effect of porosity on the elastic modulus, and the cubic relationship is based on the effective field theory of porous media mechanics, reflecting the strong attenuation effect of pores as stress concentration sources on the overall stiffness of the material.

[0024] Based on the elastic modulus parameter, a dynamic stiffness model is constructed to calculate the normal force of each ice particle in the grid cell under test. Specifically, the value of the normal force is obtained by multiplying a calibration coefficient by the elastic modulus of the ice particle, then by the radius of the ice particle, and finally by the modulus of the normal velocity component of the ice particle. The normal velocity component is obtained by performing a dot product operation between the ice particle velocity vector and the contact point normal vector. The formula used in the above process is: in, in, Indicates the radius of the ice particle; Indicates the first The first grid cell of the first grid cell The normal force of an ice particle; Indicates the first The first grid cell of the first grid cell The normal vector of the contact point of each ice particle; Indicates the first The first grid cell of the first grid cell The instantaneous velocity of an ice particle; Indicates the first The first grid cell of the first grid cell The normal relative velocity of an ice particle; Indicates the calibration coefficient; In the above process, the dependent variable in the formula is the normal force when the ice particles come into contact with the ship's hull. The physical meaning of this formula is the instantaneous force exerted by ice particles on the hull along the normal direction of the contact surface during collision. The technical advantage lies in realizing the dynamic quantification of contact force based on the material properties, geometric characteristics, and motion state of the ice particles, providing accurate mechanical input for subsequent load distribution and fatigue damage analysis of the mesh elements; the independent variable in the formula includes the elastic modulus of the ice particles. ,radius and normal relative velocity These factors determine the mechanical response of a contact collision from three dimensions: material stiffness, geometric dimensions, and motion state. The elastic modulus reflects the ice's resistance to deformation and directly affects the contact stiffness. The radius determines the characteristic scale of the contact area. The normal velocity characterizes the intensity of the collision. The normal force increases monotonically with increasing elastic modulus, radius, or normal velocity. This conforms to the physical law that the contact force in an elastic collision is proportional to the material stiffness, contact size, and approach speed. The calibration coefficient... Used to match the consistency between discrete element models and macroscopic mechanical responses.

[0025] Furthermore, based on the normal force of each ice particle, the tangential force of the ice particles in the grid cell under test is calculated using the Coulomb friction method. Specifically, the tangential velocity vector is obtained by subtracting the normal velocity component from the velocity vector of the ice particle; the tangential velocity vector is then divided by the sum of its magnitude and a minimal near-zero constant to obtain the unit tangential vector; the tangential force is finally determined by multiplying the product of the friction coefficient and the normal force by the unit tangential vector, and then by the hyperbolic tangent function value of the ratio of the tangential velocity magnitude to the ideal reference velocity.

[0026] The formula used in the above process is: in, in, Indicates the first The first grid cell of the first grid cell The tangential force of an ice particle; Indicates the zero constant; Indicates the first The first grid cell of the first grid cell The tangential velocity component of an ice particle; Indicates the first The first grid cell of the first grid cell The unit tangent vector of an ice particle; This represents the coefficient of friction of the ice particles; This represents the ideal velocity of the ice particles; This indicates the calculation of the modulus.

[0027] In the above process, the dependent variable in the formula is the tangential force when the ice particles come into contact with the ship's hull. Physically, this refers to the frictional force exerted on the ship's hull along the tangential direction of the contact surface by ice particles during relative sliding. The technical advantage lies in the introduction of a hyperbolic tangent function to smooth the transition, achieving a continuous and stable simulation of the transition from static to dynamic friction. This effectively avoids discontinuous jumps in numerical calculations and improves the numerical stability and physical realism of the ice-ship dynamic interaction simulation. The independent variable in the formula is the normal force. coefficient of friction Unit tangent vector and tangential velocity modulus These factors together determine the direction and magnitude of the tangential frictional force: the normal force provides the basis for its magnitude through Coulomb's law of friction; the coefficient of friction reflects the roughness of the contact surface; the unit tangential vector determines the spatial direction of the force; and the tangential velocity is determined by... The function modulates the smooth variation of frictional force in both low-speed and high-speed regions; the tangential force increases with increasing normal force or friction coefficient. Simultaneously, through... The term demonstrates a nonlinear dependence on tangential velocity, increasing approximately linearly at low speeds and tending to saturate at high speeds, consistent with the actual physical behavior of friction, and avoiding zero constant. This ensures the computational stability of the unit tangent vector at zero tangential velocity.

[0028] S3: Count the number of contact points between the hull mesh cells to be tested and the ice particles, and calculate the normal load and tangential friction force of the mesh cells to be tested respectively; divide the seawater monitoring area into fluid particle groups, calculate the acceleration of each fluid particle, and update the velocity of each fluid particle in the seawater monitoring area; based on the velocity of each fluid particle at the corresponding ice particle and the velocity of the ice particle, calculate the total drag force of all neighboring fluid particles on the ice particle in the preset search area. The normal load and tangential friction force of the mesh element under test are calculated using the force region average distribution method. The specific steps are as follows: The number of contact points between the tested grid cells and ice particles on the hull was counted. The normal load and tangential friction force of the tested grid cells were calculated using the force region average distribution method. Specifically, the normal load was calculated by summing the normal forces of the ice particles corresponding to all contact points on the grid cell and then dividing the sum by the surface area of ​​the grid cell. The tangential friction force was calculated by summing the tangential forces of the ice particles corresponding to all contact points on the grid cell and then dividing the sum by the surface area of ​​the grid cell.

[0029] The formula upon which the above process is based is: in, Indicates the first The number of contact points between each grid cell and the ice particle; Indicates the first Area of ​​each grid cell; Indicates the first Normal load of each grid cell; Indicates the first Tangential friction force of each grid cell.

[0030] In the above process, the contact forces of a large number of discontinuous and discrete ice particles at the microscale are spatially integrated and averaged according to the number of contact points and the area of ​​the grid cells by the "force region average distribution method". This generates a standard mechanical load input that can be used for structural stress and fatigue analysis at the macroscale. This not only reflects the non-uniform distribution characteristics of ice load on the hull surface, but also overcomes the computational complexity and numerical instability caused by directly processing a large number of discrete contact points.

[0031] The velocity of each fluid particle is updated using the Leapfrog integral method. The specific steps are as follows: The cubic array method is used to divide the seawater monitoring area into fluid particle swarms. Specifically, a three-dimensional monitoring cubic domain is defined according to the spatial range of the seawater monitoring area. It is regarded as a regular hexahedral computation space. Within the cubic domain, the entire water body is equally cut into a series of regularly arranged tiny cubic grids according to a preset uniform spatial interval. The vertex position of each cubic grid is used as the spatial coordinate of an initial fluid particle. The SPH method is used to calculate the acceleration of each fluid particle within a preset neighboring fluid particle region. Specifically, it involves: first, identifying all neighboring fluid particles within the preset neighboring region except for itself; for each neighboring particle, multiplying its mass by the sum of the quotient of the current particle's pressure and the square of its own density, and the quotient of the neighboring particle's pressure and the square of its own density, then multiplying by the gradient value of the kernel function between these two particles, and summing the calculation results for all neighboring particles within the region. This negative summation, plus the gravitational acceleration, is the acceleration of the fluid particle at the current time step. The pressure of the current particle is calculated using a state equation: the state equation is based on the ratio of particle density to reference density, where the particle density is obtained by dividing the sum of the weighted contributions of neighboring particles to its kernel function by the reference density, and then raising the calculated ratio to a specified power and subtracting one, finally multiplying by the stiffness coefficient of seawater. The formula upon which the above process is based is: in, in, Indicates the first The acceleration of a fluid particle at the current time step; Indicates the adjacent first The mass of a fluid particle; Indicates the first The pressure of a fluid particle; Indicates the adjacent first The pressure of a fluid particle; Indicates the first step of the current time step The fluid particle and its neighboring... The kernel function gradient value of each fluid particle; Represents gravitational acceleration; Indicates the time step index; Indicates the first A preset fluid particle region for each fluid particle; Indicates the reference density of seawater; Indicates the stiffness coefficient of seawater; Indicates Centered on, with radius as The area; Indicates the first The fluid particle and its neighboring... The Euclidean distance of a fluid particle; Represents the kernel function weights; Indicates the first The position vector of each fluid particle; Indicates the first The position vector of each fluid particle; In the above process, an accurate numerical simulation of seawater flow field motion based on the smooth particle fluid dynamics method was achieved. The core technical effect lies in the transformation of the continuous medium equations describing fluid motion, such as the Navier-Stokes equations, into mechanical equations on a discrete particle system through the particle approximation method of kernel function weighted summation. This allows for the meshless calculation of the acceleration of each fluid particle. The state equation dynamically correlates the local density of the particles, obtained by weighting the contributions of neighboring particles through kernel functions, with pressure, fully incorporating the compressibility effect of the fluid. The acceleration calculation term accurately simulates complex flow states such as pressure propagation and turbulent diffusion in seawater by including the particle interaction forces of pressure gradient and gravity, providing a high-fidelity time-varying flow field dynamics foundation for subsequent SPH-DEM coupled calculations. By constructing local particle densities using kernel function weighting and substituting them into the equation of state to calculate pressure, the density-pressure constitutive relationship of the macroscopic continuous medium is essentially discretized to the particle system, thus characterizing the compressibility of the fluid from a Lagrange perspective. The acceleration calculation formula is based on the particle method discretization of the pressure gradient and gravity terms in the Navier-Stokes equations. The kernel function gradient is used to establish an interaction force model between neighboring particles, ensuring momentum conservation and accurately simulating key fluid dynamic phenomena such as pressure propagation and vortex generation. This calculation method overcomes the mesh distortion problem of traditional mesh methods when dealing with free surfaces, large deformations, and fluid-structure interaction interfaces.

[0032] The velocity of each fluid particle is updated using the Leapfrog integral method. Specifically, the velocity value of each fluid particle in the previous half-time step is added to the product of the acceleration value of the fluid particle in the current integer time step and the preset integration time step, so as to obtain the velocity value of the fluid particle in the new half-time step.

[0033] The formula used in the above process is: in, This indicates the preset integration time step; Indicates the first The velocity of a fluid particle in the first half of the time step; Indicates the first The velocity of the fluid particles in the new half-time step; Indicates the first A fluid particle at the current integer time step The acceleration.

[0034] In the above process, the Leapfrog integral method cleverly avoids the numerical damping or energy drift problems that may exist in the direct Euler method by interleaving the calculation of velocity and position in half time steps. Its core is to use the acceleration of the current integer time step to update the velocity in the half time step. This "velocity lags behind position" update strategy can effectively maintain the long-term stability of the system's kinetic and potential energy. It is particularly suitable for simulating long-term, large-scale fluid dynamic processes such as seawater, ensuring the numerical robustness and physical realism of SPH fluid simulation in time progression. By using the Leapfrog integral scheme, the acceleration-velocity update and position update are staggered by half a step in time. This effectively suppresses the numerical energy dissipation or spurious growth problems common in the traditional explicit Euler method while ensuring second-order time accuracy, thus ensuring the stability of the system's kinetic energy during long-term simulations.

[0035] The specific steps for calculating the drag force of fluid particles on ice particles are as follows: Based on the velocity of each fluid particle at its corresponding ice particle and the velocity of the ice particle itself, the SPH-DEM bidirectional coupling method is used to calculate the total drag force exerted by all neighboring fluid particles on the ice particle within a preset search area. Specifically, for each neighboring fluid particle within the preset search area, the vector difference between the fluid particle's velocity and the ice particle's velocity, as well as its modulus, are first calculated. Then, half of this difference is multiplied by the drag coefficient, then by the seawater density, then by the projected area of ​​the ice particle, then by the modulus, and finally by the vector difference to obtain the drag force contribution of that single fluid particle on the ice particle. Finally, the drag force contributions of all neighboring fluid particles within the search area are vector-summed to obtain the total drag force exerted by all fluid particles on the ice particle.

[0036] The formula used in the above process is: in, This represents the drag coefficient, determined based on expert scoring. Represents the area of ​​the ice particles; Indicates the density of seawater; Indicates the velocity of the ice particles; Indicates the first The velocity of each fluid particle; Indicates the first The first grid cell of the first grid cell The total drag force of all neighboring fluid particles on the ice particle within the preset search area; This indicates all neighboring fluid particle regions within the preset search area.

[0037] In the above process, by dynamically linking the velocity field of SPH fluid particles with the motion state of DEM ice particles, and by locally vector superimposing the contribution of each neighboring fluid particle according to Newton's law of drag, the distributed drag force exerted by the seawater flow field on irregular ice particles is accurately calculated in a meshless particle framework. This calculation reflects the nonlinear characteristic that the drag force is proportional to the square of the relative velocity, and by summing all fluid particles in the search area, it naturally considers the influence of the non-uniformity of the flow field on the force on ice particles. This provides key hydrodynamic input for updating the acceleration and trajectory of ice particles in real ocean currents, thus ensuring the physical reality and numerical stability of ice-water coupled motion. The dependent variable in the formula is the total drag force of the fluid on the ice particles. The physical meaning of this formula is the resultant hydrodynamic drag force experienced by ice particles in flowing seawater. The technical advantage lies in dynamically converting the velocity field of continuous fluid particles into a distributed drag force acting on discrete ice particles through the SPH-DEM bidirectional coupling mechanism, thereby achieving a refined quantification of the driving effect of water flow on ice particle motion. The independent variable of this formula includes the drag coefficient. Seawater density Projected area of ​​ice particles and fluid-ice relative velocity These factors together determine the magnitude and direction of the drag force: the drag coefficient characterizes the influence of particle shape and surface roughness on resistance; seawater density and projected area together constitute the reference scale of resistance; the relative velocity term is based on the law in fluid mechanics that resistance is proportional to the square of velocity, and at the same time, it determines the direction of the force through a vector form; the total drag force increases monotonically with the increase of the drag coefficient, the increase of seawater density, the expansion of the projected area of ​​ice particles, or the increase of fluid-ice relative velocity. The formula, through the quadratic term of the relative velocity vector, ensures that the direction of the drag force is always consistent with the direction of the relative velocity, which conforms to the physical essence of Newton's law of drag.

[0038] S4: Based on the drag force of fluid particles on ice particles, calculate the acceleration of ice particles and update the velocity of ice particles. Based on the updated velocity of ice particles, recalculate the normal force and tangential friction force of ice particles. Simultaneously update the normal load and tangential friction force of the mesh element under test, and complete the update of the hull and ice ridge load status at the current time step. The specific steps for updating the normal load and tangential friction force of the mesh element under test are as follows: Based on the drag force of fluid particles on ice particles, the acceleration of ice particles is calculated using Newton's second law. Specifically, the total drag force on the ice particle is vector-added with the gravity obtained by multiplying the mass of the ice particle by the acceleration due to gravity. The total force on the ice particle is then divided by the mass of the ice particle, and the result is the acceleration of the ice particle. The formula used in the above process is: in, Indicates the first The first grid cell of the first grid cell The acceleration of an ice particle; Indicates the first The first grid cell of the first grid cell The mass of each ice particle; Indicates the first The first grid cell of the first grid cell The total force on each ice particle; In the above process, the vector sum of fluid drag force and gravity is converted into the instantaneous acceleration of ice particles according to Newton's second law, thereby accurately driving the motion evolution of each ice particle within the discrete element framework. This calculation ensures the physical self-consistency of the dynamic response of ice particles, providing accurate and real-time kinematic input for subsequent velocity updates, position movements, and a new round of contact force calculations. It is a key dynamic bridge connecting the hydrodynamic environment and the motion of the ice particles themselves, thus influencing the dynamic changes of ice-ship contact load. The dependent variable in the formula is the updated ice particle acceleration. The physical meaning of this formula is the rate of change of the instantaneous motion state of an ice particle under the combined action of fluid drag and gravity. The technical advantage lies in directly converting force into a motion response based on Newton's second law, providing a core dynamic driver for the dynamic updating of the ice particle's velocity and position, and ensuring the physical correctness and numerical stability of the ice particle motion simulation. The independent variable of this formula includes the total force acting on the ice particle. and its quality The two are directly related: the total force is the vector sum of the fluid drag force and gravity, which determines the net external force applied to the particle; the mass characterizes the particle's inertia and affects its response to external forces. Acceleration is quantified by the quotient of the total force and the mass. This formula reflects the positive correlation between the dependent variable and the total force, i.e., the acceleration increases with the increase of the total force, and the negative correlation with the mass, i.e., the acceleration decreases with the increase of the mass. It strictly follows the basic mechanical laws of Newton's second law, ensuring that the dynamic update of ice particles is physically self-consistent.

[0039] At the same time, the velocity of the ice particles is updated by vector addition of the current velocity of the ice particles and the calculated acceleration of the ice particles multiplied by the preset integration time step. The sum is the updated velocity of the ice particles. The formula upon which the above process is based is: in, Indicates the updated number The first grid cell of the first grid cell The speed of each ice particle; In the above process, by multiplying the newly calculated acceleration by the preset time step and superimposing it with the current velocity vector, the momentum accumulation and state evolution of the ice particle within a time step are completed. This update mechanism provides the ice particle with a continuous motion trajectory that conforms to the laws of Newtonian mechanics while ensuring computational efficiency. It ensures the self-consistency and synchronization of its position, velocity and force state in the time series. It is a key time-series advancement link connecting the dynamic response of the ice particle with subsequent contact collisions, load transfer and spatial position updates.

[0040] Based on the updated ice particle velocities, the normal force and tangential friction force of the ice particles are recalculated. Specifically, the updated ice particle velocity is subtracted from the ship's velocity, and the resulting relative velocity vector is projected onto the normal direction of the contact point between the ice particle and the ship. The magnitude and direction of the normal relative velocity are determined by calculating the cosine of the angle between the projection and the normal vector. The normal force of the ice particle consists of two parts: the first part is calculated by multiplying the calibration coefficient by the elastic modulus of the ice particle, then by the radius of the ice particle, and finally by the modulus of the normal relative velocity; the second part is calculated by multiplying the normal damping coefficient by the normal relative velocity. The updated normal force of the ice particle is obtained by adding these two parts together. The formula used in the above process is: in, Indicates ship speed; This indicates the angle between the velocity direction of the ice particles and the normal direction of the ship's hull; Indicates the first The first grid cell of the first grid cell The relative velocity of an ice particle to the ship's hull in the normal direction; Indicates the normal damping coefficient; Indicates the updated number The first grid cell of the first grid cell The normal force of an ice particle; In the above process, by introducing a linear damping term related to the normal relative velocity, the physical characterization of the energy dissipation effect during the collision is added to the original stiffness force model based on elastic modulus, particle size and collision velocity. This more realistically simulates the inelastic behaviors that may occur during ice-ship contact, such as plastic deformation, internal friction and local fragmentation, and improves the physical completeness and numerical convergence of the contact force calculation. The dependent variable in this formula is the updated ice particle normal force. The physical meaning of this formula is the collision force between ice particles and the ship's hull along the normal direction under the updated motion state. The technical advantage lies in achieving a more accurate simulation of energy transfer and dissipation during the collision process by combining contributions from both stiffness and damping, thus improving the physical realism and numerical stability of the ice-ship contact force calculation. The independent variable in this formula includes the elastic modulus. ice particle radius Normal relative velocity and normal damping coefficient These factors collectively determine the magnitude of the normal force: the elastic modulus and radius reflect the contact stiffness through their product; the normal velocity affects the contribution of the rigid force through its modulus and determines the direction of the damping force through its vector direction; the damping coefficient controls the intensity of velocity-related dissipation. This formula reflects the positive correlation between the dependent variable and the elastic modulus, radius, and normal velocity modulus, while also demonstrating the relationship through an independent damping term. A linear positive correlation with the normal velocity vector was achieved, thus simultaneously characterizing the elastic recovery and energy dissipation mechanisms of the contact in the model.

[0041] The relative velocity obtained by subtracting the ship's speed from the updated ice particle's velocity is further subtracted by its projection onto the normal direction at the contact point between the ice particle and the ship, yielding the tangential velocity component. The tangential force is calculated by multiplying four parts consecutively: the first part is the product of the friction coefficient and the updated normal force; the second part is the quotient obtained by dividing the tangential velocity component by its modulus and the sum of a minimal near-zero constant; the third part is the hyperbolic tangent function value of the ratio of the tangential velocity modulus to an ideal reference velocity; the product of these three parts is then multiplied by a unit tangential vector used to determine the direction, finally yielding the updated tangential frictional force of the ice particle. The formula used in the above process is: in, in, Indicates the updated number The first grid cell of the first grid cell The tangential friction of an ice particle; Indicates the updated number The first grid cell of the first grid cell The tangential velocity component of an ice particle; In the above process, by constructing a friction model that includes velocity normalization and hyperbolic tangent function modulation, a smooth transition from static friction to dynamic friction with physical continuity and numerical stability was achieved within the Coulomb friction framework. This calculation not only accurately reflects the proportional relationship between tangential friction and normal load and its correlation with the friction coefficient, but also simulates the real nonlinear behavior of friction force increasing rapidly with increasing tangential velocity and then gradually saturating through the hyperbolic tangent function. At the same time, the use of zero constant and normalization treatment effectively avoids numerical singularities at low or zero speeds, thus ensuring the numerical robustness and physical fidelity of the ice-ship tangential interaction simulation. The dependent variable in this formula is the updated tangential frictional force of the ice particles. The physical meaning of this formula is the sliding friction force between ice particles and the ship hull along the tangential direction under the updated motion state. The technical advantage lies in achieving a smooth transition and numerically stable calculation of static-dynamic friction by introducing hyperbolic tangential function modulation and unit tangential vector normalization. This avoids numerical singularities that may occur when the tangential velocity is too small or zero, thereby improving the continuity and physical realism of the ice-ship tangential interaction simulation. The independent variable of this formula includes the updated normal force. coefficient of friction Updated tangential velocity and ideal reference speed These factors together determine the magnitude and direction of the tangential force: the product of the normal force and the friction coefficient forms the basis of Coulomb friction; the tangential velocity is normalized to obtain the unit tangential vector to determine its direction, and the hyperbolic tangent function is used to achieve a smooth transition in the low-speed region and saturation characteristics in the high-speed region. This formula reflects the positive correlation between the dependent variable and the normal force and the friction coefficient, while also using the hyperbolic tangent function... It achieves nonlinear dependence on tangential velocity: it increases approximately linearly at low speeds and tends to saturate at high speeds, preventing zero constant. Numerical stability at zero tangential velocity was ensured.

[0042] Based on the normal force and tangential friction of the ice particles, the updated normal load and tangential friction of the measured mesh element are further calculated. Specifically, for each mesh element, the updated normal forces of the ice particles corresponding to all its contact points are accumulated, and the sum is divided by the surface area of ​​the mesh element to obtain the updated normal load; at the same time, the updated tangential friction of the ice particles corresponding to all its contact points are accumulated, and the sum is divided by the surface area of ​​the mesh element to obtain the updated tangential friction.

[0043] The formula used in the above process is: in, Indicates the updated number Normal load of each grid cell; Indicates the updated number Tangential friction force of each grid cell.

[0044] In the above process, the "force region average distribution method" is used to perform spatial integration and area averaging on all updated normal forces and tangential friction forces in the mesh element. Thus, based on the dynamic evolution of the ice particle motion state, an equivalent surface load that strictly corresponds to the ice-ship interaction state at the current moment and can be used for structural mechanics analysis is generated. This not only reflects the transient characteristics of ice load changing with time and particle motion, but also ensures the continuity and applicability of load data in finite element analysis through spatial averaging.

[0045] S5: For each grid cell to be tested, the stress components are solved based on the updated normal load and tangential friction force. The hot spot stress is calculated and corrected to obtain the equivalent alternating stress amplitude. The fatigue damage coefficient is calculated based on the equivalent alternating stress amplitude and the number of contacts between the grid cell to be tested and the ice particles. The fatigue damage coefficient is classified according to the preset threshold.

[0046] The equivalent alternating stress amplitude is calculated using the following steps: Based on the updated normal load and tangential friction of the grid cells to be tested, the hot spot stress of each grid cell is calculated using the Von Mises criterion method. Specifically, the normal load is squared, and then three times the square of the tangential friction is added. The sum is then squared, and the result is the hot spot stress value of the grid cell. The formula used in the above process is: in, Indicates the first Hot spot stress in each grid cell; In the above process, the updated normal load and tangential friction force, two mutually perpendicular stress components, are synthesized using the Von Mises yield criterion to calculate the hot spot equivalent stress, which reflects the yielding tendency of the material under complex stress conditions. This calculation simplifies the multiaxial stress state into a scalar value that can be directly compared with the material fatigue performance parameters, providing a unified and physically meaningful stress input for subsequent fatigue damage assessment based on the equivalent stress amplitude. It is a key mechanical conversion link connecting local ice loads and fatigue life prediction of ship structures.

[0047] In the above embodiments, 20 sets of data on normal load and corresponding hot spot stress are given to reflect the change of hot spot stress with the change of normal load, as shown in Table 1: Table 1: Relationship between normal load and corresponding hot spot stress In Table 1 above, given a tangential friction force of 10 MPa for a given mesh element, changing the normal load on the mesh element shows that the larger the value of the normal load, the higher the hot spot stress of the mesh element.

[0048] Furthermore, Goodman correction is applied to the hot spot stress to obtain the equivalent alternating stress amplitude. Specifically, the hot spot stress is divided by the quotient of the normal load and the ultimate tensile strength of the marine steel, and the result is the equivalent alternating stress amplitude of the mesh element.

[0049] The formula used in the above process is: in, Indicates the ultimate tensile strength of marine steel; Indicates the first The equivalent alternating stress amplitude of each grid cell.

[0050] In the above process, the Goodman model is used to convert the actual alternating stress containing static average stress components into the equivalent alternating stress amplitude under the zero average stress reference, thereby eliminating the influence of static stress caused by normal load on the material fatigue life assessment. This correction ensures that the subsequent fatigue damage calculation based on the SN curve can accurately reflect the fatigue characteristics of the material under pure alternating load, and significantly improves the accuracy and reliability of fatigue life prediction of ship structures under complex ice load. The dependent variable in this formula is the equivalent alternating stress amplitude. Its physical meaning is the equivalent stress amplitude, after mean stress correction, applicable to fatigue life assessment. The technical effect lies in eliminating the influence of static mean stress on material fatigue strength through the Goodman correction model, thereby transforming the complex stress state under actual working conditions into an effective alternating stress that can be directly compared with the standard SN curve, significantly improving the accuracy of fatigue damage prediction; the independent variables of this formula include hot spot stress. Normal load and the ultimate tensile strength of materials These factors collectively determine the corrected stress amplitude: the hot spot stress is the original alternating stress to be corrected; the normal load characterizes the average stress component in the stress cycle; and the ultimate tensile strength provides a benchmark for the material's ability to resist static tensile failure, and is a key material constant in the Goodman correction. The formula reflects the positive correlation between the dependent variable and the hot spot stress, and the negative correlation with the normal load, i.e., the equivalent amplitude decreases as the average stress increases, as reflected in the denominator. The reduction and correction of the model conforms to the classical fatigue mechanics law that increased average stress reduces fatigue strength.

[0051] Based on the threshold results, the corresponding alarm is issued. The specific steps are as follows: The number of contacts between the grid cell under test and the ice particles is counted. Based on the equivalent alternating stress amplitude, the hull fatigue damage coefficient of the grid cell under test is calculated. Specifically, the total number of contacts between the grid cell under test and the ice particles is counted. Based on the equivalent alternating stress amplitude generated by each contact, the hull fatigue damage coefficient of the grid cell is calculated. In the calculation, the equivalent alternating stress amplitude corresponding to each contact of the grid cell is taken as the specified power of the slope of the SN curve. Then, the power of the results of all contact numbers is accumulated. Finally, the accumulated sum is divided by the SN curve constant of the hull material. The quotient is the fatigue damage coefficient of the grid cell. In the above process, the Palmgren-Miner linear cumulative damage theory is used to convert the different equivalent alternating stress amplitudes experienced by each grid element in multiple ice collisions into linearly superimposed damage contributions based on the fatigue performance characterized by the material's SN curve. Finally, the cumulative fatigue damage coefficient of the element is obtained by summing these contributions. This calculation process transforms discrete and random ice load events into continuous and quantifiable structural damage indicators, realizing a reliable mapping from local dynamic stress history to macroscopic fatigue damage degree. This provides a direct and quantitative scientific basis for fatigue life prediction and maintenance decisions of ship structures in polar environments. The SN curve constants for hull materials typically refer to the fatigue strength coefficient and slope, obtained through standardized material fatigue tests. Specifically, a series of fatigue tests are conducted on steel specimens representing the hull structure under different constant stress amplitudes. The number of cycles the specimens undergo until failure is recorded. Subsequently, a linear fit is performed on the stress amplitude and the number of cycles until failure in a double logarithmic coordinate system. The stress amplitude value at the position corresponding to the specified number of cycles on the fitted line can be converted to obtain the fatigue strength coefficient. The Palmgren-Miner linear cumulative damage theory is a widely used cumulative damage hypothesis in engineering fatigue analysis. This theory assumes that: fatigue damage to materials under variable amplitude loads can be linearly accumulated; the damage caused by each stress cycle is related to the stress amplitude of that cycle and is independent of other cycles; fatigue failure occurs when the cumulative damage value reaches 1; each contact with ice particles is regarded as an independent stress cycle, and its equivalent alternating stress amplitude is converted into the theoretical failure cycle number through the SN curve, thereby realizing the linear cumulative quantification of fatigue damage under random ice loads. The fatigue damage coefficient of the hull is classified into thresholds, and corresponding alarms are issued based on the classification results: when the calculated fatigue damage coefficient is greater than or equal to 0.3, it indicates that the damage accumulation of the grid cell has entered the accelerated stage, and the system will issue a reminder to prompt for enhanced monitoring; when the fatigue damage coefficient is greater than or equal to 0.7, it indicates that the remaining life of the grid cell has been significantly reduced, and an immediate overhaul should be arranged; when the fatigue damage coefficient is greater than or equal to 1.0, it indicates that the grid cell has entered a state of fatigue failure risk, and it must be immediately stopped and overhauled.

[0052] In the above process, under the Palmgren-Miner linear cumulative damage theory framework, when the damage coefficient reaches around 0.3, the micro-cracks inside the material often begin to enter a stable propagation stage, and the damage accumulation rate accelerates significantly. At this point, early warning can achieve early intervention during the damage acceleration period. When the damage coefficient reaches around 0.7, most of the remaining life of the structure has been consumed, and crack propagation enters an unstable stage. Maintenance needs to be arranged within a limited period to avoid functional failure in the short term. This threshold corresponds to the critical point of the time window for planned maintenance. The threshold of 1.0 directly corresponds to the theoretical fatigue failure limit. This classification not only conforms to the three-stage physical law of fatigue damage from initiation to failure, but also takes into account the progressive decision-making process of monitoring, early warning, and maintenance in ship operation, achieving an engineering balance between safety and economy.

[0053] Please see Figure 3 The present invention also provides a discrete element method-based ice load calculation device for ships interacting with ice ridges, the device being used to perform the above-described calculation method, comprising: The ice particle construction module is used to construct a 3D model of the ship under test and establish a global coordinate system. Under the global coordinate system, the seawater monitoring area is delineated, the ship under test is meshed using finite element methods, and the meshes in contact with the ice ridges are marked. The CT scan results of the ice ridges in contact with the ship are obtained, and the discrete element ice particles consistent with the geometric boundaries are generated using a fractal algorithm. The force calculation module is used to collect temperature, salinity, and porosity parameters for each ice particle according to the grid position and calculate the elastic modulus parameter of each ice particle; a dynamic stiffness model is established based on the elastic modulus parameter and the normal force of each ice particle is calculated using the dynamic stiffness model; then, the tangential force of each ice particle is calculated using the Coulomb friction method based on the normal force of each ice particle. The fluid force calculation module is used to count the number of contact points between the hull mesh cells to be tested and the ice particles, and to calculate the normal load and tangential friction force of the mesh cells to be tested; to divide the seawater monitoring area into fluid particle swarms, calculate the acceleration of each fluid particle, and update the velocity of each fluid particle within the seawater monitoring area; based on the velocity of each fluid particle at the corresponding ice particle and the velocity of the ice particle, to calculate the total drag force of all neighboring fluid particles on the ice particle within the preset search area. The parameter update module is used to calculate the acceleration of ice particles and update the velocity of ice particles based on the drag force of fluid particles on ice particles. Based on the updated velocity of ice particles, the normal force and tangential friction force of ice particles are recalculated, and the normal load and tangential friction force of the mesh element under test are updated synchronously to complete the update of the hull and ice ridge load status at the current time step. The classification alarm module is used to solve the stress components for each grid cell under test based on the updated normal load and tangential friction force, calculate the hot spot stress and correct it to obtain the equivalent alternating stress amplitude, calculate the fatigue damage coefficient based on the equivalent alternating stress amplitude and the number of contactes between the grid cell under test and the ice particles, and classify the fatigue damage coefficient according to the preset threshold.

[0054] The above formulas are all dimensionless calculations. The formulas are derived from software simulations based on a large amount of collected data to obtain the most recent real-world results. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.

[0055] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented in software, the above embodiments can be implemented, in whole or in part, as a computer program product. Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented by electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution.

[0056] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment, depending on actual needs.

[0057] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A method for calculating ice loads under the action of ships and ice ridges based on discrete element method, characterized in that, include: S1: Construct a 3D model of the ship to be tested and establish a global coordinate system. Delineate the seawater monitoring area under the global coordinate system. Divide the ship to be tested into a finite element mesh and mark the mesh in contact with the ice ridge. Obtain the CT scan results of the ice ridge in contact with the ship and use a fractal algorithm to generate discrete element ice particles consistent with the geometric boundary. S2: Collect temperature, salinity, and porosity parameters for each ice particle according to the grid position and calculate the elastic modulus parameter of each ice particle; establish a dynamic stiffness model based on the elastic modulus parameter and calculate the normal force of each ice particle using the dynamic stiffness model; then calculate the tangential force of each ice particle based on the normal force of each ice particle using the Coulomb friction method. S3: Count the number of contact points between the hull mesh cells to be tested and the ice particles, and calculate the normal load and tangential friction force of the mesh cells to be tested respectively; divide the seawater monitoring area into fluid particle groups, calculate the acceleration of each fluid particle, and update the velocity of each fluid particle in the seawater monitoring area; based on the velocity of each fluid particle at the corresponding ice particle and the velocity of the ice particle, calculate the total drag force of all neighboring fluid particles on the ice particle in the preset search area. S4: Based on the drag force of fluid particles on ice particles, calculate the acceleration of ice particles and update the velocity of ice particles. Based on the updated velocity of ice particles, recalculate the normal force and tangential friction force of ice particles. Simultaneously update the normal load and tangential friction force of the mesh element under test, and complete the update of the hull and ice ridge load status at the current time step. S5: For each grid cell to be tested, solve the stress components based on the updated normal load and tangential friction force, calculate the hot spot stress and correct it to obtain the equivalent alternating stress amplitude, calculate the fatigue damage coefficient based on the equivalent alternating stress amplitude and the number of contactes between the grid cell to be tested and the ice particles, and classify the fatigue damage coefficient according to the preset threshold. The velocity of each fluid particle is updated using the Leapfrog integral method. The specific steps are as follows: The cubic array method is used to divide the seawater monitoring area into fluid particle swarms. Specifically, a three-dimensional monitoring cubic domain is defined according to the spatial range of the seawater monitoring area. It is regarded as a regular hexahedral computation space. Within the cubic domain, the entire water body is equally cut into a series of regularly arranged tiny cubic grids according to a preset uniform spatial interval. The vertex position of each cubic grid is used as the spatial coordinate of an initial fluid particle. The SPH method is used to calculate the acceleration of each fluid particle within a preset neighboring fluid particle region. Specifically, it involves: first, identifying all neighboring fluid particles within the preset neighboring region except for itself; for each neighboring particle, multiplying its mass by the sum of the quotient of the current particle's pressure and the square of its own density, and the quotient of the neighboring particle's pressure and the square of its own density, then multiplying by the gradient value of the kernel function between these two particles, and summing the calculation results for all neighboring particles within the region. This negative summation, plus the gravitational acceleration, is the acceleration of the fluid particle at the current time step. The pressure of the current particle is calculated using a state equation: the state equation is based on the ratio of particle density to reference density, where the particle density is obtained by dividing the sum of the weighted contributions of neighboring particles to its kernel function by the reference density, and then raising the calculated ratio to a specified power and subtracting one, finally multiplying by the stiffness coefficient of seawater. The velocity of each fluid particle is updated using the Leapfrog integral method. Specifically, the velocity value of each fluid particle in the previous half-time step is added to the product of the acceleration value of the fluid particle in the current integer time step and the preset integration time step, so as to obtain the velocity value of the fluid particle in the new half-time step.

2. The method for calculating ice loads under the action of ships and ice ridges based on discrete element method according to claim 1, characterized in that, The ice particles that generate discrete elements are as follows: A three-dimensional model of the ship to be tested is constructed and a global coordinate system is established. The positive direction of the Z-axis of the global coordinate system is defined as vertically upward. The X-axis is set to point towards the bow and stern of the ship, and the Y-axis points towards the port and starboard sides of the ship. The seawater monitoring area is delineated under the global coordinate system. The monitoring area extends forward from the foremost point of the bow by 1.0 to 1.5 times the ship length, backward from the stern by 0.5 to 1.0 times the ship length, downward from the bottom of the ship by 2.5 to 3.0 times the ship's draft, and downward to 0.5 to 1.0 times the ship's width above the deck. The total width of the ship, including both port and starboard sides, is 2.0 to 3.0 times the ship's width. The hull under test was meshed using finite element methods, and the meshes in contact with the ice ridge were marked. Specifically, the outer surface geometric data were extracted based on the hull CAD model, and meshes were generated for the bow, weld seams, deck, and stern. A denser mesh was used at the bow and weld seams, while a sparser mesh was used at the deck and stern. All meshes were then filtered and marked to identify those meshes in contact with the ice ridge. The CT scan results of the ice ridges in contact with the ship's hull are obtained, and discrete element ice particles are generated using a fractal algorithm. Specifically, the CT scan results are denoised, the porosity distribution characteristic parameters of the ice ridge CT scan results are extracted, and multiple fractal surfaces are generated in three-dimensional space using the random midpoint displacement method. These surfaces divide the space into ice matrix regions and porosity regions. The Voronoi cutting algorithm is used to generate discrete element ice particles of polyhedra in the ice matrix regions.

3. The method for calculating ice loads under the action of ships and ice ridges based on discrete element method according to claim 2, characterized in that, The specific steps for calculating the normal and tangential forces on the ice particles are as follows: The temperature, salinity, and porosity parameters of ice particles at each grid point are obtained, and the elastic modulus parameter of each ice particle is calculated. Specifically, based on 10.2, an exponential function value with the natural constant e as the base is multiplied. The exponent of the exponential function value is calculated by adding 273 to the temperature value of the ice particle, dividing by 120, and taking the negative value to obtain the first calculation term. The first calculation term is then multiplied by a salinity correction factor obtained by subtracting 0.15 from 1 and multiplying by the square root of salinity to obtain the second calculation term. The first and second calculation terms are multiplied together, and the product is then multiplied by a porosity correction factor obtained by subtracting the cube of the porosity from 1, thereby obtaining the elastic modulus of the ice particles in that grid cell. Based on the elastic modulus parameter, a dynamic stiffness model is constructed to calculate the normal force of each ice particle in the grid element under test. Specifically, the value of the normal force is obtained by multiplying a calibration coefficient by the elastic modulus of the ice particle, then by the radius of the ice particle, and finally by the modulus of the normal velocity component of the ice particle. The normal velocity component is obtained by performing a dot product operation between the ice particle velocity vector and the contact point normal vector. Furthermore, based on the normal force of each ice particle, the Coulomb friction method is used to calculate the tangential force of the ice particles in the grid cell under test. Specifically, the tangential velocity vector is obtained by subtracting the normal velocity component from the velocity vector of the ice particle; the tangential velocity vector is then divided by the sum of its magnitude and a minimal zero constant to obtain the unit tangential vector; the tangential force is finally determined by multiplying the product of the friction coefficient and the normal force by the unit tangential vector, and then by the hyperbolic tangent function value of the ratio of the tangential velocity magnitude to the ideal reference velocity.

4. The method for calculating ice loads under the action of ships and ice ridges based on discrete element method according to claim 3, characterized in that, The normal load and tangential friction force of the mesh element under test are calculated using the force region average distribution method. The specific steps are as follows: The number of contact points between the tested grid cells and ice particles on the hull was counted. The normal load and tangential friction force of the tested grid cells were calculated using the force region average distribution method. Specifically, the normal load was calculated by summing the normal forces of the ice particles corresponding to all contact points on the grid cell and then dividing the sum by the surface area of ​​the grid cell. The tangential friction force was calculated by summing the tangential forces of the ice particles corresponding to all contact points on the grid cell and then dividing the sum by the surface area of ​​the grid cell.

5. The method for calculating ice loads under the action of ships and ice ridges based on discrete element method according to claim 1, characterized in that, The specific steps for calculating the drag force of fluid particles on ice particles are as follows: Based on the velocity of each fluid particle at its corresponding ice particle and the velocity of the ice particle itself, the SPH-DEM bidirectional coupling method is used to calculate the total drag force exerted by all neighboring fluid particles on the ice particle within a preset search area. Specifically, for each neighboring fluid particle within the preset search area, the vector difference between the fluid particle's velocity and the ice particle's velocity, as well as its modulus, are first calculated. Then, half of this difference is multiplied by the drag coefficient, then by the seawater density, then by the projected area of ​​the ice particle, then by the modulus, and finally by the vector difference to obtain the drag force contribution of that single fluid particle on the ice particle. Finally, the drag force contributions of all neighboring fluid particles within the search area are vector-summed to obtain the total drag force exerted by all fluid particles on the ice particle.

6. The method for calculating ice loads under the action of ships and ice ridges based on discrete element method according to claim 5, characterized in that, The specific steps for updating the normal load and tangential friction force of the mesh element under test are as follows: Based on the drag force of fluid particles on ice particles, the acceleration of ice particles is calculated using Newton's second law. Specifically, the total drag force on the ice particle is vector-added with the gravity obtained by multiplying the mass of the ice particle by the acceleration due to gravity. The total force on the ice particle is then divided by the mass of the ice particle, and the result is the acceleration of the ice particle. At the same time, the velocity of the ice particles is updated by vector addition of the current velocity of the ice particles and the calculated acceleration of the ice particles multiplied by the preset integration time step. The sum is the updated velocity of the ice particles. Based on the updated ice particle velocities, the normal force and tangential friction force of the ice particles are recalculated. Specifically, the updated ice particle velocity is subtracted from the ship's velocity, and the resulting relative velocity vector is projected onto the normal direction of the contact point between the ice particle and the ship. The magnitude and direction of the normal relative velocity are determined by calculating the cosine of the angle between the projection and the normal vector. The normal force of the ice particle consists of two parts: the first part is calculated by multiplying the calibration coefficient by the elastic modulus of the ice particle, then by the radius of the ice particle, and finally by the modulus of the normal relative velocity; the second part is calculated by multiplying the normal damping coefficient by the normal relative velocity. The updated normal force of the ice particle is obtained by adding these two parts together. The relative velocity obtained by subtracting the ship's speed from the updated ice particle's velocity is further subtracted by its projection onto the normal direction at the contact point between the ice particle and the ship, yielding the tangential velocity component. The tangential force is calculated by multiplying four parts consecutively: the first part is the product of the friction coefficient and the updated normal force; the second part is the quotient obtained by dividing the tangential velocity component by its modulus and the sum of a minimal near-zero constant; the third part is the hyperbolic tangent function value of the ratio of the tangential velocity modulus to an ideal reference velocity; the product of these three parts is then multiplied by a unit tangential vector used to determine the direction, finally yielding the updated tangential frictional force of the ice particle. Based on the normal force and tangential friction of the ice particles, the updated normal load and tangential friction of the measured mesh element are further calculated. Specifically, for each mesh element, the updated normal forces of the ice particles corresponding to all its contact points are accumulated, and the sum is divided by the surface area of ​​the mesh element to obtain the updated normal load; at the same time, the updated tangential friction of the ice particles corresponding to all its contact points are accumulated, and the sum is divided by the surface area of ​​the mesh element to obtain the updated tangential friction.

7. The method for calculating ice loads under the action of ships and ice ridges based on discrete element method according to claim 6, characterized in that, The equivalent alternating stress amplitude is calculated using the following steps: Based on the updated normal load and tangential friction of the grid cells to be tested, the hot spot stress of each grid cell is calculated using the Von Mises criterion method. Specifically, the normal load is squared, and then three times the square of the tangential friction is added. The sum is then squared, and the result is the hot spot stress value of the grid cell. Furthermore, Goodman correction is applied to the hot spot stress to obtain the equivalent alternating stress amplitude. Specifically, the hot spot stress is divided by the quotient of the normal load and the ultimate tensile strength of the marine steel, and the result is the equivalent alternating stress amplitude of the mesh element.

8. The method for calculating ice loads under the action of ships and ice ridges based on discrete element method according to claim 7, characterized in that, Based on the threshold results, the corresponding alarm is issued. The specific steps are as follows: The number of contacts between the grid cell under test and the ice particles is counted. Based on the equivalent alternating stress amplitude, the hull fatigue damage coefficient of the grid cell under test is calculated. Specifically, the total number of contacts between the grid cell under test and the ice particles is counted. Based on the equivalent alternating stress amplitude generated by each contact, the hull fatigue damage coefficient of the grid cell is calculated. In the calculation, the equivalent alternating stress amplitude corresponding to each contact of the grid cell is taken as the specified power of the slope of the SN curve. Then, the power of the results of all contact numbers is accumulated. Finally, the accumulated sum is divided by the SN curve constant of the hull material. The quotient is the fatigue damage coefficient of the grid cell. The fatigue damage coefficient of the hull is classified into thresholds, and corresponding alarms are issued based on the classification results: when the calculated fatigue damage coefficient is greater than or equal to 0.3, it indicates that the damage accumulation of the grid cell has entered the accelerated stage, and the system will issue a reminder to prompt for enhanced monitoring; when the fatigue damage coefficient is greater than or equal to 0.7, it indicates that the remaining life of the grid cell has been significantly reduced, and an immediate overhaul should be arranged; when the fatigue damage coefficient is greater than or equal to 1.0, it indicates that the grid cell has entered a state of fatigue failure risk, and it must be immediately stopped and overhauled.

9. A discrete element method-based ice load calculation device for ships interacting with ice ridges, characterized in that: The computing device is used to execute the computing method according to any one of claims 1-8, including: The ice particle construction module is used to construct a 3D model of the ship under test and establish a global coordinate system. Under the global coordinate system, the seawater monitoring area is delineated, the ship under test is meshed using finite element methods, and the meshes in contact with the ice ridges are marked. The CT scan results of the ice ridges in contact with the ship are obtained, and the ice particles with discrete elements consistent with the geometric boundaries are generated using a fractal algorithm. The force calculation module is used to collect temperature, salinity, and porosity parameters for each ice particle according to the grid position and calculate the elastic modulus parameter of each ice particle; a dynamic stiffness model is established based on the elastic modulus parameter and the normal force of each ice particle is calculated using the dynamic stiffness model; then, the tangential force of each ice particle is calculated using the Coulomb friction method based on the normal force of each ice particle. The fluid force calculation module is used to count the number of contact points between the hull mesh cells to be tested and the ice particles, and to calculate the normal load and tangential friction force of the mesh cells to be tested; to divide the seawater monitoring area into fluid particle swarms, calculate the acceleration of each fluid particle, and update the velocity of each fluid particle within the seawater monitoring area; based on the velocity of each fluid particle at the corresponding ice particle and the velocity of the ice particle, to calculate the total drag force of all neighboring fluid particles on the ice particle within the preset search area. The parameter update module is used to calculate the acceleration of ice particles and update the velocity of ice particles based on the drag force of fluid particles on ice particles. Based on the updated velocity of ice particles, the normal force and tangential friction force of ice particles are recalculated, and the normal load and tangential friction force of the mesh element under test are updated synchronously to complete the update of the hull and ice ridge load status at the current time step. The classification alarm module is used to solve the stress components for each grid cell under test based on the updated normal load and tangential friction force, calculate the hot spot stress and correct it to obtain the equivalent alternating stress amplitude, calculate the fatigue damage coefficient based on the equivalent alternating stress amplitude and the number of contactes between the grid cell under test and the ice particles, and classify the fatigue damage coefficient according to the preset threshold.

Images