High-precision numerical method capable of simulating ice breaking process and ice resistance of ship in ice area, program, equipment and storage medium

Abstract

According to the method, a sea ice model is constructed in a three-dimensional fluid calculation domain by adopting a gravity falling method, a ship model is dispersed into triangular units, calculation of a parallel bonding model between sea ice particles, calculation of ice loads between the sea ice particles and the triangular units and calculation of fluid action on the sea ice model and the ship model are considered, and ship icebreaking process simulation is executed. According to the invention, randomly arranged discrete units are adopted to simulate a sea ice structure, mesoscopic model parameters are calibrated according to mechanical characteristics of on-site sea ice, and isotropic mechanical characteristics of the sea ice in a horizontal plane are realized; a discrete element method is adopted to simulate numerical sea ice, the discreteness of the sea ice is simulated on the microscale, and the mechanical property of the sea ice is represented by presenting the crushing characteristic of the sea ice on the macroscopic scale; the GPU parallel acceleration technology is adopted for numerical simulation, the calculation cost is greatly reduced, and high-precision numerical simulation of the working conditions of icebreaking, rotation and tail icebreaking of a hectometer-scale ship can be achieved.

Description

Technical Field

[0001] The present invention belongs to the technical field of numerical simulation of ship icebreaking process, and in particular relates to a high-precision numerical method, program, device and storage medium capable of simulating the ship icebreaking process and ice resistance in ice areas. Background Art

[0002] The polar regions have extreme environments, and numerical simulation of sea ice is a key technical support for polar development. However, a series of difficulties, such as difficulty in field testing in the polar regions, have greatly restricted the development and utilization of the polar regions.

[0003] Computer numerical simulation is a comprehensive application technology with significant value in teaching, scientific research, design, production, management, and decision-making. Penetration and explosion tests are extremely expensive and dangerous, but numerical simulation offers significant economic benefits and accelerates the progress of theoretical and experimental research.

[0004] In the field of numerical sea ice simulation, it is common to treat continuous sea ice as regularly arranged discrete units. This approach simulates the mechanical properties of sea ice and addresses its complex behaviors, such as dynamic processes, fracture, and accumulation. Real physical sea ice is not a homogeneous medium, possessing a complex multi-scale structure. Its internal ice crystal structure and porosity exhibit diverse and random characteristics, resulting in complex mechanical properties. Therefore, traditional continuum medium methods may struggle to address the random fracture and collision of sea ice under stress. Discrete element methods can simulate these processes in detail, while traditional continuum methods may only provide averaged results. Summary of the Invention

[0005] The purpose of the present invention is to provide a high-precision numerical method, program, device and storage medium that can simulate the icebreaking process and ice resistance of ships in ice areas.

[0006] A high-precision numerical method for simulating the icebreaking process and ice resistance of ships in ice areas includes the following steps:

[0007] Construct a three-dimensional fluid calculation domain, and set the waterplane and ship model in the three-dimensional fluid calculation domain;

[0008] The sea ice model is constructed in the computational domain using the gravity drop method. The sea ice model includes spherical sea ice particles of varying diameters. The sea ice model is immersed to a certain depth below the waterline, and the microscopic parameters of the sea ice model are calibrated.

[0009] The ship model is discretized into triangular units, the time step is set, and the ship icebreaking process simulation begins. Initially, there is no collision between the ship model and the sea ice model. When a triangular unit comes into contact with the sea ice particles, the velocity vector of the ship model and the position of each triangular unit are obtained. The parallel bonding model between the sea ice particles, the ice load between the sea ice particles and the triangular unit, and the fluid action on the sea ice model and the ship model are calculated.

[0010] For sea ice particles, the parallel bonding force, ice load force, fluid buoyancy and self-weight are combined into a resultant force, and the parallel bonding moment and ice load moment are combined into a resultant moment. Based on the resultant force and moment, the velocity and position of each sea ice particle are updated.

[0011] For the ship model, the reverse force of the ice load force on all sea ice particles, the fluid buoyancy force on all triangular elements, and the fluid drag force on the ship model are combined into the resistance force on the ship model. The reverse moment of the ice load moment on all sea ice particles and the fluid drag moment on the ship model are combined into the resistance moment on the ship model.

[0012] Repeat the above process until the simulation of the ship's icebreaking process is completed.

[0013] Furthermore, the gravity drop method is used to construct a sea ice model in the computational domain. Specifically, a random function is used to generate multiple groups of spherical sea ice particles and hollow troughs of varying diameters. All of the spherical sea ice particles are released simultaneously at a certain height. The spherical sea ice particles accumulate to form a sea ice model that is the same size as the hollow troughs and arranged randomly and irregularly.

[0014] Get the mass m of each sea ice particle i , moment of inertia I i , initial position x i (0); i = 1, 2, ..., N1, N1 is the total number of sea ice particles; calibrate the microscopic parameters of the sea ice model, including the normal contact stiffness k between sea ice and sea ice n , tangential contact stiffness k between sea ice and sea ice s , the normal bonding strength between sea ice and sea ice σ n , tangential bonding strength between sea ice and sea ice σ s , normal contact stiffness k between sea ice and ship wn , tangential contact stiffness k between sea ice and ship ws , the maximum friction coefficient between sea ice and ship μ w , the mesoscopic parameters of all sea ice particles are the same and remain unchanged during the simulated ship icebreaking process.

[0015] Furthermore, the water plane position and the draft of the ship model remain unchanged during the simulation of ship breaking ice, and the length L of the ship model and the volume V of the immersed fluid of the ship model are obtained. sub , the density of the fluid ρ w , velocity vector V w , drag coefficient C α , drag torque coefficient C β ; Discretize the ship model into triangular units, N2 is the total number of triangular units, and the speed of all triangular units is the same as the movement speed of the ship model;

[0016] Initially, there is no collision between the ship model and the sea ice model, and the velocity vector V of all sea ice particles is i (0) are all zero vectors, and the velocity vector V of the ship model is obtained. s (0) and the angular velocity vector ω s (0); When there is a triangle element in contact with the sea ice particles, d = 1, and the ice load calculation between the ship model and the sea ice model begins.

[0017] Furthermore, the parallel bonding model calculation between the sea ice particles is specifically as follows:

[0018] Step 1.1: For the i-th sea ice particle, obtain the index of the sea ice particles that it contacts and construct the set A i (d);

[0019] Step 1.2: For set A i In (d), the parallel bonding model is used for the mth sea ice particle. An elastic bonding disk is set between the two bonded sea ice particles to transmit force and torque. The i-th sea ice particle is subjected to the parallel bonding force F from the m-th sea ice particle. im (d) is:

[0020] F im (d) = F imn (d)+F ims (d)

