Liquid film oscillation mechanism and fluctuation stability calculation method and device based on particle method

Abstract

The invention provides a liquid film oscillation mechanism and fluctuation stability calculation method and device based on a particle method. The method provided by the invention comprises the following steps: firstly, constructing a calculation model by arranging five types of particles, namely, sole, fluid, wall, dummy and ghost; secondly, applying a net heat flow boundary considering wall enthalpy correction, radiation heat dissipation and a gas injection factor based on a virtual particle method to realize heat transfer and phase change calculation; then solving a viscosity item and a shear force by adopting a fully implicit method, and updating a particle motion state by combining a non-Newton flow power law model; then, a potent model is introduced to calculate surface tension and overload influence; solving a pressure field through a pressure Poisson equation and iteratively updating particle positions; and finally, dynamically marking ghost particles at the outlet to optimize the calculation efficiency. According to the method, the oscillation phenomenon of the liquid film can be accurately simulated, the generation mechanism and stability influence factors of the liquid film can be revealed, theoretical support and technical means are provided for thermal protection design of the silicon-based material, and the technical problem that the three-dimensional dynamic ablation morphology cannot be analyzed through a traditional method is solved.

Description

TECHNICAL FIELD

[0001] The present document relates to the technical field of liquid film oscillation mechanism and wave stability calculation, and particularly relates to a liquid film oscillation mechanism and wave stability calculation method and device based on a particle method. BACKGROUND

[0002] The thermal protection design of intercontinental missiles and reentry vehicles has always been a core challenge in the field of national defense science and technology. The ablative thermal protection technology has become a key solution due to its reliability and maturity. As the first generation of reentry warhead thermal protection materials, silicon-based materials are still widely used in aircraft, solid rocket engine nozzles, and antenna windows / shields of strategic missiles due to their low thermal conductivity, excellent ablative and thermal insulation performance. However, with the development of new weapons such as strategic glide missiles, the non-ablation requirements in long-time, low-heat-flow environments have put higher demands on the ablation behavior prediction of silicon-based materials (such as SiO2 transparent materials). A small error in surface chemical reaction mass loss may accumulate in long-time aerodynamic heating, significantly affecting the wave transmission performance of the antenna cover and the guidance accuracy, so it is urgent to study the ablation heat transfer mechanism of silicon-based materials.

[0003] The simulation method of silicon-based material high-temperature molten product loss started from Adams' study on the flow of pure quartz. He solved the ablation velocity by simplifying the boundary layer equation and gave an analytical solution for the liquid layer velocity at the stagnation point. Hidalgo further extended to the non-stagnation point condition, simplified the mass conservation equation by using the high viscosity characteristics of SiO2 melt, and formed an algebraic relationship. Domestic research on the ablation mechanism of silicon-based materials was carried out based on Adams' theory and applied to the thermal protection design of the type. Subsequently, under the support of national projects, the liquid-solid coupling simulation method was proposed by aerospace agencies, which calculated the liquid layer thickness through energy conservation, realized the unified solution of liquid flow and solid heat transfer, and improved the prediction accuracy of thermal response.

[0004] However, there are still unresolved key problems in existing research. Experiments have found that high-temperature molten liquid layer will form a corrugated flow under aerodynamic heating. This dynamic morphology may significantly affect the ablation prediction and wave transmission performance through the environment-flow coupling effect. The current understanding of the mechanism of SiO2 melt corrugated flow is still insufficient, and there is a lack of relevant model and experimental data support, which makes it difficult to accurately evaluate its influence on the performance of the thermal protection system. Therefore, it is urgent to develop new analysis methods or technologies to reveal the melt corrugated flow law and optimize the ablation prediction model of silicon-based materials in complex thermodynamic environments. SUMMARY

[0005] Embodiments of the present application provide a liquid film oscillation mechanism and wave stability calculation method and device based on a particle method, which aims to solve the above technical problems.

[0006] According to the embodiment of the present application, a liquid film oscillation mechanism and wave stability calculation method based on a particle method are provided, comprising:

[0007] S1, a calculation model is constructed by using a particle method, different types of particles are arranged to simulate physical boundary conditions, and aerodynamic heating, shear force loading and particle flow control are realized;

