A numerical simulation method for underwater ice impact in polar environment based on riemann-sph

Abstract

This invention relates to the field of numerical simulation technology for explosive fluid-structure interaction, specifically to a numerical simulation method for underwater impact icebreaking in polar environments based on the Riemann-SPH model. The method includes establishing an axisymmetric multi-medium particle model of water, explosives, detonation products, air, and ice; establishing a local normal Riemann problem for adjacent particles; constructing an equivalent interface pressure on the ice side based on the total normal traction for the interface containing intact ice particles and solving for the Riemann star region state; updating the ice layer stress and damage state based on a material model that separates the volumetric response and deviatoric stress response of the ice layer; retaining the mass, position, density, and velocity of ice particles after failure, and converting the deviatoric stress state into a volumetric response state; and using a predictive-correction time-margin to output the simulation results of the icebreaking process. This invention is used to achieve accurate load transfer at the ice-fluid interface and to continuously simulate the fragmented state of ice after damage.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation technology for explosive fluid-structure interaction, and specifically to a numerical simulation method for underwater impact icebreaking in polar environments based on Riemann-SPH. Background Technology

[0002] Underwater icebreaking processes involve multiple physical processes, including the propagation of the explosion shock wave, the expansion and contraction of the detonation bubble, water movement, ice deformation and destruction, and ice fragmentation. This is a typical high-impact fluid-structure interaction problem. Under explosive loads, the ice layer is subjected not only to instantaneous impact pressure but also to continuous loads generated by bubble movement, water flow, and changes in the free surface. Therefore, continuous numerical simulations of the explosion load propagation process, the ice damage evolution process, and the ice fragmentation process are necessary.

[0003] Existing numerical simulation methods for underwater explosive icebreaking mainly employ the finite element method, the arbitrary Lagrange-Eulerian method (ALE), and traditional fluid-structure interaction (FSI) methods. These methods calculate the explosive load and ice response by establishing discrete models of the fluid and ice regions. However, during strong impacts, large deformations, and ice crack propagation, these methods rely on computational meshes to describe material motion, which easily leads to mesh distortion, interface penetration, and complex handling of FSI boundaries. Furthermore, after ice failure, the crack and ice fragment regions are typically handled by deleting failed elements, resulting in the loss of mass and momentum information in the failed region, preventing the ice fragments from participating in subsequent water movement and bubble loading processes.

[0004] Meshless particle methods, due to their independence from fixed meshes, are suitable for large deformation and free interface problems, and are therefore increasingly applied in numerical simulations of explosion impacts and fluid-structure interactions. Among these, the Smoothed Particle Hydrodynamics (SPH) method can describe material motion processes through interparticle interactions. However, in handling interfaces between different media such as water and ice, existing SPH methods typically use pressure averaging or artificial viscosity to transfer interfacial forces, failing to fully consider the influence of deviatoric stresses within the ice layer on the normal traction of the interface. This makes it difficult to transfer solid-side loads and fluid-side pressures within a unified flux framework.

[0005] On the other hand, ice is a brittle material with obvious nonlinear deformation characteristics, and under explosive impact, it simultaneously undergoes compressive plastic deformation, tensile cracking, and damage propagation. Existing ice material models are usually difficult to describe circumferential deformation, deviatoric stress evolution, and the transition of ice fragmentation state after damage failure under axisymmetric conditions within a meshless computational framework.

[0006] Furthermore, the underwater explosion process includes a detonation stage, a bubble expansion stage, and a subsequent load propagation stage. The calculation process involves rapid changes in material state and strong nonlinear response issues. Traditional explicit time-progression methods are easily affected by abrupt changes in state, resulting in insufficient computational stability. Summary of the Invention

[0007] This invention aims to solve the above-mentioned technical problems by providing a numerical simulation method for underwater impact icebreaking in polar environments based on Riemann-SPH. This method describes the interaction between water, explosives, detonation products, air, and ice in a unified computational model, realizes accurate load transfer at the ice-fluid interface, and continuously simulates the fragmented state of ice after ice damage.

[0008] The objective of this invention is achieved through the following technical solution: A numerical simulation method for underwater impact icebreaking in polar environments based on Riemann-SPH includes: S1. Establish an axisymmetric multi-medium particle model and discretize water, explosives, detonation products, air and ice into meshless particles. S2. Determine the particle pairs between adjacent particles based on the axisymmetric multi-medium particle model, and establish a local normal Riemann problem for the particle pairs; S3. When the particle pair contains complete ice particles, construct the equivalent interface pressure on the ice side based on the total normal traction of the complete ice particles, and solve the Riemann region state by combining the fluid side pressure and the acoustic impedance on both sides to obtain the particle flux. S4. Update the stress state of ice particles based on the material model that separates the volume response and deviatoric stress response of the ice layer, and determine the damage state based on the plastic and tensile states of the ice particles. S5. When the ice particles reach the failure state, retain the mass, position, density and velocity of the ice particles, and convert the deviatoric stress state of the ice particles into the volume response state. S6. The axisymmetric multi-medium particle model is updated using a prediction-correction time-progression method, and the damage state, equivalent plastic strain, equivalent stress, pressure, velocity and position data of ice particles are output.