[0021] Where m∈A i (d), F imn (d) is F im (d) The normal component, F ims (d) is F im (d) the tangential component;

[0022] F imn (d) = F imn (d-1)+ΔF imn (d), F ims (d) = F ims (d-1)+ΔF ims (d)

[0023] ΔF imn (d)=(-k n A im ΔU imn (d))n im , ΔF ims (d) = -k s A im ΔU im (d)

[0024] Among them, A im is the cross-sectional area of ​​the elastically bonded disk between the ith and mth sea ice particles, R i is the radius of the i-th sea ice particle; n im is the vector from the center of the ith sea ice particle to the center of the mth sea ice particle, n im =n imi -n imm , n imi is the vector from the center of the ith sea ice particle to its contact point with the mth sea ice particle, n imm is the vector from the mth sea ice particle to its contact point with the ith sea ice particle; ΔU imn (d) = ΔU im (d)n im , ΔU im (d) = V im (d-1)Δt; V im (d-1) is the relative velocity between the ith sea ice particle and the mth sea ice particle in the d-1th calculation, V im (d-1)=V i (d-1)+ω i (d-1)×n imi -V m (d-1)-ω m (d-1)×n imm , V i (d-1) and ω i (d-1) is the velocity vector and angular velocity vector of the i-th sea ice particle in the d-1-th calculation;

[0025] The i-th sea ice particle is subjected to the parallel bonding moment M from the m-th sea ice particle. im (d) is:

[0026] M im (d) = M imn (d)+M ims (d)

[0027] Among them, M imn (d) is Mim (d) The normal component, M ims (d) is M im (d) the tangential component;

[0028] M imn (d) = M imn (d-1)+ΔM imn (d), M ims (d) = M ims (d-1)+ΔM ims (d)

[0029] ΔM imn (d)=(-k s J im Δθ imn (d))n im , ΔM ims (d) = -k n I im Δθ im (d)

[0030] Among them, J im is the polar moment of inertia of the elastically bonded disk between the ith and mth sea ice particles, I im is the moment of inertia of the elastically bonded disk between the ith and mth sea ice particles, Δθ imn (d) = Δθ im (d)n im , Δθ im (d)=(ω i (d-1)-ω m (d-1))Δt;

[0031] Step 1.3: Calculate the maximum tensile strength σ of the i-th sea ice particle from the m-th sea ice particle immax (d) and maximum shear strength τ immax (d);

[0032]

[0033]

[0034] Step 1.4: If σ immax (d)<σ n And τ immax (d)<σ s , it is determined that the i-th sea ice particle and the m-th sea ice particle are still in a bonded state and no fracture has occurred; otherwise, it is determined that the i-th sea ice particle and the m-th sea ice particle are fractured;

[0035] Step 1.5: Traverse Ai For all sea ice particles in (d), obtain the bonding state between the i-th sea ice particle and other sea ice particles through steps 1.2 to 1.4, obtain the indexes of all sea ice particles in a bonding state with the i-th sea ice particle, and construct the set B i (d);

[0036] Step 1.6: Aggregate the parallel bonding forces F of all the sea ice particles that are bonded to the i-th sea ice particle. im (d) and moment M im (d) The parallel bonding force on the i-th sea ice particle is obtained and torque

[0037]

[0038] Furthermore, the ice load calculation between the sea ice particles and the triangular element is specifically as follows:

[0039] Step 2.1: For the i-th sea ice particle, obtain the index of the triangular unit that contacts it and construct the set C i (d);

[0040] The method for determining whether the i-th sea ice particle is in contact with the j-th triangular unit is as follows:

[0041] Step 2.1.1: Determine the projection point of the center of the ith sea ice particle on the plane where the jth triangular unit is located. If the distance from the center of the ith sea ice particle to the projection point is less than its radius R i , then it is determined that the i-th sea ice particle overlaps with the plane where the j-th triangular unit is located, and step 2.1.2 is executed; otherwise, it is determined that the i-th sea ice particle does not contact the j-th triangular unit;

[0042] Step 2.1.2: Get the three vertex coordinates h of the jth triangle unit j1 (d), h j2 (d), h j3 (d) Obtain the projection coordinates P of the center of the sphere of the i-th sea ice particle on the plane where the j-th triangular unit is located. ij (d);

[0043] Let a = h j2 (d)-h j1 (d), b=h j3 (d)-h j1 (d), c=P ij (d)-h j1 (d)

[0044] If u≥0, v≥0 and u+v≤1 are satisfied, then the projection point of the i-th sea ice particle on the plane where the j-th triangular unit is located is inside the j-th triangular unit, and the i-th sea ice particle is in contact with the j-th triangular unit; otherwise, it is determined that the i-th sea ice particle is not in contact with the j-th triangular unit;

[0045] Step 2.2: For set C i In (d), determine the contact point P between the jth triangular element and the i-th sea ice particle. Cij (d) Calculate the embedding amount ΔL ij (d);

[0046] ΔL ij (d)=|P Cij (d)P ij (d)|-R i

[0047] Step 2.3: Calculate the relative displacement Δx between the i-th sea ice particle and the j-th triangular element ij (d), and Δx ij (d) Decompose into the normal component Δx ijn (d) and the tangential component Δx ijs (d);

[0048] Δx ij (d)=(V i (d-1)-V s (d-1))Δt

[0049] Δx ijn (d)=(n wij (d)Δx ij (d))n wij (d), Δx ijs (d) = Δx ij (d)-Δx ijn (d)

[0050] Among them, V s (d-1) is the velocity vector of the ship model in the d-1th calculation, that is, the velocity vector of all triangular elements;

[0051] Step 2.4: Calculate the ice load force F of the jth triangular element on the i-th sea ice particle fij (d);

[0052] F fij (d)=[F fijn (d),F fijs (d)]

[0053] F fijn (d) = k wn ΔLij (d)n wij (d)

[0054]

[0055] Step 2.5: Traverse C i For all triangular elements in (d), obtain the ice load force F between the i-th sea ice particle and each triangular element through steps 2.2 to 2.4. fij (d) The ice load force F acting on the i-th sea ice particle from all the triangular elements in contact with it fij (d) The ice load force on the i-th sea ice particle is obtained

[0056]

[0057] Calculate the ice load moment on the i-th sea ice particle

[0058]

[0059] Furthermore, the calculation of the sea ice model and the ship model under the action of fluid is specifically as follows:

[0060] Get the immersed volume V of each triangular element subsj , calculate the buoyancy F on each triangular element bsj ;

[0061] F bsj =-ρ w gV subsj

[0062] Where g is the gravitational acceleration vector; the immersed volume V of the triangular element subsj The calculation method is:

[0063] The intersection of the normal vector of the center of the triangular unit and the waterline plane is taken as the fourth vertex; if the j-th triangular unit is completely immersed in the fluid, the fourth vertex and the lines connecting the vertices of the triangular unit form a tetrahedron, and the volume of the tetrahedron is taken as the immersed volume V of the triangular unit. subsj If the jth triangular unit is partially immersed in the fluid, then connect the fourth vertex with the vertex of the triangular unit immersed in the fluid and the intersection of the triangular unit and the waterline to form a tetrahedron, and take the volume of the tetrahedron as the immersed volume V of the triangular unit. subsj ;

[0064] Get the submerged volume V of each sea ice particle subi (d) Calculate the buoyancy F on each sea ice particle bi (d);

[0065] Fbi (d)=-ρ w gV subi (d)

[0066] Calculate the drag force F on the ship model under the action of fluid e (d) and moment M e (d);

[0067] F e (d) = -C α ρ w V sub (V ω -V s (d))

[0068] M e (d) = -C β ρ w V sub L 2 ω s (d).

[0069] Furthermore, the resultant force F acting on the sea ice particles is Zi (d) and moment M Zi The calculation method of (d) is:

[0070]

[0071]

[0072] The resistance F of the ship model Zs (d) and moment M Zs The calculation method of (d) is:

[0073]

[0074]

[0075] The velocity vector V of each sea ice particle is updated i (d), angular velocity vector ω i (d) and position vector X i (d) Specifically:

[0076]

[0077]

[0078] X i (d) = X i (d-1)+V i (d)·Δt.

[0079] A computer device / equipment / system comprising a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the steps of the above-mentioned high-precision numerical method for simulating the icebreaking process and ice resistance of ships in ice areas.

[0080] A computer-readable storage medium stores a computer program / instruction thereon, which, when executed by a processor, implements the steps of the above-mentioned high-precision numerical method for simulating the icebreaking process and ice resistance of ships in ice areas.

[0081] A computer program product includes a computer program / instructions, which, when executed by a processor, implements the steps of the above-mentioned high-precision numerical method for simulating the icebreaking process and ice resistance of ships in ice areas.

[0082] The beneficial effects of the present invention are:

[0083] The present invention simulates sea ice structure using randomly arranged discrete elements. The microscopic model parameters are calibrated based on the mechanical properties of the in-situ sea ice, achieving isotropic mechanical properties of sea ice in the horizontal plane. The discrete element method is used to numerically simulate sea ice, simulating the discreteness of sea ice at the microscopic scale and characterizing its mechanical properties by presenting its breakup characteristics at the macroscopic scale. GPU parallel acceleration technology is used for numerical simulation, significantly reducing computational costs and enabling high-precision numerical simulation of icebreaking, turning, and stern-breaking conditions for ships operating at a scale of 100 meters. This invention is primarily used to address large-scale simulations of sea ice breakage and accumulation processes in polar engineering, as well as ice-structure interactions (e.g., ship navigation). BRIEF DESCRIPTION OF THE DRAWINGS

[0084] Figure 1 Schematic diagram of the sea ice model constructed using the gravity drop method.

[0085] Figure 2 Schematic diagram of the loading model for uniaxial compression and three-point bending tests.

[0086] Figure 3 Schematic diagram of the parallel bonding model between spherical sea ice particles.

[0087] Figure 4 Schematic diagram of the projection of the center of a spherical sea ice particle onto the surface where the triangular unit is located.

[0088] Figure 5 Schematic diagram of the contact type between spherical sea ice particles and triangular elements.

[0089] Figure 6 A schematic diagram for determining whether a projection point is inside a triangle unit.

[0090] Figure 7 Schematic diagram of the edge contact between spherical sea ice particles and triangular units.

[0091] Figure 8 Schematic diagram of the contact between spherical sea ice particles and the vertices of triangular units.

[0092] Figure 9 Schematic diagram of the calculation of the submerged volume when the triangular element is completely immersed in the fluid.

[0093] Figure 10 Schematic diagram of the calculation of the submerged volume when the triangular element is partially immersed in the fluid.

[0094] Figure 11 Schematic diagram of searching between spherical sea ice particle units using GPU parallel acceleration technology.

[0095] Figure 12 This is a flow chart of the overall architecture of the present invention.

[0096] Figure 13 Schematic diagram of a uniaxial compression sea ice sample and a three-point bending sea ice sample generated by the gravity drop method in an embodiment of the present invention.

[0097] Figure 14 Schematic diagram of microscopic parameter calibration of the sea ice model in an embodiment of the present invention.

[0098] Figure 15 Schematic diagram of a hull model containing triangular units in an embodiment of the present invention.

[0099] Figure 16 Schematic diagram of a simulated ship traveling through an ice zone (viewed from above) according to an embodiment of the present invention.

[0100] Figure 17 Schematic diagram of ice resistance changing with sailing distance in an embodiment of the present invention. DETAILED DESCRIPTION

[0101] The present invention will be further described below with reference to the accompanying drawings.

[0102] A high-precision numerical method for simulating the icebreaking process and ice resistance of ships in ice areas includes the following steps:

[0103] Step 1: Construct a 3D fluid calculation domain, set the waterplane and ship model in the 3D fluid calculation domain. The position of the waterplane and the draft of the ship model remain unchanged during the simulation of the ship breaking ice. Obtain the length L of the ship model and the volume V of the ship model immersed in the fluid. sub , the density of the fluid ρ w , velocity vector V w , drag coefficient C α , drag torque coefficient C β ;

[0104] Step 2: Use the gravity drop method to construct a sea ice model in the computational domain. The sea ice model includes spherical sea ice particles of different diameters. The sea ice model is immersed in the water plane to a certain depth, and the mass m of each sea ice particle is obtained. i , moment of inertia I i , initial position x i (0); i = 1, 2, ..., N1, N1 is the total number of sea ice particles;

[0105] Step 3: Calibrate the microscopic parameters of the sea ice model, including the Poisson's ratio λ of sea ice and the density ρ of sea ice s , the tangential modulus G of sea ice, the normal contact stiffness k between sea ice and sea ice n , tangential contact stiffness k between sea ice and sea ice s , the normal bonding strength between sea ice and sea ice σ n , tangential bonding strength between sea ice and sea ice σ s , normal contact stiffness k between sea ice and ship wn , tangential contact stiffness k between sea ice and ship ws , the maximum friction coefficient between sea ice and ship μ w , the mesoscopic parameters of all sea ice particles are the same and remain unchanged during the simulated ship icebreaking process;