[0008] S2, a net heat flow boundary condition is applied to the model, the net heat flow boundary condition considers the influence of wall enthalpy correction, radiation heat dissipation and gas injection factor, heat transfer calculation is performed, and the phase change state of the particle is determined according to the temperature: when the particle temperature reaches the melting point, it is converted into a fluid particle, otherwise it remains a solid particle, wherein the viscosity coefficient of the fluid particle is calculated by a dynamic calculation method varying with temperature;

[0009] S3, based on the phase change determination result of S2, the viscosity term and the change of particle velocity and position caused by shear force are calculated for the fluid particle to obtain the first update of particle velocity and position;

[0010] S4, according to the updated particle distribution, the acceleration caused by surface tension and overload is calculated, and the temporary velocity and position of the particle are updated;

[0011] S5, the pressure field is solved based on the temporary velocity field, the pressure gradient force is calculated, and the second update of the particle velocity and position is completed;

[0012] S6, the outlet boundary treatment is performed, when the particle flows out of the calculation domain boundary, it is marked as a ghost particle and excluded in the subsequent calculation;

[0013] S7, S2 to S6 are executed in a loop with a preset time step as a unit, and particle state data is output every time step until a set calculation time is reached;

[0014] S8, by loading different types of boundary conditions, the dominant influencing factors and generation mechanism of liquid film oscillation are analyzed.

[0015] According to the embodiment of the present application, a liquid film oscillation mechanism and wave stability calculation device based on a particle method are provided, comprising:

[0016] The model construction module constructs a calculation model by using a particle method, simulates physical boundary conditions by arranging different types of particles, and realizes aerodynamic heating, shear force loading and particle flow control;

[0017] The heat flow and phase change module applies a net heat flow boundary condition to the model, the net heat flow boundary condition considers the influence of wall enthalpy correction, radiation heat dissipation and gas injection factor, performs heat transfer calculation, and determines the phase change state of the particle according to the temperature: when the particle temperature reaches the melting point, it is converted into a fluid particle, otherwise it remains a solid particle, wherein the viscosity coefficient of the fluid particle is calculated by a dynamic calculation method varying with temperature;

[0018] a viscosity force calculation module, configured to calculate the viscosity term and the particle velocity and position changes caused by the shear force based on the phase change determination result of the heat flow and phase change module, to obtain the first update of the particle velocity and position;

[0019] a tension and overload influence module, configured to calculate the acceleration caused by the surface tension and the overload according to the updated particle distribution, and update the particle temporary velocity and position;

[0020] a pressure influence module, configured to solve the pressure field based on the temporary velocity field, calculate the pressure gradient force, and complete the second update of the particle velocity and position;

[0021] an exclusion module, configured to perform the outlet boundary processing, and mark the particles flowing out of the calculation domain boundary as ghost particles and exclude them in the subsequent calculation;

[0022] an iteration execution module, configured to be executed in a preset time step unit, and the heat flow and phase change module, the viscosity force calculation module, the tension and overload influence module, the pressure influence module and the exclusion module are executed in each time step, and the particle state data is outputted in each time step until the set calculation time is reached;

[0023] an analysis module, configured to analyze the dominant influencing factors and generation mechanism of the liquid film oscillation by loading different types of boundary conditions.

[0024] By adopting the embodiment of the present application, the heat transfer, phase change and molten liquid layer flow process of the silicon-based material under the condition of aerodynamic heating are accurately simulated through the virtual particle technology and the multi-physical field coupling modeling. The method innovatively combines the key factors such as surface tension, non-Newtonian fluid effect, shear force and dynamic phase change, adopts full-implicit numerical solution and intelligent particle type switching, and effectively solves the problem that the traditional method cannot depict the three-dimensional dynamic topography. The advantage is that the oscillation behavior and stability of the molten liquid layer are accurately predicted, which provides an efficient and reliable theoretical tool for the thermal protection design of the aircraft, the wind tunnel test analysis and the performance optimization of the wave-transparent material, and significantly improves the simulation capability of the ablation behavior of the silicon-based material under the complex thermodynamic environment. BRIEF DESCRIPTION OF DRAWINGS

[0025] In order to more clearly illustrate the technical solutions in the one or more embodiments of the present application or the prior art, the drawings needed to be used in the embodiment or prior art description will be briefly introduced. Obviously, the drawings in the following description are only some embodiments described in the present application, and other drawings can be obtained by those skilled in the art without creative labor.