[0009] Furthermore, step S1 includes: Establish using radial coordinates and axial coordinates The axially symmetric coordinate system is constructed, and the meshless particle is set as a torus formed by rotating around the axis of symmetry. The mass of the meshless particle is determined according to the volume of the torus. The initial mass of the meshless particles satisfies: in, Indicates the first The mass of a meshless particle Indicates the first The density of a meshless particle Indicates the first The radial coordinates of a meshless particle. Indicates the radial discrete length. Indicates the discrete length along the axial direction; Set the position, velocity, mass, density, pressure, specific internal energy, sound velocity, smooth length, and material type for meshless particles; set the deviatoric stress, circumferential deviatoric stress, equivalent plastic strain, and damage variables for ice particles. and damage variables The initial value is set to 0; The equivalent toroidal volume is determined based on the mass and density of the meshless particles. The equivalent toroidal volume satisfies the following: in, Indicates the first The equivalent toroidal volume of a meshless particle.

[0010] Furthermore, step S2 includes: A grid bucket is used to perform a neighborhood search on gridless particles to determine neighboring particles within the support domain, and particle pairs are formed based on the neighboring particles. For particle pairs corresponding to liquids of the same material type and particle pairs consisting of complete ice particles, the midpoint state of the particle pair is obtained by linear reconstruction using MUSCL with a limiter. For particle pairs consisting of meshless particles of different material types, particle pairs corresponding to air, and particle pairs corresponding to detonation products, a local normal Riemann problem is established using the particle center state. The unit normal is determined based on the positions of the two meshless particles in the particle pair. Unit legal direction satisfy: in, Represents the first particle in a pair The position of a gridless particle Represents the first particle in a pair The position of a gridless particle; Project the velocities on both sides of the local normal Riemann problem onto the unit normal. The normal velocity on the left side is obtained. and right-side normal velocity left-side normal velocity and right-side normal velocity They respectively satisfy: in, Indicates the left-hand velocity of the local normal Riemann problem. This represents the right-hand velocity of the local normal Riemann problem.

[0011] Furthermore, step S3 includes: Obtain the material type and damage state of the particle pair, and determine whether the particle pair contains intact ice particles based on the material type and damage state; When the particle pair contains intact ice particles, according to the pressure and deviatoric stress tensor of the intact ice particles... and unit legal direction Construct the equivalent interface pressure on the ice side, which satisfies the following: in, This represents the equivalent interfacial pressure on the ice side. This represents the pressure of a complete ice particle. The deviatoric stress tensor represents the complete ice particle. Indicates unit normal transpose; The ice-side equivalent interface pressure and the fluid-side pressure are respectively used as the left and right reconstruction pressures of the local normal Riemann problem. and ; The acoustic impedances on both sides are determined based on the density and sound velocity on both sides of the local normal Riemann problem, and the acoustic impedances on both sides satisfy: in, Indicates the acoustic impedance on the left side. Indicates the acoustic impedance on the right side. Indicates the density on the left side. Indicates the density on the right side. Indicates the speed of sound on the left. Indicates the speed of sound on the right; Based on the left and right reconstruction pressure , Left Acoustic Impedance right-side acoustic impedance Left normal velocity and the normal velocity on the right side Solving for sector pressure and star sector normal velocity ,in: According to star sector pressure and star sector normal velocity The state of the Riemann region is determined, and the density flux, momentum flux, and energy flux of particle pairs are formed based on the state of the Riemann region to obtain the particle flux.

[0012] Furthermore, step S4 includes: The volumetric pressure of ice particles is determined based on their density, and the volumetric pressure satisfies the following conditions: Based on elastic modulus Compared to Poisson Determine the bulk modulus bulk modulus satisfy: Based on bulk modulus and reference density Determine the sound velocity of ice particles used in the acoustic impedance calculation on both sides. speed of sound satisfy: in, This represents the pressure of ice particles. Indicates reference pressure. This represents the density of ice particles. Indicates the reference density; Based on elastic modulus Compared to Poisson Determine the shear modulus shear modulus satisfy: .

[0013] shear modulus Used to update the deviatoric stress state of ice particles.

[0014] Furthermore, step S4 also includes: The radial strain rate, axial strain rate, shear strain rate, and circumferential strain rate are determined based on the velocity gradient of ice particles, where the circumferential strain rate satisfies: At positions near the axis of symmetry of ice particles, the circumferential strain rate is determined based on the radial velocity gradient, and the circumferential strain rate satisfies: The deviatoric stress state of ice particles is updated using a deviatoric stress update method with Jaumann objective rate, while satisfying zero trace constraints: According to radial deviatoric stress Axial deviatoric stress Circumferential stress and shear deviatoric stress Determine the second deviatoric stress invariant The second deviatoric stress invariant satisfy: According to the second deviatoric stress invariant Determine equivalent stress Equivalent stress satisfy: Based on the initial yield strength Plastic modulus and equivalent plastic strain Determine the yield strength Yield strength satisfy: When equivalent stress Greater than yield strength At that time, the deviatoric stress state was subjected to Radial regression, and cumulative equivalent plastic strain Among them, the increase in plasticity satisfy: Furthermore, step S4 also includes: Set the equivalent plastic strain for failure. and tensile strength ; Based on equivalent plastic strain Equivalent plastic strain to failure The relationship between these factors forms the plastic failure criterion, which satisfies the following: Based on tensile cut-off strength Determine the tensile cut-off pressure Tensile cut-off pressure satisfy: According to pressure With tensile cut-off pressure The relationship between the two conditions forms the tensile failure criterion, which satisfies the following: When the plastic failure criterion is met, the damage variable will be... Set to 1; When the tensile failure criterion is met, the damage variable will be... Set to 1; Based on damage variables Determine the damage state and damage variables. A value of 1 indicates that the ice particles have reached a failure state.

