A numerical simulation method and system for explosion and shock fluid-solid coupling based on Riemann-ULSPH

Abstract

This invention provides a numerical simulation method and system for explosion and impact fluid-structure interaction based on Riemann-ULSPH, belonging to the field of computational mechanics. This invention employs an updated Lagrange SPH method, uniformly discretizing the fluid, solid, and interface as particles and establishing a coordinated solution mechanism. The solid deviatoric stress normal projection is converted into equivalent pressure for input into the Riemann solver. MUSCL reconstruction and PVRS approximation solvers are used to calculate interface flux. On the solid side, the stress state is updated based on objective stress rates, and the material failure process is described by combining yield criteria and asymptotic damage evolution. Cylindrical coordinate geometric additions and axisymmetric safety treatment are added to the axisymmetric treatment, and adaptive time step progression is used to achieve stable integration. This invention improves the stability and accuracy of calculations involving strong discontinuities, large deformations, and multi-material coupling, and is suitable for high-speed collisions, explosive impacts, and structural blast resistance assessment.

Description

Technical Field

[0001] This invention relates to the field of computational mechanics and numerical simulation technology, and more particularly to a numerical simulation method and system for explosion and impact fluid-structure interaction based on Riemann-ULSPH. Background Technology

[0002] Fluid-structure interaction problems under extreme dynamic loads are often accompanied by complex physical processes such as rapid propagation of shock waves, drastic changes in the medium interface, large fluid movement, large structural deformation, and material yielding, damage and failure. They are widely present in engineering scenarios such as high-speed collisions, penetrating impacts, explosive loads, underwater explosions and fluid impacts on structures.

[0003] Existing numerical simulation methods mainly include the finite element method, finite volume method, arbitrary Lagrange-Euler method, coupled Euler-Lagrange method, and smoothed particle hydrodynamics method. Among them, traditional mesh-based methods have a mature theoretical foundation for handling small deformations of structures and continuous media problems. However, in the process of large deformations, material fracture, interface separation, and topological changes caused by extreme loads, they are prone to problems such as severe mesh distortion, negative element volume, and frequent remapping, which affect computational stability and solution efficiency. Although Euler-based methods can avoid mesh distortion due to material motion, they usually require a large fixed computational domain and suffer from problems such as numerical diffusion at material interfaces and insufficient interface identification accuracy. Smoothed particle hydrodynamics method uses a meshless Lagrange particle discretization form, which can naturally track free surfaces, material interfaces, and large deformation processes, and is therefore gradually being used for extreme dynamic fluid-structure interaction calculations.

[0004] However, traditional SPH methods often employ symmetric discretization based on kernel function weighted averaging, which can easily lead to non-physical pressure oscillations, velocity noise, excessive numerical dissipation, and interfacial particle penetration near shock waves, contact discontinuities, and high-density fluid-solid interfaces. Furthermore, the commonly used artificial viscosity parameters are highly empirical; if the parameters are too small, it is difficult to fully suppress shock oscillations, while if the parameters are too large, it will cause excessive broadening of the shock wave and loss of physical details.

[0005] Some existing Riemann-SPH methods improve shock wave capture capabilities by introducing local Riemann problems between particle pairs. However, most studies focus on single fluids or weakly coupled problems. They lack a method that considers kinematic compatibility, traction continuity, and numerical stability in unifying the incorporation of fluid pressure and the intact stress state of the solid into the interface solver. Particularly in cases where the solid exhibits significant deviatoric stress, plastic yielding, and damage failure, directly matching the fluid pressure with the solid hydrostatic pressure can easily neglect the influence of the solid's normal deviatoric stress on the interface load, leading to errors in interface pressure transmission and structural response calculations. Furthermore, existing axisymmetric SPH methods also need to address issues such as cylindrical coordinate geometric additions, circumferential stress, singularities near the axis, and objective stress updates. Some methods only perform axisymmetric corrections to the fluid conservation equations, failing to integrate solid deviatoric stress, circumferential strain rate, material strength, and damage evolution into the solution process. Therefore, they struggle to accurately describe the elastoplastic large deformation and failure behavior of axisymmetric structures under explosions and high-speed impacts.

[0006] Meanwhile, fluid flux calculation, solid constitutive updating, interface coupling, geometric source term correction, and time integration are often handled independently in existing programs, which can easily lead to inconsistent variable definitions, discontinuous interface information transmission, and insufficient computational stability between different physical modules. Therefore, there is an urgent need to propose a unified axisymmetric meshless fluid-structure interaction (FSI) calculation method applicable to strongly discontinuous, large deformation, and multi-material conditions, in order to achieve coordinated solutions for fluid state, solid stress, and cross-medium interface flux, and improve the stability, accuracy, and engineering applicability of numerical simulations of extreme dynamic loads. Summary of the Invention

[0007] To address the technical problems of existing numerical methods in handling extreme dynamic load fluid-structure interaction problems, such as the tendency of traditional mesh methods to suffer from mesh distortion and frequent remapping under large deformations and topological changes, the tendency of smoothed particle hydrodynamic methods to generate non-physical pressure oscillations and particle penetration at strong discontinuities and cross-medium interfaces, the lack of solid normal deviatoric stress contribution at the fluid-structure interface leading to discontinuous load transfer, and the difficulty in coordinating the calculation of solid elastoplastic response and geometric source terms in axisymmetric problems, this paper proposes a numerical simulation method and system for explosion and impact fluid-structure interaction based on Riemann-ULSPH.

[0008] The technical means employed in this invention are as follows: A numerical simulation method for explosion and impact fluid-structure interaction based on Riemann-ULSPH includes: S1. Establish a computational model and generate initial particles: Determine the fluid region, solid structure, explosive medium, boundary region and symmetry axis position according to the problem to be simulated, input geometric parameters and material parameters, and generate fluid particles, solid particles and boundary particles in the radial-axial plane according to the preset particle spacing. Assign each particle mass, position, density, velocity, pressure, specific internal energy, sound speed, material type, initial stress, temperature, plastic strain and damage variables. S2. Establish neighborhood relationships and calculate material state: Construct a linked list search unit based on the current particle position, determine the neighboring particles in the kernel function support domain, and call the corresponding state equations according to the material type of each particle to obtain pressure and sound velocity. For solid particles, calculate hydrostatic pressure and deviatoric stress separately. S3. Identify particle pair categories and construct interface states: Based on the material types on both sides of the particle pair, classify the particle pairs into fluid-fluid, solid-solid, and fluid-solid particle pairs. For particle pairs of the same material, use spatial gradient and limiter to reconstruct the left and right states of the interface. For fluid-solid particle pairs, calculate the projection of solid deviatoric stress on the normal direction of the particle pair and construct the equivalent normal pressure of the solid. S4. Solve the local Riemann problem: Input the left and right pressures, normal velocities, densities, and sound speeds of the interface into the approximate Riemann solver to obtain the interface pressures and velocities. Calculate the density change rate, velocity change rate, and specific internal energy change rate of each particle from the interface state. S5. Update solid stress and damage: Calculate strain rate, rotation rate and circumferential strain rate based on the neighborhood velocity field of solid particles, update the deviatoric stress test state, determine whether the material has undergone plastic deformation according to the yield criterion, and perform plastic correction, temperature update and damage accumulation when yielding. Particles that reach the failure condition reduce or lose their deviatoric stress bearing capacity. S6. Apply axisymmetric correction and advance particle state: Calculate axisymmetric mass, momentum and energy additions based on particle radial position, use axisymmetric safety treatment for particles close to the axis of symmetry, and update particle position, velocity, density, specific internal energy, pressure, stress, temperature and damage variables using prediction-correction time integration method, while determining the next time step based on stability conditions. S7. Loop Calculation and Result Output: Repeat steps S2 to S6 until the set termination time is reached. During the calculation, particle field data and measurement point time history data are output at preset intervals.

[0009] Further, step S1 includes: S11. Based on the high-speed collision, explosion impact, underwater explosion, or fluid impact structural problem to be simulated, determine the fluid region, solid structure, explosion medium, boundary region, and symmetry axis position. S12. Input the geometric and material parameters of each region. The geometric parameters include the computational region, structural dimensions, and symmetry axis position. The material parameters include material type, initial velocity, initial pressure, initial density, specific internal energy, boundary conditions, particle spacing, smooth length, computation termination time, and result output interval. S13. Generate fluid particles, solid particles and boundary particles in the radial-axial plane according to the preset particle spacing, and assign each particle mass, position, density, velocity, pressure, specific internal energy, sound speed and material type. Solid particles are also assigned initial stress, temperature, plastic strain and damage variables. S14. For axisymmetric problems, particles are arranged in the radial-axial computational plane. Each particle represents a ring-shaped material unit formed by rotating around the axis of symmetry. During the calculation, the motion and deformation of the actual three-dimensional axisymmetric medium are reflected by updating the radial and axial coordinates of the particles, without the need to establish a complete three-dimensional computational model.

[0010] Further, step S2 includes: S21. Construct a linked list search unit based on the current particle position, divide the computation area into multiple search units, and store each particle in the corresponding search unit according to its spatial position. S22. For any particle to be calculated, only search for particles in its own cell and adjacent cells, and determine whether they constitute a valid particle pair based on the distance between particles and the support domain of the kernel function. S23. At each time step or at a set update time, re-establish the adjacency relationship based on the updated position of the particles to adapt to changes in the relative position of particles after large fluid movement, large structural deformation, interface movement, and material breakage. S24. Based on the material type of each particle, call the corresponding equation of state to obtain the pressure and speed of sound; S25. For fluid particles, thermodynamic pressure is obtained based on the fluid state equation. S26. For metallic solid particles, the stress tensor is decomposed into hydrostatic pressure and deviatoric stress components, where hydrostatic pressure is calculated by the solid state equation and deviatoric stress is updated by subsequent solid strength and damage calculations. S27. For explosion products or gaseous particles, determine their pressure and sound velocity according to the gas state equation or the detonation product state equation. S28. Each material's equation of state is set up with an independent interface and called through the material type identifier, enabling the system to extend material models for different fluids, metals, gases, and explosive media.

