Elastoplastic smooth particle state updating method and device

Abstract

The application relates to the field of simulation technology, in particular to an elastic-plastic smooth particle state updating method and device; a particle swarm variable array is initialized according to simulation parameters; neighborhood searching is performed on the initialized variable array according to a preset particle support domain radius to obtain adjacent particle list data; a GPU parallel computing state is determined according to the adjacent particle list data; the GPU parallel computing state is generated according to the adjacent particle list data; thread allocation is performed on each particle in the initialized variable array; the initialized variable array is calculated according to the thread allocation result to obtain state parameters; and state updating is performed on all particles in the initialized variable array according to the state parameters; the scheme reads simulation parameters, initializes a particle swarm, performs neighborhood searching on the particle support domain radius, combines a GPU state, uses a 'one particle one thread' parallel architecture, updates a particle array, realizes multi-time step simulation, considers generalization and expansion, and supports multi-field engineering simulation.

Description

Technical Field

[0001] This invention relates to the field of simulation technology, and in particular to a method and apparatus for updating the state of elastoplastic smooth particles. Background Technology

[0002] Existing smooth particle dynamics simulation techniques suffer from process defects. The generation of the particle swarm variable array relies on manual adjustment based on experience, and the initialization process lacks a standardized procedure, often resulting in uneven particle distribution. The neighborhood search stage uses a serial traversal approach, and its efficiency drops sharply as the number of particles increases. The computation process does not fully utilize the parallel advantages of GPUs, lacks an effective parallel computing state management mechanism, and thread allocation lacks dynamic optimization, leading to insufficient hardware resource utilization. State parameter calculation relies on single-threaded, particle-by-particle processing, which is too time-consuming in large-scale particle scenarios. Particle state updates are not fully automated, requiring manual intervention to adjust the iteration step size, which not only reduces simulation efficiency but also easily introduces errors due to human operation, making it difficult to meet the high-precision and high-efficiency requirements of engineering applications. Summary of the Invention

[0003] In order to overcome the shortcomings of the prior art, the purpose of this invention is to provide a method and apparatus for updating the state of elastoplastic smooth particles.

[0004] A method for updating the state of elastoplastic smooth particles includes: acquiring simulation parameters; generating a particle swarm variable array based on a preset simulation model; initializing the particle swarm variable array according to the simulation parameters to obtain an initialized variable array; performing a neighborhood search on the initialized variable array based on a preset particle support domain radius to obtain a list of neighboring particles; determining the GPU parallel computing state based on the neighboring particle list data; assigning threads to each particle in the initialized variable array under the adjusted GPU parallel computing state to obtain thread assignment results; calculating the initialized variable array based on the thread assignment results to obtain state parameters; and updating the state of all particles in the initialized variable array based on the state parameters to obtain an updated particle variable array.

[0005] Further, the step of calculating the initial variable array according to the thread allocation result to obtain state parameters includes: traversing and calculating the initial variable array according to a preset density change rate calculation formula, a preset velocity gradient calculation formula, a preset torsional rate calculation formula, and the thread allocation result to obtain density change rate correction terms, velocity gradient correction terms, and torsional rate correction terms; traversing and calculating the initial variable array according to a preset strain rate calculation formula, a preset artificial viscosity correction term calculation formula, and the thread allocation result to obtain strain rate correction terms and artificial viscosity correction terms; traversing and calculating the initial variable array according to a preset virtual stress calculation formula, a preset dynamic correction formula, and the thread allocation result to obtain virtual stress correction terms and extended smooth particle dynamics correction terms; and traversing and calculating the initial variable array according to a preset normal plane... The initialization variable array is calculated using the contact stiffness coefficient, preset tangential stiffness coefficient, preset primary contact surface, and preset secondary contact surface to obtain the normal contact force set and tangential contact force set. The initialization variable array is then calculated using artificial viscosity correction terms, virtual stress correction terms, strain rate correction terms, and torsional rate correction terms to obtain the updated particle displacement set and deviatoric stress change rate. The initialization variable array is then calculated using the deviatoric stress change rate and preset yield function to obtain the deviatoric stress tensor. Finally, state parameters are generated using density change rate correction terms, velocity gradient correction terms, torsional rate correction terms, strain rate correction terms, artificial viscosity correction terms, virtual stress correction terms, extended smooth particle dynamics correction terms, normal contact force set, tangential contact force set, updated particle displacement set, deviatoric stress change rate, and deviatoric stress tensor.

[0006] Furthermore, the step of performing a neighborhood search on the initialization variable array based on the preset particle support domain radius to obtain the neighboring particle list data includes: obtaining particle position information from the initialization variable array; analyzing the particle position information to obtain the particle distribution range; dividing the particle distribution range into a grid based on the particle support domain radius to obtain the grid division result; numbering the initialization variable array according to the grid division result to obtain the particle unit number set; sorting the particle unit number set to obtain the unit particle number list; and performing a neighborhood search on the unit particle number list based on the particle position information to obtain the neighboring particle list data.

[0007] Further, the step of traversing and calculating the initial variable array according to the preset density change rate calculation formula, the preset velocity gradient calculation formula, the preset torsion rate calculation formula, and the thread allocation result to obtain the density change rate correction term, the velocity gradient correction term, and the torsion rate correction term includes: traversing and calculating the initial variable array according to the density change rate calculation formula, the neighboring particle list data, and the thread allocation result to obtain a density change rate set; accumulating the density change rate set to obtain a density change rate correction term; traversing and calculating the initial variable array according to the velocity gradient calculation formula and the thread allocation result to obtain a velocity gradient set; accumulating the velocity gradient set to obtain a velocity gradient correction term; traversing and calculating the initial variable array according to the torsion rate calculation formula and the thread allocation result to obtain a torsion rate set; accumulating the torsion rate set to obtain a torsion rate correction term.

[0008] Further, the step of traversing and calculating the initial variable array according to the preset strain rate calculation formula, the preset artificial viscosity correction term calculation formula, and the thread allocation result to obtain the strain rate correction term and the artificial viscosity correction term includes: traversing and calculating the initial variable array according to the strain rate calculation formula and the thread allocation result to obtain a strain rate set; accumulating the strain rate set to obtain the strain rate correction term; traversing and calculating the initial variable array according to the artificial viscosity correction term calculation formula and the thread allocation result to obtain an artificial viscosity correction parameter set; and accumulating the artificial viscosity correction parameter set to obtain the artificial viscosity correction term.

[0009] Further, the step of traversing and calculating the initialization variable array according to the preset virtual stress calculation formula, the preset dynamic correction formula, and the thread allocation result to obtain the virtual stress correction term and the extended smooth particle dynamics correction term includes: obtaining the current stress state of the particle group and analyzing the current stress state of the particle group; if the current stress state of the particle group is the principal stress state, then obtaining the first principal stress set from the initialization variable array according to the current stress state of the particle group; determining whether there are particles with principal stress greater than zero in the initialization variable array according to the first principal stress set; when there are particles with principal stress greater than zero in the initialization variable array, obtaining the principal stresses corresponding to all particles with principal stress greater than zero to obtain the second principal stress set; traversing and calculating the initialization variable array according to the virtual stress calculation formula, the second principal stress set, and the thread allocation result to obtain the virtual stress set; accumulating and calculating the virtual stress set to obtain the virtual stress correction term; traversing and calculating the initialization variable array according to the dynamic correction formula and the thread allocation result to obtain the velocity correction term set; accumulating and calculating the velocity correction term set to obtain the extended smooth particle dynamics correction term.