[0015] Furthermore, step S5 includes: When damage variables When ice particles reach a failure state, the radial deviatoric stress of the ice particles is used to characterize this. Axial deviatoric stress Circumferential stress Shear deviatoric stress And the work done by the deviatoric stress is reset to zero; The mass, position, density, pressure, and velocity of ice particles that have reached a failure state are preserved. Based on the material model that separates the volumetric response and deviatoric stress response of the ice layer, the volumetric response of the ice particles that have reached the failure state is retained, and the deviatoric stress state of the ice particles is converted into the volumetric response state based on the volumetric response. Ice particles that have reached a failure state will continue to participate as meshless particles in neighborhood search, local normal Riemann problem solving, and particle flux calculation.

[0016] Furthermore, step S6 includes: Save the axisymmetric multi-medium particle model at the current moment; Based on the axisymmetric multi-medium particle model at the current moment, the following steps are executed sequentially: state equation calculation, detonation transformation, mirror boundary processing, neighborhood search, solving the local normal Riemann problem for particle pairs, updating the stress state and damage state of ice particles, and calculating the axisymmetric geometric source term. The predicted state is obtained in half a time step. Detonation conversion includes: setting the explosive to an unreacted explosive state; when the distance from the detonation point to the corresponding meshless particle is not greater than the product of the detonation velocity and the time after detonation, converting the material state of the corresponding meshless particle into detonation products; when the absolute value of the exponent term of the JWL equation of state of the detonation products is lower than a preset threshold, converting the detonation products into an ideal gas branch. Based on the predicted state, the state equation calculation, detonation transformation, mirror boundary processing, neighborhood search, local normal Riemann problem solution for particle pairs, stress state and damage state update of ice particles, and axisymmetric geometric source term calculation are performed again to obtain the predicted state change rate. The axisymmetric multi-medium particle model for the next time step is obtained based on the axisymmetric multi-medium particle model saved at the current time step and the predicted state change rate corresponding to a time step.

[0017] Furthermore, step S6 also includes: The time step is determined based on the sound speed-velocity constraint, acceleration constraint, and density decrease rate constraint. Time step satisfy: in, This represents the sound speed-velocity constraint coefficient. Indicates the acceleration constraint coefficient. This represents the density reduction rate constraint coefficient. Indicates the minimum smooth length. Indicates the maximum speed of sound. Indicates the maximum speed. Indicates the maximum acceleration. This represents the density of meshless particles. This represents the rate of change of density of meshless particles; When the density of meshless particles decreases, the density decrease rate constraint is used to determine the time step. ; After completing the prediction-correction time advance method for one time step, particle position correction is performed on the meshless particles corresponding to water, air and detonation products, while the position of intact ice particles is not corrected. The updated axisymmetric multi-medium particle model outputs the damage state, equivalent plastic strain, equivalent stress, pressure, velocity, and position data of ice particles.

[0018] Compared with the prior art, the present invention has the following beneficial effects: 1. This invention establishes an axisymmetric multi-medium particle model that includes water, explosives, detonation products, air, and ice, and calculates the interaction process between different media based on the Riemann-SPH method. For the interface interaction between ice and fluid, the equivalent interface pressure on the ice side is constructed based on the total normal traction of ice particles, and the Riemann star region state is solved in combination with the fluid side pressure, so that the ice load state participates in the particle flux calculation.

[0019] 2. This invention establishes an ice layer material model that separates volumetric response and deviatoric stress response, updates the stress state of ice layer particles based on the J2 elastoplastic model containing circumferential components, and determines the damage state of ice layer particles according to plastic and tensile states. After ice layer particles reach the failure state, the mass, position, density, and velocity of ice layer particles are retained, and the deviatoric stress state of ice layer particles is converted into the volumetric response state, so that the failed ice particles can continue to participate in the subsequent particle calculation process.

[0020] 3. This invention employs an explosion state conversion and prediction-correction time-progression method to continuously calculate the explosive detonation process, the expansion process of detonation products, and the subsequent bubble movement process. Based on the updated particle state, it outputs ice layer damage state, equivalent plastic strain, equivalent stress, pressure, velocity, and location data to obtain information on the damaged area and ice fragmentation area during the ice-breaking process. Attached Figure Description

[0021] Figure 1 This is a flowchart illustrating the overall process of a numerical simulation method for underwater impact icebreaking in polar environments based on Riemann-SPH, according to the present invention. Figure 2 This is a flowchart illustrating the Riemann coupling process of the water-ice interface normal traction in this invention. Figure 3 This is a flowchart illustrating the axisymmetric J2 elastoplasticity, dual-criteria damage, and failure state transition of ice particles according to the present invention. Detailed Implementation

[0022] 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.

[0023] like Figures 1 to 3 The numerical simulation method for underwater impact icebreaking in polar environments based on Riemann-SPH includes the following steps: Step S1 includes: An axisymmetric coordinate system consisting of radial coordinate r and axial coordinate z is used to establish water, air, explosive, and ice regions within the (r, z) parent plane. Water, explosive, air, and ice are discretized into meshless particles. For axisymmetric calculations, each meshless particle represents a torus formed by rotating the particle around the axis of symmetry, ensuring that the particle mass in the two-dimensional parent plane corresponds to the material volume in the three-dimensional axisymmetric space.