[0011] Further, step S3 includes: S31. Based on the material types on both sides of the particle pair, the particle pairs are classified into fluid-fluid particle pairs, solid-solid particle pairs and fluid-solid particle pairs. Different types of particle pairs adopt different state construction methods. S32. For a pair of particles of the same material, the reconstructed state at the particle pair interface is obtained through a restricted linear reconstruction method based on the local pressure, density, velocity, specific internal energy and spatial gradient of the particles. S33. Limit the reconstruction slope and make a physical validity judgment on the reconstructed density, pressure, specific internal energy and sound speed. When the reconstruction state exceeds the preset physical range, the corresponding particle pair will be degenerated into a low-order reconstruction or the local original state of the particle will be adopted. S34. For fluid-solid cross-material particle pairs, instead of directly performing cross-medium high-order reconstruction, the Riemann initial value problem is constructed using the local states of particles on both sides of the interface. S35. For a fluid-solid particle pair, determine the direction of interaction between the two particles based on their position difference, and use this direction as the normal direction for the local Riemann problem. The normal direction is:

[0012] in, This represents the unit normal vector of the particle to the Riemann solution direction. Indicates the right-hand particle The position vector, Indicates the particle on the left The position vector; S36. The fluid-solid interface in the continuous sense satisfies the normal kinematic compatibility condition and the dynamic traction continuity condition, as follows:

[0013]

[0014] in, Indicates the velocity at the fluid interface. Indicates the solid-solid interface velocity. Represents the unit normal vector of the interface. This represents the Cauchy stress tensor of the fluid. Represents the Cauchy stress tensor of a solid. S37. For inviscid fluids, their Cauchy stress tensor degenerates into... ,in, Represents the thermodynamic pressure of a fluid. Represents the unit tensor; the fluid-side normal traction is... Solid side normal traction is The solid side normal component is ,in, Represents the hydrostatic pressure of a solid. Represents the deviatoric stress tensor of a solid; S38, Based on traction continuity conditions The equivalent normal pressure entering the Riemann solver from the solid side is defined as:

[0015] S39. Because the direction of the particle's Riemann solution is located... In the calculation plane, i.e. ,exist The plane unfolds as follows:

[0016] in, Represents the radial deviatoric stress component. Represents the shear deviatoric stress component. Indicates the axial deviatoric stress component. Indicates the normal direction In the unit component of the radial direction, Indicates the normal direction Unit component in the axial direction; S310. Based on the equivalent pressure definition, the left and right normal pressure states of the fluid-solid particle pair are transformed into a standard form one-dimensional Riemann initial value problem. For any fluid-solid particle pair, the normal pressure entering the Riemann solver is defined as:

[0017]

[0018] in, This indicates the normal pressure entering the Riemann solver from the left. This indicates the normal pressure entering the Riemann solver from the right side. This represents the hydrostatic pressure of the particle on the left. This represents the hydrostatic pressure of the particle on the right. Represents the fluid domain. Represents a solid domain. This represents the deviatoric stress tensor of the solid particles on the left. This represents the deviatoric stress tensor of the solid particles on the right. S311. Calculate the normal velocities on the left and right sides of the fluid-solid interface, using the following formula:

[0019]

[0020] in, Indicates the normal velocity on the left side. Indicates the normal velocity on the right side. Indicates the particle on the left The velocity vector, Indicates the right-hand particle The velocity vector.

[0021] Further, step S4 includes: S41. Input the left and right interface pressures, normal velocities, densities, and sound velocities into the PVRS approximation Riemann solver based on acoustic impedance weighting to calculate the contact pressures and normal velocities, as shown in the following formulas:

[0022]

[0023] in, Indicates the acoustic impedance on the left side. , This represents the particle density on the left. Indicates the sound speed of the particle on the left; Indicates the acoustic impedance on the right side. , This represents the particle density on the right. Indicates the sound speed of the particle on the right; S42. The acoustic impedance of the solid-side PVRS is constructed using the local volume sound velocity given by the Mie-Grüneisen equation of state, as shown in the following formula:

[0024] in, Represents solid particles The speed of sound, Indicates the volume compression ratio. , , All are Hugoniot polynomial coefficients. Indicates the reference density. The reference sound velocity is indicated; the contribution of deviatoric stress is not incorporated into the sound velocity, but is explicitly calculated through the constitutive update of the solid and the deviatoric stress flux in the momentum equation. S43. Calculate the rate of change of density, rate of change of velocity, and rate of change of specific internal energy of each particle based on the interface state. The normal pressure between the fluid and solid intermaterial particle pairs is incorporated into the momentum equation in the following form:

[0025] in, Represents particles The rate of change of velocity, Represents particles density, Indicates the interface contact pressure. Representing kernel function Regarding particles Gradient of position, Representing neighboring particles Volume; S44. The intact stress flux within a solid enters the momentum equation in the following form:

[0026] in, This indicates interfacial deviatoric stress.

[0027] Further, step S5 includes: S51. Calculate the radial strain rate, axial strain rate, shear strain rate, circumferential strain rate and rotation rate based on the neighborhood velocity field of solid particles, and obtain the experimental value of solid deviatoric stress using the Jaumann objective stress update relationship. S52. Calculate the equivalent stress based on the deviatoric stress tensor and compare it with the current material yield strength; when the equivalent stress does not exceed the material yield strength, the solid particles remain in an elastic state, and the deviatoric stress components do not need correction; when the equivalent stress exceeds the material yield strength, the particles... Upon entering the plastic yield state, the deviatoric stress components must proportionally recede back to the yield surface. The deviatoric stress components are corrected as follows:

[0028] in, This represents the corrected deviatoric stress tensor; Represents the deviatoric stress test tensor; The yield strength of the material is represented by the Johnson-Cook yield model. , Indicates static yield strength. and All of these represent strain hardening parameters. Indicates the strain rate strengthening parameter. Indicates the thermal softening parameter. For equivalent plastic strain, Represents the equivalent plastic strain rate. Indicates the reference strain rate. Represents dimensionless temperature. , Melting point Room temperature; Indicates the equivalent Mises stress; This represents the second invariant of the deviatoric stress tensor. , Represents the matrix trace operation. This represents the components of the deviatoric stress tensor in cylindrical coordinates. S53. The adiabatic temperature rise is calculated using the Taylor-Quinney coefficient and plastic work, as shown in the following formula:

[0029]

[0030] in, This represents the Taylor-Quinney coefficient; This represents the equivalent plastic strain increment within a single time step, i.e., the amount of plastic deformation that occurs in the material during this step, based on the von Mises yield criterion and Plasticity theory , This represents the yield strength at the current moment. Indicates shear modulus; Indicates constant volume specific heat capacity; This represents the temperature rise within a single time step; S54. Update temperature using a prediction-correction time integral format:

[0031]

[0032] The starting point is indicated by a superscript. The midpoint is indicated by a superscript. This indicates that the new time step is indicated by a superscript. express; Indicates based on The trial temperature rise calculated based on the instantaneous state and half-step size; Indicates based on High-precision flow field and full-step flow at any time The calculated actual temperature rise.

[0033] S55. The yield strength of the material is determined based on at least one of the state variables: plastic strain, strain rate, and temperature. Part of the plastic work generated by plastic deformation is converted into heat energy, which is used to update the temperature of solid particles and feeds the updated temperature back to the material thermal softening calculation. S56. Damage variables are accumulated based on the ratio of the single-step plastic strain increment to the current failure strain. When the damage variable reaches the preset failure threshold, the corresponding particle is marked as a failed particle. In subsequent time steps, the deviatoric stress and material strength contribution of the particle are reduced or eliminated, but at least a portion of its mass, inertia, and pressure state are retained to simulate the motion process after material failure, fragmentation, and loss of load-bearing capacity.

[0034] Further, step S6 includes: S61. Add cylindrical coordinate geometric additions related to the radial position of the particle to the mass conservation, momentum conservation and energy conservation relationships, and add circumferential strain rate and circumferential deviatoric stress to the solid stress update, so that the radial stress, axial stress, shear stress and circumferential stress are updated in a coordinated manner within the same time step. S62. Set an axial safety zone. When the particle is outside the axial safety zone, the circumferential strain rate is calculated using the ratio of radial velocity to radial coordinate. When the particle is inside the axial safety zone, the circumferential strain rate is calculated using an axial approximation method based on symmetry conditions, velocity gradient, or limit values ​​to avoid numerical divergence near the axis. S63. In the prediction phase, based on the particle position, velocity, density, specific internal energy, pressure, and stress change rate at the current moment, the predicted state at the intermediate moment is obtained. S64. Re-perform neighborhood search, material state calculation, interface identification, Riemann flux solution, and solid constitutive update using the predicted state; S65. During the correction phase, the particle position, velocity, density, specific internal energy, pressure, stress, temperature, and damage variables for the next moment are obtained based on the recalculated rate of change. S66. The time step is limited based on the particle spacing, smooth length, sound speed, particle velocity, particle acceleration, and the degree of material state change. When high sound speed, high velocity, strong pressure gradient, or severe structural deformation occurs in the computational domain, the time step is reduced. When the particle state change is relatively gentle, a larger time step is allowed.

[0035] Further, step S7 includes: S71. Output particle position, velocity, density, pressure, specific internal energy, temperature, material type, stress components, equivalent stress, equivalent plastic strain, yield state and damage variables at preset time intervals. S72. Set up virtual measuring points at specified spatial locations and record the pressure, velocity, density, and structural response near the measuring points in chronological order. S73. Extract data on shock wave propagation, bubble boundary motion, free surface changes, structural deformation, and material failure processes based on the spatial positions of fluid particles, gas particles, solid particles, or material interface particles.