[0106] Step 4: Discretize the ship model into triangular units, N2 is the total number of triangular units, and the speed of all triangular units is the same as the movement speed of the ship model; set the time step Δt, R min is the minimum radius of the sea ice particle; the ship breaking process simulation starts, and initially the ship model and the sea ice model do not collide, and the velocity vector V of all sea ice particles i (0) are all zero vectors; when there is a triangle unit in contact with the sea ice particles, the velocity vector V of the ship model is obtained. s (0), d = 1;

[0107] Step 5: Get the velocity vector V of the ship model s (d) and the angular velocity vector ω s (d) The position of each triangular element, using GPU parallel technology to simultaneously process the ice load calculation between the ship model and the sea ice model, including the parallel bonding model calculation between sea ice particles, the ice load calculation between sea ice particles and triangular elements, and the calculation of the sea ice model and the ship model under the action of fluid;

[0108] Step 5.1: Calculation of the parallel bonding model between sea ice particles;

[0109] Step 5.1.1: For the i-th sea ice particle, obtain the index of the sea ice particles that it contacts and construct the set A i(d);

[0110] Step 5.1.2: For set A i In (d), the parallel bonding model is used for the mth sea ice particle. An elastic bonding disk is set between the two bonded sea ice particles to transmit force and torque. The i-th sea ice particle is subjected to the parallel bonding force F from the m-th sea ice particle. im (d) is:

[0111] F im (d) = F imn (d)+F ims (d)

[0112] Where m∈A i (d), F imn (d) is F im (d) The normal component, F ims (d) is F im (d) the tangential component;

[0113] F imn (d) = F imn (d-1)+ΔF imn (d), F ims (d) = F ims (d-1)+ΔF ims (d)

[0114] ΔF imn (d)=(-k n A im ΔU imn (d))n im , ΔF ims (d) = -k s A im ΔU im (d)

[0115] Among them, A im is the cross-sectional area of ​​the elastically bonded disk between the ith and mth sea ice particles, R i is the radius of the i-th sea ice particle; n im is the vector from the center of the ith sea ice particle to the center of the mth sea ice particle, n im =n imi -n imm , n imi is the vector from the center of the ith sea ice particle to its contact point with the mth sea ice particle, n imm is the vector from the mth sea ice particle to its contact point with the ith sea ice particle; ΔU imn (d) = ΔU im (d)nim , ΔU im (d) = V im (d-1)Δt; V im (d-1) is the relative velocity between the ith sea ice particle and the mth sea ice particle in the d-1th calculation, V im (d-1)=V i (d-1)+ω i (d-1)×n imi -V m (d-1)-ω m (d-1)×n imm , V i (d-1) and ω i (d-1) is the velocity vector and angular velocity vector of the i-th sea ice particle in the d-1-th calculation;

[0116] The i-th sea ice particle is subjected to the parallel bonding moment M from the m-th sea ice particle. im (d) is:

[0117] M im (d) = M imn (d)+M ims (d)

[0118] Among them, M imn (d) is M im (d) The normal component, M ims (d) is M im (d) the tangential component;

[0119] M imn (d) = M imn (d-1)+ΔM imn (d), M ims (d) = M ims (d-1)+ΔM ims (d)

[0120] ΔM imn (d)=(-k s J im Δθ imn (d))n im , ΔM ims (d) = -k n I im Δθ im (d)

[0121] Among them, J im is the polar moment of inertia of the elastically bonded disk between the ith and mth sea ice particles, I im is the moment of inertia of the elastically bonded disk between the ith and mth sea ice particles, Δθ imn (d) = Δθ im (d)n im , Δθ im (d)=(ω i (d-1)-ω m (d-1))Δt;

[0122] Step 5.1.3: Calculate the maximum tensile strength σ of the i-th sea ice particle from the m-th sea ice particle immax (d) and maximum shear strength τ immax (d);

[0123]

[0124]

[0125] Step 5.1.4: If σ immax (d)<σ n And τ immax (d)<σ s , it is determined that the i-th sea ice particle and the m-th sea ice particle are still in a bonded state and no fracture has occurred; otherwise, it is determined that the i-th sea ice particle and the m-th sea ice particle are fractured;

[0126] Step 5.1.5: Traverse A i For all sea ice particles in (d), obtain the bonding state between the i-th sea ice particle and other sea ice particles through steps 5.1.2 to 5.1.4, obtain the indexes of all sea ice particles that are in a bonding state with the i-th sea ice particle, and construct the set B i (d);

[0127] Step 5.1.6: Aggregate the parallel bonding forces F of all the sea ice particles that are bonded to the i-th sea ice particle. im (d) and moment M im (d) The parallel bonding force on the i-th sea ice particle is obtained and torque

[0128]

[0129] Step 5.2: Calculation of ice loads on sea ice particles and triangular elements;

[0130] Step 5.2.1: For the i-th sea ice particle, obtain the index of the triangular unit that contacts it and construct the set C i (d);

[0131] Step 5.2.1.1: Determine whether the i-th sea ice particle overlaps with the plane where the j-th triangular unit is located. If so, proceed to step 5.2.1.2; otherwise, determine that the i-th sea ice particle does not contact the j-th triangular unit.

[0132] Determine the projection point of the center of the ith sea ice particle on the plane where the jth triangular unit is located. If the distance from the center of the ith sea ice particle to the projection point is less than its radius R i , then it is determined that the i-th sea ice particle and the plane where the j-th triangular unit is located have an overlap;

[0133] Step 5.2.1.2: Use the centroid method to determine the projection point of the sphere center of the i-th sea ice particle on the plane where the j-th triangular unit is located. If so, it is inside the j-th triangular unit. If so, the i-th sea ice particle is determined to be in contact with the j-th triangular unit. Otherwise, the i-th sea ice particle is determined to be not in contact with the j-th triangular unit.

[0134] Get the three vertex coordinates h of the jth triangle unit j1 (d), h j2 (d), h j3 (d) Obtain the projection coordinates P of the center of the sphere of the i-th sea ice particle on the plane where the j-th triangular unit is located. ij (d) Let a = h j2 (d)-h j1 (d), b=h j3 (d)-h j1 (d), c=P ij (d)-h j1 (d) If u≥0, v≥0 and u+v≤1 are satisfied, then the projection point of the i-th sea ice particle on the plane where the j-th triangular unit is located is determined to be inside the j-th triangular unit;

[0135] Step 5.2.2: For set C i In (d), determine the contact point P between the jth triangular element and the i-th sea ice particle. Cij (d) Calculate the embedding amount ΔL ij (d);

[0136] ΔL ij (d)=|P Cij (d)P ij (d)|-R i

[0137] Step 5.2.3: Calculate the relative displacement Δx between the i-th sea ice particle and the j-th triangular element ij (d), and Δx ij (d) Decompose into the normal component Δx ijn(d) and the tangential component Δx ijs (d);

[0138] Δx ij (d)=(V i (d-1)-V s (d-1))Δt