[0010] Further, the step of calculating the set of normal contact force and the set of tangential contact force based on the preset normal contact stiffness coefficient, the preset tangential stiffness coefficient, the preset main contact surface, and the preset secondary contact surface to initialize the variable array includes: obtaining the set of secondary surface particles and the set of main surface particles from the initialized variable array based on the main contact surface and the secondary contact surface; calculating the set of secondary surface particles and the set of main surface particles according to the preset normal vector calculation formula to obtain the set of contact normal vectors; calculating the set of secondary surface particles and the set of main surface particles according to the preset distance calculation formula to obtain the set of contact distances; analyzing the set of contact distances according to the preset contact threshold to obtain the contact analysis result; obtaining the set of contact particle pairs from the initialized variable array based on the contact analysis result; and calculating the set of contact particle pairs based on the normal contact stiffness coefficient, the set of contact normal vectors, and the tangential stiffness coefficient to obtain the set of normal contact force and the set of tangential contact force.

[0011] Further, the step of calculating the initial variable array based on the artificial viscosity correction term, the virtual stress correction term, the strain rate correction term, and the torsional rate correction term to obtain the updated particle displacement set and the deviatoric stress change rate includes: calculating the initial variable array based on the preset equilibrium equation, the artificial viscosity correction term, and the virtual stress correction term to obtain the particle acceleration set; obtaining the particle displacement set from the initial variable array; updating the particle displacement set based on the particle acceleration set to obtain the updated particle displacement set; and calculating the initial variable array based on the strain rate correction term, the torsional rate correction term, and the preset elastic constitutive relation calculation formula to obtain the deviatoric stress change rate.

[0012] Furthermore, the step of calculating the initial variable array based on the deviatoric stress change rate and the preset yield function to obtain the deviatoric stress tensor includes: obtaining a material constitutive model and analyzing the material constitutive model to obtain the material constitutive model type; determining whether the material constitutive model type is an inelastic model; and when the material constitutive model type is an inelastic model, iteratively calculating the initial variable array based on the deviatoric stress change rate and the yield function to obtain the deviatoric stress tensor.

[0013] Furthermore, an elastoplastic smooth particle state update device includes: a parameter acquisition module for acquiring simulation parameters; an array generation module for generating a particle swarm variable array according to a preset simulation model; an initialization module for initializing the particle swarm variable array according to the simulation parameters to obtain an initialization variable array; a neighborhood search module for performing a neighborhood search on the initialization variable array according to a preset particle support domain radius to obtain a list of neighboring particles; a state determination module for determining the GPU parallel computing state based on the list of neighboring particles; a thread allocation module for allocating threads to each particle in the initialization variable array under the adjusted GPU parallel computing state to obtain a thread allocation result; a state parameter calculation module for calculating the initialization variable array based on the thread allocation result to obtain state parameters; and a state update module for updating the state of all particles in the initialization variable array according to the state parameters to obtain an updated particle variable array.

[0014] In the technical solution of this invention, simulation parameters are first read to provide the input basis for calculation; by generating and initializing a particle swarm variable array, the initial state of the particles is ensured to be highly matched with the simulation target; the core advantage lies in adopting a neighborhood search strategy based on the particle support domain radius, combined with GPU radix sorting technology, to reduce invalid calculations and improve the generation efficiency of the neighboring particle list data; through a "one particle, one thread" GPU parallel computing architecture, an independent thread is allocated to each particle to fully utilize hardware computing power and improve the simulation efficiency of large-scale particle systems; by updating the particle variable array, multi-time-step continuous simulation is achieved to capture the elastoplastic evolution process of materials under dynamic loads and adapt to the needs of complex mechanical behavior analysis; the overall solution takes into account both versatility and scalability, supports multi-domain simulation parameter configuration, and the modular process facilitates functional expansion, providing efficient and reliable numerical simulation support for structural strength analysis, material failure prediction, and other scenarios in engineering design. Attached Figure Description

[0015] The above and / or additional aspects and advantages of the present invention will become apparent and readily understood from the description of the embodiments taken in conjunction with the following drawings, in which:

[0016] Figure 1 This is a first flowchart of a method for updating the state of elastoplastic smooth particles provided in an embodiment of the present invention;

[0017] Figure 2 This is a second flowchart of a method for updating the state of elastoplastic smooth particles provided in an embodiment of the present invention;

[0018] Figure 3 This is a third flowchart of a method for updating the state of elastoplastic smooth particles provided in an embodiment of the present invention;

[0019] Figure 4 This is a fourth flowchart of a method for updating the state of elastoplastic smooth particles provided in an embodiment of the present invention;

[0020] Figure 5 A fifth flowchart of a method for updating the state of elastoplastic smooth particles provided in an embodiment of the present invention;

[0021] Figure 6 The sixth flowchart of a method for updating the state of elastoplastic smooth particles provided in an embodiment of the present invention;

[0022] Figure 7 The seventh flowchart of a method for updating the state of elastoplastic smooth particles provided in an embodiment of the present invention;

[0023] Figure 8 The eighth flowchart of a method for updating the state of elastoplastic smooth particles provided in an embodiment of the present invention;

[0024] Figure 9 A ninth flowchart of a method for updating the state of elastoplastic smooth particles provided in an embodiment of the present invention;

[0025] Figure 10 This is a schematic diagram of an elastoplastic smooth particle state update device provided in an embodiment of the present invention. Detailed Implementation

[0026] The terms "first," "second," "third," "fourth," etc. (if present) in the specification, claims, and accompanying drawings of this invention 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 can be interchanged where appropriate so that the embodiments described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "comprising" or "having" and any variations thereof are intended to cover a non-exclusive inclusion; for example, a process, method, system, product, or apparatus that comprises a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0027] For ease of understanding, the specific process of the embodiments of the present invention is described below. Please refer to [link / reference]. Figure 1 One embodiment of the method for updating the state of elastoplastic smooth particles according to the present invention includes:

[0028] 101. Obtain simulation parameters;

[0029] In this embodiment, simulation parameters are read from the preprocessing file. The simulation parameters include mechanical parameters such as material elastic modulus, Poisson's ratio, and density, particle smooth radius, calculation dimension, contact stiffness coefficient, friction coefficient, contact threshold, geometry, boundary conditions, material properties, and contact properties.

[0030] 102. Generate a particle swarm variable array based on the preset simulation model;

[0031] 103. Initialize the particle swarm variable array according to the simulation parameters to obtain the initial variable array;

[0032] 104. Perform a neighborhood search on the initial variable array based on the preset particle support domain radius to obtain the list of neighboring particles;

[0033] In this embodiment, neighborhood search is performed based on the particle support domain radius. This search strategy reduces invalid computation and, in conjunction with GPU radix sorting technology, improves the efficiency of generating neighboring particle list data.

[0034] 105. Determine the GPU parallel computing status based on the neighboring particle list data;

[0035] 106. Under the adjusted GPU parallel computing state, thread allocation is performed on each particle in the initialization variable array to obtain the thread allocation result;

[0036] In this embodiment, in order to achieve parallel computing on the GPU, a separate thread is allocated to each particle to fully utilize the hardware computing power and improve the simulation efficiency of large-scale particle systems (such as models containing tens of thousands to millions of particles).

[0037] 107. Calculate the initialization variable array based on the thread allocation result to obtain the state parameters;

[0038] 108. Update the state of all particles in the initial variable array according to the state parameters to obtain the updated particle variable array;