[0036] This invention also provides a Riemann-ULSPH-based numerical simulation system for explosion and impact fluid-structure interaction, implemented based on the aforementioned Riemann-ULSPH-based numerical simulation method. Deployed in a computing device, the computing device includes a processor, memory, data input interface, and result output interface, comprising: a model input and particle initialization module, a particle neighborhood search module, a material state calculation module, a spatial state reconstruction module, a fluid-structure interface processing module, a Riemann flux solution module, a solid strength and damage calculation module, an axisymmetric geometric correction module, a time integration module, and a result output module, wherein: The model input and particle initialization module is used to receive parameters such as computational region, structural size, material type, initial velocity, initial pressure, initial density, specific internal energy, boundary conditions, particle spacing, smooth length, computation termination time, and result output interval. It discretizes the computational region into Lagrange particles carrying mass, position, velocity, density, pressure, specific internal energy, sound speed, stress tensor, temperature, equivalent plastic strain, and damage variables. Each particle is assigned a material type identifier to distinguish between fluid particles, solid particles, detonation product particles, boundary particles, and failure particles. For axisymmetric problems, the particles are arranged in the radial-axial computational plane, and each particle represents a ring-shaped material unit formed by rotating around the axis of symmetry. The particle neighborhood search module is used to establish the interaction relationship between particles based on the current position of the particles, divide the computation region into multiple search units, and determine whether a valid particle pair is formed based on the distance between particles and the support domain of the kernel function. The adjacency relationship is re-established at each time step or at a set update time. The material state calculation module is used to call the corresponding state relationship based on the material type of the particle to calculate the pressure, sound velocity and temperature state variables from the particle density and specific internal energy. For solid particles, hydrostatic pressure and deviatoric stress are calculated separately. The spatial state reconstruction module is used to construct the left and right states on both sides of the particle interface between particle pairs of the same material. The fluid-solid interface processing module is used to identify cross-medium particle pairs formed by fluid particles and solid particles, and to convert the projection of solid deviatoric stress onto the normal direction of the particle pair into an equivalent normal pressure. This allows the fluid pressure and the total solid stress to be processed in the same local Riemann problem. For particle pairs of the same material, the reconstructed state at the interface is obtained through a restricted linear reconstruction method. For fluid-solid cross-material particle pairs, the local state of the particles on both sides of the interface is used to construct the Riemann initial value problem. The Riemann flux solving module is used to establish a local one-dimensional Riemann initial value problem between each effective particle pair. It calculates the interface star region state based on the pressure, velocity, density and sound speed on both sides of the particle pair. It weights the pressure and normal velocity according to the acoustic impedance on both sides to obtain the interface pressure and interface normal velocity, which act on both sides of the particle pair in pairs and opposite directions to update the mass, momentum and energy change rate of the particles. The solid strength and damage calculation module is used to update the deviatoric stress, equivalent stress, yield state, equivalent plastic strain, temperature and damage variables of solid particles. It calculates the strain rate based on the neighborhood velocity field of solid particles and obtains the deviatoric stress trial value using the Jaumann objective stress update relationship. It performs plastic correction based on the von Mises yield criterion and the Johnson-Cook constitutive model, calculates the adiabatic temperature rise using Taylor-Quinney coefficients, and accumulates the damage variables based on the plastic strain increment. The axisymmetric geometry correction module is used to add cylindrical coordinate geometry additions related to the radial position of the particle to the mass conservation, momentum conservation and energy conservation relationships, add circumferential strain rate and circumferential deviatoric stress to the solid stress update, and adopt axisymmetric safety treatment for particles close to the axis of symmetry. The time integration module is used to update the state variables of each particle according to the explicit prediction-correction method. In the prediction stage, it obtains the predicted state at the intermediate time and recalculates each rate of change. In the correction stage, it obtains the particle state at the next time and determines the next time step based on the stability condition. The results output module is used to output particle state data at preset time intervals, set virtual measuring points at specified spatial locations and record time history data, and extract data on shock wave propagation, structural deformation and material failure processes. The data connections between the model input and particle initialization module, particle neighborhood search module, material state calculation module, spatial state reconstruction module, fluid-solid interface processing module, Riemann flux solution module, solid strength and damage calculation module, axisymmetric geometry correction module, time integration module, and result output module are as follows: The data output terminal of the model input and particle initialization module is connected to the data input terminals of the particle neighborhood search module, the material state calculation module and the time integration module, and is used to provide initial particle states and calculation control parameters to each calculation module; The output of the particle neighborhood search module is connected to the fluid-solid interface processing module and the solid strength and damage calculation module, respectively, and is used to provide the neighbor list, inter-particle distance, interaction direction and kernel function gradient of each particle. The output of the material state calculation module is connected to the fluid-solid interface processing module and the Riemann flux solving module, and is used to provide particle pressure, density, specific internal energy and sound velocity. The output of the spatial state reconstruction module is connected to the Riemann flux solving module to provide the left and right reconstruction states of the same material particle pairs. The output of the fluid-solid interface processing module is connected to the Riemann flux solving module, and is used to provide the left and right reconstruction states of the same material particle pairs, as well as the fluid normal pressure state and solid equivalent normal pressure state across material particle pairs. The output of the Riemann flux solving module is connected to the time integration module to provide particle mass, momentum, and energy change rate; the outputs of the solid strength and damage calculation module and the axisymmetric geometry correction module are also connected to the time integration module to provide solid stress change rate, plastic state, damage state, and axisymmetric additional terms. The updated particle state from the time integration module is re-input into the particle neighborhood search module and the material state calculation module, thereby forming a data closed loop for time-step cyclic solution; when the output condition is met, the time integration module transmits the current particle state to the result output module.

[0037] Furthermore, the system also includes a parallel computing module, which is used to distribute tasks such as particle initialization, neighborhood search, particle pair interaction, material state update, solid constitutive update, and result statistics to multiple processor threads; independent particle loops are executed in parallel using threads, and particle pair calculations that may be written to the same particle variable at the same time are protected by using thread-private variables, local accumulation, atomic operations, or partitioned calculations to reduce data contention in parallel computing.

[0038] Compared with the prior art, the present invention has the following advantages: 1. The present invention provides a numerical simulation method for explosion and impact fluid-structure interaction based on Riemann-ULSPH. By converting the projection of solid deviatoric stress on the interface normal direction into equivalent normal pressure, the fluid pressure and the total stress of the solid can be handled in the same local Riemann problem, realizing accurate calculation of cross-medium interface traction continuity, momentum exchange and energy transfer.

[0039] 2. The present invention provides a numerical simulation method for explosion and shock fluid-structure interaction based on Riemann-ULSPH, which achieves a balance between strong shock wave resolution and numerical stability at the cross-medium interface by reconstructing the higher-order state of particle pairs of the same material and by treating the local state and equivalent traction of cross-material particle pairs of fluid-structure interaction.

[0040] 3. The present invention provides a numerical simulation method for explosion and impact fluid-structure interaction based on Riemann-ULSPH, which achieves the preservation of the solid's internal elastoplastic bearing capacity by retaining the complete deviatoric stress flux inside the solid and adopting a normal traction equivalent method at the fluid-structure interface, thus avoiding the inappropriate application of the solid shear response to the inviscid fluid.

[0041] 4. The present invention provides a numerical simulation method for explosion and impact fluid-structure interaction based on Riemann-ULSPH, which unifies radial, axial, shear and circumferential stresses into the solid constitutive update and adds axisymmetric geometric corrections to the calculation of mass, momentum and energy, thereby realizing the description of three-dimensional axisymmetric extreme dynamic processes using two-dimensional computational scale.

[0042] 5. The present invention provides a numerical simulation method for explosion and impact fluid-structure interaction based on Riemann-ULSPH. By setting an axial safety region and using an axial approximation method to calculate the circumferential strain rate for particles near the axis of symmetry, the risk of geometric additions being amplified when the radial coordinate approaches zero is reduced, and the stability of particle calculation near the axis of symmetry is improved.

[0043] 6. The present invention provides a numerical simulation method for explosion and impact fluid-structure interaction based on Riemann-ULSPH. By using meshless updated Lagrange particles to describe large fluid motion, large structural deformation, interface movement and material fragmentation, it can be applied to extreme dynamic problems with continuously changing topological relationships without the need for traditional mesh reconstruction.

[0044] 7. The present invention provides a numerical simulation method for explosion and impact fluid-structure interaction based on Riemann-ULSPH, which improves the consistency of multi-material fluid-structure interaction calculations by setting Riemann flux calculation, solid elastoplasticity update, material damage failure, dynamic time step and parallel neighborhood search in a unified time-progression process.

[0045] In summary, by applying the technical solution of this invention, existing technologies suffer from several drawbacks. Traditional meshing methods are prone to mesh distortion and frequent remapping under large deformations and topological changes; traditional SPH methods are prone to non-physical pressure oscillations and particle penetration at strong discontinuities and cross-medium interfaces; the lack of solid normal deviatoric stress contribution at the fluid-solid interface leads to discontinuous load transfer; and the solid elastoplastic response and geometric source terms are difficult to coordinate in axisymmetric problems. This invention establishes a unified axisymmetric Riemann-ULSPH meshless particle computation framework, unifying fluid pressure solution, solid elastoplastic response, fluid-solid interface flux exchange, axisymmetric geometric correction, material damage and failure, and time progression within this framework. Therefore, the technical solution of this invention solves the problems of insufficient computational stability, discontinuous pressure transfer at the fluid-solid interface, significant non-physical oscillations near shock waves, incomplete handling of solid deviatoric stress contribution, difficulty in coordinating axisymmetric geometric terms and material strength response, and mesh quality degradation after large structural deformations in existing technologies.

[0046] Based on the above reasons, this invention can be widely applied in fields such as high-speed collision, penetration impact, explosive load, underwater explosion, strong shock wave propagation, fluid impact structure, aerospace protection, explosion-proof design of ship and marine engineering structures, and explosion-proof and impact-resistant assessment of major engineering facilities. Attached Figure Description

[0047] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0048] Figure 1 This is a flowchart of the method of the present invention.

[0049] Figure 2 This is a flowchart illustrating the solid stress, plasticity, and damage update process of this invention.

[0050] Figure 3 This is a functional block diagram of the fluid-structure interaction numerical simulation system of the present invention. Detailed Implementation

[0051] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0052] It should be noted that the terms "comprising" and "having" and any variations thereof in the specification, claims and accompanying drawings of this invention are intended to cover non-exclusive inclusion. For example, a process, method, system, product or device that includes 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 process, method, product or device.