[0024] The initial mass of the i-th unmeshable particle is expressed as: in, For the first The mass of a meshless particle For the first The initial density of a meshless particle, For the first The radial coordinates of a meshless particle from the axis of symmetry. For radial discrete length, It represents the discrete length along the axial direction.

[0025] As the density changes during the calculation, the equivalent ring volume of the i-th meshless particle is determined by the following formula: in, For the first The equivalent toroidal volume of a meshless particle. Based on the definitions of mass and equivalent toroidal volume above, the particle mass in the subsequent continuity equation, momentum equation, and energy equation maintains a correspondence with the axisymmetric three-dimensional physical volume. For Meshless particles approaching zero are treated using axial mirroring and limiting forms. The axially symmetric term.

[0026] Each meshless particle saves its location. ,speed ,quality ,density ,pressure Internal energy Speed ​​of sound Smooth length And the type of material; ice particles simultaneously preserve the deviatoric stress tensor Circumferential stress Equivalent plastic strain and damage variables Damage variables of initially intact ice particles Set to 0. The total Cauchy stress is expressed as: in, For the total Cauchy stress, For the deviatoric stress tensor, For pressure, The stress is a unit tensor. This stress decomposition allows the volumetric pressure and deviatoric stress of the ice layer to be incorporated into the subsequent interfacial coupling and elastoplastic renewal processes, respectively.

[0027] In one implementation, a smaller particle spacing is used near the explosion source and the ice layer, while the particle spacing is gradually increased at positions further away from the explosion source and the ice layer. This results in higher spatial resolution for shock wave propagation, the water-ice interface, and the ice layer destruction area, while reducing the number of particles in the far-field region.

[0028] The interaction range between particles is determined using a smooth kernel function, preferably a two-dimensional Wendland C2 kernel function. A grid bucket approach is used for neighborhood search, identifying neighboring particles and forming particle pairs within the support domain of each unmesh particle. For liquid particles and intact ice particles of the same material type, a constraint-controlled MUSCL linear reconstruction is used to obtain the midpoint state of the particle pair. For particle pairs of different materials, particle pairs corresponding to air, and particle pairs corresponding to detonation products, the particle center state is used as the input state for the local normal Riemann problem to reduce state overshoot caused by cross-material interface reconstruction.

[0029] Step S2 includes: Pressure and sound velocity were calculated separately for each type of material without meshed particles. The Tait equation of state was used for water; the JWL equation of state was used for detonation products; and the ideal gas equation of state was used for air and later rarefied detonation products. For ice, volumetric response and deviatoric stress response were calculated separately, allowing ice pressure and shear deformation to be handled in different computational stages.

[0030] The volume pressure of ice particles is expressed as: The bulk modulus K is determined by the following formula: The velocity of sound of ice particles used in Riemann acoustic impedance calculations is determined by the following formula: in, For reference pressure, This represents the current density of ice particles. For reference density, Bulk modulus For elastic modulus, It is Poisson's ratio.

[0031] The shear modulus of ice is expressed as: The sound velocity of ice particles participating in Riemann acoustic impedance calculations is determined by the bulk modulus. Determine the shear modulus. This is used for deviatoric stress updating, so that the volume response and shear response are respectively incorporated into the acoustic impedance solution and deviatoric stress updating process.

[0032] For explosives, during initialization, explosives that have not yet detonated are set to the unreacted explosive state. Based on the location of the detonation point, the detonation velocity, and the current time after detonation, it is determined whether the detonation wave has propagated to the location of the corresponding explosive particles; when the distance from the detonation point to the explosive particle is not greater than the product of the detonation velocity and the time after detonation, the explosive is converted into JWL detonation products.

[0033] The JWL detonation products continuously update their states during the expansion process. When the absolute value of the relevant exponential term in the JWL equation of state decreases below a preset threshold, the corresponding detonation products are converted to the ideal gas branch, and the corresponding state parameters are reset based on the pressure and specific internal energy at the time of conversion. This state conversion allows for continuous calculation of the early high-pressure detonation stage and the later bubble expansion and contraction stages within the same particle frame.

[0034] Step S3 includes: For any pair of particles within the support domain, the unit normal n pointing from particle i to particle j is determined based on the particle positions: in, and These represent the positions of the two ungrid-free particles in the particle pair.

[0035] Projecting the velocities on the left and right sides of the local normal Riemann problem onto the unit normal n respectively, we obtain the left normal velocity. and right-side normal velocity : in, and These represent the velocities on the left and right sides of the local normal Riemann problem, respectively. The normal projection transforms the two-dimensional axisymmetric inter-particle interaction into a one-dimensional Riemann problem along the normal of the particle pair, and the tangential velocity is preserved in the subsequent combination of star sector vector velocities.

[0036] Read the material type and damage variable D on both sides of the particle pair. When the damage variable D of the ice particle is less than the failure threshold, the ice particle is identified as an intact ice particle; when the damage variable D reaches the failure threshold, the ice particle participates in the interface calculation as if it only has a volume response.

[0037] When one side of the particle pair is a complete ice particle, the thermodynamic pressure of the complete ice particle is not directly used as the input for the Riemann solution. Instead, the equivalent interface pressure on the ice side is constructed based on the pressure of the complete ice particle, the deviatoric stress tensor, and the unit normal of the particle pair. in, For the ice-side equivalent interface pressure, The pressure of an intact ice particle, For the deviatoric stress tensor of a complete ice particle, For unit normal direction The transpose of . This relation converts the hydrostatic pressure component and the normal deviatoric stress component of the intact ice particle into the normal pressure state in the Riemann problem.