[0139] Δx ijn (d)=(n wij (d)Δx ij (d))n wij (d), Δx ijs (d) = Δx ij (d)-Δx ijn (d)

[0140] Among them, V s (d-1) is the velocity vector of the ship model in the d-1th calculation, that is, the velocity vector of all triangular elements;

[0141] Step 5.2.4: Calculate the ice load force F of the jth triangular element on the i-th sea ice particle fij (d);

[0142] F fij (d)=[F fijn (d),F fijs (d)]

[0143] F fijn (d) = k wn ΔL ij (d)n wij (d)

[0144]

[0145] Step 5.2.5: Traverse C i For all triangular elements in (d), obtain the ice load force F between the i-th sea ice particle and each triangular element through steps 5.2.2 to 5.2.4. fij (d) The ice load force F acting on the i-th sea ice particle from all the triangular elements in contact with it fij (d) The ice load force on the i-th sea ice particle is obtained

[0146]

[0147] Step 5.2.6: Calculate the ice load moment on the i-th sea ice particle

[0148]

[0149] Step 5.3: Calculation of the sea ice model and ship model under the action of fluid:

[0150] Step 5.3.1: Since the waterplane position and the draft of the ship model remain unchanged during the simulation of ship breaking ice, obtain the submerged volume V of each triangular element subsj , the submerged volume V of each sea ice particle subi (d);

[0151] Take the intersection of the normal vector of the face center of the triangular element and the waterline plane as the fourth vertex;

[0152] If the jth triangular unit is completely immersed in the fluid, the fourth vertex and the lines connecting the vertices of the triangular unit form a tetrahedron, and the volume of the tetrahedron is taken as the immersed volume V of the triangular unit. subsj ;

[0153] If the jth triangular unit is partially immersed in the fluid, the fourth vertex is connected with the vertex of the triangular unit immersed in the fluid and the intersection of the triangular unit and the waterline to form a tetrahedron. The volume of the tetrahedron is taken as the immersed volume V of the triangular unit. subsj ;

[0154] Step 5.3.2: Calculate the buoyancy F on each triangular element bsj , the buoyancy F of each sea ice particle bi (d);

[0155] F bsj =-ρ w gV subsj

[0156] F bi (d)=-ρ w gV subi (d)

[0157] Where g is the gravitational acceleration vector;

[0158] Step 5.3.3: Calculate the drag force F on the ship model e (d) and moment M e (d);

[0159] F e (d) = -C α ρ w V sub (V ω -V s (d))

[0160] M e (d) = -C β ρ w Vsub L 2 ω s (d)

[0161] Step 6: Calculate the net force F on each sea ice particle Zi (d) and moment M Zi (d) Calculate the resistance F of the ship model Zs (d) and moment M Zs (d);

[0162]

[0163]

[0164]

[0165]

[0166] Step 7: Update the velocity vector V of each sea ice particle i (d), angular velocity vector ω i (d) and position vector X i (d);

[0167]

[0168]

[0169] X i (d) = X i (d-1)+V i (d)·Δt

[0170] Step 8: If the simulation of the ship's icebreaking process is not completed, set d=d+1 and return to step 5.

[0171] Example 1:

[0172] like Figure 12 As shown, the present invention generates a random irregular sequence through the gravity drop method to simulate the structure of natural sea ice, and then calibrates the microscopic parameters of the generated sea ice through uniaxial compression and three-point bending to give it mechanical properties consistent with real sea ice; combined with a hull model divided into triangular units, high-precision calculation of the ship's ice load is performed; at the same time, GPU parallel technology is integrated to greatly improve computing efficiency and reduce simulation time, thereby accurately and efficiently simulating the icebreaking process of ships in ice areas and predicting the law of ice resistance.

[0173] 1. Generating Randomly Arranged Numerical Sea Ice. Due to the anisotropy of sea ice, this method uses a gravity drop method to generate a sea ice model. First, a random function, rand, is used to generate spherical sea ice particles of varying diameters and a hollow trough of the desired sample size. These particles are then dropped from a certain height and released simultaneously to fill the hollow trough. These spherical sea ice particles accumulate to form a sea ice sample of the same size as the hollow trough and arranged in a random and irregular pattern. This method generates a random and irregular arrangement of sea ice particles, closely resembling the structure of real, natural sea ice. Figure 1 This is a sea ice sample generated by the gravity drop method.

[0174] 2. Calibration of sea ice micro-parameters: After the sea ice model is generated, several key sea ice parameters need to be calibrated to determine the selection of micro-parameters for the macro-intensity of polar sea ice.

[0175] like Figure 2 As shown in Figure 2, in order to determine these mesoscopic parameters, it is generally necessary to generate specimens for uniaxial compression and three-point bending tests, and to determine the mesoscopic parameter values ​​for the macroscopic compressive strength and bending strength values ​​of typical polar sea ice by numerically simulating the uniaxial compression and three-point bending tests.

[0176] 3. Input the structure containing triangular elements and sea ice information containing microscopic parameters and calculate the ice load between the ice zone and the structure. Among them, the structure is mainly divided into triangular elements using the starccm software. The ice load calculation between the ice zone and the structure is mainly divided into three force calculation modules, namely the parallel bonding model calculation of sea ice particles, the ice load calculation between sea ice and the structure, and the calculation of the fluid action on sea ice and the structure. The sea ice numerical model uses a discrete unit parallel bonding model to calculate the dynamic action and crushing process of sea ice; the Hertz model is used to calculate the ice load between sea ice and the structure; at the same time, the influence of seawater buoyancy and ocean current drag on sea ice and structures is added to realize the ice load calculation between the ice zone and the structure.

[0177] 1. Calculation of parallel bonding model of sea ice particles. There are two common bonding models for particle units: contact bonding and parallel bonding. Contact bonding only occurs at the contact points between particles and can only transmit contact forces between particles; while parallel bonding glues two spheres together, which can transmit both force and torque. This invention adopts Figure 3 The parallel bonding model shown is used to simulate the bonding between sea ice units.

[0178] 2. The Hertz model is used to calculate ice loads between sea ice and structures. In discrete element methods, complex structures are typically divided into a collection of triangular elements, each connected by vertices. Therefore, the problem is simplified to determining the contact between spherical sea ice particles and triangular elements, and then calculating the force exerted by the sea ice on the entire structure through traversal summation.

[0179] To determine whether a triangular boundary element f contacts a spherical particle p, we first need to determine whether the spherical particle overlaps with the plane where the triangular boundary element is located. Define the coordinates of the center P of the spherical particle p as X P , the particle radius is R ball , the projection point of the sphere center P on the plane where the triangle boundary unit is located is Q, as shown Figure 4 As shown, if |PQ|<R ball , it indicates that the spherical particle overlaps with the plane of the triangular boundary element. The spherical particle and the triangular boundary element may be in contact, and further contact determination is required. Otherwise, the spherical particle and the triangular element will not be in contact. This is called the first determination.