[0053] like Figure 1 As shown, this invention provides a numerical simulation method for explosion and impact fluid-structure interaction based on Riemann-ULSPH, including: S1. Establish the computational model: Determine the fluid region, solid structure, explosion medium, boundary region, and axis of symmetry position according to the problem to be simulated, and input the geometric parameters and material parameters; S2. Generate initial particles: Generate fluid particles, solid particles and boundary particles in the radial-axial plane according to the preset particle spacing, and assign each particle mass, position, density, velocity, pressure, specific internal energy, sound speed, material type, initial stress, temperature, plastic strain and damage variables. S3. Establish neighborhood relationships: Construct a linked list search unit based on the current particle position to determine neighboring particles within the kernel function support domain; as fluid motion and structural deformation occur, the relative positions of the particles change, and the neighborhood relationships are updated accordingly in subsequent time steps; S4. Calculate the material state: According to the material type of each particle, call the corresponding state equation to obtain the pressure and sound velocity. For solid particles, calculate the hydrostatic pressure and deviatoric stress separately to provide state variables for subsequent interface flux and material strength calculations. S5. Identify particle pair categories: Based on the material types on both sides of the particle pair, the particle pairs are classified into fluid-fluid, solid-solid, and fluid-solid particle pairs. Different categories of particle pairs adopt different state construction methods. S6. Constructing the left and right states of the interface: For particle pairs of the same material, the left and right states of the interface are reconstructed using spatial gradient and constraint. For fluid-solid particles, cross-material high-order reconstruction is not performed. Instead, the projection of the solid deviatoric stress onto the normal direction of the particle pair is calculated, and the equivalent normal pressure of the solid is constructed. S7. Solve the local Riemann problem: Input the left and right pressures, normal velocities, densities, and sound speeds of the interface into the approximate Riemann solver to obtain the interface pressures and velocities. Calculate the density change rate, velocity change rate, and specific internal energy change rate of each particle from the interface state. S8. Update solid stress and damage: Calculate strain rate, rotation rate and circumferential strain rate based on the neighborhood velocity field of solid particles, update the deviatoric stress test state, determine whether the material has undergone plastic deformation according to the yield criterion, and perform plastic correction, temperature update and damage accumulation when yielding. Particles that reach the failure condition reduce or lose the deviatoric stress bearing capacity. S9. Apply axisymmetric correction: Calculate the axisymmetric mass, momentum, and energy additions based on the particle's radial position. Apply axisymmetric safety treatment to particles close to the axis of symmetry to avoid numerical anomalies caused by the radial coordinates approaching zero. S10, Propulsion Particle State: The particle position, velocity, density, specific internal energy, pressure, stress, temperature and damage variables are updated using a prediction-correction time integration method, while the next time step is determined based on stability conditions; During the above update process, fluid particles migrate with the velocity field, and shock waves propagate between particles through pressure and density changes; the position of the fluid-solid interface is determined by the real-time coordinates of particles near the interface and changes automatically with fluid motion and solid deformation; the shape changes of the solid structure, such as indentation, bending, stretching, perforation and fragmentation, are represented by the overall evolution of the solid particle positions, without the need for re-meshing.

[0054] S11. Loop Calculation and Result Output: Repeat steps S2 to S10 until the set termination time is reached. During the calculation, particle field data and measurement point time history data are output at preset intervals.

[0055] In a specific implementation, as a preferred embodiment of the present invention, step S1 includes: S11. Based on the high-speed collision, explosion impact, underwater explosion, or fluid impact structural problem to be simulated, determine the fluid region, solid structure, explosion medium, boundary region, and symmetry axis position. S12. Input the geometric and material parameters of each region. The geometric parameters include the calculation region, structural dimensions, and symmetry axis position. The material parameters include material type, initial velocity, initial pressure, initial density, specific internal energy, boundary conditions, particle spacing, smooth length, calculation termination time, and result output interval.

[0056] In a specific implementation, as a preferred embodiment of the present invention, step S2 includes: S21. Generate fluid particles, solid particles and boundary particles in the radial-axial plane according to the preset particle spacing, and assign each particle mass, position, density, velocity, pressure, specific internal energy, sound speed and material type. Solid particles are also assigned initial stress, temperature, plastic strain and damage variables. S22. For axisymmetric problems, particles are arranged in a radial-axial computational plane. Each particle represents a ring-shaped material unit formed by rotating around the axis of symmetry. During the calculation, the motion and deformation of the actual three-dimensional axisymmetric medium are reflected by updating the radial and axial coordinates of the particles, without the need to establish a complete three-dimensional computational model.

[0057] In a specific implementation, as a preferred embodiment of the present invention, step S3 includes: S31. Construct a linked list search unit based on the current particle position, divide the computation area into multiple search units, and store each particle in the corresponding search unit according to its spatial position; S32. For any particle to be calculated, only search for particles in its own cell and adjacent cells, and determine whether they constitute a valid particle pair based on the distance between particles and the support domain of the kernel function. S33. At each time step or at a set update time, the adjacency relationship is re-established based on the updated position of the particles to adapt to changes in the relative position of particles after large fluid movement, large structural deformation, interface movement, and material breakage. In a specific implementation, as a preferred embodiment of the present invention, step S4 includes: S41. Based on the material type of each particle, call the corresponding equation of state to obtain the pressure and speed of sound; S42. For fluid particles, thermodynamic pressure is obtained based on the fluid state equation. S43. For metallic solid particles, the stress tensor is decomposed into hydrostatic pressure and deviatoric stress components, where hydrostatic pressure is calculated by the solid state equation and deviatoric stress is updated by subsequent solid strength and damage calculations. S44. For explosion products or gaseous particles, determine their pressure and sound velocity according to the gas state equation or the detonation product state equation. S45. Each material's equation of state is set up with an independent interface and called through the material type identifier, enabling the system to extend material models for different fluids, metals, gases, and explosive media.

[0058] In a specific implementation, as a preferred embodiment of the present invention, step S5 includes: Based on the material types on both sides of the particle pair, the particle pairs are classified into fluid-fluid particle pairs, solid-solid particle pairs, and fluid-solid particle pairs. Different types of particle pairs adopt different state construction methods. In a specific implementation, as a preferred embodiment of the present invention, step S6 includes: S61. For a pair of particles of the same material, the reconstructed state at the particle pair interface is obtained through a restricted linear reconstruction method based on the local pressure, density, velocity, specific internal energy and spatial gradient of the particles. S62. Limit the reconstruction slope and make a physical validity judgment on the reconstructed density, pressure, specific internal energy and sound speed. When the reconstruction state exceeds the preset physical range, the corresponding particle pair will be degenerated into a low-order reconstruction or the local original state of the particle will be adopted. S63. For fluid-solid cross-material particle pairs, instead of directly performing high-order reconstruction across the medium, the Riemann initial value problem is constructed using the local states of particles on both sides of the interface. S64. For a fluid-solid particle pair, determine the direction of interaction between the two particles based on their position difference, and use this direction as the normal direction for the local Riemann problem. The normal direction is:

[0059] in, This represents the unit normal vector of the particle to the Riemann solution direction. Indicates the right-hand particle The position vector, Indicates the particle on the left The position vector; S65. The fluid-solid interface in the continuous sense satisfies the normal kinematic compatibility condition and the dynamic traction continuity condition, as follows:

[0060]

[0061] in, Indicates the velocity at the fluid interface. Indicates the solid-solid interface velocity. Represents the unit normal vector of the interface. This represents the Cauchy stress tensor of the fluid. Represents the Cauchy stress tensor of a solid. S66. For inviscid fluids, their Cauchy stress tensor degenerates into... ,in, Represents the thermodynamic pressure of a fluid. Represents the unit tensor; the fluid-side normal traction is... Solid side normal traction is The solid side normal component is ,in, Represents the hydrostatic pressure of a solid. Represents the deviatoric stress tensor of a solid; S67, Based on traction continuity conditions The equivalent normal pressure entering the Riemann solver from the solid side is defined as:

[0062] S68. Because the particle's Riemann solution direction is located... In the calculation plane, i.e. ,therefore It does not directly enter the particle pair normal projection, but it still affects the solid response through axisymmetric geometric source terms and constitutive updates. The plane unfolds as follows:

[0063] in, Represents the radial deviatoric stress component. Represents the shear deviatoric stress component. Indicates the axial deviatoric stress component. Indicates the normal direction In the unit component of the radial direction, Indicates the normal direction Unit component in the axial direction; S69. Based on the equivalent pressure definition, the left and right normal pressure states of the fluid-solid particle pair are transformed into a standard form one-dimensional Riemann initial value problem. For any fluid-solid particle pair, the normal pressure entering the Riemann solver is defined as:

[0064]

[0065] in, This indicates the normal pressure entering the Riemann solver from the left. This indicates the normal pressure entering the Riemann solver from the right side. This represents the hydrostatic pressure of the particle on the left. This represents the hydrostatic pressure of the particle on the right. Represents the fluid domain. Represents a solid domain. This represents the deviatoric stress tensor of the solid particles on the left. This represents the deviatoric stress tensor of the solid particles on the right. S610. Calculate the normal velocities on the left and right sides of the fluid-solid interface, using the following formula:

[0066]

[0067] in, Indicates the normal velocity on the left side. Indicates the normal velocity on the right side. Indicates the particle on the left The velocity vector, Indicates the right-hand particle The velocity vector.

[0068] In a specific implementation, as a preferred embodiment of the present invention, step S7 includes: S71. At the cross-medium interface, to avoid non-physical oscillations caused by abrupt changes in density, sound velocity, and the equation of state, cross-material MUSCL reconstruction is not performed. Instead, the local first-order states of the left and right particles are used to enter the Riemann solver. For example... Figure 2 As shown, the contact pressure, normal velocity, density, and sound velocity at the interface are input into the PVRS approximation Riemann solver based on acoustic impedance weighting to calculate the contact pressure and normal velocity, as shown in the following formula:

[0069]

