Ice load calculation method for ship-ice ridge interaction based on discrete element method

Abstract

The present application relates to a kind of ice load calculation method based on discrete element under the action of ship and ice ridge, according to the sea area where ice ridge is located, the sea ice discrete element unit model of particle unit with corresponding cohesion effect of mass and size is established;Based on the buoyancy method, the discrete element model of ice ridge is generated by combining sea ice discrete element unit;The entire ship surface is discretized into a series of triangular units using grid division preprocessing, and the contact model between the discrete triangular units of the ship body and the spherical sea ice discrete element units is established;Dynamic calculation and analysis of ice load on the interaction between spherical sea ice discrete element units and discrete triangular units of the ship body;Parallel programming based on GPU is realized on CUDA platform using discrete element method, which realizes fast modeling of ice ridge model and ice load calculation under the dynamic action of ship-ice ridge.It greatly reduces the modeling difficulty of ice ridge, and improves the reusability and modifiability of ice ridge model.

Description

Technical Field

[0001] This invention relates to a ship hull design technology, and more particularly to a method for calculating ice loads on ships under the action of ice ridges based on discrete element method. Background Technology

[0002] Marine structures and ships operating or navigating in icy marine environments are significantly affected by ice loads. Traditionally, research on ice loads has relied heavily on prototype observations and model tests. However, in recent years, theoretical analysis methods, which offer lower costs and higher repeatability, have become increasingly widely used. Currently, research on ice loads faces several major challenges:

[0003] (1) Prototype observation and model test require high costs and are limited by the data obtained by monitoring technology.

[0004] (2) Due to the inherent characteristics of real sea ice, which is composed of individual ice crystal particles, the traditional finite element method uses a unit model composed of a continuum. During the interaction between the ship and the ice, the sea ice will undergo a transformation from a continuum to a discrete block. The traditional finite element method is difficult to accurately simulate the destruction mode of sea ice.

[0005] (3) Current research hotspots on ice loads mainly focus on the interaction between the ship hull and flat ice. However, the actual polar ice conditions are complex, and icebreakers have a high probability of encountering ice ridges when sailing in the polar regions. There are very few current research methods on ice loads under the interaction between the ship hull and ice ridges.

[0006] The latest relevant technologies are as follows:

[0007] Currently, there are two main technologies for predicting ice loads under structure-ice interaction:

[0008] Technology 1: Dalian University of Technology has developed a discrete element method (DEM) calculation program for ice load calculation. This program can analyze the microscopic parameters of the sea ice discrete element model to reflect its macroscopic mechanical properties, and can then calculate the ice load on offshore platforms from both smooth ice and broken ice.

[0009] This technology cannot construct ice ridge models to calculate ice loads under the interaction of ice ridges and structures based on discrete element method; it cannot accurately simulate the dynamic process of icebreakers hitting sea ice, which will result in large errors in the calculation of ice loads under ship-ice interaction.

[0010] Technology 2: Harbin Engineering University has developed a method for calculating ice loads by embedding a subroutine into the commercial software LS-dyna for finite element analysis. This subroutine can control the sea ice model to simulate the motion state under the action of the water flow field.

[0011] This technique requires too much time to build the ice ridge finite element model in the early stage, and it is inconvenient to modify it, which is not conducive to the reuse of the model. It is difficult to intuitively and accurately simulate the actual damage of sea ice using the traditional finite element method. If the ice ridge is simulated with the smallest possible mesh size, the demand for computing resources will increase significantly, and the calculation time will also increase significantly. Summary of the Invention

[0012] To address the aforementioned issues, a method for calculating ice loads under the interaction of ships and ice ridges based on the discrete element method (DEM) is proposed. This method employs GPU-based parallel programming on the CUDA platform using the DEM, enabling rapid modeling of ice ridges and calculation of ice loads under the dynamic interaction between the ship and the ice ridge. This solves the problems of Technique 1, which cannot simulate the dynamic process of icebreakers impacting sea ice and cannot establish a discrete element model of the ice ridge. It also reduces the time-consuming process of establishing ice ridge models, low model reusability, and lack of ability to accurately simulate the failure modes of ice ridges associated with the traditional finite element method (Technique 2).