[0026] Figure 1A flow chart of the liquid film oscillation mechanism and wave stability calculation method based on the particle method according to the embodiment of the present application is shown in FIG. 1.

[0027] Figure 2 A specific implementation flow chart of the embodiment of the present application is shown in FIG. 2.

[0028] Figure 3 A force schematic diagram of a flat plate model under aerodynamic heating according to the embodiment of the present application is shown in FIG. 3.

[0029] Figure 4 A pressure action schematic diagram of a silicon-based material flat plate model after surface ablation morphology appears according to the embodiment of the present application is shown in FIG. 4.

[0030] Figure 5 A wave condition comparison schematic diagram of liquid layer particles under certain aerodynamic heating conditions according to the embodiment of the present application is shown in FIG. 5.

[0031] Figure 6 A schematic diagram of the liquid film oscillation mechanism and wave stability calculation device based on the particle method according to the embodiment of the present application is shown in FIG. 6. DETAILED DESCRIPTION

[0032] In order to enable the person skilled in the art to better understand the technical solutions in one or more embodiments of the present specification, the technical solutions in one or more embodiments of the present specification will be described clearly and completely in conjunction with the drawings in one or more embodiments of the present specification. Obviously, the described embodiments are only a part of the embodiments of the present specification, not all. Based on one or more embodiments of the present specification, all other embodiments obtained by the person skilled in the art without creative labor should belong to the protection scope of the present document.

[0033] Method embodiment

[0034] According to the embodiment of the present application, a liquid film oscillation mechanism and wave stability calculation method based on the particle method are provided, Figure 1 A flow chart of the liquid film oscillation mechanism and wave stability calculation method based on the particle method according to the embodiment of the present application is shown in FIG. 1. Figure 1 The liquid film oscillation mechanism and wave stability calculation method based on the particle method according to the embodiment of the present application specifically includes:

[0035] S1, a calculation model is constructed by using the particle method, different types of particles are arranged to simulate physical boundary conditions, and aerodynamic heating, shear force loading and particle flow control are realized;

[0036] The types of the particles include solid particles, fluid particles, wall particles, dummy particles and ghost particles, the solid particles are used to form a model body, the wall particles are arranged around the solid particles except the aerodynamic heating surface and the outlet, and are used to simulate the solid boundary of the actual wind tunnel test model, the dummy particles are arranged at the outlet, and are used to ignore the collision between the liquid particles and the outlet particles.

[0037] S2, a net heat flow boundary condition is applied to the model, the net heat flow boundary condition considers the wall enthalpy correction, radiation heat dissipation and gas entrainment factor influence, heat transfer calculation is performed, and the phase change state of the particle is determined according to the temperature: when the particle temperature reaches the melting point, the particle is converted into a fluid particle, otherwise the particle remains a solid particle, wherein the viscosity coefficient of the fluid particle is calculated by using a dynamic calculation method changing with the temperature;

[0038] By referring to the virtual particle idea in the MPS method, a method of applying a net heat flow boundary considering the wall enthalpy correction, radiation heat dissipation and gas entrainment factor influence to the surface is proposed, heat transfer, phase change and viscosity coefficient changing with the temperature are calculated, when the temperature reaches the melting point, the particle type is fluid, otherwise it is solid, the viscosity coefficient of the fluid is calculated by using the viscosity coefficient expression of the silicon-based material changing with the temperature, and the specific expression is as follows:

[0039] The aerodynamic heating and heat transfer calculation are performed:

[0040]

[0041] In the formula, T is the temperature of the particle, Tm is the melting point of the particle, and μ is the viscosity coefficient of the particle. <n> i and <n' i are the particle number density and the particle number density considering virtual particles, respectively, and their calculation formulas are n' i i , for surface particles, n' 0 i i i i Thus, the aerodynamic heating is realized only on the surface by using the particle number density.

[0042] q n (i) is the net heat flux into the surface, and its expression is

[0043] q n (i) = q or · ψ · (1 - h w (i) / h r ) - εσ[(T w (i)) 4 - T envir 4 ] Formula 2;

[0044] where q or is the cold-wall heat flux, h r is the total enthalpy of the gas flow, h w is the wall enthalpy, ε is the radiation coefficient, σ is the Stefan-Boltzmann constant, and ψ is the entrainment factor, whose expression can be obtained from the derivation process in the silicon-based ablation engineering solution method:

[0045]

[0046]

[0047] P v = 10 5 · exp[18.48 - (57780 / T w )] is the vapor pressure, P g is the external gas pressure. M av is the ratio of the air molecular weight to the SiO2 molecular weight, which is taken as 0.72. S f is equal to 0.62 and 0.2 for laminar flow and turbulent flow, respectively;

[0048] Phase change calculation: taking the solid melting temperature as the boundary, when the temperature is greater than T melt , the particle type is changed to fluid, and the flow calculation is started; when the temperature is less than T melt , the particle type is solid and fixed. ​​​​​

[0049] S3, based on the phase change determination result of S2, calculate the changes of particle velocity and position caused by viscosity and shear force, and obtain the first update of particle velocity and position;

[0050] The changes of particle velocity and position caused by viscosity and shear force are calculated, the application method of surface shear force is proposed by referring to the virtual particle idea in MPS method, and the full implicit method is used for the calculation of viscosity to be suitable for high viscosity fluid. The influence of power law model of non-Newtonian flow is considered in the calculation process, which is specifically:

[0051] (3.1) calculate the changes of particle velocity and position caused by viscosity and shear by using implicit method;

[0052]

[0053] In the formula The shear force in x direction is calculated by using engineering calculation method:

[0054]

[0055] The velocity of the outer edge of the boundary layer is

[0056] (3.2) consider the non-Newtonian flow effect by using power law model;

[0057]

[0058] S4, according to the updated particle distribution, calculate the acceleration caused by surface tension and overload, and update the temporary velocity and position of particles;

[0059] The potential-based model is used to consider the influence of surface tension, and the gravity is considered to calculate the new velocity distribution and update the particle position, which is specifically:

[0060] (4.1) calculation of acceleration caused by surface tension

[0061]

[0062]

[0063]

[0064]

[0065]

[0066] (4.2) update of temporary velocity and position of particles under the action of surface tension and overload

[0067] Surface tension and overloading induced temporary velocity update:

[0068]

[0069] where is the acceleration of surface tension, is the acceleration of overloading.

[0070] Particle temporary position calculation:

[0071]

[0072] After the temporary position is obtained, the temporary value of particle number density is calculated:

[0073]

[0074] S5, solve the pressure field based on the temporary velocity field, calculate the pressure gradient force and complete the second update of particle velocity and position;

[0075] Solve the pressure Poisson equation to calculate the pressure distribution, and calculate the particle velocity and position change caused by the pressure, specifically:

[0076] (5.1) The pressure of each particle is obtained by solving the pressure Poisson equation. The pressure Poisson equation based on the virtual particle is:

[0077]

[0078] (5.2) From the obtained pressure distribution, the velocity distribution and position distribution of the new time step are calculated:

[0079]

[0080]

[0081] S6, perform outlet boundary processing, when the particle flows out of the calculation domain boundary, mark it as a ghost particle and exclude it in subsequent calculations;

[0082] Perform particle judgment at the outlet. After the particle flows out of the boundary, change the particle type to ghost particle. Do not consider ghost particles in all calculation and output processes. Specifically:

[0083] Taking a flat plate as an example, after the particle flows out of the boundary, change the particle type to ghost particle. Do not consider the particle in all calculation processes, that is, move the particle out of the calculation domain, avoid the abnormal increase of particle number density at the outlet, and the abnormal movement of particles.

[0084] S7, performing S2 to S6 in a preset time step unit, outputting particle state data in each time step until reaching a set calculation time; specifically comprising: outputting the coordinates, pressure, temperature and the like of the particles every certain time in the calculation process, and automatically exiting the program after reaching the calculation time, so that the calculation results of each time can be viewed, and the change in the calculation process can be observed by making an animation.

[0085] S8, analyzing the dominant influencing factors and generation mechanism of the liquid film oscillation by loading different types of boundary conditions.

[0086] Figure 2 A specific implementation flowchart of the embodiment of the application, Figure 3 A force schematic diagram under the aerodynamic heating environment of the flat plate model, Figure 4 A pressure action schematic diagram after the surface ablation morphology of the silicon-based material flat plate model appears, Figure 5 Under certain aerodynamic heating conditions, the wave conditions of the liquid layer particles are compared, combined with Figures 2-5 The specific steps of a specific implementation of the embodiment of the application are described as follows:

[0087] (1) For the flat plate model of the wind tunnel test, the model is arranged by the particle method, including solid particles, wall particles and dummy particles.

[0088] (2) Set the initial conditions, including the initial spacing of the particles, the initial temperature, the specific heat capacity, the density, the thermal conductivity of the material, the surface tension coefficient, the surface radiation coefficient and the like.

[0089] (3) By referring to the virtual particle idea in the MPS method, the net heat flow boundary considering the wall enthalpy correction, radiation heat dissipation and gas injection factor is applied to the surface, and the heat transfer, phase change and viscosity coefficient changing with temperature are calculated, when the temperature reaches the melting point, the particle type is fluid, otherwise it is solid, the viscosity coefficient of the fluid is calculated by using the viscosity coefficient expression of the silicon-based material changing with temperature, such as the viscosity coefficient of high-silicon oxide μ = 0.01exp(68890 / T w -20).

[0090] (2.1) Set the heat flow, total enthalpy of the gas flow and total pressure conditions;

[0091] (2.2) Phase change calculation: taking the solid melting temperature as the boundary, when the temperature is greater than T melt , the particle type is changed to fluid, and the flow calculation is started, when the temperature is less than T melt , the particle type is solid and fixed;

[0092] (4) Calculate the change of particle velocity and position caused by viscosity and shear force, the calculation of surface shear force is improved by the virtual particle idea in MPS method, the calculation of viscosity term uses full implicit method, which is suitable for high viscosity fluid, and the influence of power law model of non-Newtonian flow is considered in the calculation process.

[0093] (5) Use potential based model to consider the influence of surface tension, and consider the effect of gravity, calculate the new velocity distribution, and update the particle position.

[0094] (6) Solve the pressure Poisson equation to calculate the pressure distribution, and calculate the particle velocity and position change caused by pressure.

[0095] (7) Perform the exit particle judgment, change the particle type to ghost particle after the particle flows out of the boundary, and do not consider the ghost particle in all calculation and output processes, specifically:

[0096] Take a flat plate as an example, change the particle type to ghost particle after the particle flows out of the boundary, and do not consider the particle in all calculation processes, i.e. remove the particle from the calculation domain, avoid the abnormal increase of particle number density at the outlet and the abnormal movement of particles.

[0097] (8) Output the coordinates, pressure, temperature, particle number density, etc. of the particles every certain time step in the calculation process, and the program automatically exits after the calculation time, so that the calculation results at each time can be viewed, and the change in the calculation process can be observed by making an animation.

[0098] (9) Load the gas flow pressure, gradient distribution of gas flow pressure, y-direction pressure gradient and gradient distribution of y-direction pressure gradient boundary conditions respectively, and discuss the generation mechanism of liquid film oscillation phenomenon.

[0099] The gas flow pressure is applied by formula 18, P g is the gas flow pressure; the gradient distribution of gas flow pressure is applied by formula 19, in which gradient_P g = P g / x i ; the y-direction pressure gradient is applied by formulas 20, 21 and 23; and the gradient distribution of y-direction pressure is applied by formulas 20, 22 and 23.

[0100]

[0101]

[0102]

[0103]

[0104]

[0105]

[0106] With the embodiment of the present application, the following beneficial effects are achieved: through virtual particle technology and multi-physical field coupling modeling, the heat transfer, phase change and molten liquid laminar flow process of silicon-based materials under aerodynamic heating conditions are accurately simulated. This method innovatively combines key factors such as surface tension, non-Newtonian fluid effect, shear force and dynamic phase change, adopts full-implicit numerical solution and intelligent particle type switching, effectively solving the problem that traditional methods cannot depict three-dimensional dynamic topography. The advantage is that it can accurately predict the oscillation behavior and stability of the molten liquid layer, providing an efficient and reliable theoretical tool for aircraft thermal protection design, wind tunnel test analysis and wave-transparent material performance optimization, significantly improving the simulation capability of silicon-based material ablation behavior under complex thermodynamic environment.

[0107] Device embodiment

[0108] According to the embodiment of the present application, a liquid film oscillation mechanism and wave stability calculation device based on particle method is provided, Figure 2 The schematic diagram of the liquid film oscillation mechanism and wave stability calculation device based on particle method of the embodiment of the present application is shown in Figure 2 The liquid film oscillation mechanism and wave stability calculation device based on particle method of the embodiment of the present application specifically includes:

[0109] The model construction module 60 adopts particle method to construct a calculation model, simulates physical boundary conditions by arranging different types of particles, and realizes aerodynamic heating, shear force loading and particle flow control;

[0110] The heat flow and phase change module 61 applies a net heat flow boundary condition to the model, which considers the influence of wall enthalpy correction, radiation heat dissipation and gas injection factor, performs heat transfer calculation and determines the phase change state of the particles according to the temperature: when the particle temperature reaches the melting point, it is converted into a fluid particle, otherwise it remains a solid particle, wherein the viscosity coefficient of the fluid particle is calculated dynamically with temperature;

[0111] The viscous force calculation module 62 calculates the change of particle velocity and position caused by viscous term and shear force based on the phase change determination result of the heat flow and phase change module, to obtain the first update of particle velocity and position;

[0112] The tension and overload influence module 63 is used to calculate the acceleration caused by surface tension and overload according to the updated particle distribution, and update the temporary velocity and position of the particles;

[0113] a pressure influence module 64 for solving pressure field based on the temporary velocity field, calculating pressure gradient force and completing the second update of particle velocity and position;

[0114] an exclusion module 65 for performing outlet boundary processing, marking the particles flowing out of the calculation domain boundary as ghost particles and excluding them in subsequent calculation;

[0115] an iteration execution module 66 for cyclically executing the heat flow and phase change module, the viscous force calculation module, the tension and overload influence module, the pressure influence module and the exclusion module in preset time steps, and outputting particle state data in each time step until reaching the set calculation time;

[0116] an analysis module 67 for analyzing the dominant influence factors and generation mechanism of liquid film oscillation by loading different types of boundary conditions.

[0117] Finally, it should be noted that: the above embodiments are only used to illustrate the technical solutions of the present application, and not to limit them; although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that: it can still modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacement for part or all of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the present application.< / n>

Claims

1. A method for calculating a liquid membrane oscillation mechanism and wave stability based on a particle method, characterized by The method comprises the following steps: S1, a particle method is used to construct a calculation model, different types of particles are arranged to simulate physical boundary conditions, and aerodynamic heating, shear force loading and particle flow control are realized; S2, a net heat flow boundary condition is applied to the model, the net heat flow boundary condition considers the influence of wall enthalpy correction, radiation heat dissipation and gas injection factor, heat transfer calculation is performed, and the phase change state of the particles is determined according to the temperature: when the particle temperature reaches the melting point, the particle is converted into a fluid particle, otherwise the particle remains a solid particle, wherein the viscosity coefficient of the fluid particle is calculated by a dynamic calculation method varying with temperature; S3, based on the phase change determination result of S2, the viscosity term and the change of particle velocity and position caused by shear force are calculated for the fluid particle to obtain the first update of the particle velocity and position; S4, according to the updated particle distribution, the acceleration caused by surface tension and overload is calculated, and the temporary velocity and position of the particle are updated; S5, the pressure field is solved based on the temporary velocity field, the pressure gradient force is calculated, and the second update of the particle velocity and position is completed; S6, the outlet boundary treatment is performed, when the particle flows out of the calculation domain boundary, it is marked as a ghost particle and excluded in the subsequent calculation; S7, steps S2 to S6 are executed in a preset time step, and particle state data is output every time step until the set calculation time is reached; S8, different types of boundary conditions are loaded to analyze the dominant influence factors and generation mechanism of liquid film oscillation.

2. The method of claim 1, wherein, The types of particles include solid particles, fluid particles, wall particles, dummy particles and ghost particles, the solid particles are used to form the main body of the model, the wall particles are arranged around the solid particles except the aerodynamic heating surface and the outlet, and are used to simulate the solid boundary of the actual wind tunnel test model, the dummy particles are arranged at the outlet, and the collision between the liquid particles and the outlet particles is ignored.