[0070] in, Indicates the acoustic impedance on the left side. , This represents the particle density on the left. Indicates the sound speed of the particle on the left; Indicates the acoustic impedance on the right side. , This represents the particle density on the right. Indicates the sound speed of the particle on the right; S72. The acoustic impedance of the solid-side PVRS is constructed using the local volume sound velocity given by the Mie-Grüneisen equation of state, as shown in the following formula:

[0071] in, Represents solid particles The speed of sound, Indicates the volume compression ratio. , , All are Hugoniot polynomial coefficients. Indicates the reference density. The reference sound velocity is indicated; the contribution of deviatoric stress is not incorporated into the sound velocity, but is explicitly calculated through the constitutive update of the solid and the deviatoric stress flux in the momentum equation. S73. Calculate the rate of change of density, rate of change of velocity, and rate of change of specific internal energy of each particle based on the interface state. The normal pressure between the fluid and solid intermaterial particle pairs is incorporated into the momentum equation in the following form:

[0072] in, Represents particles The rate of change of velocity, Represents particles density, Indicates the interface contact pressure. Representing kernel function Regarding particles Gradient of position, Representing neighboring particles Volume; S74. Unlike fluid-solid intermedium interactions, the complete elastoplastic deviatoric stress contribution is preserved between particle pairs within a solid. The complete stress flux within the solid enters the momentum equation in the following form:

[0073] in, This indicates interfacial deviatoric stress.

[0074] Through the above processing, the pulling effect of the total stress of the solid on the normal direction of the particle pair is transformed into an equivalent pressure state, so that the influence of the solid deviatoric stress on the normal load of the fluid-solid interface can be included in the interface flux calculation, instead of just matching the hydrostatic pressure of the solid with the fluid pressure.

[0075] In a specific implementation, as a preferred embodiment of the present invention, step S8 includes: S81. Calculate the radial strain rate, axial strain rate, shear strain rate, circumferential strain rate and rotation rate based on the neighborhood velocity field of solid particles, and obtain the experimental value of solid deviatoric stress using the Jaumann objective stress update relationship. S82. Calculate the equivalent stress based on the deviatoric stress tensor and compare it with the current material yield strength; when the equivalent stress does not exceed the material yield strength, the solid particles remain in an elastic state, and the deviatoric stress components do not need correction; when the equivalent stress exceeds the material yield strength, the particles... Upon entering the plastic yield state, the deviatoric stress components must proportionally recede back to the yield surface. The deviatoric stress components are corrected as follows:

[0076] in, This represents the corrected deviatoric stress tensor; Represents the deviatoric stress test tensor; The yield strength of the material is represented by the Johnson-Cook yield model. , Indicates static yield strength. and All of these represent strain hardening parameters. Indicates the strain rate strengthening parameter. Indicates the thermal softening parameter. For equivalent plastic strain, Represents the equivalent plastic strain rate. Indicates the reference strain rate. Represents dimensionless temperature. , Melting point Room temperature; Indicates the equivalent Mises stress; This represents the second invariant of the deviatoric stress tensor. , Represents the matrix trace operation. This represents the components of the deviatoric stress tensor in cylindrical coordinates. S83. High-speed impact or explosive impact processes are considered adiabatic processes, neglecting heat conduction and radiation. The adiabatic temperature rise is calculated using the Taylor-Quinney coefficient and plastic work, as shown in the following formula:

[0077]

[0078] in, This represents the Taylor-Quinney coefficient, which is set to 0.9. This represents the equivalent plastic strain increment within a single time step, i.e., the amount of plastic deformation that occurs in the material during this step, based on the von Mises yield criterion and Plasticity theory , This represents the yield strength at the current moment. Indicates shear modulus; Indicates constant volume specific heat capacity; This represents the temperature rise within a single time step; S84. Update the temperature using the predictor-corrector time integral format:

[0079]

[0080] The starting point is indicated by a superscript. The midpoint is indicated by a superscript. This indicates that the new time step is indicated by a superscript. express; Indicates based on The trial temperature rise calculated based on the instantaneous state and half-step size; Indicates based on High-precision flow field and full-step flow at any time The calculated actual temperature rise.

[0081] S85. The yield strength of the material is determined based on at least one of the state variables: plastic strain, strain rate, and temperature. Part of the plastic work generated by plastic deformation is converted into heat energy, which is used to update the temperature of solid particles and feeds the updated temperature back to the material thermal softening calculation. S86. Damage variables are accumulated based on the ratio of the single-step plastic strain increment to the current failure strain. When the damage variable reaches the preset failure threshold, the corresponding particle is marked as a failed particle. In subsequent time steps, the deviatoric stress and material strength contribution of the particle are reduced or eliminated, but at least a portion of its mass, inertia, and pressure state are retained to simulate the motion process after material failure, fragmentation, and loss of load-bearing capacity.

[0082] In a specific implementation, as a preferred embodiment of the present invention, step S9 includes: S91. Add cylindrical coordinate geometric additions related to the radial position of the particle to the mass conservation, momentum conservation and energy conservation relationships, and add circumferential strain rate and circumferential deviatoric stress to the solid stress update, so that the radial stress, axial stress, shear stress and circumferential stress are updated in a coordinated manner within the same time step. S92. Set an axial safety zone. When the particle is outside the axial safety zone, the circumferential strain rate is calculated by the ratio of radial velocity to radial coordinate. When the particle is inside the axial safety zone, the circumferential strain rate is calculated by an axial approximation method based on symmetry conditions, velocity gradient, or limit values ​​to avoid numerical divergence near the axis.

[0083] In a specific implementation, as a preferred embodiment of the present invention, step S10 includes: S101. In the prediction phase, based on the particle position, velocity, density, specific internal energy, pressure, and stress change rate at the current moment, the predicted state at the intermediate moment is obtained. S102. Re-perform neighborhood search, material state calculation, interface identification, Riemann flux solution, and solid constitutive update using the predicted state; S103. During the correction phase, the particle position, velocity, density, specific internal energy, pressure, stress, temperature, and damage variables for the next moment are obtained based on the recalculated rate of change. S104. The time step is limited based on the particle spacing, smooth length, sound speed, particle velocity, particle acceleration, and the degree of material state change. When high sound speed, high velocity, strong pressure gradient, or severe structural deformation occurs in the computational domain, the time step is reduced. When the particle state change is relatively gentle, a larger time step is allowed.

[0084] In a specific implementation, as a preferred embodiment of the present invention, step S11 includes: S111, Output particle position, velocity, density, pressure, specific internal energy, temperature, material type, stress components, equivalent stress, equivalent plastic strain, yield state and damage variables at preset time intervals; S112. Set up virtual measuring points at specified spatial locations and record the pressure, velocity, density, and structural response near the measuring points in chronological order. S113. Extract data on shock wave propagation, bubble boundary motion, free surface changes, structural deformation, and material failure processes based on the spatial positions of fluid particles, gas particles, solid particles, or material interface particles.

[0085] This invention also provides a Riemann-ULSPH-based numerical simulation system for explosion and impact fluid-structure interaction, implemented using the aforementioned Riemann-ULSPH-based numerical simulation method. This system is deployed in a computing device, which includes a processor, memory, a data input interface, and a result output interface. Figure 3 As shown, it includes: a model input and particle initialization module, a particle neighborhood search module, a material state calculation module, a spatial state reconstruction module, a fluid-solid interface processing module, a Riemann flux solution module, a solid strength and damage calculation module, an axisymmetric geometric correction module, a time integration module, and a result output module, wherein: The model input and particle initialization module receives parameters such as computational region, structural dimensions, material type, initial velocity, initial pressure, initial density, specific internal energy, boundary conditions, particle spacing, smoothing length, computation termination time, and result output interval. It discretizes the computational region into Lagrange particles carrying mass, position, velocity, density, pressure, specific internal energy, sound velocity, stress tensor, temperature, equivalent plastic strain, and damage variables. Each particle is assigned a material type identifier to distinguish between fluid particles, solid particles, detonation product particles, boundary particles, and failure particles. For axisymmetric problems, the particles are arranged in a radial-axial computational plane, with each particle representing a toroidal material unit formed by rotation around the axis of symmetry. During the calculation, the motion and deformation of the actual three-dimensional axisymmetric medium are reflected by updating the radial and axial coordinates of the particles, without the need to establish a complete three-dimensional computational model.

[0086] The particle neighborhood search module is used to establish inter-particle interaction relationships based on the current position of the particles, divide the computational region into multiple search units, and determine whether a valid particle pair is formed based on the inter-particle distance and the kernel function support domain. Adjacency relationships are re-established at each time step or a set update time. In this embodiment, the module divides the computational region into multiple search units and stores each particle in its corresponding search unit according to its spatial position. For any particle to be calculated, only particles within its own unit and adjacent units are searched, and a valid particle pair is determined based on the inter-particle distance and the kernel function support domain. This avoids global pairwise traversal of all particles, reducing the computational load of neighborhood search. At each time step or a set update time, the neighborhood search module re-establishes adjacency relationships based on the updated particle position, thus adapting to situations where the relative positions of particles change due to large fluid movements, large structural deformations, interface movement, and material breakage.

[0087] The material state calculation module is used to calculate pressure, sound velocity, and temperature state variables based on particle density and specific internal energy by invoking the corresponding state relationship according to the particle's material category. For solid particles, hydrostatic pressure and deviatoric stress are calculated separately. In this embodiment, for fluid particles, the material state calculation module obtains the thermodynamic pressure based on the fluid state equation. For metallic solid particles, the stress tensor is decomposed into hydrostatic pressure and deviatoric stress components, where hydrostatic pressure is calculated by the solid state equation, and deviatoric stress is updated by the solid strength and damage calculation module. For explosion products or gaseous particles, their pressure and sound velocity can be determined according to the corresponding gas state equation or detonation product state equation. Each material state equation is set up as an independent interface and called through material type identifiers, enabling the system to extend material models for different fluids, metals, gases, and explosion media.