[0180] If the spherical particle overlaps with the plane where the triangular boundary unit is located, it is necessary to further determine whether the spherical particle is in contact with the triangular unit itself. The present invention uses the centroid method to determine whether the projection point of the sphere center is inside the triangular unit, which is called the second judgment. The three vertices of the triangular unit are located on the same plane. If one of the vertices is selected, the other two vertices can be regarded as translation displacements relative to the point, such as Figure 6 As shown in the figure, if vertex A is selected as the starting point, then vertex B is equivalent to moving a distance in the AB direction, and vertex C is equivalent to moving a distance in the AC direction.

[0181] If spherical particles may come into contact with triangular elements, their contact type needs to be determined in order to correctly calculate the contact force between the particles and the triangular elements. There are three types of contact between spherical particles and triangular boundary elements: sphere-surface contact, sphere-edge contact, and sphere-vertex contact, hereinafter referred to as surface contact, edge contact, and angular contact. Surface contact means that the projection point of the center of the spherical particle is within the triangular surface; edge contact means that the spherical particle is in contact with the three edges of the triangular boundary element; angular contact means that the spherical particle is in contact with the three vertices of the triangular boundary element. The contact types are as follows: Figure 5 shown.

[0182] The condition for edge contact with the triangular boundary element is that the second judgment is not satisfied on the basis of satisfying the first judgment. If the projection Q of the particle center P on the triangle edge is inside the triangle edge and the distance from the center to the projection point is less than the particle radius R ball When, such as Figure 7As shown, it indicates that edge contact occurs, which is called the third judgment. If the distance between the particle center P and the vertices A, B, and C is less than the particle radius R ball , then angular contact occurs, such as Figure 8 shown.

[0183] 3. In the discrete element model, the structure is composed of several triangular plates. The following method is used to calculate the volume of the structure immersed in water: take point p as the vertex and make a tetrahedron with the triangle below the waterline as the base. The sum of the volumes of all tetrahedrons is the volume of the structure immersed in water, as shown in the following example: Figure 9 When calculating the volume, special consideration should be given to the triangular unit part exposed to the water surface after the structure is discretized. Figure 10 For the triangular plate shown in the figure that is not completely immersed in water, the part exposed above the water needs to be cut off and only the volume of the underwater part needs to be calculated.

[0184] 4. The time step Δt of the entire calculation cycle can be set as:

[0185]

[0186] Where R min is the minimum radius of the sea ice unit, λ is the Poisson's ratio, ρ s represents the density of sea ice, and G represents the tangential modulus of sea ice

[0187] 5. GPU parallel acceleration technology. Due to the massive computational scale, involving millions of sea ice cells and tens of thousands of structural triangles, to reduce computational cost and time, we developed discrete element numerical calculations based on GPU parallel acceleration technology. This technology rapidly generates a list of contact relationships between a sea ice cell and other sea ice cells and structural triangles, determines whether there is contact and interaction force between them, and accelerates the updating of the sea ice cell's velocity and position.

[0188] GPU parallel acceleration technology mainly accelerates the search module and force transfer module between computing units.

[0189] 1. Such as Figure 11 As shown in the figure, the search between particle units is accelerated. The computational domain is divided into grids of equal size and the grids are sequentially numbered to serve as the background for the search particle units. The particle units are reordered and numbered based on the grid, so that they have similar sequence numbers and form a neighbor list. CUDA is used to run the Thrust library for related calculations, achieving the goal of accelerating the search.

[0190] 2. Accelerated force transfer between units. Accelerated force calculation in discrete element methods is primarily reflected in the update of tangential forces. When a pair of contacting particle units maintains contact in both the previous and current time steps, the tangential force is updated using an incremental overlay method. This means that the tangential force from the previous time step is transferred to the current time step for incremental overlay.

[0191] like Figure 13 As shown, in this embodiment, the gravity drop method is first used to generate uniaxial compression sea ice samples and three-point bending sea ice samples. The uniaxial compression sample has a size of 20 cm × 20 cm × 50 cm, and the three-point bending sample has a size of 140 cm × 20 cm × 20 cm.

[0192] Numerical simulations of uniaxial compression and three-point bending were carried out, stress-strain curves were drawn, and contour maps characterizing macroscopic strength and microscopic model parameters were established. By changing the microscopic normal bonding strength and tangential normal strength ratio between particles, numerical simulations of a total of 56 uniaxial compression tests and 56 three-point bending tests were carried out at 7 microscopic normal strengths of 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, and 3.5 MPa and 8 tangential normal strength ratios of 0.2, 0.25, 0.33, 0.5, 1.0, 2.0, 3.0, and 4.0. The compression and bending strength values ​​obtained under different microscopic parameters were plotted as contour cloud maps, as shown in the figure below. Figure 14 shown.

[0193] Figure 14 The horizontal axis is the microscopic normal bonding strength between particles, and the vertical axis is the microscopic strength relationship after the tangential normal strength ratio is converted through the relationship.

[0194] When the typical polar sea ice compressive strength is 1.63MPa and the bending strength is 0.55MPa, the intersection of these two strength contour lines is found in the compressive strength and bending strength contour cloud map. The microscopic parameter value corresponding to this intersection is the calibration value. The microscopic parameter values ​​calibrated for typical polar sea ice in this example are: microscopic normal bonding strength 0.787MPa, tangential normal strength ratio 7.991, and friction coefficient 1.0. Considering the issues of calculation accuracy and efficiency, the elastic modulus is selected as 0.5GPa. These microscopic values ​​are used to set the mechanical parameters of sea ice. At the same time, the hull model is a simplified wedge-shaped hull model established by Caeses software based on a certain ship type as the parent model, and the triangular unit grid is divided by starccm software. Finally, the sea ice and hull model information is input into the program, and the numerical simulation of the icebreaking process of ships in ice areas is carried out.

[0195] In the numerical simulation, an ice sheet 200m long, 100m wide and 1m thick was built and placed in the open sea. The model ship sailed at a constant speed of 1m / s in the ice area.

[0196] In order to better observe the resistance mechanism of broken ice on the hull after the ice layer is broken and summarize the movement law of broken ice on the bottom of the ship, we output the calculation results in an upward-looking manner as follows Figure 16 shown. Figure 17 It shows how the ice resistance of the ship changes with the increase of sailing distance when the ship travels through the ice layer.

[0197] The foregoing description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Those skilled in the art will readily appreciate that various modifications and variations of the present invention are possible. Any modifications, equivalent substitutions, or improvements made within the spirit and principles of the present invention are intended to be within the scope of protection of the present invention.