[0038] The fluid-side pressure and the ice-side equivalent interface pressure are respectively taken as the reconstruction pressures on both sides of the local normal Riemann problem. The acoustic impedances on the left and right sides are expressed as: in, and The densities are on the left and right sides, respectively. and These represent the sound speeds on the left and right sides, respectively.

[0039] The star region pressure p* is calculated based on the reconstruction pressure, acoustic impedance, and normal velocity on both sides: And calculate the star region normal velocity u*: The star region vector velocity is obtained by combining the normal component corresponding to the star region normal velocity with the tangential velocities on both sides of the particle pair, and density flux, momentum flux, and energy flux are formed based on the obtained Riemann star region state. The normal deviatoric stress of the intact ice particle enters the same Riemann star region solution process as the fluid side through the equivalent interface pressure on the ice side, thereby exchanging the fluid side pressure and the total normal traction on the ice side in the same particle flux.

[0040] like Figure 2 As shown, the process first reads the material type, damage variable D, pressure, density, sound velocity, velocity, and deviatoric stress tensor of the ice particles from adjacent particle pairs, and determines whether the particle pair contains intact ice particles with damage variable D less than the failure threshold. If the particle pair does not contain intact ice particles, thermodynamic pressure is applied to both sides; if it contains intact ice particles, the reconstructed pressure on the corresponding side is replaced with the equivalent interface pressure on the ice side, depending on whether the intact ice particles are located on the left or right side of the Riemann problem. After the branches converge, the acoustic impedance and normal velocity on both sides are calculated uniformly, and then the star region pressure, star region normal velocity, and particle flux are obtained.

[0041] Step S4 includes: calculating the radial strain rate, axial strain rate, shear strain rate, and circumferential strain rate based on the velocity gradient of the ice particles. The circumferential strain rate is expressed as: The circumferential strain rate is determined by using the radial velocity gradient at positions near the axis of symmetry of ice particles. An axisymmetric strain rate state of the ice layer is constructed from the strain rates in each direction, and a deviatoric stress trial state is obtained by updating the deviatoric stress using a method incorporating the Jaumann objective law based on the elastic increment and rotation tensor. The axisymmetric model preserves the circumferential deviatoric stress. This allows the calculation of the two-dimensional (r, z) parent plane to include the stress response in the axial direction.

[0042] According to radial deviatoric stress Axial deviatoric stress Circumferential stress and shear deviatoric stress Calculate the second deviatoric stress invariant J2: Calculate the equivalent stress based on the second deviatoric stress invariant J2: The yield strength of ice is expressed as: in, The initial yield strength, Plastic modulus, To accumulate equivalent plastic strain, when equivalent stress Greater than yield strength At that time, the deviatoric stress trial state was performed. Radial regression and cumulative equivalent plastic strain.

[0043] In the implementation of ideal plasticity, Taking 0, the plasticity increment is expressed as: Set the equivalent plastic strain for failure. and tensile strength The criterion for plastic failure is: The tensile cutting pressure is determined by the following formula: The corresponding tensile failure criterion is: When either the plastic failure criterion or the tensile failure criterion is met, the damage variable will be... Setting it to 1 causes the ice particles to enter a failure state. The plastic failure criterion is used to characterize the plastic accumulation formed during impact compression and shear processes, while the tensile failure criterion is used to characterize the tensile failure formed during reflective tension, bending, and subsequent bubble rise.

[0044] like Figure 3 As shown, for ice particles that have not yet reached the failure state, the axisymmetric strain rate, the deviatoric stress trial state with Jaumann objective law, the second deviatoric stress invariant J2, and the equivalent stress are calculated sequentially. and yield strength When the equivalent stress is greater than the yield strength, J2 radial regression is performed and the equivalent plastic strain is updated. Then, plastic failure determination and tensile failure determination are performed respectively.

[0045] Step S5 includes: when the damage variable When an ice particle reaches a failure state, the radial deviatoric stress, axial deviatoric stress, circumferential deviatoric stress, shear deviatoric stress, and corresponding deviatoric stress work of the ice particle are reset to zero, while the mass, position, density, pressure, and velocity of the ice particle are retained.

[0046] Ice particles that have reached a failure state continue to update the pressure according to the volumetric response relationship of the ice layer, and continue to participate in neighborhood search, local normal Riemann problem solving, and particle flux calculation. Thus, the ice particles change from a state of bearing solid capable of bearing both volumetric and deviatoric stress responses to a state that retains only the volumetric response, while their mass and momentum remain within the computational domain.

[0047] Under the subsequent loading of bubbles, water mounds, and jets, the ice particles in the failed state can continue to move and generate feedback to the surrounding flow field. Therefore, the cracked zone and the ice fragmentation zone will not disappear from the computational domain due to the numerical deletion of the failed particles. The initial damage of the shock wave and the subsequent bubble load propagation damage proceed continuously in the same calculation process.