[0013] The technical solution of this invention is: a method for calculating ice loads under the interaction of ships and ice ridges based on discrete element method, specifically including the following steps:

[0014] Step 1: Establish a sea ice material failure model: Based on the sea area where the ice ridge is located, establish a discrete element model of sea ice composed of particle units with corresponding bonding effects, based on their mass and size.

[0015] Step 2: Based on the buoyancy method, combine the discrete element units of sea ice to generate a discrete element model of the ice ridge;

[0016] Step 3: Establish the discrete element model of the hull: The entire hull surface is discretized into a series of triangular elements using mesh generation preprocessing, and a contact model is established between the discrete triangular elements of the hull and the discrete element elements of the spherical sea ice.

[0017] Step 4: Based on steps 2 and 3, perform dynamic calculation and analysis on the ice load of the interaction between the discrete element unit of the spherical sea ice and the discrete triangular element unit of the ship hull, obtain the force acting on each particle, and update the velocity and coordinates of each particle at the same time.

[0018] Step 5: Based on Step 4, implement GPU-based parallel programming on the CUDA platform to perform hull-ice ridge coupling calculation and analysis, and realize ice load analysis and calculation.

[0019] Furthermore, the bonding action in step one adopts a particle bonding model, which takes into account the freezing effect between units. The two spheres are bonded together in parallel, which can not only transmit force but also torque. In the parallel bonding model, an elastic bonding disk is set between the two bonding particle units, and the disk transmits the force and torque between the two units.