[0039] In this embodiment, the time is increased to solve the next time step, and it is determined whether the termination time has been reached. If not, the time is updated and the process proceeds to step 104 to solve the next time step. If the termination time has been reached, the calculation ends. By dynamically updating the particle variable array, multi-time step continuous simulation is achieved until the preset termination time is reached. This iterative mechanism can simulate the elastoplastic evolution process of materials under dynamic loads and is suitable for the analysis of complex mechanical behavior.

[0040] In this embodiment, simulation parameters are first read to provide the input basis for calculation. By generating and initializing a particle swarm variable array, the initial state of the particles is ensured to be highly matched with the simulation target. The core advantage lies in adopting a neighborhood search strategy based on the particle support domain radius, combined with GPU radix sorting technology, to reduce invalid calculations and improve the generation efficiency of the neighboring particle list data. Through a GPU parallel computing architecture of "one particle, one thread", an independent thread is allocated to each particle to give full play to the hardware computing power and improve the simulation efficiency of large-scale particle systems. By updating the particle variable array, continuous simulation at multiple time steps is achieved to capture the elastoplastic evolution process of materials under dynamic loads and adapt to the needs of complex mechanical behavior analysis. The overall solution takes into account both versatility and scalability, supports the configuration of simulation parameters in multiple fields, and the modular process facilitates functional expansion, providing efficient and reliable numerical simulation support for structural strength analysis, material failure prediction and other scenarios in engineering design.

[0041] Please see Figure 2 In a second embodiment of the elastoplastic smooth particle state update method of the present invention, step 107 includes:

[0042] 201. Based on the preset density change rate calculation formula, preset velocity gradient calculation formula, preset torsion rate calculation formula and thread allocation results, the initial variable array is traversed and calculated to obtain the density change rate correction term, velocity gradient correction term and torsion rate correction term.

[0043] In this embodiment, three types of basic mechanical correction terms are generated through traversal calculations to accurately capture the density evolution, velocity spatial changes, and rotational effects between particles, providing underlying data support for subsequent mechanical response analysis; combined with thread allocation to achieve parallel computing, the efficiency and accuracy of basic physical quantity calculations in large-scale particle systems are ensured, laying the foundation for simulation realism;

[0044] 202. Based on the preset strain rate calculation formula, the preset artificial viscosity correction term calculation formula, and the thread allocation results, the initial variable array is traversed and calculated to obtain the strain rate correction term and the artificial viscosity correction term.

[0045] In this embodiment, the strain rate correction term quantifies the particle deformation rate and is the core input for solving the elastoplastic constitutive relation; the artificial viscosity correction term compensates for numerical oscillations through empirical formulas, effectively suppressing the tensile instability commonly found in the Smooth Particle Dynamics (SPH) method, improving the numerical stability of dynamic simulation, and avoiding interference from non-physical phenomena.

[0046] 203. Based on the preset virtual stress calculation formula, the preset dynamic correction formula, and the thread allocation results, the initial variable array is traversed and calculated to obtain the virtual stress correction term and the extended smooth particle dynamics correction term;

[0047] In this embodiment, the virtual stress correction term optimizes the stress field distribution for the principal stress state, enhancing the calculation accuracy under complex stress scenarios; the extended smooth particle dynamics correction term stabilizes particle motion trajectories through velocity smoothing, reduces abnormal behaviors such as particle penetration, and ensures physical consistency in the simulation process.

[0048] 204. Calculate the set of normal contact force and the set of tangential contact force based on the preset normal contact stiffness coefficient, preset tangential stiffness coefficient, preset primary contact surface, and preset secondary contact surface to initialize the variable array;

[0049] In this embodiment, based on the main contact surface, the sub-contact surface division, and the stiffness coefficient, the contact force set is accurately solved, clearly depicting the normal compression and tangential friction between particles; effective contact pairs are screened by the contact threshold criterion to avoid invalid contact calculations, thus ensuring the authenticity of the contact mechanical behavior and improving computational efficiency.

[0050] 205. Calculate the initial variable array based on the artificial viscosity correction term, virtual stress correction term, strain rate correction term, and torsional rate correction term to obtain the updated particle displacement set and deviatoric stress change rate;

[0051] In this embodiment, multiple correction terms are integrated to solve for particle acceleration and displacement updates, enabling dynamic tracking of the deformation process; the deviatoric stress change rate is calculated through elastic constitutive relations to quantify the internal stress evolution trend of the material, providing a key basis for judging the elastoplastic state and ensuring the coherence of mechanical response analysis.

[0052] 206. Calculate the initial variable array based on the deviatoric stress change rate and the preset yield function to obtain the deviatoric stress tensor;

[0053] In this embodiment, the yield function is used to iteratively solve the deviatoric stress tensor for inelastic materials, accurately capturing the transition behavior of the material from elastic to plastic; for elastic materials, the deviatoric stress is directly reused to update, achieving optimized allocation of computational resources and taking into account both the adaptability of constitutive models for different materials and computational efficiency.

[0054] 207. Generate state parameters based on the density change rate correction term, velocity gradient correction term, torsional rate correction term, strain rate correction term, artificial viscosity correction term, virtual stress correction term, extended smooth particle dynamics correction term, normal contact force set, tangential contact force set, updated particle displacement set, deviatoric stress change rate and deviatoric stress tensor.

[0055] In this embodiment, multi-dimensional correction terms, force sets, and mechanical tensors are integrated to form a complete particle state parameter system, which comprehensively reflects the particle's motion, deformation, stress, and contact state. This parameter set provides a unified input for particle state updates, ensuring the systematic nature and accuracy of multi-physics coupling analysis and supporting the efficient advancement of the entire simulation process.

[0056] In this embodiment, the efficiency and accuracy of elastoplastic particle simulation are improved through multi-dimensional correction term calculation and parallel optimization. First, basic correction terms such as density change rate and velocity gradient are calculated by traversal calculation, and parallel computation is achieved through thread allocation, providing accurate underlying data for mechanical analysis. Strain rate and artificial viscosity correction terms are introduced to support constitutive relation solving and suppress numerical oscillations. Virtual stress and dynamic correction terms optimize the stress field distribution and stabilize particle trajectory. The contact force set is solved based on master-slave contact surface partitioning, balancing the realism of contact behavior with computational efficiency. Multiple correction terms are integrated to update particle displacement and deviatoric stress change rate, ensuring coherence in mechanical response analysis. Deviatoric stress tensor solutions are adapted for different materials, optimizing computational resource allocation. Finally, a complete state parameter system is formed, comprehensively reflecting the particle state, ensuring the systematicity and accuracy of multi-physics coupling analysis, and providing efficient and reliable technical support for the simulation of complex elastoplastic problems.

[0057] Please see Figure 3 In a third embodiment of the elastoplastic smooth particle state update method of the present invention, step 104 includes:

[0058] 301. Obtain particle position information from the initialized variable array;

[0059] In this embodiment, the maximum and minimum values ​​of particle positions are used as the boundary range;

[0060] 302. Analyze the particle position information to obtain the particle distribution range;

[0061] 303. Divide the particle distribution range into a grid based on the particle support domain radius to obtain the grid division result;

[0062] 304. Number the initialization variable array according to the mesh division results to obtain the particle unit number set;

[0063] 305. Sort the particle unit number set to obtain a list of unit particle numbers;

[0064] In this embodiment, the particle support domain search radius is determined based on the smooth particle dynamics simulation parameters in the preprocessing file, and this radius is used as the cell size. Within the mesh boundary, cells are divided according to the cell size and assigned cell numbers. All particles are traversed, and the cell corresponding to the particle is determined based on the particle position.

[0065] 306. Perform a neighborhood search on the list of unit particle numbers based on the particle position information to obtain the list of neighboring particles;