[0048] Furthermore, to suppress numerical instability caused by particle aggregation in intact ice particles under tension, the total Cauchy stress of the ice particles is decomposed into principal stresses. An artificial stress tensor is constructed from the positive principal stresses based on the artificial stress coefficient, and weighted according to the power of the kernel function ratio. The weighted artificial stress term is then incorporated into the interaction between ice particle pairs. Preferably, the artificial stress coefficient is 0.33, and the power is... Take 4. The artificial stress term only acts on short-range particle interactions under tension and does not replace the ice-side equivalent interface pressure and Riemann region solution in step S3.

[0049] Step S6 includes: exist Mirror particles are placed outside the axis of symmetry. Pressure, density, and specific internal energy are mapped in an even-symmetric manner; radial velocity, gradient components containing a single radial index, and deviatoric stress are mapped in an odd-symmetric manner; axial velocity, normal stress, and circumferential stress are mapped in an even-symmetric manner, so that particles near the axis of symmetry satisfy axisymmetric boundary relations and participate in neighborhood search and particle flux calculation.

[0050] Time progression employs a predictive-corrective approach. First, the time points are saved. All particle states, based on time The particle states are sequentially processed through state equations and detonation state transitions, axis mirror boundary treatment, neighborhood search, particle state reconstruction, local normal Riemann problem solving, particle flux calculation, and ice layer elastoplasticity and damage state updates, based on half a time step. The predicted state is obtained.

[0051] After obtaining the predicted state, the state equations and detonation state transition, axis mirror boundary processing, neighborhood search, particle state reconstruction, local normal Riemann problem solving, particle flux calculation, and ice layer elastoplasticity and damage state updates are re-executed under the predicted state to obtain the rate of change corresponding to the predicted state; then, at time... Save the particle state plus a time step The corresponding predicted state change amount is obtained at time [time]. The particle states. This process allows the pressure, sound velocity, detonation state, yield state, and damage state, which have changed after the half-time step, to re-participate in the second physical update.

[0052] Meshless particles corresponding to water, air, and detonation products undergo particle position correction after calibration; intact ice particles do not undergo this particle position correction in order to maintain the geometric relationship of the intact ice body.

[0053] The time step is determined by the sound speed-velocity constraint, acceleration constraint, and density decrease rate constraint. in, For time steps, This represents the density of meshless particles. This represents the rate of change of density for meshless particles. The density decrease rate constraint is enabled only when the density of meshless particles decreases. In one implementation, , , .

[0054] After the preset simulation time is completed, output the state data of the meshless particles according to the material type. For ice particles, output position, velocity, pressure, density, and damage variables. Equivalent plastic strain And equivalent stress. Based on damage variables Identify ice particles that have reached a failure state, and determine the areas formed by adjacent failed ice particles as the damage area and the ice fragmentation area based on spatial adjacency.

[0055] Furthermore, the maximum radial damage range is statistically determined based on the spatial distribution of the failed ice layer particles, the penetration damage state is determined based on the damage distribution along the ice layer thickness, and the ice fragmentation characteristics are statistically determined based on the velocity of the failed ice layer particles; the peak value of the shock wave, the arrival time of the shock wave, and the bubble boundary are obtained based on the pressure time history and spatial distribution of the meshless particles corresponding to water and detonation products.

[0056] In the specific calculations, the water depth, ice thickness, explosion source location, explosive quantity, explosive type, material state equation parameters, and ice mechanical parameters are input according to the icebreaking conditions to be simulated. Poisson's ratio Initial yield strength The tensile breaking strength and failure equivalent plastic strain were determined. For natural ice, the above parameters were calibrated based on uniaxial compression tests, splitting tests, three-point bending tests, and literature data. The simulation results were verified using pressure time history, bubble size, ice-breaking radius, and crack image data.

[0057] The above calculations were implemented using a Fortran program and parallelized using OpenMP. Other programming languages ​​and parallel computing methods with the same computational logic can also be used. The program executes in a loop following the sequence of particle initialization, material state update, neighborhood search, local normal Riemann solution, ice constitutive update, damage state transition, prediction-correction time progression, and result output, until the preset simulation termination time is reached.

[0058] Through the above implementation methods, the axisymmetric multi-medium particle model, the total normal traction of the complete ice particle, the equivalent interface pressure on the ice side, the acoustic impedance on both sides, the Riemann star region state, the material model that separates the volume response and deviatoric stress response of the ice layer, the damage state, the volume response state after failure, and the prediction-correction time progression method are all corresponding to clear input states, calculation relationships, judgment conditions, and state update processes, and there is a continuous data transfer relationship between each step.

[0059] It should be noted that the terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit the exemplary embodiments according to this application. As used herein, the singular form is intended to include the plural form as well, unless the context clearly indicates otherwise. Furthermore, it should be understood that when the terms "comprising" and / or "including" are used in this specification, they indicate the presence of features, steps, operations, devices, components, and / or combinations thereof.

[0060] It should be noted that the terms "first," "second," etc., in the specification, claims, and accompanying drawings of this application are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data are interchangeable where appropriate so that the embodiments of this application described herein can be implemented in sequences other than those illustrated or described herein.

[0061] The present invention has been further described above with reference to specific embodiments. However, it should be understood that the specific description herein should not be construed as limiting the nature and scope of the present invention. Various modifications made to the above embodiments by those skilled in the art after reading this specification are all within the scope of protection of the present invention.

Claims