[0020] Furthermore, the discrete element model establishment method for the ice ridge in step two is as follows: First, the number of ice fragments N is set within a given computational region. This number is also the number of seed points in the segmentation algorithm. The N seed points are arranged in a close-packed hexagonal arrangement within the computational region. The coordinates of these seed points are {c i}, i = 1, ..., N; the polygons generated in this way are regular hexagons, and the initial coordinates of the seed points are {c i The following processing is performed: x i =c i +σξ i , Where c i x is the initial coordinate value of the i-th seed; i Based on the initial coordinates c i A series of coordinate points are calculated, and these points form a polygon; ξ i Divided into The coordinates are obtained by participating in the calculation separately; among them, the noise term and The value of follows a standard normal distribution; σ is a disorder factor, whose value determines the irregularity of the polygon; the spherical elements are filled into the generated polygon, and the final ice ridge shape is generated under the action of buoyancy, thus completing the modeling.

[0021] Furthermore, the ice ridge morphology is based on statistical analysis of the morphology of common ice ridges, including sail height, freezing layer, keel, and flat ice from top to bottom.

[0022] Furthermore, the method for determining the contact type between the discrete triangular unit of the hull and the discrete element of the spherical sea ice is as follows: First, it is necessary to determine whether the spherical particle overlaps with the plane where the triangular unit is located.

[0023] Then, if there is an overlap between the spherical particle and the plane containing the triangular unit, it is necessary to further determine whether the spherical particle is in contact with the triangular unit itself. There are three types of contact between the spherical particle and the triangular unit: sphere-triangle surface contact, sphere-triangle edge contact, and sphere-triangle vertex contact, hereinafter referred to as surface contact, edge contact, and angular contact. Surface contact means that the projection point Q of the center of the spherical particle is inside the triangular surface; edge contact means that the spherical particle is in contact with the three edges of the triangular unit; angular contact means that the spherical particle is in contact with the three vertices ABC of the triangular unit. The centroid method is used to determine whether the projection point of the center of the sphere is inside the triangular unit.

[0024] Finally, determine whether the spherical particle comes into contact with the three sides of the triangular unit, and whether the spherical particle comes into angular contact with the triangular boundary unit.

[0025] Furthermore, the contact force calculation method in step four, the dynamic calculation and analysis of ice load, is as follows:

[0026] Calculate the embedding amount ΔL of the spherical element in the boundary element.

[0027] ΔL=|P C P|-R ball

[0028] In the formula, point P C The contact point between the spherical particle and the boundary element is represented by: The normal vector between the contact point of the particle and the triangular element is expressed as:

[0029]

[0030] The relative displacement Δx is expressed as:

[0031] Δx=(v p -v w )Δt

[0032] Among them, v p Let v be the velocity of particle P. w Contact point P on the triangular unit C The velocity is denoted by Δt, which is the displacement time.

[0033] The relative displacement along the contact surface is decomposed into a normal component Δx. n and tangential component Δx s :

[0034] Δx n =(n w Δx)n w

[0035] Δx s =Δx - Δx n

[0036] The interaction forces between particles and boundary elements are calculated using a linear contact model, where the normal force is expressed as:

[0037]

[0038] The tangential force is calculated incrementally and obeys the Mohr-Coulomb law:

[0039] If |F s |<|μ w F n |

[0040] F s ←F s (μ w F n / |Fs |), if |F s |≥|μ w F n |

[0041] In the formula, μ represents the normal and tangential contact stiffness between the particles and the boundary elements. w This represents the maximum friction coefficient between the particle and the boundary unit.

[0042] A method for calculating ice load under dynamic interaction between a ship and an ice ridge is proposed. This method uses a discrete element method to calculate ice load under the interaction between a ship and an ice ridge, and inputs different real-time bollard thrusts to calculate the ship's ice resistance.

[0043] A computer device includes a memory and a processor, the memory storing a computer program, characterized in that the processor executes the computer program to implement the aforementioned method for calculating ice loads under the action of ships and ice ridges based on discrete element method.

[0044] A computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the aforementioned method for calculating ice loads under the action of ships and ice ridges based on discrete element method.

[0045] A computer program product includes a computer program that, when executed by a processor, implements the aforementioned method for calculating ice loads under the action of a ship and an ice ridge based on discrete element method.

[0046] The beneficial effects of this invention are as follows: The method for calculating ice loads under the interaction between a ship and an ice ridge based on the discrete element method greatly reduces the difficulty of modeling ice ridges, improves the reusability and modifiability of ice ridge models, realizes dynamic calculations similar to those of the traditional finite element method for simulating active icebreaking of the ship, fully utilizes the obvious advantages of the discrete element method in calculating ship ice loads and accurately simulating sea ice destruction modes, and solves the problem of the lack of effective and relatively accurate research methods under the interaction between the ship and ice ridge. Attached Figure Description

[0047] Figure 1 This is a model diagram of the bonding between particles in this invention;

[0048] Figure 2 This is a geometric diagram of the ice ridge of the present invention;

[0049] Figure 3 This is a diagram of the ice ridge morphology generated by the present invention;

[0050] Figure 4 This is a schematic diagram of the discrete element model of the ship hull of the present invention;

[0051] Figure 5 This is a volume diagram of the tetrahedron of the present invention;

[0052] Figure 6 This is a schematic diagram of the triangular plate of the present invention that is not fully submerged in water;

[0053] Figure 7 This is a schematic diagram of the projection point of the sphere center P of the present invention onto the plane containing the triangular unit;

[0054] Figure 8 This is a diagram showing the contact types between the spherical particles and the triangular units of the present invention.

[0055] Figure 9 This is a schematic diagram illustrating the method for determining whether a projection point is inside a triangular unit according to the present invention.

[0056] Figures 10A-10B This is a schematic diagram showing the contact between the spherical particle and the triangular side of the present invention;

[0057] Figure 11 This is a schematic diagram showing the contact between the spherical particle and the vertex of the triangle in this invention;

[0058] Figure 12 This is a model diagram of an ice ridge of the present invention;

[0059] Figures 13A-13C This is a top view of the ship-ice ridge interaction process simulated by the discrete element method of this invention;

[0060] Figures 14A-14C This is a side view of the ship-ice ridge interaction process simulated by the discrete element method of this invention;

[0061] Figures 15A-15B For the present invention T B Schematic diagram of ship ice resistance at 2400kN;

[0062] Figures 16A-16B For the present invention T B Schematic diagram of ship ice resistance at 3000kN. Detailed Implementation

[0063] The present invention will now be described in detail with reference to the accompanying drawings and specific embodiments. These embodiments are based on the technical solution of the present invention and provide detailed implementation methods and specific operating procedures. However, the scope of protection of the present invention is not limited to the following embodiments.

[0064] Step 1: Establish a sea ice material failure model

[0065] A discrete element model of sea ice is established, composed of particle units of a certain mass and size with corresponding binding effect. For example... Figure 1 The particle bonding model shown considers the freezing effect between units, bonding two spheres together with parallel adhesive. This not only transmits forces but also torques. The coordinates of the element are denoted by . In the parallel bonding model, an elastic bonding disk is set between two bonding particle elements. The disk can transmit the forces and moments between the two elements, namely axial force, shear force, bending moment, and torque. Furthermore, both forces and moments can be represented by their respective normal and tangential components. These are the normal and tangential components of the force and torque, respectively, n. i Let R be the normal vector connecting the centers of the two spheres, and R be the radius of the elastically bonded disk.

[0066] Step 2, Ice Ridge Morphology

[0067] The geometric shape of the ice ridge is as follows Figure 2 As shown, this is based on statistical analysis of the morphology of common ice ridges. A: sail height; B: frozen layer; C: keel; D: flat ice. Discrete element models of ice ridges are generated by combining sea ice discrete element units using the buoyancy method. The Voronoi segmentation algorithm is used to generate ice fragments of different sizes and shapes within the ice ridge. First, the number of ice fragments N is set within a given computational region. This number is also the number of seed points in the segmentation algorithm. The N seed points are arranged in a close-packed hexagonal arrangement (two-dimensional) within the computational region. The coordinates of these seed points are {c i}, i = 1, ..., N. The polygons generated in this way are regular hexagons, with the initial coordinates of the seed point {c i The following processing is performed: x i =c i +σξ i , Where c i x is the initial coordinate value of the i-th seed; i Based on the initial coordinates c i A series of coordinate points are calculated, and these points form a polygon; ξ i Divided into The coordinates are obtained by participating in the calculation separately.

[0068] Among them, noise item and The values ​​of σ follow a standard normal distribution; σ is a disorder factor, and its value determines the irregularity of the polygon. Spherical elements are filled into the generated polygon, and under the action of buoyancy, the final ice ridge shape is generated, as shown in the figure. Figure 3 As shown.

[0069] Step 3: Establish the discrete element model of the ship hull.

[0070] The entire hull surface was discretized into a series of triangular elements using mesh generation preprocessing software, and a contact model was established between the triangular elements and the discrete element elements of the spherical sea ice. For example... Figure 4This is a schematic diagram of the discrete element model of the ship's hull.

[0071] Step 4: Calculation of hull buoyancy and towing force

[0072] Considering the effects of seawater buoyancy and ocean current drag on the ship's motion, the buoyancy is calculated using Archimedes' principle, and the formula is as follows:

[0073] F b =-ρ w gV sub

[0074] In the formula, ρ w Let g be the density of seawater, g be the acceleration due to gravity, and V be the acceleration due to gravity. sub This refers to the volume of the ship's hull submerged in water.

[0075] In the discrete element model, the hull is composed of several triangular plates. The following method is used to calculate the volume of the hull submerged in water: Take any point P on the water surface as a vertex. All triangular plates are calculated with P as the vertex. Construct tetrahedra with triangles V1, V2, and V3 below the waterline as their bases. The sum of the volumes of all tetrahedra is the volume of the hull submerged in water. The volume V of any one tetrahedron (e.g., V1V2V3) is the volume of the hull submerged in water. Figure 5 (As shown) can be calculated using the following formula:

[0076]

[0077] In the formula, a = V2 - V1, b = V3 - V1, r = p - V1, where V1, V2, and V3 are the coordinates of the three vertices of the triangular unit; p is the coordinate of any point on the water surface; a, b, and r are pointers; and the body center vector of the tetrahedron is:

[0078]

[0079] Therefore, the volume of the ship's hull submerged in water can be expressed as the sum of the volumes of all tetrahedrons with the triangles below the waterline as their bases:

[0080]

[0081] The body center vector of the submerged portion can be represented as:

[0082]

[0083] V i c i Let be the volume and body center vector of the i-th tetrahedron submerged in water, respectively.

[0084] The hull discretized from triangular elements will contain triangular plates that protrude above the waterline, requiring special consideration during volume calculations. For Figure 6The triangular plate shown is not fully submerged in water. The part that is above the water surface needs to be cut off, and only the volume of the underwater part needs to be calculated.

[0085] Thus, the volume V of the ship submerged in water was obtained. sub This allows us to calculate the buoyancy of the hull in the water. The drag force on a floating structure in water consists of two parts: the drag force F when the structure undergoes translational motion. d And the drag torque M during rotation d The two can be represented as:

[0086] F d =-C α ρ w V sub (v w -v s ))

[0087] M d =-C β ρ w V sub L 2 ω s

[0088] In the formula, C α C β These are the drag coefficients corresponding to the drag force and drag moment, respectively, and v. w For water flow velocity, v s Let ω be the speed of the ship's translation. s Let L be the angular velocity of the ship's rotation and L be the length of the floating body.

[0089] Step 5, Calculation of hull thrust

[0090] Propeller thrust X T The calculation is performed according to the following formula. In the formula, u is the longitudinal velocity of the ship in the hull coordinate system, u ow T is the maximum speed of the ship in open water. B This refers to the bollard tension (the tension on the mooring line when the propeller is running, as the ship is moored to a shore bollard).

[0091] Step 6, Contact Judgment

[0092] To determine whether a triangular element f is in contact with a spherical particle p, we first need to determine whether the spherical particle overlaps with the plane containing the triangular element. Let the coordinates of the center P of the spherical particle p be X. P The particle radius is R ball The projection point of the center of the sphere P onto the plane ABC containing the triangular unit is Q, as shown below. Figure 7 As shown.

[0093] The following vectors can be obtained: AP, BP, CP, BA, AB, CA, AC, CB, BC, PQ, if |PQ|≤R ball This indicates that the spherical particle overlaps with the plane containing the triangular unit ABC, suggesting that the spherical particle and the triangular unit may be in contact, requiring further determination of the contact relationship; for example, |PQ|>R. ball If the spherical particle does not make contact with the triangular unit, then the spherical particle will not come into contact with the triangular unit. This is called the first determination.

[0094] If the spherical particle overlaps with the plane containing the triangular element, it is necessary to further determine whether the spherical particle is in contact with the triangular element itself. If contact between the spherical particle and the triangular element is possible, it is necessary to determine their contact type in order to correctly calculate the contact force between the particle and the triangular element. There are three types of contact between the spherical particle and the triangular element: sphere-triangle surface contact, sphere-triangle edge contact, and sphere-triangle vertex contact, hereinafter referred to as surface contact, edge contact, and angular contact. Surface contact means that the projection point Q of the center of the spherical particle lies within the triangular surface; edge contact means that the spherical particle is in contact with the three edges of the triangular element; angular contact means that the spherical particle is in contact with the three vertices A, B, and C of the triangular element. The contact types between the spherical particle and the triangular element are as follows: Figure 8 As shown.

[0095] This method uses the centroid method to determine whether the projection point of the sphere's center lies inside the triangular element; this is called the second determination. The three vertices of the triangular element lie on the same plane. If one vertex is selected, the other two vertices can be considered as translations relative to that point, such as... Figure 9 The diagram shown illustrates how to determine whether a projection point is inside a triangular unit. If vertex A is chosen as the starting point, then vertex B is obtained by moving a certain distance in the AB direction, and vertex C is obtained by moving a certain distance in the AC direction.

[0096] For any point P in the plane, it can be represented as:

[0097] P = A + u(CA) + v(BA)

[0098] If the movement coefficients u or v of vertices C and B are negative, it is equivalent to point P moving in the opposite direction, i.e., the BA or CA direction; if the projection point P is located inside triangle ABC, u and v must simultaneously satisfy the following three conditions:

[0099] u≥0, v≥0, u+v≤1.

[0100] If the movement coefficients u and v of vertices C and B satisfy the above formula, it indicates that the projection point Q of the sphere center P within the triangular unit is inside the triangular unit, indicating that the spherical particle is in contact with the surface of the triangular unit.

[0101] The next step is to determine whether the spherical particle makes contact with the three sides of the triangular unit. Contact with the triangular unit occurs when the first condition is met but the second condition is not. If the projection of the particle's center P onto the side of the triangular unit lies inside the triangle and the distance from the center to the projection point is less than the particle radius R... ball At times, such as Figure 10A , 10B The diagram shown illustrates the contact between a spherical particle and the edge of a triangle, indicating that edge contact has occurred, which we refer to as the third determination.

[0102] Taking side AB as an example, the projection of point P onto side AB is Q. AB ,but

[0103]

[0104] If (AP·BA)(BP·BA)≤0, it indicates that the projection point Q AB On side AB, if (AP·BA)(BP·BA)>0, it indicates that the projection point Q is... AB It is not on side AB.

[0105] Therefore, the criterion for determining whether a spherical particle makes contact with side AB of a triangular unit is: when (AP·BA)(BP·BA)≤0 and |Q AB P|≤R ball The particle makes edge contact with the triangular unit; when (AP·BA)(BP·BA)>0 or |Q AB P|>R ball The particles do not make edge contact with the triangular units. Each triangular unit has three sides, so three judgments are required.

[0106] If the contact relationship between the spherical particle and the triangular unit satisfies the first criterion but not the second and third, it indicates that the particle may be in contact with the vertices of the triangular unit. This is further confirmed if the distance between the particle's center P and vertices A, B, and C is less than the particle's radius R. ball Angular contact occurs, such as Figure 11 The diagram shows the contact between the spherical particle and the vertex of the triangle.

[0107] The criterion for determining angular contact between spherical particles and triangular boundary elements is:

[0108] |AP|≤R ball It makes angular contact with vertex A.

[0109] |BP|≤R ball It makes angular contact with vertex B.

[0110] |CP|≤R ball It makes angular contact with vertex C.

[0111] Step 7, Contact Force Calculation

[0112] After determining the contact type between the spherical particles and the triangular elements, the embedding amount ΔL of the spherical elements in the boundary elements can be calculated.

[0113] ΔL=|P C P|-R ball

[0114] In the formula, point P C The contact points between the spherical particle and the boundary element are represented by points Q (surface contact, edge contact), A, B, and C (vertex contact). The normal vector between the contact points of the particle and the triangular element can be expressed as:

[0115]

[0116] The relative displacement Δx can be expressed as:

[0117] Δx=(v p -v w )Δt

[0118] Among them, v p Let v be the velocity of particle P. w Contact point P on the triangular unit C The velocity is given by Δt, where Δt is the displacement time.

[0119] The relative displacement along the contact surface is decomposed into a normal component Δx. n and tangential component Δx s :

[0120] Δx n =(n w Δx)n w

[0121] Δx s =Δx - Δx n

[0122] The interaction force between particles and boundary elements is calculated using a linear contact model, where the normal force can be expressed as:

[0123]

[0124] The tangential force is calculated incrementally and obeys the Mohr-Coulomb law:

[0125] If |F s |<|μ w F n |

[0126] Fs ←F s (μ w F n / |F s |), if |F s |≥|μ w F n |

[0127] In the formula, μ represents the normal and tangential contact stiffness between the particles and the boundary elements. w This represents the maximum friction coefficient between the particle and the boundary element. Once the force acting on each particle is determined, the velocity and coordinates of each particle can be updated.

[0128] Step 8, Ice Load Coupling Calculation and Analysis

[0129] The discrete element model of ice ridge-ship interaction is implemented on the CUDA platform using GPU-based parallel programming, defining computational conditions, and performing ice load prediction analysis.

[0130] Application examples:

[0131] 1. Based on steps 1 and 2, perform rapid modeling of the ice ridges using the ice ridge parameters in Table 1. The ice ridge model is as follows: Figure 12 As shown

[0132] Table 1 Ice Ridge Parameters

[0133]

[0134]

[0135] 2. Based on the actual ship parameters, establish the discrete element model of the icebreaker in step 3. It is composed of 4272 triangular boundary elements.

[0136] 3. Based on steps 4 to 7, write a dynamic analysis program for calculating ice loads.

[0137] 4. Based on step 8, implement GPU-based parallel programming on the CUDA platform to perform hull-ice ridge coupling calculation and analysis, and realize ice load analysis and calculation.

[0138] Ship-ice ridge interaction processes such as Figures 13A-13C Top view of the ship-ice ridge interaction process simulated by the discrete element method. Figures 14A-14C The figure shows a side view of the ship-ice ridge interaction process simulated by the discrete element method. The ship's ice resistance under different bollard thrusts is as follows: Figures 15A-15B T shown B Schematic diagram of ship ice resistance at 2400kN, shown in 15A~15B. B Schematic diagram of ship ice resistance at 3000kN.

[0139] Depend on Figures 15A-15B and Figures 16A-16B It can be seen that during the ship-ice ridge interaction process, different bollard thrusts have a significant impact on the ice resistance in the X and Z directions. When the bollard thrust is 2400 kN, the ice resistance in both the X and Z directions reaches its maximum value simultaneously. At approximately 80 s, the ice resistance in the X direction is about 2000 kN, and the ice resistance in the Z direction is about 3500 kN. After 100 s, the ice resistance in the X direction tends to stabilize at about 200 kN. When the bollard thrust is 3000 kN, the ice resistance in both the X and Z directions reaches its maximum value simultaneously. At approximately 70 s, the ice resistance in the X direction is about 2000 kN, and the ice resistance in the Z direction is about 3800 kN. After 85 s, the ice resistance in the X direction tends to stabilize at about 200 kN. It is evident that the two bollard thrusts have a significant difference in their impact on the ice resistance in the Z direction; the greater the bollard thrust, the greater the ice resistance in the Z direction. Furthermore, during the ship-ice ridge interaction process, the unloading rate of ice resistance in the X direction is significantly faster than that in the Z direction.

[0140] a. The discrete element method (DEM) program for calculating ice loads developed by Dalian University of Technology can only calculate and analyze the loads exerted on structures by smooth ice and broken ice. This invention has a wider range of applications and can calculate the ice loads acting between structures and ice ridges.

[0141] Solution: Based on the discrete element method, this invention creatively establishes a discrete element model of the ice ridge by using the buoyancy method, thereby realizing the calculation of ice load under the action of structure-ice ridge.

[0142] b. The discrete element method (DEM) program for calculating ice loads developed by Dalian University of Technology primarily enables mutual calculation between the sea ice DEM model and static structures. Compared to this program, this invention can calculate ice loads in the dynamic simulation process of a structure impacting sea ice.

[0143] Solution: Implement GPU-based parallel programming on the CUDA platform to simulate the dynamic process of structures impacting sea ice.

[0144] c. The ice load calculation method developed by Harbin Engineering University requires a significant amount of time to construct the finite element model of the ice ridge. In contrast, this invention, based on the discrete element method, can quickly and automatically generate the ice ridge model by setting the geometric parameters of the ice ridge, thus greatly shortening the preliminary calculation preparation time.

[0145] Solution: Based on the discrete element method, this invention creatively uses the buoyancy method to quickly generate a discrete element model of the ice ridge based on its geometric parameters.

[0146] d. The ice load calculation method developed by Harbin Engineering University uses the traditional finite element method to simulate sea ice. This method cannot accurately reflect the sea ice failure process, leading to significant errors in the calculation results. To achieve better simulation results, the sea ice element size needs to be as small as possible, which greatly increases the computational resource requirements and is time-consuming. In contrast, this invention uses the discrete element method, which has significant advantages in calculating ice loads and accurately simulating sea ice failure modes.

[0147] Solution: Based on the discrete element method, implement GPU-based parallel programming on the CUDA platform to calculate the ice load under the structure-ice interaction.

[0148] The embodiments described above are merely illustrative of several implementations of the present invention, and while the descriptions are relatively specific and detailed, they should not be construed as limiting the scope of the invention patent. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the present invention, and these all fall within the protection scope of the present invention. Therefore, the protection scope of this invention patent should be determined by the appended claims.

Claims

1. A method for calculating ice loads under the interaction of a ship and an ice ridge based on the discrete element method, characterized in that, Specifically, the steps include the following: Step 1: Establish a sea ice material failure model: Based on the sea area where the ice ridge is located, establish a discrete element model of sea ice composed of particle units with corresponding bonding effects, based on their mass and size. Step 2: Based on the buoyancy method, the discrete element units of sea ice are combined to generate a discrete element model of the ice ridge; the method for establishing the discrete element model of the ice ridge is as follows: First, the number N of ice fragments is set in a given computational region. This number is also the number of seed points in the segmentation algorithm. The N seed points are arranged in a close-packed hexagonal arrangement in the computational region. The coordinate values ​​of these seed points are... , The polygons generated in this way are regular hexagons, and the initial coordinates of the seed points are set... Perform the following processing: ,in c i It is the first i The initial coordinates of each seed; x i Based on initial coordinates c i A series of coordinate points are calculated, and these points form a polygon. Divided into , , , The coordinates are obtained by participating in the calculation separately; among them, the noise term and The values ​​of follow a standard normal distribution; The disorder factor determines the degree of irregularity of the polygon; spherical elements are filled into the generated polygon, and under the action of buoyancy, the final ice ridge shape is generated, completing the modeling; Step 3: Establish the discrete element model of the hull: The entire hull surface is discretized into a series of triangular elements using mesh generation preprocessing, and a contact model is established between the discrete triangular elements of the hull and the discrete element elements of the spherical sea ice. Step 4: Based on steps 2 and 3, perform dynamic calculation and analysis on the ice load of the interaction between the discrete element unit of the spherical sea ice and the discrete triangular element unit of the ship hull, obtain the force acting on each particle, and update the velocity and coordinates of each particle at the same time. Step 5: Based on Step 4, implement GPU-based parallel programming on the CUDA platform to perform hull-ice ridge coupling calculation and analysis, and realize ice load analysis and calculation.

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 ridge morphology is based on statistical analysis of the morphology of common ice ridges, including sail height, freezing layer, keel and flat ice from top to bottom.

3. 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 bonding action in step one adopts a particle bonding model, which takes into account the freezing effect between units. The two spheres are bonded together in parallel, which can not only transmit force but also torque. In the parallel bonding model, an elastic bonding disk is set between the two bonding particle units. The disk transmits the force and torque between the two units.

4. 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 method for determining the contact type between the discrete triangular element of the hull and the discrete element of the spherical sea ice is as follows: First, it is necessary to determine whether the spherical particle overlaps with the plane of the triangular element. Then, if there is an overlap between the spherical particle and the plane containing the triangular unit, it is necessary to further determine whether the spherical particle is in contact with the triangular unit itself. There are three types of contact between the spherical particle and the triangular unit: sphere-triangle surface contact, sphere-triangle edge contact, and sphere-triangle vertex contact, hereinafter referred to as surface contact, edge contact, and angular contact. Surface contact means that the projection point Q of the center of the spherical particle is inside the triangular surface; edge contact means that the spherical particle is in contact with the three edges of the triangular unit; angular contact means that the spherical particle is in contact with the three vertices ABC of the triangular unit. The centroid method is used to determine whether the projection point of the center of the sphere is inside the triangular unit. Finally, determine whether the spherical particle comes into contact with the three sides of the triangular unit, and whether the spherical particle comes into angular contact with the triangular boundary unit.

5. The method for calculating ice loads under the action of ships and ice ridges based on discrete element method according to claim 4, characterized in that, The contact force calculation method in step four, the dynamic calculation analysis of ice load, is as follows: Calculate the embedding amount of the spherical element in the boundary element. , In the formula, point P C The contact point between the spherical particle and the boundary element is represented by: The normal vector between the contact point of the particle and the triangular element is expressed as: relative displacement Represented as: in, Let P be the velocity of the particle. Contact point on the triangular unit P C speed, For displacement time; the relative displacement along the contact surface is decomposed into normal components. and tangential components : The interaction forces between particles and boundary elements are calculated using a linear contact model, where the normal force is expressed as: The tangential force is calculated incrementally and obeys the Mohr-Coulomb law: In the formula, , The normal and tangential contact stiffness between particles and boundary elements. This represents the maximum friction coefficient between the particle and the boundary element.

6. A method for calculating ice load under dynamic interaction between a ship's hull and an ice ridge, characterized in that, Using the ice load calculation method based on discrete element method described in claim 1 for ship-ice ridge interaction, the ship's ice resistance is calculated by inputting real-time bollard thrusts of different values.

7. A computer device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that, When the processor executes the computer program, it implements the method of any one of claims 1 to 5.

8. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the method of any one of claims 1 to 5.

9. A computer program product, comprising a computer program, characterized in that, When the computer program is executed by a processor, it implements the method of any one of claims 1 to 5.

Images