[0066] In this embodiment, the efficient radix sorting of the GPU is used to determine the list of particle labels in the cell. The particles are sorted according to the cell number corresponding to the particle to obtain the rearranged particle list and the start and end index of each cell. The unit particle number list can be obtained based on the rearranged particle list and the start and end index of each cell. Since the particle support domain radius is the same as the grid cell size, a neighborhood search is performed based on the particle position information and the unit particle list. The search is only performed in the cell where the particle is located and the adjacent cells to determine whether the particle pair distance is less than the support domain radius. The neighboring particle list data is obtained, including particle pair distance, kernel function and its derivative value, etc.

[0067] In this embodiment, the boundary range is determined by the maximum and minimum values ​​of particle positions, and the mesh is divided and numbered in combination with the particle support domain radius to form a standardized set of particle unit numbers. The unit particle number list is sorted using GPU efficient cardinality sorting, and particles are located according to cell indexes. The neighborhood search range is limited to the cell where the particle is located and the adjacent cells, reducing the search range and avoiding inefficient global search, thus reducing the proportion of search time. At the same time, the particle support domain radius is matched with the mesh cell size, and the neighboring particle data is obtained by combining distance judgment. This not only ensures search accuracy, but also improves hardware resource utilization through GPU parallel acceleration, effectively solving the problem of low computational efficiency in large-scale particle scenes and laying the foundation for high efficiency in the entire simulation process.

[0068] Please see Figure 4 In a fourth embodiment of the elastoplastic smooth particle state update method of the present invention, step 201 includes:

[0069] 401. Based on the density change rate calculation formula, the neighboring particle list data, and the thread allocation results, the initial variable array is traversed and calculated to obtain the density change rate set;

[0070] In this embodiment, the density change rate of all particles in the initialization variable array is calculated. The expression for the density change rate calculation formula is as follows:

[0071] (1)

[0072] In the formula, Let i be the rate of change of particle density. Let i be the density of particle i. For particle i, there is a list of neighboring particles. Let j be the mass of the particle. Let the particle density be j. The relative velocity between particle i and particle j is directional components, Let i be the value of the kernel function at position j of particle i. For particle i kernel function The directional derivative takes its value at particle position j; the formula incorporates particle mass, density, relative velocity components, and a kernel function (including the derivative) to quantify the impact of mass exchange between particles on density, and then obtains a correction term through accumulating the set, accurately capturing the spatiotemporal changes of the density field during dynamic processes. Representing particle i in Directional coordinates, representing the partial derivative of the kernel function W(x) in that direction. ;

[0073] 402. Accumulate the set of density change rates to obtain the density change rate correction term;

[0074] 403. Based on the velocity gradient calculation formula and thread allocation results, traverse and calculate the initial variable array to obtain the velocity gradient set;

[0075] In this embodiment, the expression for the velocity gradient calculation formula is as follows:

[0076] (2)

[0077] In the formula, The relative velocity between particle i and other particles Directional components in Directional derivative, The relative velocity between particle i and particle j is Directional components; by using the relative velocity components between particles and their spatial derivatives, the velocity spatial change rate of each particle is solved by iterating through the velocity gradient formula. This calculation reflects the non-uniformity of the velocity field and is the basis for subsequent strain rate and stress field analysis. The correction term formed after accumulating the set can accurately describe the gradient characteristics of the velocity field.

[0078] 404. Accumulate the velocity gradient set to obtain the velocity gradient correction term;

[0079] 405. Based on the torsion ratio calculation formula and thread allocation results, traverse and calculate the initial variable array to obtain the torsion ratio set;

[0080] In this embodiment, the expression for the torsion ratio calculation formula is as follows:

[0081] (3)

[0082] In the formula, Let i be the torsion rate of particle i. The rotational effect of the particle is calculated by iterating through the torsion rate formula (based on the antisymmetric component of the velocity gradient). By quantifying the degree of torsion of the particle around its own axis, the correction term obtained by accumulating the set can supplement the description of the rotational characteristics of the particle motion and improve the characterization of complex deformation behavior.

[0083] 406. Accumulate the torsion rate set to obtain the torsion rate correction term;

[0084] In this embodiment, the accuracy and efficiency of particle dynamics simulation are improved by calculating the density change rate, velocity gradient, and torsional rate correction term. The density change rate calculation combines parameters such as particle mass and density with kernel functions to quantify the impact of mass exchange on density and capture the spatiotemporal changes of the density field. The velocity gradient calculation reflects the non-uniformity of the velocity field, laying the foundation for subsequent strain rate and stress field analysis. The torsional rate calculation quantifies the particle rotation effect and improves the characterization of complex deformation behavior. All three calculations rely on thread allocation to achieve parallel computation, and the correction term is formed by accumulating sets. While ensuring the integrity of capturing key mechanical characteristics of particle motion, the computational efficiency of large-scale particle systems is improved, providing high-precision and stable underlying data support for subsequent multiphysics coupling analysis.

[0085] Please see Figure 5 In the fifth embodiment of the elastoplastic smooth particle state update method of the present invention, step 202 includes:

[0086] 501. Based on the strain rate calculation formula and thread allocation results, the initial variable array is traversed and calculated to obtain the strain rate set;

[0087] In this embodiment, the expression for the strain rate calculation formula is as follows:

[0088] (4)

[0089] In the formula, Let be the strain rate of particle i, and represent the time derivative of the strain of particle i. Representing particle i in Direction coordinates, representing velocity overall. In this direction of partial derivative ;

[0090] 502. Perform cumulative calculations on the strain rate set to obtain the strain rate correction term;

[0091] In this embodiment, the strain rate calculation formula is based on the relative motion state between particles, quantifies the deformation rate of the particle system, and is a key physical quantity describing the flow and deformation characteristics of materials; the strain rate correction term is directly related to the subsequent constitutive relation solution, providing an important basis for analyzing the transition of materials from elastic to plastic; the calculation process combines thread allocation to achieve parallel processing, ensuring efficient acquisition of strain rate information of each particle in a large-scale particle system;

[0092] 503. Based on the calculation formula of the artificial viscosity correction term and the thread allocation results, the initial variable array is traversed and calculated to obtain the set of artificial viscosity correction parameters;

[0093] In this embodiment, the expression for the formula for calculating the artificial viscosity correction parameter is as follows:

[0094] (5)

[0095] (6)

[0096] In the formula, It serves as an artificial viscosity correction parameter between particle i and particle j, and is used in the calculation of particle acceleration. The first artificial viscosity coefficient, This is the second artificial viscosity correction parameter. Let i be the relative velocity between particle i and particle j. Let i be the relative position between particle i and particle j. Let i be the distance between particle i and particle j. For the smooth length of the particle, This is the artificial viscosity coefficient, and this empirical value is usually taken as 0.01; , Let be the sound speeds at particle i and particle j, respectively. , The values ​​are the densities at particle i and particle j, respectively, both obtained from material properties. The average speed of sound, For intermediate parameters, Average density;

[0097] 503. Accumulate and calculate the set of artificial viscosity correction parameters to obtain the artificial viscosity correction term;

[0098] In this embodiment, the artificial viscosity correction parameter calculation formula determines the mutual motion trend of particles by relative velocity and relative position, adjusts the range of action by combining particle spacing and smooth length, uses sound speed and density parameters to correlate material properties, and controls the correction intensity by artificial viscosity coefficient. This design can effectively simulate the viscosity effect in physical systems and compensate for the oscillation problem that may occur in numerical calculations. It also relies on thread allocation to achieve parallel computing and ensures computational efficiency in large-scale particle scenarios.