1. A numerical simulation method for underwater impact icebreaking in polar environments based on Riemann-SPH, characterized in that, include: S1. Establish an axisymmetric multi-medium particle model and discretize water, explosives, detonation products, air and ice into meshless particles. S2. Determine the particle pairs between adjacent particles based on the axisymmetric multi-medium particle model, and establish a local normal Riemann problem for the particle pairs; S3. When the particle pair contains complete ice particles, construct the equivalent interface pressure on the ice side based on the total normal traction of the complete ice particles, and solve the Riemann star region state by combining the fluid side pressure and the acoustic impedance on both sides to obtain the particle flux. S4. Update the stress state of ice particles based on the material model that separates the volume response and deviatoric stress response of the ice layer, and determine the damage state based on the plastic and tensile states of the ice particles. S5. When the ice particles reach the failure state, retain the mass, position, density and velocity of the ice particles, and convert the deviatoric stress state of the ice particles into the volume response state. S6. Update the axisymmetric multi-medium particle model using a prediction-correction time-progression method, and output the damage state, equivalent plastic strain, equivalent stress, pressure, velocity, and position data of the ice particles.

2. The numerical simulation method for underwater impact icebreaking in polar environments based on Riemann-SPH as described in claim 1, characterized in that, Step S1 includes: Establish using radial coordinates and axial coordinates An axisymmetric coordinate system is constructed, the meshless particle is set as a toroidal body formed by rotating around the axis of symmetry, and the mass of the meshless particle is determined according to the volume of the toroidal body; The initial mass of the meshless particles satisfies: in, Indicates the first The mass of the aforementioned meshless particle Indicates the first The density of the aforementioned meshless particles, Indicates the first The radial coordinates of the aforementioned meshless particles, Indicates the radial discrete length. Indicates the discrete length along the axial direction; The meshless particles are assigned position, velocity, mass, density, pressure, specific internal energy, sound speed, smooth length, and material type. The ice particles are assigned deviatoric stress, circumferential deviatoric stress, equivalent plastic strain, and damage variables. and the damage variable The initial value is set to 0; The equivalent ring volume is determined based on the mass and density of the meshless particles, and the equivalent ring volume satisfies: in, Indicates the first The equivalent ring volume of the meshless particles.

3. The numerical simulation method for underwater impact icebreaking in polar environments based on Riemann-SPH as described in claim 1, characterized in that, Step S2 includes: A neighborhood search is performed on the unmesh-free particles using a grid bucket to determine the neighboring particles within the support domain, and the particle pairs are formed based on the neighboring particles. For the particle pairs corresponding to liquids of the same material type and the particle pairs composed of the complete ice particles, the midpoint state of the particle pair is obtained by linear reconstruction using MUSCL with a limiter. For the particle pairs composed of the meshless particles of different material types, the particle pairs corresponding to the air, and the particle pairs corresponding to the detonation products, the local normal Riemann problem is established using the particle center state. The unit normal is determined based on the positions of the two meshless particles in the particle pair. The unit normal satisfy: in, Indicates the first particle in the particle pair The positions of the aforementioned meshless particles Indicates the first particle in the particle pair The positions of the aforementioned meshless particles; Project the velocities on both sides of the local normal Riemann problem onto the unit normal. The normal velocity on the left side is obtained. and right-side normal velocity The left-side normal velocity and the right-side normal velocity They respectively satisfy: in, This represents the left-hand velocity of the local normal Riemann problem. This represents the right-hand velocity of the local normal Riemann problem.

4. The numerical simulation method for underwater impact icebreaking in polar environments based on Riemann-SPH as described in claim 3, characterized in that, Step S3 includes: Obtain the material type and damage state of the particle pair, and determine whether the particle pair contains the intact ice particle based on the material type and damage state; When the particle pair contains the intact ice particle, according to the pressure and deviatoric stress tensor of the intact ice particle... and the unit normal The ice-side equivalent interface pressure is constructed, and the ice-side equivalent interface pressure satisfies: in, This represents the equivalent interfacial pressure on the ice side. This represents the pressure of the intact ice particle. This represents the deviatoric stress tensor of the complete ice particle. Indicates the unit normal transpose; The ice-side equivalent interface pressure and the fluid-side pressure are respectively used as the left and right reconstruction pressures of the local normal Riemann problem. and ; The acoustic impedances on both sides are determined based on the density and sound velocity on both sides of the local normal Riemann problem, and the acoustic impedances on both sides satisfy: in, Indicates the acoustic impedance on the left side. Indicates the acoustic impedance on the right side. Indicates the density on the left side. Indicates the density on the right side. Indicates the speed of sound on the left. Indicates the speed of sound on the right; According to the left and right reconstruction pressure , The left acoustic impedance The right-side acoustic impedance The left-side normal velocity and the right-side normal velocity Solving for sector pressure and star sector normal velocity ,in: According to the star region pressure and the normal velocity of the star region The state of the Riemann region is determined, and the density flux, momentum flux, and energy flux of the particle pair are formed based on the state of the Riemann region to obtain the particle flux.

5. The numerical simulation method for underwater impact icebreaking in polar environments based on Riemann-SPH as described in claim 1, characterized in that, Step S4 includes: The volumetric pressure of the ice particles is determined based on their density, and the volumetric pressure satisfies: Based on elastic modulus Compared to Poisson Determine the bulk modulus The bulk modulus satisfy: According to the bulk modulus and reference density Determine the sound velocity of the ice particles used in the acoustic impedance calculation on both sides. The speed of sound satisfy: in, This indicates the pressure of the ice particles. Indicates reference pressure. This represents the density of the ice particles. This represents the reference density; According to the elastic modulus and the Poisson ratio Determine the shear modulus The shear modulus satisfy: The shear modulus Used to update the deviatoric stress state of the ice particles.