Claims

1. A high-precision numerical method for simulating the icebreaking process and ice resistance of ships in ice areas, characterized by: Construct a three-dimensional fluid calculation domain, and set the waterplane and ship model in the three-dimensional fluid calculation domain; The sea ice model is constructed in the computational domain using the gravity drop method. The sea ice model includes spherical sea ice particles of varying diameters. The sea ice model is immersed to a certain depth below the waterline, and the microscopic parameters of the sea ice model are calibrated. The ship model is discretized into triangular units, the time step is set, and the ship icebreaking process simulation begins. Initially, there is no collision between the ship model and the sea ice model. When a triangular unit comes into contact with the sea ice particles, the velocity vector of the ship model and the position of each triangular unit are obtained. The parallel bonding model between the sea ice particles, the ice load between the sea ice particles and the triangular unit, and the fluid action on the sea ice model and the ship model are calculated. For sea ice particles, the parallel bonding force, ice load force, fluid buoyancy and self-weight are combined into a resultant force, and the parallel bonding moment and ice load moment are combined into a resultant moment. Based on the resultant force and moment, the velocity and position of each sea ice particle are updated. For the ship model, the reverse force of the ice load force on all sea ice particles, the fluid buoyancy force on all triangular elements, and the fluid drag force on the ship model are combined into the resistance force on the ship model. The reverse moment of the ice load moment on all sea ice particles and the fluid drag moment on the ship model are combined into the resistance moment on the ship model. Repeat the above process until the simulation of the ship's icebreaking process is completed.

2. A high-precision numerical method for simulating the icebreaking process and ice resistance of ships in ice areas according to claim 1, characterized in that: The gravity drop method is used to construct a sea ice model in the computational domain. Specifically, a random function is used to generate multiple groups of spherical sea ice particles and hollow troughs of varying diameters. All of the spherical sea ice particles are released simultaneously at a certain height. The spherical sea ice particles accumulate to form a sea ice model that is the same size as the hollow troughs and arranged randomly and irregularly. Get the mass m of each sea ice particle i , moment of inertia I i , initial position x i (0); i = 1, 2, ..., N1, N1 is the total number of sea ice particles; calibrate the microscopic parameters of the sea ice model, including the normal contact stiffness k between sea ice and sea ice n , tangential contact stiffness k between sea ice and sea ice s , the normal bonding strength between sea ice and sea ice σ n , tangential bonding strength between sea ice and sea ice σ s , normal contact stiffness k between sea ice and ship wn , tangential contact stiffness k between sea ice and ship ws , the maximum friction coefficient between sea ice and ship μ w , the mesoscopic parameters of all sea ice particles are the same and remain unchanged during the simulated ship icebreaking process.

3. The high-precision numerical method for simulating the icebreaking process and ice resistance of ships in ice areas according to claim 2 is characterized by: The water plane position and the draft of the ship model remain unchanged during the simulation of ship breaking ice. The length L of the ship model and the volume V of the immersed fluid of the ship model are obtained. sub , the density of the fluid ρ w , velocity vector V w , drag coefficient C α , drag torque coefficient C β ; Discretize the ship model into triangular units, N2 is the total number of triangular units, and the speed of all triangular units is the same as the movement speed of the ship model; Initially, there is no collision between the ship model and the sea ice model, and the velocity vector V of all sea ice particles is i (0) are all zero vectors, and the velocity vector V of the ship model is obtained. s (0) and the angular velocity vector ω s (0); When there is a triangle element in contact with the sea ice particles, d = 1, and the ice load calculation between the ship model and the sea ice model begins.