[0088] The spatial state reconstruction module is used to construct left and right states on both sides of the particle interface between particle pairs of the same material. Based on the local pressure, density, velocity, specific internal energy, and spatial gradient of the particles, this module obtains the reconstructed state at the particle pair interface through a constrained linear reconstruction method. To avoid new non-physical solutions arising from higher-order reconstruction in extreme regions, the spatial state reconstruction module restricts the reconstruction slope and performs a physical validity assessment on the reconstructed density, pressure, specific internal energy, and sound velocity. When the reconstructed state exceeds a preset physical range, the corresponding particle pair is degenerated into a lower-order reconstruction or the original local state of the particles is adopted. For fluid-solid cross-material particle pairs, instead of directly performing high-order cross-medium reconstruction, the local states of the particles on both sides of the interface are used to construct a Riemann initial value problem to reduce non-physical oscillations caused by abrupt changes in material parameters, equations of state, and acoustic impedance.

[0089] The fluid-solid interface processing module is used to identify cross-medium particle pairs formed by fluid and solid particles, and to convert the projection of the solid deviatoric stress onto the normal direction of the particle pair into an equivalent normal pressure. This allows the fluid pressure and the total solid stress to be processed in the same local Riemann problem. For particle pairs of the same material, the reconstructed state at the interface is obtained through a restricted linear reconstruction method. For fluid-solid cross-material particle pairs, the local states of the particles on both sides of the interface are used to construct the initial Riemann problem. In this embodiment, the normal velocities on the fluid side and the solid side are obtained by projecting the corresponding particle velocities onto the unit direction vector. The fluid-solid interface processing module transmits the equivalent pressure, normal velocity, density, and acoustic impedance to the Riemann flux solving module and receives the solved interface pressure and interface velocity.

[0090] The Riemann flux solving module is used to establish a local one-dimensional Riemann initial value problem between each effective particle pair. It calculates the interface star region state based on the pressure, velocity, density, and sound velocity on both sides of the particle pair. For fluid-solid cross-material particle pairs, it uses the fluid-side pressure state and the equivalent normal pressure state on the solid side. For particle pairs within the solid, the hydrostatic pressure is obtained through Riemann calculation, while the deviatoric stress is used as a complete stress flux in the particle interaction calculation. The Riemann flux solving module weights the pressure and normal velocity based on the acoustic impedance on both sides to obtain the interface pressure and interface normal velocity, which act on both sides of the particle pair in pairs and opposite directions to update the particle's mass, momentum, and rate of energy change. The solid strength and damage calculation module is used to update the deviatoric stress, equivalent stress, yield state, equivalent plastic strain, temperature, and damage variables of solid particles. It calculates the strain rate based on the neighborhood velocity field of the solid particles and obtains the deviatoric stress trial value using the Jaumann objective stress update relationship. Plastic correction is performed according to the von Mises yield criterion and the Johnson-Cook constitutive model. The adiabatic temperature rise is calculated using Taylor-Quinney coefficients, and damage variables are accumulated based on the plastic strain increment. In this embodiment, the equivalent stress is calculated based on the deviatoric stress tensor and compared with the current material yield strength. When the equivalent stress does not exceed the material yield strength, the solid particles remain in an elastic state. When the equivalent stress exceeds the material yield strength, the deviatoric stress state is corrected to the yield surface using a radial return method, and the equivalent plastic strain and equivalent plastic strain rate are updated simultaneously. The elastic limit is determined by the von Mises yield condition, and the yield strength is given by the Johnson-Cook constitutive model. The von Mises yield criterion is used to determine the particle's yield strength when updating the stress. Is the current state a plastic yield state? If the Mises stress is less than the yield strength... So, particles It remains in an elastic state, and the deviatoric stress components do not require correction. If the Mises stress exceeds the yield strength... Then the particle Upon reaching the plastic yield state, the deviatoric stress components must proportionally recede back to the yield surface, and the deviatoric stress components are corrected to... ; The axisymmetric geometry correction module is used to convert the radial-axial two-dimensional particle calculation results into a three-dimensional equivalent response of an axisymmetric physical problem. This module adds cylindrical coordinate geometric terms related to the particle's radial position to the mass, momentum, and energy conservation relationships, and incorporates circumferential strain rate and circumferential deviatoric stress into the solid stress update, ensuring coordinated updates of radial stress, axial stress, shear stress, and circumferential stress within the same time step. To address the potential numerical amplification caused by dividing the geometric terms by the radial coordinate when the particle's radial coordinate approaches zero, the axisymmetric geometry correction module establishes an axisymmetric safety region. When the particle is outside the axisymmetric safety region, the circumferential strain rate is calculated using the ratio of radial velocity to radial coordinate; when the particle is inside the axisymmetric safety region, the circumferential strain rate is calculated using an axisymmetric approximation based on symmetry conditions, velocity gradients, or constraint values ​​to avoid numerical divergence near the axis.

[0091] The time integration module is used to update the state variables of each particle according to an explicit prediction-correction method. In the prediction phase, based on the particle position, velocity, density, specific internal energy, pressure, and stress change rate at the current moment, the predicted state at an intermediate moment is obtained. The predicted state is then used to re-perform neighborhood search, material state calculation, interface identification, Riemann flux solution, and solid constitutive update. In the correction phase, the particle position, velocity, density, specific internal energy, pressure, stress, temperature, and damage variables at the next moment are obtained based on the recalculated change rate. The time step is limited by the particle spacing, smooth length, sound velocity, particle velocity, particle acceleration, and the degree of material state change. The time step is reduced when high sound velocity, high velocity, strong pressure gradients, or severe structural deformation occur in the computational domain; a larger time step is allowed when particle state changes are relatively gradual.

[0092] The results output module outputs particle position, velocity, density, pressure, specific internal energy, temperature, material type, stress components, equivalent stress, equivalent plastic strain, yield state, and damage variables at preset time intervals. The module can also set virtual measuring points at specified spatial locations and record pressure, velocity, density, and structural response near these points in chronological order. It can extract process data such as shock wave propagation, bubble boundary motion, free surface changes, structural deformation, and material failure based on the spatial location of fluid particles, gas particles, solid particles, or material interface particles. The output data is transmitted to a visualization post-processing device via the results output interface to generate pressure contour maps, velocity contour maps, stress contour maps, damage contour maps, structural deformation maps, and time history curves of characteristic physical quantities.

[0093] The data connections between the model input and particle initialization module, particle neighborhood search module, material state calculation module, fluid-solid interface processing module, Riemann flux solution module, solid strength and damage calculation module, axisymmetric geometry correction module, time integration module, and result output module are as follows: The data output terminal of the model input and particle initialization module is connected to the data input terminals of the particle neighborhood search module, the material state calculation module and the time integration module, and is used to provide initial particle states and calculation control parameters to each calculation module; The output of the particle neighborhood search module is connected to the fluid-solid interface processing module and the solid strength and damage calculation module, respectively, and is used to provide the neighbor list, inter-particle distance, interaction direction and kernel function gradient of each particle. The output of the material state calculation module is connected to the fluid-solid interface processing module and the Riemann flux solving module, and is used to provide particle pressure, density, specific internal energy and sound velocity. The output of the fluid-solid interface processing module is connected to the Riemann flux solving module, and is used to provide the left and right reconstruction states of the same material particle pairs, as well as the fluid normal pressure state and solid equivalent normal pressure state across material particle pairs. The output of the Riemann flux solving module is connected to the time integration module to provide particle mass, momentum, and energy change rate; the outputs of the solid strength and damage calculation module and the axisymmetric geometry correction module are also connected to the time integration module to provide solid stress change rate, plastic state, damage state, and axisymmetric additional terms. The updated particle state from the time integration module is re-input into the particle neighborhood search module and the material state calculation module, thereby forming a data closed loop for time-step cyclic solution; when the output condition is met, the time integration module transmits the current particle state to the result output module.

[0094] In a specific implementation, as a preferred embodiment of the present invention, the system of the present invention further includes a parallel computing module, which is used to distribute the tasks of particle initialization, neighborhood search, particle pair interaction, material state update, solid constitutive update and result statistics to multiple processor threads; the independent particle loops are executed in parallel by threads, and the calculation of particle pairs that may be written to the same particle variable at the same time is protected by using thread-private variables, local accumulation, atomic operations or partitioned calculation methods to reduce data competition in parallel computing.

[0095] The hardware processor, memory, and data interface in this invention do not undergo mechanical position changes during operation; the dynamic changes of the system are mainly manifested in the changes of the position, velocity, and material state of various Lagrange particles in the calculation model over time.

[0096] Fluid particles undergo displacement driven by pressure gradients and interfacial forces. The outer boundary of the particle set represents a free surface, a shock wave region, a bubble boundary, or a fluid-solid interface. Solid particles undergo displacement under fluid pressure, impact loads, and internal stresses. The shape of a solid structure is determined by the coordinate set of the solid particles. When a solid is in the elastic stage, it can partially or completely recover after unloading. After entering the plastic stage, it undergoes irreversible deformation. When damage accumulates to the failure condition, local particles lose or reduce their load-bearing capacity and can continue to move due to inertia, thus forming numerical forms such as material cracking, perforation, breakage, or fragmentation.

[0097] The interface between fluid and solid does not require a fixed mesh connection; instead, it is determined by the real-time positions and material types of adjacent fluid and solid particles. When the interface moves, bends, separates, or re-engages, the system automatically adjusts the interacting objects by updating the particle neighborhood relationships, thereby maintaining continuous calculation of the fluid-solid load transfer process.

