Ablation morphology simulation method for quartz material axial symmetry model

Abstract

An ablation morphology simulation method for an axial symmetry model of a quartz material comprises the following steps: (1) constructing an axial symmetry hemisphere model by using a sold particle and a dummy particle, and forming a computational domain; coordinates of all particles of the axisymmetric hemisphere model are collected, and a coordinate set is obtained; (2) collecting the heat conductivity, the density, the specific heat capacity and the initial temperature of each particle based on a preset pneumatic heating condition; (3) calculating the temperature of each particle according to the collected thermal conductivity, density, specific heat capacity and initial temperature of the particles; (4) calculating the speed and coordinates of each particle in the computational domain according to a preset shear force parameter, a preset post-based model, the temperature of each particle, the coordinate set and a preset time node; and (5) outputting the surface morphology of the axisymmetric hemisphere model according to the temperature, speed and coordinate of each particle, and completing simulation of the ablation morphology of the quartz material axisymmetric model. According to the method, the defect of research on the three-dimensional ablation morphology of the quartz material in the prior art is overcome.

Description

Technical Field

[0001] The invention relates to a method for simulating ablation morphology of an axisymmetric model of a quartz material, and belongs to the technical field of pneumatic thermal protection. Background Art

[0002] Thermal protection has always been one of the core technologies for hypersonic flight and the development of various aircraft. For aircraft that re-enter the atmosphere at high Mach numbers or cruise in the atmosphere for a long time, the thermal protection of the structure and its impact on the aerodynamic performance of the aircraft are always a key issue. As the head of a blunt-nosed aircraft, the end cap has the worst aerodynamic heating environment and is also one of the key areas of thermal protection for the aircraft.

[0003] Traditional ablation calculation is a one-dimensional calculation. According to the heat flux density and the physical and chemical properties of the material, the ablation amount in the normal direction of the vertical structure surface, that is, the ablation retreat distance, is obtained.

[0004] However, different parts of the aircraft have different protective layers and different degrees of ablation. The one-dimensional ablation prediction of the traditional algorithm can only capture the stagnation point ablation of axisymmetric models such as blunt-headed bodies, and it is difficult to capture the surface morphology of quartz composite materials at different positions during ablation. In addition, the momentum equation in the traditional silicon-based material solution method only considers the viscosity and pressure gradient terms, and cannot consider the unsteady effects of flow, surface tension effects, overload effects, asymmetric ablation effects, etc., resulting in the actual ablation retreat of the aircraft not matching the calculated ablation retreat, thereby affecting the accuracy of the ablation calculation. Summary of the invention

[0005] The purpose of the present invention is to overcome the shortcomings of the prior art and propose a method for simulating the ablation morphology of an axisymmetric model of a quartz material, thereby solving the shortcomings of the prior art in studying the three-dimensional ablation morphology of quartz materials.

[0006] The technical solution adopted by the present invention is: A method for simulating ablation morphology of an axisymmetric model of quartz material comprises the following steps: (1) Use solid particles and dummy particles to construct an axisymmetric hemispherical model to form a computational domain; and collect the coordinates of each particle in the axisymmetric hemispherical model to obtain a coordinate set; (2) Based on the preset aerodynamic heating conditions, the thermal conductivity, density, specific heat capacity and initial temperature of each particle are collected; (3) Calculate the temperature of each particle based on the collected thermal conductivity, density, specific heat capacity and initial temperature of the particles; (4) calculating the velocity and coordinates of each particle in the calculation domain according to a preset shear force parameter, a preset potential based model, the temperature of each particle, the coordinate set and a preset time node; (5) According to the temperature, velocity and coordinates of each particle, the surface morphology of the axisymmetric hemispherical model is output, and the simulation of the ablation morphology of the axisymmetric model of quartz material is completed.

[0007] Furthermore, the solid particles represent solid particles, simulating the state of the axisymmetric hemispherical model before heating; the dummy particles represent virtual particles, which are used to ensure that aerodynamic heating conditions are only applied on the hemispherical surface.

[0008] Furthermore, when creating an axisymmetric hemispherical model, solid particles are used to arrange the shape of the axisymmetric hemispherical model, and dummy particles are arranged behind the cross section of the hemispherical model to prevent the solid particles at the cross section of the hemispherical model from being affected by aerodynamic heating during the calculation process; based on the established axisymmetric hemispherical model, the coordinates of each solid particle are collected to obtain a coordinate set.

[0009] Furthermore, the thermal conductivity, density, specific heat capacity and initial temperature of each particle are collected, specifically: Aerodynamic heating is applied to the heating surface of the axisymmetric hemispherical model in the computational domain, and a fixed temperature is set. During the aerodynamic heating process, the temperature of the solid particles gradually rises, and the thermal conductivity, density, specific heat capacity and initial temperature of each solid particle are collected over a period of time.

[0010] Furthermore, the temperature of each particle is calculated based on the collected thermal conductivity, density, specific heat capacity and initial temperature of the particles, specifically:

[0011] in, represents the temperature of the particle in the i-th axisymmetric hemispherical model at the k+1th step; represents the temperature of the i-th particle at the k-th step; the j-th particle represents all particles around the i-th particle and different from the i-th particle; represents the temperature of the jth particle at the kth step; represents the time step; represents the thermal conductivity of the ith particle; represents the density of the i-th particle; represents the specific heat capacity of the ith particle; Represents dimension; , are the position vectors of the jth and ith particles respectively; then Represents the distance between particles j and i ; is the effective radius, take ; represents the standard particle number density at the initial moment; represents the weight function; represents the position vector of particle j; represents the position vector of particle i; represents the net heat flow into the surface; represents the spatial step length, i.e., the distance between particles at the initial moment; represents the particle number density after considering virtual particles; represents the particle number density, ; For surface particles, , internal particles , the particle number density was used to achieve aerodynamic heating only on the surface; It represents the average value of the square of the distance between particles.