3. The method of claim 1, wherein, In step S2, the heat transfer calculation step comprises: Aerodynamic heating and heat transfer calculation is performed by formula 1: In the formulae <n> i and (n i are the particle number density and the particle number density considering virtual particles, respectively, and the calculation formulas are n i i , for surface particles, n 0 i > <n> i , internal particle <n' i = <n> i ；< / n> < / n> ​​< / n> q n (i) the net heat flow into the surface, expressed by: q n (i) = q or • ψ • (1 - h w (i) / h r )- εσ[(T w (i)) 4 - T envir 4 ] Equation 2; where q or is the cold wall heat flux, h r is the total enthalpy of the gas flow, h w is the wall enthalpy, ε is the emissivity, σ is the Stefan-Boltzmann constant, and ψ is the entrainment factor, whose expression can be obtained from the derivation process in the silicon-based ablation engineering solution method as follows: where P v = 10 5 exp[18.48-(57780 / T w )] is the vapor pressure, P g is the external gas pressure, M av is the ratio of the air molecular weight to the SiO2 molecular weight, taken as 0.72, and S f is equal to 0.62 and 0.2 for laminar and turbulent flow, respectively.

4. The method of claim 1, wherein, The phase change calculation specifically includes: with solid melting temperature as a boundary, when temperature is greater than T melt , particle type is changed to fluid, and flow calculation is started, when temperature is less than T melt , particle type is solid, and is fixed, and T melt is melting point.

5. The method of claim 1, wherein, The calculation of the viscosity term adopts a fully implicit method, which is suitable for high viscosity fluid, and the influence of the power law model of non-Newtonian flow is considered in the calculation process. S4 specifically comprises:

6. The method of claim 1, wherein, First, the acceleration caused by surface tension is calculated, the interaction force between surface particles is considered based on the potential model, then the comprehensive acceleration under the joint action of surface tension and overload is calculated combined with the given overload condition, and finally the temporary velocity and position of the particle are updated according to the acceleration, and the temporary value of the updated particle number density is calculated. S5 specifically comprises:

7. The method of claim 1, wherein, First, the pressure Poisson equation based on the virtual particle method is constructed, the pressure distribution of each particle in the calculation domain is solved, then the pressure gradient force is calculated according to the pressure gradient, and the velocity and position of the particle are updated based on the force field to complete the second update of the particle state. S6 specifically comprises:

8. The method of claim 1, wherein, ​ Real-time monitoring of particle position, when the particle flow out of the calculation domain boundary, immediately marked as ghost particles; In all subsequent calculations automatically ignore ghost particles, avoid export boundary particle number density anomalies and calculation distortion.

9. The method of claim 1, wherein, The S8 specifically includes: Four different types of boundary conditions are loaded in turn: uniform air flow pressure, gradient distribution air flow pressure, y pressure gradient and gradient distribution y pressure gradient; By comparing and analyzing the oscillation characteristics of liquid film under different boundary conditions, the key factors affecting the stability of liquid film and their action mechanism are determined.

10. A device for calculating a liquid membrane oscillation mechanism and wave stability based on a particle method, characterized by Including: The model construction module adopts the particle method to construct the calculation model, simulates the physical boundary conditions by arranging different types of particles, and realizes the aerodynamic heating, shear force loading and particle flow control; The heat flow and phase change module applies the net heat flow boundary condition to the model, which considers the influence of wall enthalpy correction, radiation heat dissipation and gas injection factor, carries out heat transfer calculation and determines the phase change state of particles according to temperature: when the particle temperature reaches the melting point, it is converted into fluid particles, otherwise it remains as solid particles, wherein the viscosity coefficient of fluid particles is calculated by a dynamic calculation method varying with temperature; The viscous force calculation module calculates the viscous term and the change of particle velocity and position caused by shear force based on the phase change determination result of the heat flow and phase change module, to obtain the first update of particle velocity and position; The tension and overload influence module is used to calculate the acceleration caused by surface tension and overload according to the updated particle distribution, and update the temporary velocity and position of particles; The pressure influence module is used to solve the pressure field based on the temporary velocity field, calculate the pressure gradient force and complete the second update of particle velocity and position; The exclusion module is used for outlet boundary processing, when the particle flows out of the calculation domain boundary, it is marked as ghost particle and excluded in subsequent calculation; The iteration execution module is used to execute in a preset time step unit, the heat flow and phase change module, the viscous force calculation module, the tension and overload influence module, the pressure influence module and the exclusion module, and the particle state data is output every time step until the set calculation time is reached; The analysis module is used to analyze the dominant influence factors and generation mechanism of liquid film oscillation by loading different types of boundary conditions.

Images