[0099] In this embodiment, the scheme improves the accuracy, stability, and efficiency of particle simulation by calculating strain rate and artificial viscosity correction term. Strain rate calculation is based on the deformation rate quantified by the relative motion of particles, providing a key basis for elastoplastic constitutive analysis. Combined with thread allocation, parallel processing is achieved, efficiently acquiring strain information of large-scale particle systems. Artificial viscosity correction adjusts the range of action based on relative motion trends, particle spacing, etc., correlates with material properties, and controls the correction intensity, effectively simulating viscous effects and compensating for numerical oscillations. Both rely on parallel computing to ensure efficiency, accurately supporting the analysis of material transition from elastic to plasticity, suppressing non-physical phenomena, adapting to complex deformation scenarios, and providing high-quality data for subsequent simulations.

[0100] Please see Figure 6 In the sixth embodiment of the elastoplastic smooth particle state update method of the present invention, step 203 includes:

[0101] 601. Obtain the current stress state of the particle group and analyze the current stress state of the particle group;

[0102] 602. If the current stress state of the particle group is the principal stress state, then the first principal stress set is obtained from the initialization variable array according to the current stress state of the particle group.

[0103] In this embodiment, the current stress state of the particle group is analyzed first, and only the particles in the principal stress state are processed. The particles with principal stress greater than zero (i.e., tensile particles) are selected through the first principal stress set to form the second principal stress set, thus avoiding ineffective correction of the compressive particles.

[0104] 603. Determine whether there are particles with principal stresses greater than zero in the initialization variable array based on the first principal stress set;

[0105] 604. If there are particles with principal stresses greater than zero in the initialization variable array, obtain the principal stresses corresponding to all particles with principal stresses greater than zero to obtain the second principal stress set;

[0106] In this embodiment, the virtual stress is calculated in a targeted manner by using principal stress screening to avoid wasting resources; the thread parallel mechanism is adapted to large-scale particle systems, which improves computational efficiency and shortens simulation time while ensuring accuracy.

[0107] 605. Based on the virtual stress calculation formula, the second principal stress set, and the thread allocation results, the initial variable array is traversed and calculated to obtain the virtual stress set;

[0108] In this embodiment, all particles are traversed. Under principal stress conditions, if the stress in any principal stress direction of a particle is greater than zero, the virtual stress in the corresponding direction is obtained based on the stress value, particle density, and virtual stress coefficient. This virtual stress is precisely applied to the tensile stress region, thus overcoming the shortcomings of the traditional SPH method (Smoothed Particle Hydrodynamics, a meshless numerical simulation technique initially used to solve fluid motion problems in astrophysics and later widely applied to complex dynamic simulations in fluid mechanics, solid mechanics, explosion mechanics, and other fields) in describing the tensile state. The virtual stress is then transformed back to the original coordinate system to ensure the spatial consistency of the stress field. Taking a planar problem as an example, the formula for calculating the virtual stress in the principal stress direction is:

[0109] (7)

[0110] In the formula, , These are virtual stresses in different principal stress directions, which are used in the calculation of particle acceleration. This is the virtual stress coefficient. , These are the principal stresses in different principal stress directions;

[0111] 606. Perform cumulative calculations on the virtual stress set to obtain the virtual stress correction term;

[0112] 607. Based on the dynamic correction formula and thread allocation results, traverse and calculate the initial variable array to obtain the set of velocity correction terms;

[0113] In this embodiment, the formula for calculating the ExtraSmoothed Particle Hydrodynamics (XSPH) term is as follows:

[0114] (8)

[0115] In the formula, For XSPH correction speed, For speed, XSPH correction factor Let j be the mass of particle j. Let i be the average density of particle i and particle j. Let be the relative velocity between particle i and particle j. The introduced XSPH-corrected velocity can stabilize particle motion. Except for the velocity update step, all other steps that use velocity as input in the program use XSPH-corrected velocity to ensure the continuity of particle motion trajectory. XSPH correction effectively reduces non-physical phenomena such as particle penetration and aggregation through velocity smoothing, making particle motion more in line with actual mechanical laws, reducing numerical oscillations in dynamic simulation, and improving the reliability of results.

[0116] 608. Accumulate the set of velocity correction terms to obtain the extended smooth particle dynamics correction terms;

[0117] In this embodiment, virtual stress correction focuses on the principal stress state, selects tensile particles with principal stress greater than zero to form a second principal stress set, calculates virtual stress in the corresponding direction and transforms it back to the original coordinate system, accurately makes up for the problem of insufficient description of tensile state by traditional SPH, and avoids ineffective correction and waste of resources; XSPH correction uses velocity smoothing to generate correction velocity based on parameters such as particle mass and average density. Except for the update step, the correction velocity is used to reduce non-physical phenomena such as particle penetration and aggregation, and ensure the continuity of motion trajectory. The thread parallel mechanism is adapted to large-scale particle systems. While improving the calculation accuracy of tensile stress scene and reducing numerical oscillation, it optimizes resource allocation and shortens simulation time, providing efficient and reliable support for the simulation of complex mechanical behavior.

[0118] Please see Figure 7 In the seventh embodiment of the elastoplastic smooth particle state update method of the present invention, step 204 includes:

[0119] 701. Obtain the set of particles on the secondary surface and the set of particles on the primary surface from the initialized variable array based on the primary and secondary contact surfaces;

[0120] In this embodiment, the main contact surface and the secondary contact surface are determined according to the rigidity of the two contact surfaces. Generally, the surface with greater rigidity is set as the main contact surface and the other contact surface is set as the secondary contact surface. The particle sets of the main surface and the particle sets of the secondary surface are extracted respectively. This division is based on the physical characteristics to determine the reference benchmark for contact calculation. Since the main contact surface has small deformation and high stability, the reliability of contact judgment can be improved.

[0121] 702. Calculate the contact normal vector set by calculating the set of particles on the outer surface and the set of particles on the main surface according to the preset normal vector calculation formula;

[0122] In this embodiment, the expression for the normal vector calculation formula is as follows:

[0123] (9)

[0124] In the formula, For particles from the contact surface The contact normal vector to the main contact surface provides a basis for determining the direction of subsequent forces. The set of particles located within the support domain of particle i on the main contact surface. Let be the relative displacement vector between particle i and particle j. Let be the distance between particle i and particle j;

[0125] 703. Calculate the contact distance set by calculating the distance between the particle set on the outer surface and the particle set on the main surface according to the preset distance calculation formula;

[0126] In this embodiment, the distance calculation formula is:

[0127] (10)

[0128] In the formula, The contact distance is from particle i on the contact surface to the main contact surface.

[0129] 704. Analyze the set of contact distances according to the preset contact threshold to obtain contact analysis results;

[0130] In this embodiment, the particles in contact with the main contact surface from the contact surface are determined based on the contact distance and the contact threshold. The calculation formula is as follows: (11)

[0131] For contact threshold, For the contact criterion, when Contact occurs at times. Contact distance;

[0132] 705. Obtain the set of contact particle pairs from the initialization variable array based on the contact analysis results;

[0133] In this embodiment, the spatial distance from the surface particle to the main contact surface is quantified by the distance formula, and effective contact particle pairs are screened by combining the contact threshold, avoiding invalid calculations for non-contact particles and reducing the amount of calculation.

[0134] 706. Calculate the set of contact particles based on the normal contact stiffness coefficient, the set of contact normal vectors, and the tangential stiffness coefficient to obtain the set of normal contact forces and the set of tangential contact forces.

[0135] In this embodiment, the formula for calculating the normal contact force is: (12)

[0136] In the formula, For normal contact force, This is the normal contact stiffness coefficient;