[0012] Furthermore, the net heat flux into the particle surface Calculate using the following method: ; in, is the cold wall heat flow; is the total enthalpy of the air flow; is the wall enthalpy; is the radiation coefficient; is the Stefan Boltzmann constant; is the elicitation factor; represents the temperature of the particle; Indicates the ambient temperature, is the loss rate, is the evaporation rate, is the heat loss effect, is the evaporation heat effect, is the coefficient of Lees distribution of Lees spherical heat flux, is the spherical center angle; is the coefficient of Lees distribution of Lees spherical heat flux, is the sphere center angle, The calculation is done using the following formula:

[0013] .

[0014] in is the specific heat ratio of the gas, is the incoming flow Mach number; Emission Factor Calculate using the following method:

[0015] in, is the mass flow rate of the ejected gas generated by carbon ablation and pyrolysis of the original material, is the pyrolysis gas mass flow rate; The mass flow rate of carbon monoxide gas produced by carbon oxidation; Generates gas mass flow rate for carbon sublimation; is the vapor pressure; Pe is the external gas pressure; is the ratio of the molecular weight of air to the molecular weight of SiO2; for laminar and turbulent flow, Equal to 0.62 and 0.2 respectively.

[0016] Furthermore, when the particle temperature reaches the melting point, the particle type changes from solid to fluid; when the particle temperature does not reach the melting point, the particle type is solid.

[0017] Setting the solid melting temperature is the temperature threshold. When the particle temperature is greater than or equal to When the particle type is changed from solid particle to fluid particle, flow calculation begins. , the particle type is solid particle, fixed and motionless.

[0018] Furthermore, the calculation of the velocity and coordinates of each particle in the calculation domain according to the preset shear force parameter, the preset potential based model, the temperature of each particle, the coordinate set and the preset time node specifically includes: (4.1) Based on the particle’s current velocity and coordinate set, obtain the Laplacian operator of velocity; (4.2) According to the temperature of each particle, collect the viscosity coefficient of each particle; (4.3) According to the viscosity coefficient of each particle, the Laplace operator and the preset shear force parameter, the first temporary velocity and the first temporary coordinate generated by the viscosity term are obtained; (4.4) According to the preset potential based model, the surface tension of each particle is obtained; (4.5) obtaining a second temporary velocity and a second temporary coordinate of each particle according to the surface tension, the preset time node, the preset gravity, the first temporary velocity and the first temporary coordinate; (4.6) Calculate the pressure of each particle based on the second temporary velocity and the second temporary coordinate of each particle; (4.7) Based on the pressure of each particle, the velocity and coordinates of each particle are calculated.

[0019] Furthermore, The Laplace operator of the velocity is obtained based on the current velocity and coordinate set of the particle, specifically:

[0020] represents the velocity vector of particle i at the kth time step; represents the velocity vector of particle j at the kth time step; represents the position vector of particle j; represents the position vector of particle i; represents the average value of the square of the distance between particles; represents the weight function; Represents the distance between particles j and i ; is the effective radius, take ; represents the standard particle number density at the initial moment; The viscosity coefficient of each particle ;

[0021] A, B, C are constants and have different values ​​for different quartz materials. Represents the temperature of the particle.

[0022] Furthermore, the step (4.3) obtains the first temporary velocity and the first temporary coordinate generated by the viscosity term according to the viscosity coefficient of each particle, the Laplace operator and the preset shear force parameter. The calculation method is as follows:

[0023]

[0024] in, is the coordinate component of the tangential shear force of the spherical head model; Represents the first temporary velocity vector of particle i; Represents the first temporary coordinate vector of particle i; represents the first temporary velocity vector of particle j; represents the first temporary coordinate vector of particle j; represents the velocity vector of particle i at the kth time step; represents the velocity vector of particle j at the kth time step; represents the coordinate vector of particle i at the kth time step; represents the time step; μ represents the dynamic viscosity coefficient; ρ represents the density of the particle; Indicates dimension, for three-dimensional problems, d=3; represents the average value of the square of the distance between particles i and j; represents the position vector of particle j; represents the position vector of particle i; Represents the distance between particles j and i ; is the effective radius, take ; represents the standard particle number density at the initial moment; represents the set of all particles around particle i excluding ghost particles; represents the weight function; represents the particle number density after considering virtual particles; represents the particle number density; represents the spatial step length, i.e., the distance between particles at the initial moment; The shear force is calculated using engineering calculation method, and the calculation method is as follows: ; in, is the velocity at the outer edge of the boundary layer represents the elicitation factor; represents the cold wall heat flux; represents the total enthalpy of air flow; represents the Prandtl number; Represents the angle between the tangent direction of the spherical head model and the incoming flow direction.

[0025] Furthermore, the step (4.4) obtains the surface tension of each particle according to the preset potential based model, specifically: The acceleration due to surface tension is calculated as follows: ; Where A represents the set of i; B represents the set of particles around particle i; represents surface tension; is the Stefan Boltzmann constant; represents the distance between particles i and j; represents the effective radius of the force between particles, here we take ; Indicates the minimum distance between particles at the initial moment; ; in, represents the potential energy between particles i and j; ; in, represents the sum of the potential energy of all particles around particle i acting on it; ; in, represents the resultant force of particles around particle i on it; represents the weight function; ; in, Represents the particle acceleration obtained based on the potential based model; Represents the mass of the particle.