[0098] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A numerical simulation method for explosion and impact fluid-structure interaction based on Riemann-ULSPH, characterized in that, include: S1. Establish a computational model and generate initial particles: Determine the fluid region, solid structure, explosive medium, boundary region and axis of symmetry position according to the problem to be simulated, input geometric parameters and material parameters, and generate fluid particles, solid particles and boundary particles in the radial-axial plane according to the preset particle spacing. Assign each particle mass, position, density, velocity, pressure, specific internal energy, sound speed, material type, initial stress, temperature, plastic strain and damage variables. S2. Establish neighborhood relationships and calculate material state: Construct a linked list search unit based on the current particle position, determine the neighboring particles in the kernel function support domain, and call the corresponding state equations according to the material type of each particle to obtain pressure and sound velocity. For solid particles, calculate hydrostatic pressure and deviatoric stress separately. S3. Identify particle pair categories and construct interface states: Based on the material types on both sides of the particle pair, classify the particle pairs into fluid-fluid, solid-solid, and fluid-solid particle pairs. For particle pairs of the same material, reconstruct the left and right states of the interface using spatial gradient and constraint. For fluid-solid particle pairs, calculate the projection of the solid deviatoric stress onto the normal direction of the particle pair and construct the equivalent normal pressure of the solid. S4. Solve the local Riemann problem: Input the left and right pressures, normal velocities, densities, and sound speeds of the interface into the approximate Riemann solver to obtain the interface pressures and velocities. Calculate the density change rate, velocity change rate, and specific internal energy change rate of each particle from the interface state. S5. Update solid stress and damage: Calculate strain rate, rotation rate and circumferential strain rate based on the neighborhood velocity field of solid particles, update the deviatoric stress test state, determine whether the material has undergone plastic deformation according to the yield criterion, and perform plastic correction, temperature update and damage accumulation when yielding. Particles that reach the failure condition reduce or lose their deviatoric stress bearing capacity. S6. Apply axisymmetric correction and advance particle state: Calculate axisymmetric mass, momentum and energy additions based on particle radial position, use axisymmetric safety treatment for particles close to the axis of symmetry, and update particle position, velocity, density, specific internal energy, pressure, stress, temperature and damage variables using prediction-correction time integration method, while determining the next time step based on stability conditions. S7. Loop Calculation and Result Output: Repeat steps S2 to S6 until the set termination time is reached. During the calculation, particle field data and measurement point time history data are output at preset intervals.

2. The numerical simulation method for explosion and impact fluid-structure interaction based on Riemann-ULSPH according to claim 1, characterized in that, Step S1 includes: S11. Based on the high-speed collision, explosion impact, underwater explosion, or fluid impact structural problem to be simulated, determine the fluid region, solid structure, explosion medium, boundary region, and symmetry axis position. S12. Input the geometric and material parameters of each region. The geometric parameters include the computational region, structural dimensions, and symmetry axis position. The material parameters include material type, initial velocity, initial pressure, initial density, specific internal energy, boundary conditions, particle spacing, smooth length, computation termination time, and result output interval. S13. Generate fluid particles, solid particles and boundary particles in the radial-axial plane according to the preset particle spacing, and assign each particle mass, position, density, velocity, pressure, specific internal energy, sound speed and material type. Solid particles are also assigned initial stress, temperature, plastic strain and damage variables. S14. For axisymmetric problems, particles are arranged in the radial-axial computational plane. Each particle represents a ring-shaped material unit formed by rotating around the axis of symmetry. During the calculation, the motion and deformation of the actual three-dimensional axisymmetric medium are reflected by updating the radial and axial coordinates of the particles, without the need to establish a complete three-dimensional computational model.

3. The numerical simulation method for explosion and impact fluid-structure interaction based on Riemann-ULSPH according to claim 1, characterized in that, Step S2 includes: S21. Construct a linked list search unit based on the current particle position, divide the computation area into multiple search units, and store each particle in the corresponding search unit according to its spatial position. S22. For any particle to be calculated, only search for particles in its own cell and adjacent cells, and determine whether they constitute a valid particle pair based on the distance between particles and the support domain of the kernel function. S23. At each time step or at a set update time, the adjacency relationship is re-established based on the updated position of the particles to adapt to changes in the relative position of particles after large fluid movement, large structural deformation, interface movement, and material breakage. S24. Based on the material type of each particle, call the corresponding equation of state to obtain the pressure and speed of sound; S25. For fluid particles, thermodynamic pressure is obtained based on the fluid state equation. S26. For metallic solid particles, the stress tensor is decomposed into hydrostatic pressure and deviatoric stress components, where hydrostatic pressure is calculated by the solid state equation and deviatoric stress is updated by subsequent solid strength and damage calculations. S27. For explosion products or gaseous particles, determine their pressure and sound velocity according to the gas state equation or the detonation product state equation. S28. Each material's equation of state is set up with an independent interface and called through the material type identifier, enabling the system to extend material models for different fluids, metals, gases, and explosive media.

4. The numerical simulation method for explosion and impact fluid-structure interaction based on Riemann-ULSPH according to claim 1, characterized in that, Step S3 includes: S31. Based on the material types on both sides of the particle pair, the particle pairs are classified into fluid-fluid particle pairs, solid-solid particle pairs and fluid-solid particle pairs. Different types of particle pairs adopt different state construction methods. S32. For a pair of particles of the same material, the reconstructed state at the particle pair interface is obtained through a restricted linear reconstruction method based on the local pressure, density, velocity, specific internal energy and spatial gradient of the particles. S33. Limit the reconstruction slope and make a physical validity judgment on the reconstructed density, pressure, specific internal energy and sound speed. When the reconstruction state exceeds the preset physical range, the corresponding particle pair will be degenerated into a low-order reconstruction or the local original state of the particle will be adopted. S34. For fluid-solid cross-material particle pairs, instead of directly performing high-order reconstruction across the medium, the Riemann initial value problem is constructed using the local states of particles on both sides of the interface. S35. For a fluid-solid particle pair, determine the direction of interaction between the two particles based on their position difference, and use this direction as the normal direction for the local Riemann problem. The normal direction is: in, This represents the unit normal vector of the particle to the Riemann solution direction. Indicates the right-hand particle The position vector, Indicates the particle on the left The position vector; S36. The fluid-solid interface in the continuous sense satisfies the normal kinematic compatibility condition and the dynamic traction continuity condition, as follows: in, Indicates the velocity at the fluid interface. Indicates the solid-solid interface velocity. Represents the unit normal vector of the interface. This represents the Cauchy stress tensor of the fluid. Represents the Cauchy stress tensor of a solid. S37. For inviscid fluids, their Cauchy stress tensor degenerates into... ,in, Represents the thermodynamic pressure of a fluid. Represents the unit tensor; the fluid-side normal traction is... Solid side normal traction is The solid side normal component is ,in, Represents the hydrostatic pressure of a solid. Represents the deviatoric stress tensor of a solid; S38, Based on traction continuity conditions The equivalent normal pressure entering the Riemann solver from the solid side is defined as: S39. Because the direction of the particle's Riemann solution is located... In the calculation plane, i.e. ,exist The plane unfolds as follows: in, Represents the radial deviatoric stress component. Represents the shear deviatoric stress component. Indicates the axial deviatoric stress component. Indicates the normal direction In the unit component of the radial direction, Indicates the normal direction Unit component in the axial direction; S310. Based on the equivalent pressure definition, the left and right normal pressure states of the fluid-solid particle pair are transformed into a standard form one-dimensional Riemann initial value problem. For any fluid-solid particle pair, the normal pressure entering the Riemann solver is defined as: in, This indicates the normal pressure entering the Riemann solver from the left. This indicates the normal pressure entering the Riemann solver from the right side. This represents the hydrostatic pressure of the particle on the left. This represents the hydrostatic pressure of the particle on the right. Represents the fluid domain. Represents a solid domain. This represents the deviatoric stress tensor of the solid particles on the left. This represents the deviatoric stress tensor of the solid particles on the right. S311. Calculate the normal velocities on the left and right sides of the fluid-solid interface, using the following formula: in, Indicates the normal velocity on the left side. Indicates the normal velocity on the right side. Indicates the particle on the left The velocity vector, Indicates the right-hand particle The velocity vector.

5. The numerical simulation method for explosion and impact fluid-structure interaction based on Riemann-ULSPH according to claim 1, characterized in that, Step S4 includes: S41. Input the left and right interface pressures, normal velocities, densities, and sound velocities into the PVRS approximation Riemann solver based on acoustic impedance weighting to calculate the contact pressures and normal velocities, as shown in the following formulas: in, Indicates the acoustic impedance on the left side. , Indicates the particle density on the left. Indicates the sound speed of the particle on the left; Indicates the acoustic impedance on the right side. , This represents the particle density on the right. Indicates the sound speed of the particle on the right; S42. The acoustic impedance of the solid-side PVRS is constructed using the local volume sound velocity given by the Mie-Grüneisen equation of state, as shown in the following formula: in, Represents solid particles The speed of sound, Indicates the volume compression ratio. , , All are Hugoniot polynomial coefficients. Indicates the reference density. The reference sound velocity is indicated; the contribution of deviatoric stress is not incorporated into the sound velocity, but is explicitly calculated through the constitutive update of the solid and the deviatoric stress flux in the momentum equation. S43. Calculate the rate of change of density, rate of change of velocity, and rate of change of specific internal energy of each particle based on the interface state. The normal pressure between the fluid and solid intermaterial particle pairs is incorporated into the momentum equation in the following form: in, Represents particles The rate of change of velocity, Represents particles density, Indicates the interface contact pressure. Representing kernel function Regarding particles Gradient of position, Representing neighboring particles Volume; S44. The intact stress flux within a solid enters the momentum equation in the following form: in, This indicates interfacial deviatoric stress.