[0137] (13)

[0138] In the formula, For tangential contact force, The tangential stiffness coefficient, The tangential contact displacement of the previous time step (when it is the first time step, the tangential contact displacement is 0. In the cutting simulation model of this SPH program, the parts do not make contact in the initial state, that is, the initial geometry has no contact or penetration, and the initial tangential contact force of the particles is 0. Therefore, the tangential contact displacement of the initial time step is 0. As the simulation progresses, in the first time step when the parts make contact, the tangential contact force of the particles in the previous time step is 0, so the tangential contact displacement of the previous time step is 0). This represents the tangential contact displacement at the current time step. The relative velocity between particle i and particle j at the current time step. It is the tangential vector. For time step; if Then update it to To satisfy the sliding friction condition, the normal contact force reflects the compressive strength, and the tangential contact force introduces friction constraints. Both characterize the normal compression and tangential friction behavior between particles, making the contact force calculation more consistent with the mechanical nature of material contact and improving the credibility of the simulation results.

[0139] In this embodiment, the scheme improves the reliability and efficiency of particle contact simulation by dividing the contact surfaces into master and slave surfaces and calculating contact forces. The master and slave contact surfaces are divided based on rigidity differences, with the master surface serving as the contact calculation benchmark due to its small deformation and high stability, thus improving the reliability of the judgment. Targeted processing of the slave and master particle sets ensures that the contact analysis focuses on key areas of action. The contact direction is accurately determined using the normal vector formula, and the spatial distance is quantified using the distance formula. Effective contact pairs are selected using contact thresholds to avoid invalid calculations of non-contact particles, reducing the computational load. The normal contact force reflects the compressive strength, while the tangential contact force introduces friction constraints and is limited by sliding conditions. Both characterize the contact mechanical behavior, closely reflecting the essence of material contact. The overall scheme optimizes computational efficiency while ensuring the accuracy of contact judgment and the reliability of simulation results, adapting to different material contact scenarios and providing efficient and reliable support for dynamic simulation of complex contact surfaces.

[0140] Please see Figure 8 In the eighth embodiment of the elastoplastic smooth particle state update method of the present invention, step 205 includes:

[0141] 801. Calculate the initial variable array according to the preset equilibrium equation, artificial viscosity correction term and virtual stress correction term to obtain the particle acceleration set;

[0142] In this embodiment, the acceleration of each particle is calculated according to the equilibrium equation, the expression of which is as follows:

[0143] (14)

[0144] In the formula, Let i be the acceleration of particle i. Let i be the stress tensor of particle i. Let j be the stress tensor of particle j. It is the sum of the artificial viscosity correction term and the virtual stress correction term. Let be the value of the total differential of the kernel function of particle i at position j. The external force field is used to ensure that the acceleration calculations closely match the actual dynamic laws;

[0145] 802. Obtain the particle displacement set from the initialization variable array;

[0146] 803. Update the particle displacement set based on the particle acceleration set to obtain the updated particle displacement set;

[0147] In this embodiment, the dynamic evolution of the particle's trajectory is achieved through the integration process from acceleration (second derivative) to displacement, ensuring the consistency between position changes and mechanical driving forces.

[0148] The initial variable array is calculated according to the preset state equation to obtain the hydrostatic pressure of the particles;

[0149] In this embodiment, the expression for the state equation is as follows:

[0150] (15)

[0151] In the formula, Let be the hydrostatic pressure of the particle, and c be the speed of sound at the particle's location. The particle density at the current time step. The initial time step particle density is used; by calculating the particle sound velocity (material property) and the difference between the current density and the initial density, the hydrostatic pressure is solved, quantifying the pressure effect caused by the volume change of the material. This calculation directly reflects the volume deformation characteristics of compressible materials and is an indispensable parameter in the simulation of large deformation of fluids and solids.

[0152] 804. Calculate the initial variable array according to the strain rate correction term, torsion rate correction term and the preset elastic constitutive relation calculation formula to obtain the deviatoric stress change rate;

[0153] In this embodiment, the expression for the formula for calculating the elastic constitutive relation is as follows:

[0154] (16)

[0155] The rate of change of deviatoric stress, Shear modulus This is the strain rate correction term. It is the identity matrix. It is a deviatoric stress. This is the torsional rate correction term; the formula for calculating the elastic constitutive relation describes the stress evolution corresponding to shear deformation in the elastic stage, taking into account both the translational deformation and rotational effects of particles, thus improving the description of complex elastic behavior.

[0156] Update the deviatoric stress based on the rate of change of deviatoric stress:

[0157] (17)

[0158] In the formula, To update the deviatoric stress;

[0159] In this embodiment, artificial viscosity and virtual stress correction terms are integrated based on the equilibrium equation, and particle acceleration is accurately solved using the total differential of the kernel function to ensure that the dynamic calculations closely match actual laws. Acceleration drives displacement updates, and the dynamic evolution of the motion trajectory is realized through the integration process, ensuring that position changes are consistent with mechanical driving forces. The equation of state is correlated with density changes and hydrostatic pressure to quantify volume deformation effects and adapt to the simulation requirements of compressible materials. The elastic constitutive relation integrates strain rate and torsional rate correction terms to characterize the deviatoric stress changes under shear deformation and rotation effects, and provides a coherent stress state for subsequent analysis through iterative updates. The overall scheme takes into account both the realism of mechanical behavior and computational coherence, adapts to multi-material and multi-scenario simulations, and provides reliable support for complex dynamic analysis.

[0160] Please see Figure 9 In the ninth embodiment of the elastoplastic smooth particle state update method of the present invention, step 206 includes:

[0161] 901. Obtain the material constitutive model and analyze it to determine the material constitutive model type;

[0162] 902. Determine whether the constitutive model of the material is an inelastic model;

[0163] 903. When the material constitutive model is an inelastic model, the initial variable array is iteratively calculated based on the deviatoric stress change rate and yield function to obtain the deviatoric stress tensor.

[0164] In this embodiment, if the material constitutive model is elastic, it can be used directly. The deviatoric stress is updated as the deviatoric stress tensor; otherwise, the particle deviatoric stress tensor should be solved iteratively. The specific process is as follows:

[0165] The yield function that needs to be solved iteratively is:

[0166] (18)

[0167] In the formula, Let be the yield function. The deviatoric stress tensor at the current time step. The yield stress at the current time step is obtained from the Johnson-Cook plastic constitutive model. For the equivalent plastic strain at the current time step, The equivalent plastic strain rate at the current time step. For the current time step temperature, these variables are all integrals with the plastic flow factor. Related;

[0168] (19)

[0169] In the formula, For heat-to-work conversion efficiency, For specific heat capacity, The yield stress is the same as the yield stress at the previous time step (when the material is not in the plastic state, it is in the elastic state; at the first time step when it enters the plastic state, the yield stress is the same as the yield stress at the previous time step in the elastic state, the equivalent plastic strain and the equivalent plastic strain rate are all 0, there is no internal or external heat inflow, and the temperature is the initial temperature. Regarding the initial values ​​in the elastic state, this simulation program does not set the initial stress and strain field, which is in line with the common settings of existing simulation programs on the market. At the initial time step of the simulation, the equivalent plastic strain and the equivalent plastic strain rate are both 0). Shear modulus For the flow direction; use Newton's iteration method to integrate the main variable, the plastic flow factor. Perform numerical solutions until the iterative variables change. Less than the threshold; The temperature of the previous time step is used to ensure that the nonlinear evolution of the deviatoric stress tensor in the plastic stage is accurately characterized. This is the equivalent plastic strain rate at the current time step. This variable value is only related to the integral value of the main variable, the plastic flow factor. and time step Related; For the equivalent plastic strain of the previous time step, For time step size, The yield stress is the stress at the current time step; ":" indicates the double dot product of the tensor. The updated deviatoric stress tensor (i.e., the updated deviatoric stress) shown in formula (17) is a second-order tensor, therefore, The dot product of two second-order tensors results in a scalar: ,in, For a planar problem, the problem space dimension is... Regarding the space issue, ;