[0026] Furthermore, the step (4.5) obtains the second temporary velocity and the second temporary coordinate of each particle according to the surface tension, the preset time node, the preset gravity, the first temporary velocity and the first temporary coordinate, specifically: Second temporary speed Calculate using the following method: ; in, is the acceleration of surface tension, is the acceleration due to gravity.

[0027] Second temporary coordinates Calculate using the following method: .

[0028] Furthermore, (4.6) calculates the pressure of each particle according to the second temporary velocity and the second temporary coordinate of each particle, specifically: According to the second temporary coordinate, the temporary value of the particle number density is calculated , the calculation method is as follows: ; The pressure is calculated using the following Poisson's equation for pressure: ; in, represents the pressure of particle j at step k+1; represents the pressure of particle i at step k+1; represents the empirical coefficient; represents the particle number density of virtual particles; Indicates the external gas pressure; represents the standard particle number density at the initial moment; Represents the gradient.

[0029] Furthermore, (4.7) calculates the velocity and coordinates of each particle according to the pressure of each particle, specifically: After obtaining the value of the pressure at the k+1 time step, the velocity and coordinates at the k+1 time step are obtained using the following method: ; ; in, represents the velocity vector of particle i at the k+1 time step; Represents the second temporary velocity vector of particle i; represents the particle density; represents the gradient of particle pressure; in, Represents the position vector of particle i at the k+1th step.

[0030] Furthermore, the surface morphology of the axisymmetric hemispherical model is output according to the temperature, velocity and coordinates of each particle, specifically: Output the type of each particle according to the temperature of each particle; Output the position information of each particle according to the coordinates and speed of each particle; According to the particle type, flow state and the position information, the surface morphology of the axisymmetric hemispherical model is output.

[0031] Furthermore, the particle number density is calculated based on the coordinates of each particle, and the particle number density is determined based on whether it is greater than , to determine whether the position information of each particle is within the calculation domain. When the particle number density is less than When the particle flows out of the boundary, it is considered that the particle flows out of the boundary. When the particle flows out of the boundary, the particle is moved out of the calculation domain and the type is changed to a ghost particle. The ghost particle is not considered in all calculation processes and output processes. When the particle is still within the boundary, the position information of each particle is output. The ghost particle is represented as a ghost particle, which simulates the fluid particle leaving the calculation domain after being subjected to shear force.

[0032] The beneficial effects of the present invention are: (1) The present invention creates axisymmetric spherical head models of different particles in a preset calculation domain, applies ablation behavior to the particles within a certain period of time, and causes some particles to leave the calculation domain to simulate the ablation state of the protective layer during flight, thereby determining the change process of the particles in the calculation domain. At the same time, since some particles do not need to be calculated when they leave the calculation domain, it is not necessary to calculate all particles, thereby avoiding large-scale calculations, thereby making the calculation efficiency higher and also improving the calculation accuracy.

[0033] (2) The present invention is based on an axisymmetric spherical head model and outputs the temperature, velocity and coordinates of the particles. It can determine the surface morphology of the axisymmetric model during ablation, so that the calculation results match the ablation scene, thereby improving the prediction accuracy of the material ablation behavior. BRIEF DESCRIPTION OF THE DRAWINGS Figure 1 Flow chart of the method of the present invention; Figure 2 Schematic diagram of the axisymmetric particle model before ablation; Figure 3 Schematic diagram of the axisymmetric particle model after ablation. DETAILED DESCRIPTION

[0035] When an aircraft is flying, the thermal protection layer at different parts of the aircraft surface is affected by factors such as temperature, pressure and speed, and the degree of ablation of the thermal protection layer of the quartz composite material is different at different positions. The one-dimensional ablation prediction of the traditional algorithm or the calculation of the ablation retreat amount at the stagnation point position of the ball head model cannot calculate the surface tension, overload, asymmetric ablation, etc. of the molten liquid layer; the existing technology makes it difficult to simulate the ablation retreat amount of different parts in the molten state based on calculations, which affects the prediction accuracy of the thermal protection layer of the quartz composite material.

[0036] In order to overcome the above problems, the present invention proposes an ablation simulation method for an axisymmetric model of a quartz composite material under pneumatic heating conditions.

[0037] In the embodiments of the present application, the fully implicit particle method MPT is obtained based on the improvement of the MPS method. Drawing on the idea of ​​the virtual particle method, a method is proposed that can apply net heat flow and shear force boundaries on the model surface to simulate the heat transfer, phase change, and molten liquid layer flow of quartz composite materials in a ground wind tunnel state or a flight state under aerodynamic heating conditions, and obtain a two-dimensional or three-dimensional ablation morphology of the surface, as well as the ablation retreat amount at different positions.

[0038] Specifically, in the computational domain, different types of particles are used to create a particle model, and a potentialbased model is used to apply an interaction force between particles so that the distance between particles can be kept within a certain range. The particle ablation state is then simulated based on aerodynamic heating and shear force. The particle model is in the shape of a hemispherical head, including a heating surface and an outlet surface. First, solid particles are used to arrange the hemispherical surface of the axisymmetric hemispherical head model, and dummy particles are used to arrange the rear of the cross section of the spherical head model to ensure that the heat flow only acts on the hemispherical surface, such as Figure 2 shown.