4. The high-precision numerical method for simulating the icebreaking process and ice resistance of ships in ice areas according to claim 3 is characterized by: The parallel bonding model calculation between the sea ice particles is specifically as follows: Step 1.1: For the i-th sea ice particle, obtain the index of the sea ice particles that it contacts and construct the set A i (d); Step 1.2: For set A i In (d), the parallel bonding model is used for the mth sea ice particle. An elastic bonding disk is set between the two bonded sea ice particles to transmit force and torque. The i-th sea ice particle is subjected to the parallel bonding force F from the m-th sea ice particle. im (d) is: F im (d)=F imn (d)+F ims (d) Where m∈A i (d), F imn (d) is F im (d) The normal component, F ims (d) is F im (d) the tangential component; F imn (d)=F imn (d-1)+ΔF imn (d),F ims (d)=F ims (d-1)+ΔF ims (d) ΔF imn (d)=(-k n A im ΔU imn (d))n im ,ΔF ims (d)=-k s A im ΔU im (d) Among them, A im is the cross-sectional area of ​​the elastically bonded disk between the ith and mth sea ice particles, R i is the radius of the i-th sea ice particle; n im is the vector from the center of the ith sea ice particle to the center of the mth sea ice particle, n im =n imi -n imm , n imi is the vector from the center of the ith sea ice particle to its contact point with the mth sea ice particle, n imm is the vector from the mth sea ice particle to its contact point with the ith sea ice particle; ΔU imn (d) = ΔU im (d)n im , ΔU im (d) = V im (d-1)Δt; V im (d-1) is the relative velocity between the ith sea ice particle and the mth sea ice particle in the d-1th calculation, V im (d-1)=V i (d-1)+ω i (d-1)×n imi -V m (d-1)-ω m (d-1)×n imm , V i (d-1) and ω i (d-1) is the velocity vector and angular velocity vector of the i-th sea ice particle in the d-1-th calculation; The i-th sea ice particle is subjected to the parallel bonding moment M from the m-th sea ice particle. im (d) is: M im (d)=M imn (d)+M ims (d) Among them, M imn (d) is M im (d) The normal component, M ims (d) is M im (d) the tangential component; M imn (d)=M imn (d-1)+ΔM imn (d),M ims (d)=M ims (d-1)+ΔM ims (d) ΔM imn (d)=(-k s J im Δθ imn (d))n im ,ΔM ims (d)=-k n I im Δθ im (d) Among them, J im is the polar moment of inertia of the elastically bonded disk between the ith and mth sea ice particles, I im is the moment of inertia of the elastically bonded disk between the ith and mth sea ice particles, Δθ imn (d) = Δθ im (d)n im , Δθ im (d)=(ω i (d-1)-ω m (d-1))Δt; Step 1.3: Calculate the maximum tensile strength σ of the i-th sea ice particle from the m-th sea ice particle immax (d) and maximum shear strength τ immax (d); Step 1.4: If σ immax (d)<σ n And τ immax (d)<σ s , it is determined that the i-th sea ice particle and the m-th sea ice particle are still in a bonded state and no fracture has occurred; otherwise, it is determined that the i-th sea ice particle and the m-th sea ice particle are fractured; Step 1.5: Traverse A i For all sea ice particles in (d), obtain the bonding state between the i-th sea ice particle and other sea ice particles through steps 1.2 to 1.4, obtain the indexes of all sea ice particles in a bonding state with the i-th sea ice particle, and construct the set B i (d); Step 1.6: Aggregate the parallel bonding forces F of all the sea ice particles that are bonded to the i-th sea ice particle. im (d) and moment M im (d) The parallel bonding force on the i-th sea ice particle is obtained and torque 5. The high-precision numerical method for simulating the icebreaking process and ice resistance of ships in ice areas according to claim 4 is characterized by: The ice load calculation between the sea ice particles and the triangular element is specifically as follows: Step 2.1: For the i-th sea ice particle, obtain the index of the triangular unit that contacts it and construct the set C i (d); The method for determining whether the i-th sea ice particle is in contact with the j-th triangular unit is as follows: Step 2.1.1: Determine the projection point of the center of the ith sea ice particle on the plane where the jth triangular unit is located. If the distance from the center of the ith sea ice particle to the projection point is less than its radius R i , then it is determined that the i-th sea ice particle overlaps with the plane where the j-th triangular unit is located, and step 2.1.2 is executed; otherwise, it is determined that the i-th sea ice particle does not contact the j-th triangular unit; Step 2.1.2: Get the three vertex coordinates h of the jth triangle unit j1 (d), h j2 (d), h j3 (d) Obtain the projection coordinates P of the center of the sphere of the i-th sea ice particle on the plane where the j-th triangular unit is located. ij (d); Let a = h j2 (d) - h j1 (d), b = h j3 (d) - h j1 (d), c = P ij (d) - h j1 (d), If u≥0, v≥0 and u+v≤1 are satisfied, then the projection point of the i-th sea ice particle on the plane where the j-th triangular unit is located is inside the j-th triangular unit, and the i-th sea ice particle is in contact with the j-th triangular unit; otherwise, it is determined that the i-th sea ice particle is not in contact with the j-th triangular unit; Step 2.2: For set C i In (d), determine the contact point P between the jth triangular element and the i-th sea ice particle. Cij (d) Calculate the embedding amount ΔL ij (d); ΔL ij (d)=|P Cij (d)P ij (d)|-R i Step 2.3: Calculate the relative displacement Δx between the i-th sea ice particle and the j-th triangular element ij (d), and Δx ij (d) Decompose into the normal component Δx ijn (d) and the tangential component Δx ijs (d); Δx ij (d)=(V i (d-1)-V s (d-1))Δt <h2 style=";text-align:left;direction:ltr">Δx<h2 style=";text-align:left;direction:ltr"> ijn <h2 style=";text-align:left;direction:ltr"> (d)=(n<h2 style=";text-align:left;direction:ltr"> wij <h2 style=";text-align:left;direction:ltr"> (d)Δx<h2 style=";text-align:left;direction:ltr"> ij <h2 style=";text-align:left;direction:ltr"> (d))n<h2 style=";text-align:left;direction:ltr"> wij <h2 style=";text-align:left;direction:ltr"> (d),Δx<h2 style=";text-align:left;direction:ltr"> ijs <h2 style=";text-align:left;direction:ltr"> (d)=Δx<h2 style=";text-align:left;direction:ltr"> ij <h2 style=";text-align:left;direction:ltr"> (d)-Δx<h2 style=";text-align:left;direction:ltr"> ijn <h2 style=";text-align:left;direction:ltr"> (d) Among them, V s (d-1) is the velocity vector of the ship model in the d-1th calculation, that is, the velocity vector of all triangular elements; Step 2.4: Calculate the ice load force F of the jth triangular element on the i-th sea ice particle fij (d); F fij (d)=[F fijn (d),F fijs (d)] F fijn (d)=k wn ΔL ij (d)n wij (d) Step 2.5: Traverse C i For all triangular elements in (d), obtain the ice load force F between the i-th sea ice particle and each triangular element through steps 2.2 to 2.

4. fij (d) The ice load force F acting on the i-th sea ice particle from all the triangular elements in contact with it fij (d) The ice load force on the i-th sea ice particle is obtained Calculate the ice load moment on the i-th sea ice particle 6. The high-precision numerical method for simulating the icebreaking process and ice resistance of ships in ice areas according to claim 5 is characterized by: The calculation of the sea ice model and the ship model under the action of fluid is specifically as follows: Get the immersed volume V of each triangular element subsj , calculate the buoyancy F on each triangular element bsj ; F bsj =-ρ w gV subsj Where g is the gravitational acceleration vector; the immersed volume V of the triangular element subsj The calculation method is: The intersection of the normal vector of the center of the triangular unit and the waterline plane is taken as the fourth vertex; if the j-th triangular unit is completely immersed in the fluid, the fourth vertex and the lines connecting the vertices of the triangular unit form a tetrahedron, and the volume of the tetrahedron is taken as the immersed volume V of the triangular unit. subsj If the jth triangular unit is partially immersed in the fluid, then connect the fourth vertex with the vertex of the triangular unit immersed in the fluid and the intersection of the triangular unit and the waterline to form a tetrahedron, and take the volume of the tetrahedron as the immersed volume V of the triangular unit. subsj ; Get the submerged volume V of each sea ice particle subi (d) Calculate the buoyancy F on each sea ice particle bi (d); F bi (d)=-ρ w gV subi (d) Calculate the drag force F on the ship model under the action of fluid e (d) and moment M e (d); F e (d)=-C α ρ w V sub (V ω -V s (d)) M e (d)=-C β r w V sub L 2 oh s (d)。 7. The high-precision numerical method for simulating the icebreaking process and ice resistance of ships in ice areas according to claim 6 is characterized by: The resultant force F on the sea ice particles Zi (d) and moment M Zi The calculation method of (d) is: The resistance F of the ship model Zs (d) and moment M Zs The calculation method of (d) is: The velocity vector V of each sea ice particle is updated i (d), angular velocity vector ω i (d) and position vector X i (d) Specifically: X i (d)=X i (d-1)+V i (d)·Δt。 8. A computer device / apparatus / system comprising a memory, a processor, and a computer program stored in the memory, characterized in that: The processor executes the computer program to implement the steps of the method according to any one of claims 1 to 7.

9. A computer-readable storage medium having a computer program / instruction stored thereon, characterized in that: When the computer program / instructions are executed by a processor, the steps of the method according to any one of claims 1 to 7 are implemented.

10. A computer program product comprising a computer program / instructions, characterized in that: When the computer program / instructions are executed by a processor, the steps of the method according to any one of claims 1 to 7 are implemented.

Images