[0170] In this embodiment, by distinguishing the types of material constitutive models, efficient solutions for the deviatoric stress tensor are achieved. For elastic models, the deviatoric stress is directly reused for updating, simplifying the calculation. For inelastic models, the yield stress is obtained through iterative solution based on the yield function. Variables related to plastic flow factors, such as equivalent plastic strain, strain rate, and temperature, are integrated, and the nonlinear evolution of the plastic stage is accurately characterized by the Newton-Raphson iteration method, taking into account both force and thermal coupling effects. This scheme ensures both computational efficiency in elastic scenarios and improves the accuracy of stress simulation for inelastic materials. It is adaptable to multiple scenarios, and the Newton-Raphson iteration ensures numerical stability, providing reliable support for the analysis of complex mechanical behaviors and expanding the scope of engineering applications.

[0171] The above describes a method for updating the state of elastoplastic smooth particles according to an embodiment of the present invention. The following describes a device for updating the state of elastoplastic smooth particles according to an embodiment of the present invention. Please refer to [link to relevant documentation]. Figure 10 One embodiment of the elastoplastic smooth particle state update device of the present invention includes:

[0172] Parameter acquisition module 1 is used to acquire simulation parameters;

[0173] Array generation module 2 is used to generate a particle swarm variable array based on a preset simulation model;

[0174] Initialization module 3 is used to initialize the particle swarm variable array according to the simulation parameters to obtain the initial variable array;

[0175] The neighborhood search module 4 is used to perform a neighborhood search on the initial variable array based on the preset particle support domain radius to obtain the list of neighboring particles.

[0176] State determination module 5 is used to determine the GPU parallel computing state based on the neighboring particle list data;

[0177] Thread allocation module 6 is used to allocate threads to each particle in the initialization variable array under the adjusted GPU parallel computing state in order to obtain the thread allocation result;

[0178] The state parameter calculation module 7 is used to calculate the initial variable array based on the thread allocation result to obtain the state parameters;

[0179] The state update module 8 is used to update the state of all particles in the initial variable array according to the state parameters to obtain the updated particle variable array;

[0180] In this embodiment, simulation parameters are first read to provide the input basis for calculation. By generating and initializing a particle swarm variable array, the initial state of the particles is ensured to be highly matched with the simulation target. The core advantage lies in adopting a neighborhood search strategy based on the particle support domain radius, combined with GPU radix sorting technology, to reduce invalid calculations and improve the generation efficiency of the neighboring particle list data. Through a "one particle, one thread" GPU parallel computing architecture, an independent thread is allocated to each particle to fully utilize the hardware computing power and improve the simulation efficiency of large-scale particle systems. By updating the particle variable array, continuous simulation at multiple time steps is achieved, capturing the elastoplastic evolution process of materials under dynamic loads and adapting to the needs of complex mechanical behavior analysis. The overall solution takes into account both versatility and scalability, supports the configuration of simulation parameters in multiple fields, and the modular process facilitates functional expansion, providing efficient and reliable numerical simulation support for structural strength analysis, material failure prediction, and other scenarios in engineering design.

[0181] Finally, it should be noted that the above descriptions are merely preferred embodiments of the present invention and are not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. An elasto-plastic smooth particle state updating method, characterized by, The method comprises the following steps: acquiring simulation parameters; generating a particle swarm variable array according to a preset simulation model; initializing the particle swarm variable array according to the simulation parameters to obtain an initialized variable array; performing neighborhood search on the initialized variable array according to a preset particle support domain radius to obtain near particle list data; determining a GPU parallel computing state according to the near particle list data; allocating threads to each particle in the initialized variable array under the adjusted GPU parallel computing state to obtain a thread allocation result; calculating the initialized variable array according to the thread allocation result to obtain state parameters; the calculation of the initialized variable array according to the thread allocation result to obtain the state parameters comprises: performing traversal calculation on the initialized variable array according to a preset density change rate calculation formula, a preset velocity gradient calculation formula, a preset torsion rate calculation formula and the thread allocation result to obtain a density change rate correction term, a velocity gradient correction term and a torsion rate correction term; performing traversal calculation on the initialized variable array according to a preset strain rate calculation formula, a preset artificial viscosity correction term calculation formula and the thread allocation result to obtain a strain rate correction term and an artificial viscosity correction term; the expression of the artificial viscosity correction parameter calculation formula is as follows: ， ， where, is the artificial viscosity correction parameter between particle i and particle j, which participates in the calculation of particle acceleration; is the first artificial viscosity coefficient, is the second artificial viscosity correction parameter, is the relative velocity between particle i and particle j, is the relative position between particle i and particle j, is the distance between particle i and particle j, is the smoothing length of particle, is the artificial viscosity coefficient, which is an empirical value of 0.01; , is the sound speed at particle i and particle j, respectively, , is the density at particle i and particle j, respectively, which is obtained from the material properties, is the average sound speed, is the intermediate parameter, is the average density; performing traversal calculation on the initialized variable array according to a preset virtual stress calculation formula, a preset dynamics correction formula and the thread allocation result to obtain a virtual stress correction term and an extended smooth particle dynamics correction term; calculating the initialized variable array according to a preset normal contact stiffness coefficient, a preset tangent stiffness coefficient, a preset main contact surface and a preset secondary contact surface to obtain a normal contact force set and a tangent contact force set; calculating the initialized variable array according to the artificial viscosity correction term, the virtual stress correction term, the strain rate correction term and the torsion rate correction term to obtain an updated particle displacement set and a deviatoric stress change rate; calculating the initialized variable array according to the deviatoric stress change rate and a preset yield function to obtain a deviatoric stress tensor; generating state parameters according to the density change rate correction term, the velocity gradient correction term, the torsion rate correction term, the strain rate correction term, the artificial viscosity correction term, the virtual stress correction term, the extended smooth particle dynamics correction term, the normal contact force set, the tangent contact force set, the updated particle displacement set, the deviatoric stress change rate and the deviatoric stress tensor; updating the states of all particles in the initialized variable array according to the state parameters to obtain an updated particle variable array.

2. An elasto-plastic smooth particle state updating method according to claim 1, wherein, the neighborhood search on the initialized variable array according to the preset particle support domain radius to obtain the near particle list data comprises: obtaining particle position information from the initialized variable array; analyzing the particle position information to obtain a particle distribution range; performing grid division on the particle distribution range according to the particle support domain radius to obtain a grid division result; numbering the initialized variable array according to the grid division result to obtain a particle element number set; sorting the particle element number set to obtain an element particle number list; According to the particle position information, a neighborhood search is performed on a cell particle number list to obtain near particle list data.

3. An elasto-plastic smooth particle state updating method according to claim 2, wherein, The preset density change rate calculation formula, the preset velocity gradient calculation formula, the preset torsion rate calculation formula, and the thread allocation result are used to perform traversal calculation on the initialized variable array to obtain a density change rate correction term, a velocity gradient correction term, and a torsion rate correction term, including: The density change rate calculation formula, the near particle list data, and the thread allocation result are used to perform traversal calculation on the initialized variable array to obtain a density change rate set; The density change rate set is accumulated to obtain the density change rate correction term; The velocity gradient calculation formula and the thread allocation result are used to perform traversal calculation on the initialized variable array to obtain a velocity gradient set; The velocity gradient set is accumulated to obtain the velocity gradient correction term; The torsion rate calculation formula and the thread allocation result are used to perform traversal calculation on the initialized variable array to obtain a torsion rate set; The torsion rate set is accumulated to obtain the torsion rate correction term.