[0039] Specifically, by using pneumatic heating and shear force, heat flux density and airflow shear force are applied to the heating surface, and the temperature of the particles is calculated. When the temperature reaches the melting point, the particle type changes from solid particles to fluid particles; when the particles become liquid, they are affected by the airflow shear force and move tangentially along the spherical surface, and the particle type changes with the change of temperature; the viscosity coefficient and velocity of each particle are collected, and the viscosity term is calculated based on the velocity and viscosity coefficient; based on the viscosity term and shear force, the change in velocity and coordinate of each particle is calculated to obtain the first temporary velocity and the first temporary coordinate; based on the potential based model, the surface tension of each particle is obtained; using the surface tension and gravity, the second temporary velocity and temporary coordinate of each particle are calculated; the pressure of each particle is obtained by solving the pressure Poisson equation; the velocity and coordinate of each particle are calculated based on the second temporary velocity, the second temporary coordinate and the pressure.

[0040] Specifically, the particle number density of each particle is used to determine whether the particle is in the calculation domain. If so, the particle coordinates are output; if not, the particle type is changed to a ghost particle, and the particle is moved out of the calculation domain without being output.

[0041] Specifically, during the calculation process, the coordinates, velocity, temperature and other information of the particles are output at regular time steps. After the calculation time is reached, the program automatically exits, and the user can view the surface morphology of the particle model at any time.

[0042] like Figure 1 As shown, the present invention proposes an ablation simulation method for an axisymmetric model of a quartz composite material under pneumatic heating conditions, and the specific steps are: Step 1: Use solid particles and dummy particles to construct an axisymmetric hemispherical model to form a calculation domain; and collect the coordinates of each particle in the axisymmetric hemispherical model to obtain a coordinate set.

[0043] like Figure 2 As shown, when creating an axisymmetric hemispherical model, solid particles are used to arrange the shape of the axisymmetric hemispherical model, and dummy particles are arranged behind the cross section of the hemispherical model to prevent the solid particles at the cross section of the hemispherical model from being affected by aerodynamic heating during the calculation process; according to the established axisymmetric hemispherical model, the coordinates of each solid particle are collected to obtain a coordinate set.

[0044] In an embodiment of the present application, an axisymmetric hemispherical head model is created in a computational domain by using different types of particles, and the particles include solid particles, fluid particles, ghost particles, and dummy particles; wherein, solid particles represent solid particles, simulating the state of the axisymmetric hemispherical model before heating; fluid particles represent fluid particles, simulating that after the axisymmetric model is heated, as the temperature rises, the particles change from solid to fluid; ghost particles are represented as ghost particles, simulating that after the fluid particles are subjected to shear force, the particles leave the computational domain, and these particles no longer participate in the calculation; dummy particles represent virtual particles, simulating that the fluid particles in the axisymmetric model will collide with a part of the unchanged solid particles when approaching the computational domain, and these unchanged solid particles are set to virtual particles to avoid collision between these particles.

[0045] Step 2: Based on the preset pneumatic heating conditions, the thermal conductivity, density, specific heat capacity and initial temperature of each particle are collected; Aerodynamic heating is applied to the heating surface of the axisymmetric hemispherical model in the computational domain, and a fixed temperature is set. During the aerodynamic heating process, the temperature of the solid particles gradually rises, and the thermal conductivity, density, specific heat capacity and initial temperature of each solid particle are collected over a period of time.

[0046] For example: aerodynamic heating can be used to simulate the effect of temperature on the protective layer of an aircraft during flight. The aerodynamic heating conditions include the heat flow, airflow and pressure that the aircraft is subjected to within a certain period of time. When the aircraft is flying at high speed, it has a large relative speed with the airflow. When the gas molecules hit the aircraft, their speed will stagnate, and the kinetic energy of the gas will be converted into thermal energy, generating a high temperature of several thousand degrees Celsius. The high-temperature airflow will heat the surface of the aircraft, causing varying degrees of damage to different parts of the aircraft surface.

[0047] Step 3: Calculate the temperature of each particle based on the thermal conductivity, density, specific heat capacity and initial temperature of the collected particles; Calculate the temperature of each particle, specifically:

[0048] in, represents the temperature of the particle in the i-th axisymmetric hemispherical model at the k+1th step; represents the temperature of the i-th particle at the k-th step; the j-th particle represents all particles around the i-th particle and different from the i-th particle; represents the temperature of the jth particle at the kth step; represents the time step; represents the thermal conductivity of the ith particle; represents the density of the i-th particle; represents the specific heat capacity of the ith particle; Represents dimension; , are the position vectors of the jth and ith particles respectively; then Represents the distance between particles j and i ; is the effective radius, take ; represents the standard particle number density at the initial moment; represents the weight function; represents the position vector of particle j; represents the position vector of particle i; represents the net heat flow into the surface; represents the spatial step length, i.e., the distance between particles at the initial moment; represents the particle number density after considering virtual particles; represents the particle number density, ; For surface particles, , internal particles , the particle number density was used to achieve aerodynamic heating only on the surface; It represents the average value of the square of the distance between particles.

[0049] 6. The method for simulating the ablation morphology of an axisymmetric model of quartz material according to claim 5 is characterized in that: the net heat flux entering the particle surface Calculate using the following method: ; in, is the cold wall heat flow; is the total enthalpy of the air flow; is the wall enthalpy; is the radiation coefficient; is the Stefan Boltzmann constant; is the elicitation factor; represents the temperature of the particle; Indicates the ambient temperature, is the loss rate, is the evaporation rate, is the heat loss effect, is the evaporation heat effect, is the coefficient of Lees distribution of Lees spherical heat flux, is the spherical center angle; is the coefficient of Lees distribution of Lees spherical heat flux, is the sphere center angle, The calculation is done using the following formula:

[0050] .

[0051] in is the specific heat ratio of the gas, is the incoming flow Mach number; Emission Factor Calculate using the following method:

[0052] in, is the mass flow rate of the ejected gas generated by carbon ablation and pyrolysis of the original material, is the pyrolysis gas mass flow rate; The mass flow rate of carbon monoxide gas produced by carbon oxidation; Generates gas mass flow rate for carbon sublimation; is the vapor pressure; Pe is the external gas pressure; is the ratio of the molecular weight of air to the molecular weight of SiO2; for laminar and turbulent flow, Equal to 0.62 and 0.2 respectively.

[0053] Specifically, when the particle temperature reaches the melting point, the particle type changes from solid to fluid; when the particle temperature does not reach the melting point, the particle type is solid.

[0054] Setting the solid melting temperature is the temperature threshold. When the particle temperature is greater than or equal to When the particle type is changed from solid particle to fluid particle, flow calculation begins. , the particle type is solid particle, fixed and motionless.

[0055] Step 4, calculating the velocity and coordinates of each particle in the calculation domain according to a preset shear force parameter, a preset potential based model, the temperature of each particle, the coordinate set and a preset time node; In an embodiment of the present application, based on the current velocity and coordinate set of the particle, the Laplace operator of the velocity is obtained; according to the temperature of each particle, the viscosity coefficient of each particle is collected; according to the viscosity coefficient of each particle, the Laplace operator and the preset shear force parameter, the first temporary velocity and the first temporary coordinate generated by the viscosity term are obtained; according to the preset potential based model, the surface tension of each particle is obtained; according to the surface tension, the preset time node, the preset gravity, the first temporary velocity and the first temporary coordinate, the second temporary velocity and the second temporary coordinate of each particle are obtained; according to the second temporary velocity and the second temporary coordinate of each particle, the pressure of each particle is calculated; according to the pressure of each particle, the velocity and coordinate of each particle are output.

[0056] It should be noted that in actual scenarios, the speed and coordinates of each particle in the calculation domain need to be output once every cycle. Therefore, for each cycle, the current speed in this application is the speed of the particle calculated in the previous cycle. For the first cycle, the current speed is a preset initial value, usually 0.

[0057] Specifically, the current velocity is the velocity generated at each moment for the axisymmetric particle model after aerodynamic heating and shear force; the coordinate set is the change of each particle in the particle model after aerodynamic heating and shear force; according to the shear force parameter, the idea of ​​virtual particles in the MPS method is borrowed to apply a shear force to the upper surface of the axisymmetric model, and the first temporary velocity and first temporary coordinate of each particle are calculated according to the viscosity term.

[0058] The viscosity term is calculated using a fully implicit method to be suitable for high-viscosity fluids. The viscosity term is obtained from the viscosity coefficient and velocity. At the same time, the non-Newtonian flow effect will also affect the velocity of each particle during the calculation process, and the power law model is used to calculate the value of the non-Newtonian flow effect.

[0059] Specifically, according to the preset shear force parameter, the preset potential based model, the temperature of each particle, the coordinate set and the preset time node, the velocity and coordinate of each particle in the calculation domain are calculated, including the following steps: (4.1) Based on the particle’s current velocity and coordinate set, obtain the Laplacian operator of velocity;

[0060] represents the velocity vector of particle i at the kth time step; represents the velocity vector of particle j at the kth time step; represents the position vector of particle j; represents the position vector of particle i; represents the average value of the square of the distance between particles; represents the weight function; Represents the distance between particles j and i ; is the effective radius, take ; represents the standard particle number density at the initial moment; (4.2) According to the temperature of each particle, collect the viscosity coefficient of each particle ; Drawing on the idea of ​​virtual particles in the MPS method, a method is proposed to impose a net heat flux boundary on the surface that takes into account the surface wall enthalpy correction, radiation heat dissipation and gas ejection factor, and calculates heat transfer, phase change and viscosity coefficient that varies with temperature. When the temperature reaches the melting point, the particle type is fluid, otherwise it is solid. The viscosity coefficient of the fluid is calculated using the viscosity coefficient expression of silicon-based materials that varies with temperature.

[0061]

[0062] A, B, C are constants and have different values ​​for different quartz materials. Represents the temperature of the particle.

[0063] (4.3) According to the viscosity coefficient of each particle, the Laplace operator and the preset shear force parameter, the first temporary velocity and the first temporary coordinate generated by the viscosity term are obtained; The changes in particle velocity and position caused by the viscosity term and shear force are calculated. The idea of ​​virtual particles in the MPS method is used for reference, and a method for applying surface shear force is proposed. The viscosity term is calculated using a fully implicit method to be suitable for high-viscosity fluids.

[0064]

[0065]

[0066] in, is the coordinate component of the tangential shear force of the spherical head model; Represents the first temporary velocity vector of particle i; Represents the first temporary coordinate vector of particle i; represents the first temporary velocity vector of particle j; represents the first temporary coordinate vector of particle j; represents the velocity vector of particle i at the kth time step; represents the velocity vector of particle j at the kth time step; represents the coordinate vector of particle i at the kth time step; represents the time step; μ represents the dynamic viscosity coefficient; ρ represents the density of the particle; Indicates dimension, for three-dimensional problems, d=3; represents the average value of the square of the distance between particles i and j; represents the position vector of particle j; represents the position vector of particle i; Represents the distance between particles j and i ; is the effective radius, take ; represents the standard particle number density at the initial moment; represents the set of all particles around particle i excluding ghost particles; represents the weight function; represents the particle number density after considering virtual particles; represents the particle number density; represents the spatial step length, i.e., the distance between particles at the initial moment; The shear force is calculated using engineering calculation method, and the calculation method is as follows: ; in, is the velocity at the outer edge of the boundary layer represents the elicitation factor; represents the cold wall heat flux; represents the total enthalpy of air flow; represents the Prandtl number; Represents the angle between the tangent direction of the spherical head model and the incoming flow direction.

