A high-frequency thermal vibration coupled SPH simulation method

Abstract

The application discloses a high-frequency thermal vibration coupling SPH simulation method, comprising the following steps: S1, adopting a total Lagrangian SPH method to construct a solid calculation framework; S2, introducing a Riemann solver in the total Lagrangian SPH method to determine the interaction between particles, instead of an artificial viscosity term in a traditional SPH, and eliminating the influence of the artificial viscosity on numerical stability; S3, constructing a heat conduction control equation based on an internal energy form, utilizing a kernel function to disperse a thermal energy term, calculating the evolution process of enthalpy or internal energy on each particle with time, and calculating the particle temperature according to unit mass enthalpy; S4, realizing thermal force coupling through a state equation or a thermal expansion strain influence stress calculation; and S5, utilizing an explicit time integration form to update the state quantity of the particle, and cyclically iterating to complete high-frequency response calculation. The application can eliminate the influence of the artificial viscosity on numerical stability, and is suitable for high-frequency vibration response calculation under extreme conditions.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of aircraft structure dynamics simulation, in particular to a SPH simulation method for high-frequency thermal-vibration coupling. BACKGROUND

[0002] In the prior art, under hypersonic flight conditions, key structures such as thermal protection panels will encounter external excitations such as noise, impact, pulsating pressure, etc., and produce nonlinear responses such as fatigue, high-frequency vibration, etc., and the thermal-vibration coupling effect caused by aerodynamic heating will further aggravate the accumulation of fatigue damage.

[0003] Currently, mainstream commercial finite element software all have relatively complete thermal-vibration coupling functions, and are suitable for thermal expansion analysis, thermal fatigue, thermal shock, etc. However, the traditional mesh method is limited by time steps, precision, and nonlinear adaptability, etc., and there is currently no complete solution for predicting high-frequency responses under thermal-vibration coupling, such as the vibration response of a hypersonic aircraft structure under 1-20 kHz wideband loading.

[0004] The thermal-vibration coupling simulation of most commercial software has good precision and stability in the low-frequency band (1 Hz-2 kHz). In the high-frequency response analysis of structures (>5 kHz), such methods gradually reveal the following limitations:

[0005] (1) Limited frequency resolution: Commercial software can theoretically support up to several tens of kHz, but in the thermal-vibration coupling process, due to factors such as mesh size, time step precision, and modeling, conventional finite element simulation is difficult to accurately capture the structural response characteristics at higher frequencies.

[0006] (2) Mesh distortion and numerical stability problems: In high-frequency loading or large deformation scenarios, mesh-based finite element methods are prone to severe numerical dissipation and instability, leading to a significant increase in calculation errors as the frequency increases.

[0007] (3) Lack of nonlinear dynamic adaptability: High-frequency vibration is often accompanied by stress wave propagation, contact impact, and other strong nonlinear behaviors, making it difficult for standard FEM to effectively simulate.

[0008] (4) Simple coupling algorithm: Thermal-vibration coupling in commercial software is generally weakly coupled or unidirectionally coupled, with large errors in high-frequency environments.

[0009] Overall, high-frequency thermal-vibration coupling is currently limited to low-frequency coupling, and high-frequency and ultrahigh-frequency coupling are in a state of deficiency. SUMMARY

[0010] The present application aims to provide a SPH simulation method for high-frequency thermal-vibration coupling to solve the problems of insufficient prediction accuracy in the high-frequency band, introduction of additional dissipation by artificial viscosity, and limited application range of thermal-force coupling in the prior art.

[0011] To achieve the above object, the application provides a high-frequency thermal vibration coupling SPH simulation method, comprising the following steps:

[0012] S1, constructing a solid calculation framework by using a total Lagrangian SPH method;

[0013] S2, introducing a Riemann solver in the total Lagrangian SPH method to determine the interaction between particles, replacing the artificial viscosity term in the traditional SPH, and eliminating the influence of artificial viscosity on numerical stability;

[0014] S3, constructing a heat conduction control equation based on internal energy form, discretizing the heat conduction equation on the particle scale, and combining the unit mass enthalpy to update the particle temperature, so as to ensure the calculation stability under the material system with discontinuous thermal conductivity;

[0015] S4, realizing thermal coupling by affecting stress calculation through state equation or thermal expansion strain;

[0016] S5, in each time step, using the form of explicit time integration to update the displacement, stress, internal energy, temperature and other state quantities of the particle, and iteratively completing the high-frequency response calculation.