6. The numerical simulation method for explosion and impact fluid-structure interaction based on Riemann-ULSPH according to claim 1, characterized in that, Step S5 includes: S51. Calculate the radial strain rate, axial strain rate, shear strain rate, circumferential strain rate and rotation rate based on the neighborhood velocity field of solid particles, and obtain the experimental value of solid deviatoric stress using the Jaumann objective stress update relationship. S52. Calculate the equivalent stress based on the deviatoric stress tensor and compare it with the current material yield strength; when the equivalent stress does not exceed the material yield strength, the solid particles remain in an elastic state, and the deviatoric stress components do not need correction; when the equivalent stress exceeds the material yield strength, the particles... Upon entering the plastic yield state, the deviatoric stress components must proportionally recede back to the yield surface. The deviatoric stress components are corrected as follows: in, This represents the corrected deviatoric stress tensor; Represents the deviatoric stress test tensor; The yield strength of the material is represented by the Johnson-Cook yield model. , Indicates static yield strength. and All of these represent strain hardening parameters. Indicates the strain rate strengthening parameter. Indicates the thermal softening parameter. For equivalent plastic strain, Represents the equivalent plastic strain rate. Indicates the reference strain rate. Represents dimensionless temperature. , Melting point Room temperature; Indicates the equivalent Mises stress; This represents the second invariant of the deviatoric stress tensor. , Represents the matrix trace operation. This represents the components of the deviatoric stress tensor in cylindrical coordinates. S53. The adiabatic temperature rise is calculated using the Taylor-Quinney coefficient and plastic work, as shown in the following formula: in, This represents the Taylor-Quinney coefficient; This represents the equivalent plastic strain increment within a single time step, i.e., the amount of plastic deformation that occurs in the material during this step, based on the von Mises yield criterion and Plasticity theory , This represents the yield strength at the current moment. Indicates shear modulus; Indicates constant volume specific heat capacity; This represents the temperature rise within a single time step; S54. Update temperature using a prediction-correction time integral format: The starting point is indicated by a superscript. The midpoint is indicated by a superscript. This indicates that the new time step is indicated by a superscript. express; Indicates based on The trial temperature rise calculated based on the instantaneous state and half-step size; Indicates based on High-precision flow field and full-step flow at any time The calculated actual temperature rise; S55. The yield strength of the material is determined based on at least one of the state variables: plastic strain, strain rate, and temperature. Part of the plastic work generated by plastic deformation is converted into heat energy, which is used to update the temperature of solid particles and feeds the updated temperature back to the material thermal softening calculation. S56. Damage variables are accumulated based on the ratio of single-step plastic strain increment to the current failure strain. When the damage variable reaches the preset failure threshold, the corresponding particle is marked as a failed particle. In subsequent time steps, the deviatoric stress and material strength contribution of the particle are reduced or eliminated, but at least a part of its mass, inertia and pressure state are retained to simulate the motion process after material failure, fragmentation and loss of load-bearing capacity.

7. The numerical simulation method for explosion and impact fluid-structure interaction based on Riemann-ULSPH according to claim 1, characterized in that, Step S6 includes: S61. Add cylindrical coordinate geometric additions related to the radial position of the particle to the mass conservation, momentum conservation and energy conservation relationships, and add circumferential strain rate and circumferential deviatoric stress to the solid stress update, so that the radial stress, axial stress, shear stress and circumferential stress are updated in a coordinated manner within the same time step. S62. Set an axial safety zone. When the particle is outside the axial safety zone, the circumferential strain rate is calculated using the ratio of radial velocity to radial coordinate. When the particle is inside the axial safety zone, the circumferential strain rate is calculated using an axial approximation method based on symmetry conditions, velocity gradient, or limit values ​​to avoid numerical divergence near the axis. S63. In the prediction phase, based on the particle position, velocity, density, specific internal energy, pressure, and stress change rate at the current moment, the predicted state at the intermediate moment is obtained. S64. Re-perform neighborhood search, material state calculation, interface identification, Riemann flux solution, and solid constitutive update using the predicted state; S65. During the correction phase, the particle position, velocity, density, specific internal energy, pressure, stress, temperature, and damage variables for the next moment are obtained based on the recalculated rate of change. S66. The time step is limited based on the particle spacing, smooth length, sound speed, particle velocity, particle acceleration, and the degree of material state change. When high sound speed, high velocity, strong pressure gradient, or severe structural deformation occurs in the computational domain, the time step is reduced. When the particle state change is relatively gentle, a larger time step is allowed.

8. The numerical simulation method for explosion and impact fluid-structure interaction based on Riemann-ULSPH according to claim 1, characterized in that, Step S7 includes: S71. Output particle position, velocity, density, pressure, specific internal energy, temperature, material type, stress components, equivalent stress, equivalent plastic strain, yield state and damage variables at preset time intervals. S72. Set up virtual measuring points at specified spatial locations and record the pressure, velocity, density, and structural response near the measuring points in chronological order. S73. Extract data on shock wave propagation, bubble boundary motion, free surface changes, structural deformation, and material failure processes based on the spatial positions of fluid particles, gas particles, solid particles, or material interface particles.

9. A numerical simulation system for explosion and impact fluid-structure interaction based on the Riemann-ULSPH method of any one of claims 1-8, deployed in a computing device, the computing device comprising a processor, a memory, a data input interface, and a result output interface, characterized in that, include: The system comprises the following modules: model input and particle initialization, particle neighborhood search, material state calculation, spatial state reconstruction, fluid-structure interface processing, Riemann flux calculation, solid strength and damage calculation, axisymmetric geometric correction, time integration, and result output. The model input and particle initialization module is used to receive parameters such as computational region, structural size, material type, initial velocity, initial pressure, initial density, specific internal energy, boundary conditions, particle spacing, smooth length, computation termination time, and result output interval. It discretizes the computational region into Lagrange particles carrying mass, position, velocity, density, pressure, specific internal energy, sound speed, stress tensor, temperature, equivalent plastic strain, and damage variables. Each particle is assigned a material type identifier to distinguish between fluid particles, solid particles, detonation product particles, boundary particles, and failure particles. For axisymmetric problems, the particles are arranged in the radial-axial computational plane, and each particle represents a ring-shaped material unit formed by rotating around the axis of symmetry. The particle neighborhood search module is used to establish the interaction relationship between particles based on the current position of the particles, divide the computation region into multiple search units, and determine whether a valid particle pair is formed based on the distance between particles and the support domain of the kernel function. The adjacency relationship is re-established at each time step or at a set update time. The material state calculation module is used to call the corresponding state relationship based on the material type of the particle to calculate the pressure, sound velocity and temperature state variables from the particle density and specific internal energy. For solid particles, hydrostatic pressure and deviatoric stress are calculated separately. The spatial state reconstruction module is used to construct the left and right states on both sides of the particle interface between particle pairs of the same material. The fluid-solid interface processing module is used to identify cross-medium particle pairs formed by fluid particles and solid particles, and to convert the projection of solid deviatoric stress onto the normal direction of the particle pair into an equivalent normal pressure. This allows the fluid pressure and the total solid stress to be processed in the same local Riemann problem. For particle pairs of the same material, the reconstructed state at the interface is obtained through a restricted linear reconstruction method. For fluid-solid cross-material particle pairs, the local state of the particles on both sides of the interface is used to construct the Riemann initial value problem. The Riemann flux solving module is used to establish a local one-dimensional Riemann initial value problem between each effective particle pair. It calculates the interface star region state based on the pressure, velocity, density and sound speed on both sides of the particle pair. It weights the pressure and normal velocity according to the acoustic impedance on both sides to obtain the interface pressure and interface normal velocity, which act on both sides of the particle pair in pairs and opposite directions to update the mass, momentum and energy change rate of the particles. The solid strength and damage calculation module is used to update the deviatoric stress, equivalent stress, yield state, equivalent plastic strain, temperature and damage variables of solid particles. It calculates the strain rate based on the neighborhood velocity field of solid particles and obtains the deviatoric stress trial value using the Jaumann objective stress update relationship. It performs plastic correction based on the von-Mises yield criterion and Johnson-Cook constitutive model, calculates the adiabatic temperature rise using Taylor-Quinney coefficients, and accumulates the damage variables based on the plastic strain increment. The axisymmetric geometry correction module is used to add cylindrical coordinate geometry additions related to the radial position of the particle to the mass conservation, momentum conservation and energy conservation relationships, add circumferential strain rate and circumferential deviatoric stress to the solid stress update, and adopt axisymmetric safety treatment for particles close to the axis of symmetry. The time integration module is used to update the state variables of each particle according to the explicit prediction-correction method. In the prediction stage, it obtains the predicted state at the intermediate time and recalculates each rate of change. In the correction stage, it obtains the particle state at the next time and determines the next time step based on the stability condition. The results output module is used to output particle state data at preset time intervals, set virtual measuring points at specified spatial locations and record time history data, and extract data on shock wave propagation, structural deformation and material failure processes. The data connections between the model input and particle initialization module, particle neighborhood search module, material state calculation module, spatial state reconstruction module, fluid-solid interface processing module, Riemann flux solution module, solid strength and damage calculation module, axisymmetric geometry correction module, time integration module, and result output module are as follows: The data output terminal of the model input and particle initialization module is connected to the data input terminals of the particle neighborhood search module, the material state calculation module and the time integration module, and is used to provide initial particle states and calculation control parameters to each calculation module; The output of the particle neighborhood search module is connected to the fluid-solid interface processing module and the solid strength and damage calculation module, respectively, and is used to provide the neighbor list, inter-particle distance, interaction direction and kernel function gradient of each particle. The output of the material state calculation module is connected to the fluid-solid interface processing module and the Riemann flux solving module, and is used to provide particle pressure, density, specific internal energy and sound velocity. The output of the spatial state reconstruction module is connected to the Riemann flux solving module to provide the left and right reconstruction states of the same material particle pairs. The output of the fluid-solid interface processing module is connected to the Riemann flux solving module, and is used to provide the left and right reconstruction states of the same material particle pairs, as well as the fluid normal pressure state and solid equivalent normal pressure state across material particle pairs. The output of the Riemann flux solving module is connected to the time integration module to provide particle mass, momentum, and energy change rate; the outputs of the solid strength and damage calculation module and the axisymmetric geometry correction module are also connected to the time integration module to provide solid stress change rate, plastic state, damage state, and axisymmetric additional terms. The updated particle state from the time integration module is re-input into the particle neighborhood search module and the material state calculation module, thereby forming a data closed loop for time-step cyclic solution; when the output condition is met, the time integration module transmits the current particle state to the result output module.

10. The Riemann-ULSPH-based numerical simulation system for explosion and shock fluid-structure interaction according to claim 9, characterized in that, It also includes a parallel computing module, which is used to distribute tasks such as particle initialization, neighborhood search, particle pair interaction, material state update, solid constitutive update and result statistics to multiple processor threads; independent particle loops are executed in parallel using threads, and particle pair calculations that may be written to the same particle variable at the same time are protected by using thread-private variables, local accumulation, atomic operations or partitioned calculation methods.

Images