[0067] (4.4) According to the preset potential based model, the surface tension of each particle is obtained; The potential based model is used to consider the influence of surface tension and the effect of gravity, calculate the new velocity distribution, and update the particle position.

[0068] Surface tension is also called contraction force, which means that the forces on the molecules at the interface are different from those on the molecules in the bulk of the phase. The net force on a liquid molecule inside the liquid and the forces on the surrounding liquid molecules is zero, but this is not the case for liquid molecules on the surface. Because the attraction of the gas phase molecules in the upper space to it is less than the attraction of the internal liquid phase molecules to it, the net force on the molecule is not equal to zero, and the direction of the net force points vertically to the inside of the liquid, resulting in the liquid surface having a tendency to shrink automatically. There are two main ways to describe surface tension. For the grid method, the CSF surface tension model is mainly used. For the particle method of this application, this is a gridless method, and the potential based model can be used.

[0069] According to the potential based model, a certain amount of interaction force is applied to the particle model, making it difficult for the particles to disperse, presenting various models, such as the spherical head model.

[0070] The acceleration due to surface tension is calculated as follows: ; Where A represents the set of i; B represents the set of particles around particle i; represents surface tension; is the Stefan Boltzmann constant; represents the distance between particles i and j; represents the effective radius of the force between particles, here we take ; Indicates the minimum distance between particles at the initial moment; ; in, represents the potential energy between particles i and j; ; in, represents the sum of the potential energy of all particles around particle i acting on it; ; in, represents the resultant force of particles around particle i on it; represents the weight function; ; in, Represents the particle acceleration obtained based on the potential based model; Represents the mass of the particle.

[0071] (4.5) obtaining a second temporary velocity and a second temporary coordinate of each particle according to the surface tension, the preset time node, the preset gravity, the first temporary velocity and the first temporary coordinate; Second temporary speed Calculate using the following method: ; in, is the acceleration of surface tension, is the acceleration due to gravity.

[0072] Second temporary coordinates Calculate using the following method: .

[0073] (4.6) Calculate the pressure of each particle based on the second temporary velocity and the second temporary coordinate of each particle; According to the second temporary coordinate, the temporary value of the particle number density is calculated , the calculation method is as follows: ; The pressure is calculated using the following Poisson's equation for pressure: ; in, represents the pressure of particle j at step k+1; represents the pressure of particle i at step k+1; represents the empirical coefficient; represents the particle number density of virtual particles; Indicates the external gas pressure; represents the standard particle number density at the initial moment; Represents the gradient.

[0074] (4.7) Based on the pressure of each particle, the velocity and coordinates of each particle are calculated.

[0075] Specifically: After obtaining the value of the pressure at the k+1 time step, the velocity and coordinates at the k+1 time step are obtained using the following method: ; ; in, represents the velocity vector of particle i at the k+1 time step; Represents the second temporary velocity vector of particle i; represents the particle density; represents the gradient of particle pressure; in, Represents the position vector of particle i at the k+1th step.

[0076] Step 5: output the surface morphology of the axisymmetric hemispherical model according to the temperature, velocity and coordinates of each particle, and complete the simulation of the ablation morphology of the axisymmetric model of the quartz material.

[0077] include: Output the type of each particle according to the temperature of each particle; Output the position information of each particle according to the coordinates and speed of each particle; According to the particle type, flow state and the position information, the surface morphology of the axisymmetric hemispherical model is output.

[0078] Specifically, according to the coordinates of each particle, the particle number density is calculated, and the particle number density is determined based on whether it is greater than , to determine whether the position information of each particle is within the calculation domain. When the particle number density is less than When the particle flows out of the boundary, it is considered that the particle flows out of the boundary. When the particle flows out of the boundary, the particle is moved out of the calculation domain and the type is changed to ghost particle. Ghost particles are not considered in all calculation processes and output processes; when the particle is still within the boundary, the position information of each particle is output, such as Figure 3 The ghost particles are represented as ghost particles, which simulate the fluid particles leaving the calculation domain after being subjected to shear force. During the calculation process, the coordinates, pressure, temperature, etc. of the particles are output at regular time intervals. The program automatically exits after the calculation time is reached. This way, you can view the calculation results at each time, and you can also make animations to observe the changes during the calculation process.

[0079] Parts of the present invention that are not described in detail belong to common knowledge among those skilled in the art.

Claims

1. A method for simulating the ablation morphology of an axisymmetric model of quartz material, characterized in that: include: (1) Use solid particles and dummy particles to construct an axisymmetric hemispherical model to form a computational domain; and collect the coordinates of each particle in the axisymmetric hemispherical model to obtain a coordinate set; (2) Based on the preset aerodynamic heating conditions, the thermal conductivity, density, specific heat capacity and initial temperature of each particle are collected; (3) Calculate the temperature of each particle based on the collected thermal conductivity, density, specific heat capacity and initial temperature of the particles; (4) calculating the velocity and coordinates of each particle in the calculation domain according to a preset shear force parameter, a preset potential based model, the temperature of each particle, the coordinate set and a preset time node; (5) According to the temperature, velocity and coordinates of each particle, the surface morphology of the axisymmetric hemispherical model is output, and the simulation of the ablation morphology of the axisymmetric model of quartz material is completed.