[0017] Preferably, in S1, the momentum conservation equation of the solid calculation framework constructed by using the total Lagrangian SPH method is calculated as follows:

[0018] ;

[0019] The energy conservation equation of the solid calculation framework constructed by using the total Lagrangian SPH method is calculated as follows:

[0020] ;

[0021] In the formula, the subscript , are particle numbers, representing any particle i and the particles in the neighborhood of the particle ; P is the first PK stress, B 0 is a correction matrix, represents the gradient of the kernel function , wherein X is the particle coordinate, and h is the smoothing length. In the symmetric SPH method, ; , , , are the initial density, velocity, internal energy and mass of the particle, respectively; is the volume of the reference configuration particle, is the "Kronecker integral", a special form of the tensor product;

[0022] .

[0023] Preferably, in S2, the calculation formula of the Riemann solver is:

[0024] ;

[0025] ;

[0026] wherein, is the velocity approximation of the Riemann solution, is the stress approximation of the Riemann solution, is the stress on the particle.

[0027] Preferably, in S2, the calculation formula after introducing the Riemann solver in the total Lagrangian SPH method is:

[0028] ;

[0029] ;

[0030] ;

[0031] ;

[0032] wherein, is the Jacobian matrix related to the deformation gradient ;

[0033] The calculation formula of the stress approximation of the Riemann solution is: , , the subscript new represents the new value of the parameter after introducing the Riemann solver;

[0034] The particle approximation of the deformation gradient and the deformation gradient rate is:

[0035] ;

[0036] 。

[0037] Preferably, in S3, the heat conduction control equation is:

[0038] ;

[0039] ​ ;

[0040] wherein, is the thermal conductivity of the combined material, which ensures smooth heat flux between SPH calculation particles with discontinuous thermal conductivity, k is the thermal conductivity of the material; , , are the specific enthalpy, volume, kernel function of the corresponding particle respectively w derivative of the particle , distance ; , is the current temperature of the particle, is the particle spacing;

[0041] The change of is obtained by the heat conduction control equation using the explicit time integration form: :

[0042] ;

[0043] wherein, represents the heat conduction time.

[0044] Preferably, the change of the particle temperature is calculated by and the specific heat capacity , and the calculation formula is:

[0045] ;

[0046] wherein, is the change of the particle temperature.

[0047] Preferably, in S4, the influence of the temperature on the mechanical response is calculated by the thermal expansion strain, and the calculation formula is:

[0048] ;

[0049] wherein, is the thermal expansion strain, is the thermal expansion coefficient, is the initial temperature;

[0050] The thermal expansion strain is combined with the particle displacement strain to obtain the total strain , and the strain is converted into the particle coordinates by using the Lagrangian formula: , is the initial coordinate of the particle, is the characteristic length; then, the displacement gradient and velocity gradient : , , is derivative with respect to time, is inverse of ; velocity gradient is decomposed into strain rate and spin rate : , , is transpose of, then strain rate is obtained, where denotes trace of tensor ; second order unit tensor is ; partial stress is obtained from Jaumann rate form of constitutive relation : where, is shear modulus of material; partial stress is initialized to zero or assigned according to specific working condition, and subsequent partial stress is updated in the form of explicit time integration: where, is updated (n+1)th partial stress, is updated nth partial stress; stress tensor is obtained through Cauchy stress decomposition: where is given by EOS state equation.

[0051] In S4, the effect of temperature on mechanical response is calculated through state equation, and Mie-Grüneisen EOS formula is used for thermal coupling, and the stress expression is modified as:

[0052] ;

[0053] ;

[0054] ;

[0055] where, is Hugoniot linear stress, is Grüneisen coefficient, , is nonlinear correction term, denotes volume strain of material, is current density of particle, in TLSPH, mass conservation equation is written as i.e. where Jacobian matrix represents the volume ratio after material deformation and before deformation, wherein, represents determinant of represents the solid sound speed, , , are empirical coefficients;

[0056] When using the Mie-Grüneisen EOS formula, it is not necessary to separately calculate the influence of structural thermal expansion, and the enthalpy value , i.e. the change of the internal energy of the particle affects the stress state of the particle; when , i.e. the structure is in a compressed state; when , i.e. the structure is in an expanded state.

[0057] Therefore, the present application adopts the above-mentioned high-frequency thermal vibration coupling SPH simulation method, compared with the prior art, the present application has the following beneficial effects:

[0058] (1) Numerical stability is improved: the Riemann solver is introduced into the total Lagrangian SPH framework, which completely eliminates the dependence on artificial viscosity, and significantly reduces the numerical dissipation error in high-frequency vibration simulation.

[0059] (2) Thermal force coupling is more universal: the present application proposes a double-path coupling method, which can not only handle thermal expansion under low strain, but also cope with strong nonlinear scenarios such as impact and explosion, and the application scope is obviously superior to the traditional single coupling method.

[0060] (3) Strong high-frequency prediction ability: the present application method can still keep the calculation error within ±5% under 1-20kHz high-frequency load, solving the problem of accuracy decline of finite element method in high-frequency band.

[0061] (4) Good engineering applicability: the SPH method of the present application supports MPI parallel architecture, which can realize large-scale multi-node parallel in runtime, and is well adapted to supercomputers such as Tianhe and Shenwei. This function can greatly save calculation time, increase server utilization, and improve running efficiency. Unlike existing SPH methods that rely on artificial viscosity and single coupling method, the present application introduces the Riemann solver under the TLSPH framework and combines the double-path thermal-force coupling mechanism, thereby realizing higher accuracy and stronger engineering applicability in 1-20kHz high-frequency thermal vibration scenarios.

[0062] The technical solutions of the present application will be further described in detail below through the drawings and examples. DETAILED DESCRIPTION

[0063] Figure 1is a schematic diagram of steps of an embodiment of the high-frequency thermal-vibration coupled SPH simulation method of the present application;

[0064] Figure 2 is a schematic diagram of the measurement point temperature and error in the embodiment of the present application; wherein (a) is a measurement point temperature diagram, and (b) is a measurement point temperature error analysis diagram;

[0065] Figure 3 is a temperature cloud diagram obtained by SPH calculation in the embodiment of the present application; wherein (a) is at T = 0.004 s, and (b) is at T = 1 s;

[0066] Figure 4 is a free vibration displacement diagram of the free end of a cantilever beam under different particle densities in the embodiment of the present application;

[0067] Figure 5 is a high-frequency vibration-power spectral density analysis diagram of a certain hypersonic vehicle structure in the embodiment of the present application;

[0068] Figure 6 is a thermal-vibration coupled simulation result diagram of a certain nozzle in the embodiment of the present application; wherein (a) is a temperature cloud diagram, (b) is a stress vector diagram, and (c) is a mises stress cloud diagram. DETAILED DESCRIPTION

[0069] The technical solutions of the present application are further described below by means of the accompanying drawings and embodiments.

[0070] Unless otherwise defined, the technical terms or scientific terms used in the present application shall be understood as the usual meanings understood by those skilled in the art to which the present application belongs. The terms "first", "second", and similar words used in the present application do not represent any order, number, or importance, but are only used to distinguish different components. The terms "include" or "contain" and similar words mean that the elements or objects before the words cover the elements or objects listed after the words and their equivalents, and do not exclude other elements or objects. The terms "connect" or "connected" and similar words are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. The terms "up", "down", "left", "right", and the like are only used to represent relative positional relationships, and when the absolute positions of the described objects change, the relative positional relationships may also change accordingly.

[0071] Embodiment One

[0072] As shown in Figure 1 The present application provides a high-frequency thermal-vibration coupled SPH simulation method, which is applicable to a high-frequency vibration field with a frequency range of 1 kHz to 20 kHz, supports joint prediction of structural response and temperature evolution of a thermal protection panel of an aircraft under high-frequency load, and specifically includes the following steps:

[0073] S1, the total Lagrangian SPH method is used to construct a solid calculation framework, which can avoid the tensile instability in the calculation, ensure the numerical accuracy in the large deformation and nonlinear vibration scenarios, and obtain the mechanical behavior of the whole system by solving the conservation equation of the particle group and tracking the motion trajectory of each particle.

[0074] The momentum conservation equation of the total Lagrangian SPH method for constructing a solid calculation framework is calculated as follows:

[0075] ;

[0076] The energy conservation equation of the total Lagrangian SPH method for constructing a solid calculation framework is calculated as follows:

[0077] ;