4. The elasto-plastic smooth particle state updating method of claim 1, wherein, The preset strain rate calculation formula, the preset artificial viscosity correction term calculation formula, and the thread allocation result are used to perform traversal calculation on the initialized variable array to obtain a strain rate correction term and an artificial viscosity correction term, including: The strain rate calculation formula and the thread allocation result are used to perform traversal calculation on the initialized variable array to obtain a strain rate set; The strain rate set is accumulated to obtain the strain rate correction term; The artificial viscosity correction term calculation formula and the thread allocation result are used to perform traversal calculation on the initialized variable array to obtain an artificial viscosity correction parameter set; The artificial viscosity correction parameter set is accumulated to obtain the artificial viscosity correction term.

5. The elasto-plastic smooth particle state updating method of claim 1, wherein, The preset virtual stress calculation formula, the preset dynamic correction formula, and the thread allocation result are used to perform traversal calculation on the initialized variable array to obtain a virtual stress correction term and an extended smooth particle dynamics correction term, including: A current stress state of a particle group is obtained, and the current stress state of the particle group is analyzed; If the current stress state of the particle group is a principal stress state, a first principal stress set is obtained from the initialized variable array according to the current stress state of the particle group; It is judged whether there is a particle with a principal stress greater than zero in the initialized variable array according to the first principal stress set; When there is a particle with a principal stress greater than zero in the initialized variable array, a principal stress corresponding to all particles with a principal stress greater than zero is obtained to obtain a second principal stress set; The virtual stress calculation formula, the second principal stress set, and the thread allocation result are used to perform traversal calculation on the initialized variable array to obtain a virtual stress set; The virtual stress set is accumulated to obtain the virtual stress correction term; The dynamic correction formula and the thread allocation result are used to perform traversal calculation on the initialized variable array to obtain a velocity correction term set; The velocity correction term set is accumulated to obtain the extended smooth particle dynamics correction term.

6. An elasto-plastic smooth particle state updating method according to claim 1, wherein, The preset normal contact stiffness coefficient, the preset tangential stiffness coefficient, the preset main contact surface and the preset secondary contact surface are used for calculating the initialization variable array to obtain a normal contact force set and a tangential contact force set, including: A main surface particle set and a secondary surface particle set are obtained from the initialization variable array according to the main contact surface and the secondary contact surface; The main surface particle set and the secondary surface particle set are calculated according to a preset normal vector calculation formula to obtain a contact normal vector set; The main surface particle set and the secondary surface particle set are calculated according to a preset distance calculation formula to obtain a contact distance set; The contact distance set is analyzed according to a preset contact threshold to obtain a contact analysis result; A contact particle pair set is obtained from the initialization variable array according to the contact analysis result; The contact particle pair set is calculated according to the normal contact stiffness coefficient, the contact normal vector set and the tangential stiffness coefficient to obtain the normal contact force set and the tangential contact force set.

7. An elasto-plastic smooth particle state updating method according to claim 1, wherein, The initialization variable array is calculated according to an artificial viscosity correction term, a virtual stress correction term, a strain rate correction term and a torsion rate correction term to obtain an updated particle displacement set and a deviatoric stress change rate, including: The initialization variable array is calculated according to a preset balance equation, the artificial viscosity correction term and the virtual stress correction term to obtain a particle acceleration set; A particle displacement set is obtained from the initialization variable array; The particle displacement set is updated according to the particle acceleration set to obtain an updated particle displacement set; The initialization variable array is calculated according to the strain rate correction term, the torsion rate correction term and a preset elastic constitutive relation calculation formula to obtain the deviatoric stress change rate.

8. An elasto-plastic smooth particle state updating method according to claim 1, wherein, The initialization variable array is calculated according to the deviatoric stress change rate and a preset yield function to obtain a deviatoric stress tensor, including: A material constitutive model is obtained, and the material constitutive model is analyzed to obtain a material constitutive model type; It is judged whether the material constitutive model type is a non-elastic model; When the material constitutive model type is a non-elastic model, the initialization variable array is iteratively calculated according to the deviatoric stress change rate and the yield function to obtain the deviatoric stress tensor.

9. An elastoplastic smooth particle state updating device characterized by comprising: Including: A parameter acquisition module is configured to acquire simulation parameters; An array generation module is configured to generate a particle swarm variable array according to a preset simulation model; An initialization module is configured to initialize the particle swarm variable array according to the simulation parameters to obtain an initialization variable array; A neighborhood search module is configured to perform neighborhood search on the initialization variable array according to a preset particle support domain radius to obtain near particle list data; A state determination module is configured to determine a GPU parallel computing state according to the near particle list data; A thread allocation module is configured to allocate threads to each particle in the initialization variable array under the adjusted GPU parallel computing state to obtain a thread allocation result; A state parameter calculation module is configured to calculate the initialization variable array according to the thread allocation result to obtain state parameters, specifically including: According to the preset density variation rate calculation formula, the preset velocity gradient calculation formula, the preset torsion rate calculation formula and the thread allocation result, the initialized variable array is iteratively calculated to obtain a density variation rate correction term, a velocity gradient correction term and a torsion rate correction term; According to the preset strain rate calculation formula, the preset artificial viscosity correction term calculation formula and the thread allocation result, the initialized variable array is iteratively calculated to obtain a strain rate correction term and an artificial viscosity correction term; The expression of the artificial viscosity correction parameter calculation formula is as follows: ， ， where, is the artificial viscosity correction parameter between particle i and particle j, which participates in the calculation of particle acceleration; is the first artificial viscosity coefficient, is the second artificial viscosity correction parameter, is the relative velocity between particle i and particle j, is the relative position between particle i and particle j, is the distance between particle i and particle j, is the smoothing length of particle i, is the artificial viscosity coefficient, which is an empirical value of 0.01; , are the sound speeds at particle i and particle j, respectively, , are the densities at particle i and particle j, respectively, which are obtained from material properties, is the average sound speed, is the intermediate parameter, is the average density; According to the preset virtual stress calculation formula, the preset dynamics correction formula and the thread allocation result, the initialized variable array is iteratively calculated to obtain a virtual stress correction term and an extended smooth particle dynamics correction term; According to the preset normal contact stiffness coefficient, the preset tangential stiffness coefficient, the preset main contact surface and the preset secondary contact surface, the initialized variable array is calculated to obtain a normal contact force set and a tangential contact force set; According to the artificial viscosity correction term, the virtual stress correction term, the strain rate correction term and the torsion rate correction term, the initialized variable array is calculated to obtain an updated particle displacement set and a deviatoric stress change rate; According to the deviatoric stress change rate and the preset yield function, the initialized variable array is calculated to obtain a deviatoric stress tensor; According to the density variation rate correction term, the velocity gradient correction term, the torsion rate correction term, the strain rate correction term, the artificial viscosity correction term, the virtual stress correction term, the extended smooth particle dynamics correction term, the normal contact force set, the tangential contact force set, the updated particle displacement set, the deviatoric stress change rate and the deviatoric stress tensor, a state parameter is generated; The state updating module is configured to perform state updating on all particles in the initialized variable array according to the state parameter, so as to obtain an updated particle variable array.

Images