2. The method for simulating ablation morphology of an axisymmetric model of quartz material according to claim 1, characterized in that: The solid particles represent solid particles, simulating the state of the axisymmetric hemispherical model before heating; the dummy particles represent virtual particles, which are used to ensure that aerodynamic heating conditions are only applied on the hemispherical surface.

3. The method for simulating ablation morphology of an axisymmetric model of quartz material according to claim 2, characterized in that: When creating an axisymmetric hemispherical model, solid particles are used to arrange the shape of the axisymmetric hemispherical model, and dummy particles are arranged behind the cross section of the hemispherical model to prevent the solid particles at the cross section of the hemispherical model from being affected by aerodynamic heating during the calculation process; according to the established axisymmetric hemispherical model, the coordinates of each solid particle are collected to obtain a coordinate set.

4. The method for simulating ablation morphology of an axisymmetric model of quartz material according to claim 1, characterized in that: The thermal conductivity, density, specific heat capacity and initial temperature of each particle are collected as follows: Aerodynamic heating is applied to the heating surface of the axisymmetric hemispherical model in the computational domain, and a fixed temperature is set. During the aerodynamic heating process, the temperature of the solid particles gradually rises, and the thermal conductivity, density, specific heat capacity and initial temperature of each solid particle are collected over a period of time.

5. The method for simulating ablation morphology of an axisymmetric model of quartz material according to claim 4, characterized in that: The temperature of each particle is calculated based on the collected thermal conductivity, density, specific heat capacity and initial temperature of the particles, specifically: in, represents the temperature of the particle in the i-th axisymmetric hemispherical model at the k+1th step; represents the temperature of the i-th particle at the k-th step; the j-th particle represents all particles around the i-th particle and different from the i-th particle; represents the temperature of the jth particle at the kth step; represents the time step; represents the thermal conductivity of the ith particle; represents the density of the i-th particle; represents the specific heat capacity of the ith particle; Represents dimension; , are the position vectors of the jth and ith particles respectively; then Represents the distance between particles j and i ; is the effective radius, take ; represents the standard particle number density at the initial moment; represents the weight function; represents the position vector of particle j; represents the position vector of particle i; represents the net heat flow into the surface; represents the spatial step length, i.e., the distance between particles at the initial moment; represents the particle number density after considering virtual particles; represents the particle number density, ; For surface particles, , internal particles , the particle number density was used to achieve aerodynamic heating only on the surface; It represents the average value of the square of the distance between particles.

6. The method for simulating ablation morphology of an axisymmetric model of quartz material according to claim 5, characterized in that: The net heat flux into the particle surface Calculate using the following method: ； in, is the cold wall heat flow; is the total enthalpy of the air flow; is the wall enthalpy; is the radiation coefficient; is the Stefan Boltzmann constant; is the elicitation factor; represents the temperature of the particle; Indicates the ambient temperature, is the loss rate, is the evaporation rate, is the heat loss effect, is the evaporation heat effect, is the coefficient of the Lees distribution of the Lees spherical heat flux, is the spherical center angle; is the coefficient of the Lees distribution of the Lees spherical heat flux, is the sphere center angle, The calculation is done using the following formula: 。 in is the specific heat ratio of the gas, is the incoming flow Mach number; Emission Factor Calculate using the following method: in, is the mass flow rate of the ejected gas generated by carbon ablation and pyrolysis of the original material, is the pyrolysis gas mass flow rate; The mass flow rate of carbon monoxide gas produced by carbon oxidation; Generates gas mass flow rate for carbon sublimation; is the vapor pressure; Pe is the external gas pressure; is the ratio of the molecular weight of air to the molecular weight of SiO2; for laminar and turbulent flow, Equal to 0.62 and 0.2 respectively.

7. The method for simulating ablation morphology of an axisymmetric model of quartz material according to claim 5, characterized in that: When the particle temperature reaches the melting point, the particle type changes from solid to fluid; when the particle temperature does not reach the melting point, the particle type is solid. Setting the solid melting temperature is the temperature threshold. When the particle temperature is greater than or equal to When the particle type is changed from solid particle to fluid particle, flow calculation begins. , the particle type is solid particle, fixed and motionless.

8. The method for simulating ablation morphology of an axisymmetric model of quartz material according to claim 5, characterized in that: The calculating the velocity and coordinates of each particle in the calculation domain according to the preset shear force parameter, the preset potential based model, the temperature of each particle, the coordinate set and the preset time node specifically includes: (4.1) Based on the particle’s current velocity and coordinate set, obtain the Laplacian operator of velocity; (4.2) According to the temperature of each particle, collect the viscosity coefficient of each particle; (4.3) According to the viscosity coefficient of each particle, the Laplace operator and the preset shear force parameter, the first temporary velocity and the first temporary coordinate generated by the viscosity term are obtained; (4.4) According to the preset potential based model, the surface tension of each particle is obtained; (4.5) obtaining a second temporary velocity and a second temporary coordinate of each particle according to the surface tension, the preset time node, the preset gravity, the first temporary velocity and the first temporary coordinate; (4.6) Calculate the pressure of each particle based on the second temporary velocity and the second temporary coordinate of each particle; (4.7) Based on the pressure of each particle, the velocity and coordinates of each particle are calculated.