[0078] In the formula, the subscripts , are particle numbers, representing any particle and the particles in its neighborhood ; P is the first PK stress, B 0 is the correction matrix, represents the gradient of the kernel function , where X is the particle coordinate, and h is the smoothing length. In the symmetric SPH method, ; , , , are the initial density, velocity, internal energy and mass of the particle, respectively; is the volume of the reference configuration particle, is the "Kronecker integral", which is a special form of tensor product.

[0079] ;

[0080] S2, the Riemann solver is introduced into the total Lagrangian SPH method to determine the interaction between particles, instead of the artificial viscosity term in the traditional SPH, and the numerical flux between particles is directly calculated, so as to eliminate the numerical dissipation problem caused by the artificial viscosity and improve the stability and accuracy of high-frequency (>5kHz) response.

[0081] The calculation formula of the Riemann solver is as follows:

[0082] ;

[0083] ;

[0084] where, is the velocity approximation Riemann solver, is the stress approximation Riemann solver, is the stress on the particle.

[0085] The computational formula after introducing the Riemann solver in the total Lagrangian SPH method is:

[0086] ;

[0087] ;

[0088] ;

[0089] ;

[0090] where, is the Jacobian matrix related to the deformation gradient .

[0091] The computational formula of the stress approximation Riemann solver is: , , the subscript new indicates the new value of the parameter after introducing the Riemann solver.

[0092] The particle approximation of the deformation gradient and the deformation gradient rate is:

[0093] ;

[0094] ;

[0095] S3, construct the heat conduction control equation based on the internal energy form, discretize the heat conduction equation on the particle scale, and combine the unit mass enthalpy to update the particle temperature, to ensure the calculation stability under the material system with discontinuous thermal conductivity.

[0096] The heat conduction control equation is:

[0097] ;

[0098] ;

[0099] where, is the thermal conductivity coefficient of the composite material, which ensures the smooth heat flux between the SPH calculation particles with discontinuous thermal conductivity, k is the material thermal conductivity coefficient; , , respectively, the specific enthalpy, volume, kernel function of the corresponding particle w for the particle , the derivative of the distance ; , is the current temperature of the particle, is the inter-particle distance.

[0100] The change of is obtained by the heat conduction control equation using the explicit time integration form:

[0101] ;

[0102] where, denotes the heat conduction time;

[0103] The change of the particle temperature is calculated by and the specific heat capacity , and the calculation formula is:

[0104] ;

[0105] where, is the change of the particle temperature.

[0106] S4, realize the thermal-mechanical coupling by the stress calculation of the state equation or thermal expansion strain, the coupling based on thermal expansion strain is suitable for low strain amplitude scenarios, the coupling based on the Mie-Grüneisen state equation is suitable for strong nonlinear working conditions such as explosion, impact and metal large deformation, thereby expanding the applicable range of the method.

[0107] The influence of temperature on the mechanical response is calculated by the thermal expansion strain, and the calculation formula is:

[0108] ;

[0109] In the formula, is the thermal expansion strain, is the thermal expansion coefficient, is the initial temperature.

[0110] The total strain is obtained by combining the thermal expansion strain and the particle displacement strain , and the strain is converted into particle coordinates using the Lagrangian formula: , is the initial coordinate of the particle, is the characteristic length; then, the displacement gradient and the velocity gradient​ : , , is derivative with respect to time, is inverse matrix; velocity gradient decomposed into strain rate and spin rate : , , is transpose of , then partial strain rate is obtained, where denotes trace of tensor is second order identity tensor; partial stress is obtained from Jaumann rate type constitutive : where, is shear modulus of material; partial stress is initialized to zero or assigned according to specific working condition, and the value of partial stress is updated by explicit time integration: where, is updated (n+1)th partial stress, is updated nth partial stress; stress tensor is obtained by Cauchy stress decomposition: where is given by EOS state equation.

[0111] Stress at current location is calculated by linear elastic EOS formula , where, is initial stress of particle, generally set to 0, is initial density of particle, is current density of particle, reflects volume strain of material, K is bulk modulus. This method is simple and fast, and is suitable for working conditions with low strain amplitude.

[0112] For explosion, impact, and metal large deformation problems, the influence of temperature on mechanical response is calculated by state equation, and Mie-Grüneisen EOS formula is used for thermal coupling, and the stress expression is corrected as:

[0113] ;

[0114] ;

[0115] ;