6. The numerical simulation method for underwater impact icebreaking in polar environments based on Riemann-SPH as described in claim 5, characterized in that, Step S4 further includes: The radial strain rate, axial strain rate, shear strain rate, and circumferential strain rate are determined based on the velocity gradient of the ice particles, wherein the circumferential strain rate satisfies: At a position near the axis of symmetry of the ice particles, the circumferential strain rate is determined based on the radial velocity gradient, and the circumferential strain rate satisfies: The deviatoric stress state of the ice particles is updated using a deviatoric stress update method with Jaumann objective rate, while satisfying the zero-trace constraint: According to radial deviatoric stress Axial deviatoric stress Circumferential stress and shear deviatoric stress Determine the second deviatoric stress invariant The second deviatoric stress invariant satisfy: According to the second deviatoric stress invariant Determine the equivalent stress The equivalent stress satisfy: Based on the initial yield strength Plastic modulus and the equivalent plastic strain Determine the yield strength The yield strength satisfy: When the equivalent stress Greater than the yield strength At that time, the deviatoric stress state is subjected to Radial regression, and accumulation of the equivalent plastic strain. Among them, the increase in plasticity satisfy: 。 7. The numerical simulation method for underwater impact icebreaking in polar environments based on Riemann-SPH as described in claim 6, characterized in that, Step S4 further includes: Set the equivalent plastic strain for failure. and tensile strength ; According to the equivalent plastic strain Equivalent plastic strain to the failure The relationship between these factors forms the plastic failure criterion, which satisfies the following: According to the tensile cut-off strength Determine the tensile cut-off pressure The tensile cutoff pressure satisfy: According to the pressure With the tensile cut-off pressure The relationship between the two factors forms a tensile failure criterion, which satisfies the following: When the aforementioned plastic failure criterion is met, the damage variable will be... Set to 1; When the tensile failure criterion is met, the damage variable will be... Set to 1; According to the damage variable Determine the damage state, the damage variable When the value is 1, it is determined that the ice particles have reached the failure state.

8. The numerical simulation method for underwater impact icebreaking in polar environments based on Riemann-SPH as described in claim 7, characterized in that, Step S5 includes: When the damage variable When the ice particles reach the failure state, the radial deviatoric stress of the ice particles is characterized. The axial deviatoric stress The circumferential deviatoric stress The shear deviatoric stress And the work done by the deviatoric stress is reset to zero; The mass, position, density, pressure, and velocity of the ice particles that have reached the failure state are retained. Based on the material model that separates the volume response and deviatoric stress response of the ice layer, the volume response of the ice layer particles that have reached the failure state is retained, and the deviatoric stress state of the ice layer particles is converted into the volume response state based on the volume response. The ice particles that have reached the failure state will continue to participate as meshless particles in neighborhood search, local normal Riemann problem solving, and particle flux calculation.

9. The numerical simulation method for underwater impact icebreaking in polar environments based on Riemann-SPH as described in claim 1, characterized in that, Step S6 includes: Save the axisymmetric multi-medium particle model at the current moment; Based on the axisymmetric multi-medium particle model at the current moment, the following steps are executed sequentially: state equation calculation, detonation transformation, mirror boundary processing, neighborhood search, solving the local normal Riemann problem for the particle pair, updating the stress state and damage state of the ice particles, and calculating the axisymmetric geometric source term. The predicted state is obtained in half a time step. The detonation conversion includes: setting the explosive to an unreacted explosive state; when the distance from the detonation point to the corresponding meshless particle is not greater than the product of the detonation velocity and the time after detonation, converting the material state of the corresponding meshless particle to the detonation product; when the absolute value of the exponent term of the JWL equation of state of the detonation product is lower than a preset threshold, converting the detonation product to an ideal gas branch. Based on the predicted state, the state equation calculation, detonation transformation, mirror boundary processing, neighborhood search, local normal Riemann problem solution for the particle pair, stress state and damage state update of the ice particles, and axisymmetric geometric source term calculation are performed again to obtain the predicted state change rate. The axisymmetric multi-medium particle model for the next time step is obtained based on the axisymmetric multi-medium particle model saved at the current time step and the predicted state change rate corresponding to a time step.

10. A numerical simulation method for underwater impact icebreaking in polar environments based on Riemann-SPH, as described in claim 9, is characterized in that... Step S6 further includes: The time step is determined based on the sound speed-velocity constraint, acceleration constraint, and density decrease rate constraint. The time step satisfy: in, This represents the sound speed-velocity constraint coefficient. Indicates the acceleration constraint coefficient. This represents the density reduction rate constraint coefficient. Indicates the minimum smooth length. Indicates the maximum speed of sound. Indicates the maximum speed. Indicates the maximum acceleration. This represents the density of the meshless particles. This represents the rate of change of density of the meshless particles; When the density of the meshless particles decreases, the density decrease rate constraint is used to determine the time step. ; After completing the prediction-correction time advance method for one time step, particle position correction is performed on the meshless particles corresponding to the water, the air, and the detonation products, while the position of the intact ice particles is not corrected. The updated axisymmetric multi-medium particle model outputs the damage state, equivalent plastic strain, equivalent stress, pressure, velocity, and position data of the ice particles.

Images