9. The method for simulating ablation morphology of an axisymmetric model of quartz material according to claim 8, characterized in that: The Laplace operator of the velocity is obtained based on the current velocity and coordinate set of the particle, specifically: represents the velocity vector of particle i at the kth time step; represents the velocity vector of particle j at the kth time step; represents the position vector of particle j; represents the position vector of particle i; represents the average value of the square of the distance between particles; represents the weight function; Represents the distance between particles j and i ; is the effective radius, take ; represents the standard particle number density at the initial moment; The viscosity coefficient of each particle ; A, B, C are constants and have different values ​​for different quartz materials. Represents the temperature of the particle.

10. The method for simulating ablation morphology of an axisymmetric model of quartz material according to claim 9, characterized in that: The step (4.3) obtains the first temporary velocity and the first temporary coordinate generated by the viscosity term according to the viscosity coefficient of each particle, the Laplace operator and the preset shear force parameter. The calculation method is as follows: in, is the coordinate component of the tangential shear force of the spherical head model; Represents the first temporary velocity vector of particle i; Represents the first temporary coordinate vector of particle i; represents the first temporary velocity vector of particle j; represents the first temporary coordinate vector of particle j; represents the velocity vector of particle i at the kth time step; represents the velocity vector of particle j at the kth time step; represents the coordinate vector of particle i at the kth time step; represents the time step; μ represents the dynamic viscosity coefficient; ρ represents the density of the particle; Indicates dimension, for three-dimensional problems, d=3; represents the average value of the square of the distance between particles i and j; represents the position vector of particle j; represents the position vector of particle i; Represents the distance between particles j and i ; is the effective radius, take ; represents the standard particle number density at the initial moment; represents the set of all particles around particle i excluding ghost particles; represents the weight function; represents the particle number density after considering virtual particles; represents the particle number density; represents the spatial step length, i.e., the distance between particles at the initial moment; The shear force is calculated using engineering calculation method, and the calculation method is as follows: ； in, is the velocity at the outer edge of the boundary layer represents the elicitation factor; represents the cold wall heat flux; represents the total enthalpy of air flow; represents the Prandtl number; Represents the angle between the tangent direction of the spherical head model and the incoming flow direction.

11. The method for simulating ablation morphology of an axisymmetric model of quartz material according to claim 10, characterized in that: The step (4.4) obtains the surface tension of each particle according to the preset potential based model, specifically: The acceleration due to surface tension is calculated as follows: ； Where A represents the set of i; B represents the set of particles around particle i; represents surface tension; is the Stefan Boltzmann constant; represents the distance between particles i and j; represents the effective radius of the force between particles, here we take ; Indicates the minimum distance between particles at the initial moment; ； in, represents the potential energy between particles i and j; ； in, represents the sum of the potential energy of all particles around particle i acting on it; ； in, represents the resultant force of particles around particle i on it; represents the weight function; ； in, Represents the particle acceleration obtained based on the potential based model; Represents the mass of the particle.

12. The method for simulating ablation morphology of an axisymmetric model of quartz material according to claim 11, characterized in that: The step (4.5) obtains the second temporary velocity and the second temporary coordinate of each particle according to the surface tension, the preset time node, the preset gravity, the first temporary velocity and the first temporary coordinate, specifically: Second temporary speed Calculate using the following method: ； in, is the acceleration of surface tension, is the acceleration due to gravity. Second temporary coordinates Calculate using the following method: 。 13. The method for simulating ablation morphology of an axisymmetric model of quartz material according to claim 12, characterized in that: The (4.6) calculates the pressure of each particle according to the second temporary velocity and the second temporary coordinate of each particle, specifically: According to the second temporary coordinate, the temporary value of the particle number density is calculated , the calculation method is as follows: ； The pressure is calculated using the following Poisson's equation for pressure: ； in, represents the pressure of particle j at step k+1; represents the pressure of particle i at step k+1; represents the empirical coefficient; represents the particle number density of virtual particles; Indicates the external gas pressure; represents the standard particle number density at the initial moment; Represents the gradient.

14. The method for simulating ablation morphology of an axisymmetric model of quartz material according to claim 13, characterized in that: The velocity and coordinates of each particle are calculated according to the pressure of each particle in (4.7), specifically: After obtaining the value of the pressure at the k+1 time step, the velocity and coordinates at the k+1 time step are obtained using the following method: ； ； in, represents the velocity vector of particle i at the k+1 time step; Represents the second temporary velocity vector of particle i; represents the particle density; represents the gradient of particle pressure; in, Represents the position vector of particle i at the k+1th step.

15. The method for simulating ablation morphology of an axisymmetric model of quartz material according to claim 14, characterized in that: The surface morphology of the axisymmetric hemispherical model is output according to the temperature, velocity and coordinates of each particle, specifically: Output the type of each particle according to the temperature of each particle; Output the position information of each particle according to the coordinates and speed of each particle; According to the particle type, flow state and the position information, the surface morphology of the axisymmetric hemispherical model is output.

16. A method for simulating ablation morphology of an axisymmetric model of quartz material according to any one of claim 15, characterized in that: According to the coordinates of each particle, calculate the particle number density, and determine whether the particle number density is greater than , to determine whether the position information of each particle is within the calculation domain. When the particle number density is less than When the particle flows out of the boundary, it is considered that the particle flows out of the boundary. When the particle flows out of the boundary, the particle is moved out of the calculation domain and the type is changed to a ghost particle. The ghost particle is not considered in all calculation processes and output processes. When the particle is still within the boundary, the position information of each particle is output. The ghost particle is represented as a ghost particle, which simulates the fluid particle leaving the calculation domain after being subjected to shear force.

Images