[0116] where, is the Hugoniot linear stress, is the Grüneisen coefficient, , is the nonlinear correction term, denotes the volumetric strain of the material, in TLSPH, the mass conservation equation is written as i.e. where the Jacobian matrix represents the ratio of the volume after deformation to the volume before deformation, where, represents the determinant of ; denotes the solid sound speed, , , are all empirical coefficients.

[0117] When using the Mie-Grüneisen EOS formula, there is no need to separately calculate the effect of structural thermal expansion, the enthalpy value i.e. the change of the internal energy of the particles affects the stress state of the particles; when i.e. the structure is in a compressed state; when i.e. the structure is in an expanded state.

[0118] S5, in each time step, the state quantities of the particles such as displacement, stress, internal energy, temperature, etc. are updated using the explicit time integration form, and the loop iteration is completed to complete the high-frequency response calculation.

[0119] The particle state quantities (displacement, stress, temperature, etc.) are updated in the form of explicit time integration, compared with the finite element method which is limited by the grid model, the particle model constructed by the SPH method supports smaller time steps, and is more suitable for high-frequency response calculation, thereby improving the stability and accuracy of high-frequency vibration calculation.

[0120] In the embodiment, the above method is specifically applied in the process and mainly includes three parts: (1) calculating the particle temperature through the discretized heat conduction equation; (2) solving the influence of temperature on physical parameters such as strain and pressure; (3) realizing the thermal vibration coupling calculation of the hypersonic aircraft through the Rimann-TLSPH algorithm, which is specifically as follows:

[0121] (1) In order to verify the accuracy and effectiveness of the heat conduction function in the SPH simulation method, the heat transfer process on a flat plate with a given boundary heat source temperature is simulated. A flat plate with a side length (x, y direction), thickness (z direction), initial density , Young's modulus , Poisson's ratio flat plate structure, thermal conductivity coefficient , specific heat capacity . Given a thermal boundary condition: for the region of , given temperature ; other particles initial temperature is . Respectively use finite element software and the SPH simulation method of this embodiment to simulate, and carry out error analysis, as shown in Figure 2 and Figure 3 .

[0122] The temperature (energy) transfer direction on the three-dimensional flat plate is from high to low, the process is stable and reliable, and the distribution is uniform without singular points. The error analysis at the measuring point shows that the error after stable calculation is less than 0.5%, which is within the acceptable range, and it can be considered that the heat conduction function is realized in the SPH simulation method of this embodiment.

[0123] (2) Verify the vibration simulation function of the SPH simulation method. A cantilever beam structure with a length (x direction), thickness (y direction), initial density , Young's modulus . Given an initial load , wherein is a constant, is the sound speed in the solid.

[0124] ;

[0125] , calculate the vibration period under different Poisson's ratios and compare with the analytical solution, as well as the difference of the results under different particle densities. As shown in Table 1 and Figure 4 :

[0126] Table 1 Free vibration period of cantilever beam

[0127] ;

[0128] Show the vibration period when the particle density , the calculation error increases with the increase of Poisson's ratio , but all are below 5%. And also carried out grid independence verification for the cantilever beam when , respectively calculated the vibration period of the cantilever beam when , with the increase of particle density, the difference between different solutions gradually decreases.

[0129] (3) High-frequency thermal vibration coupling simulation: the high-frequency vibration analysis of a certain hypersonic vehicle structure under external load is given, as shown in Figure 5and the thermal-vibration coupling calculation of a certain nozzle, as shown in Figure 6

[0130] Therefore, the present application adopts the above-mentioned SPH simulation method of high-frequency thermal-vibration coupling, can avoid the tensile instability in the calculation, eliminate the influence of artificial viscosity on the numerical stability, overcome the difficulties in the current high-frequency thermal-vibration coupling simulation, and provide the basis and guidance for the prediction and processing of high-frequency thermal-vibration coupling in engineering practice.

[0131] Finally, it should be noted that: the above examples are only used to illustrate the technical solutions of the present application, but not to limit them, although the present application has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that: it can still modify or equivalently replace the technical solutions of the present application, and these modifications or equivalent replacements cannot make the modified technical solutions deviate from the spirit and scope of the technical solutions of the present application.​

Claims

1. A method of SPH simulation of high frequency thermal vibration coupling, characterized in that: The method comprises the following steps: S1, constructing a solid calculation framework by using a total Lagrangian SPH method; In S1, a momentum conservation equation of the solid calculation framework constructed by using the total Lagrangian SPH method is as follows: ； An energy conservation equation of the solid calculation framework constructed by using the total Lagrangian SPH method is as follows: ； where the indices , are particle indices, representing an arbitrary particle and its neighbors ; P is the first PK stress, B 0 is the correction matrix, denotes the gradient of the kernel function , where X is the particle coordinate, h is the smoothing length, and in the symmetrical SPH method, ; , , , are the initial density, velocity, internal energy, and mass of the particle, respectively; is the volume of the reference configuration particle, is the "Kronecker integral", a special form of the tensor product; ； S2, introducing a Riemann solver in the total Lagrangian SPH method to determine the interaction between particles, replacing an artificial viscosity term in SPH, and eliminating the influence of the artificial viscosity on numerical stability; In S2, the calculation formula of the Riemann solver is as follows: ； ； wherein is the velocity approximation Riemann solver, is the stress approximation Riemann solver, is the stress on the particles; The calculation formula after introducing the Riemann solver in the total Lagrangian SPH method is as follows: ； ； ； ； wherein is the Jacobian matrix associated with the deformation gradient F from the calculation of the stress-approximated Riemann solver , , the index denotes the new value of the parameter after the introduction of the Riemann solver​ deformation gradient and the rate of deformation gradient is approximated by the particle ； ； S3, constructing a heat conduction control equation based on an internal energy form, discretizing the heat conduction equation at a particle scale, and combining unit mass enthalpy to update a particle temperature, so as to ensure calculation stability in a material system with discontinuous thermal conductivity; S4, realizing thermal coupling by calculating stress through a state equation or thermal expansion strain; S5, in each time step, a state quantity of a particle is updated by using an explicit time integration form, and a high-frequency response calculation is completed through cyclic iteration.

2. The method of claim 1, wherein: In S3, the heat conduction control equation is as follows: ； ； where is the thermal conductivity of the combined material, which ensures smooth heat flux between SPH calculation particles with discontinuous thermal conductivity, is the thermal conductivity of the material; , , is the unit mass enthalpy, volume, kernel function of the corresponding particle is the derivative of the particle distance ; , is the current temperature of the particle, is the particle spacing; By the heat conduction control equation, using the explicit time integration form, the change amount of is obtained : ； wherein, represents the heat conduction time.

3. The method of claim 2, wherein: The change in particle temperature is calculated by the specific heat capacity The calculation is obtained, the calculation formula is: ； wherein, is the change in particle temperature.

4. The method of claim 2, wherein: In S4, the influence of temperature on a mechanical response is calculated through thermal expansion strain, and the calculation formula is as follows: ； wherein is the thermal expansion strain, is the thermal expansion coefficient, is the initial temperature; thermal expansion strain With particle displacement strain Combined, the total strain is obtained The strain is transformed into particle coordinates using the Lagrange formula: , The initial coordinates of the particle. The characteristic length is then determined; next, the displacement gradient is calculated. and velocity gradient : , , yes The derivative with respect to time, yes The inverse matrix; the velocity gradient is decomposed into strain rate. and rotation rate : , , for Transpose the sample and then calculate the deviatoric strain rate. ,in Tensor traces, It is a second-order unit tensor; Jaumann's rate-type constitutive Jaumann rate : where, is the shear modulus of the material; the deviatoric stress is initialized to zero or assigned a value according to the specific case, and the subsequent deviatoric stress is updated in the form of a displayed time integral: where, is the updated (n+1)th deviatoric stress, is the updated nth deviatoric stress; the stress tensor is obtained through Cauchy stress decomposition: where is given by the EOS state equation.

5. The method of claim 2, wherein: In S4, the influence of temperature on the mechanical response is calculated through the state equation, thermal coupling is realized by using a Mie-Grüneisen EOS formula, and the stress expression is corrected as follows: ； ； ； wherein is the Hugoniot linear stress, is the Grüneisen coefficient, , is the non-linear correction term, denotes the volumetric strain of the material, is the current density of the particles, in TLSPH the mass conservation equation is written as i.e. where the Jacobian matrix represents the ratio between the volume after the deformation and the volume before the deformation, wherein represents the determinant of ; denotes the solid sound speed, , , are empirical coefficients; When using the Mie-Grüneisen EOS formula, the effect of the thermal expansion of the structure does not need to be calculated separately, but is already included in the enthalpy value , i.e. the change in the internal energy of the particles , influences the stress state of the particles; when , i.e. the structure is in a state of compression; when , i.e. the structure is in a state of expansion.

Images