A method for simulating the ablation morphology of an axisymmetric model of quartz material

Abstract

A method for simulating the ablation morphology of an axisymmetric model of quartz material comprises: (1) constructing an axisymmetric hemispherical model using solid particles and dummy particles to form a calculation domain; and collecting the coordinates of each particle in the axisymmetric hemispherical model to obtain a coordinate set; (2) collecting the thermal conductivity, density, specific heat capacity and initial temperature of each particle based on a preset aerodynamic heating condition; (3) calculating the temperature of each particle based on the collected thermal conductivity, density, specific heat capacity and initial temperature of the particle; (4) calculating the velocity and coordinate of each particle in the calculation domain based on a preset shear force parameter, a preset potential-based model, the temperature of each particle, the coordinate set and a preset time node; (5) outputting the surface morphology of the axisymmetric hemispherical model based on the temperature, velocity and coordinate of each particle, thereby completing the simulation of the ablation morphology of the axisymmetric model of quartz material. The present invention solves the deficiencies in the prior art in the research on the three-dimensional ablation morphology of quartz material.

Description

Technical Field

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

[0002] Thermal protection has always been a core technology for hypersonic flight and the development of various aircraft. For aircraft undergoing high-Mach re-entry or prolonged atmospheric cruising, structural thermal protection and its impact on the vehicle's aerodynamic performance remain crucial. The nose cap, at the head of a blunt-bodied vehicle, experiences the harshest aerodynamic heating environment and is therefore a key area for thermal protection.

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

[0004] However, the protective layers of different parts of the aircraft are different, and the degree of ablation is also different. 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 deficiencies of the prior art and propose a method for simulating the ablation morphology of an axisymmetric model of a quartz material, thereby solving the deficiencies in the prior art in the study of the three-dimensional ablation morphology of quartz materials.

[0006] The technical solution adopted in the present invention is:

[0007] A method for simulating the ablation morphology of an axisymmetric model of a quartz material comprises the following steps:

[0008] (1) Using solid particles and dummy particles, an axisymmetric hemispherical model is constructed to form a computational domain; and the coordinates of each particle in the axisymmetric hemispherical model are collected to obtain a coordinate set;

[0009] (2) Based on the preset aerodynamic heating conditions, the thermal conductivity, density, specific heat capacity and initial temperature of each particle are collected;

[0010] (3) Calculate the temperature of each particle based on the collected thermal conductivity, density, specific heat capacity and initial temperature of the particles;

[0011] (4) calculating the velocity and coordinates of each particle in the computational 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;

[0012] (5) According to the temperature, velocity and coordinates of each particle, the surface morphology of the axisymmetric hemispherical model is output to complete the simulation of the ablation morphology of the axisymmetric model of quartz material.

[0013] 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 applied only on the hemispherical surface.

[0014] 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.

[0015] Furthermore, the thermal conductivity, density, specific heat capacity and initial temperature of each particle are collected as follows:

[0016] Aerodynamic heating is applied to the heating surface of an axisymmetric hemispherical model within the computational domain, and a fixed temperature is set. As the aerodynamic heating process progresses, the temperature of the solid particles gradually rises. The thermal conductivity, density, specific heat capacity, and initial temperature of each solid particle are collected over a period of time.

[0017] 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:

[0018]

[0019] 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; Δt represents the time step; k(i) represents the thermal conductivity of the i-th particle; ρ(i) represents the density of the i-th particle; c p (i) represents the specific heat capacity of the i-th particle; d represents the dimension; are the position vectors of the jth and ith particles respectively; then represents the distance r between particles j and i; r e is the effective radius, take r e=2.1l0;n 0 represents the standard particle number density at the initial moment; w represents the weight function; represents the position vector of particle j; represents the position vector of particle i; q n (i) represents the net heat flux entering the surface; l0 represents the spatial step size, that is, the distance between particles at the initial moment;

[0020] <n′> i represents the particle number density after considering virtual particles;

[0021] <n> i represents the particle number density,

[0022] <n′> i =max( <n. i ,n 0 ), for surface particles,<n′> i > <n> i , internal particles<n′> i = <n> i , the particle number density is used to achieve aerodynamic heating only on the surface;

[0023] represents the average value of the square of the distance between particles.

[0024] Furthermore, the net heat flux q entering the particle surface n (i) Calculate using the following method:

[0025]

[0026] Among them, q or is the cold wall heat flux; h r is the total enthalpy of the air flow; h w is the wall enthalpy; ε is the radiation coefficient; σ is the Stefan-Boltzmann constant; ψ is the ejection factor; T w represents the temperature of the particle; T envir Indicates the ambient temperature, v -δ is the loss velocity, v w is the evaporation rate, ΔH L is the heat loss effect, is the evaporation heat effect, is the coefficient of Lees distribution of heat flux on the Lees spherical surface, and θ is the spherical center angle;

[0027] is the coefficient of Lees distribution of heat flux on the Lees spherical surface, θ is the spherical center angle, The calculation is done using the following formula:

[0028]

[0029] Where γ represents the specific heat ratio of the gas, M ∞ is the incoming flow Mach number;

[0030] The ejection factor ψ is calculated as follows:

[0031]

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

[0033] 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.

[0034] Set the solid melting temperature T melt is the temperature threshold. When the particle temperature is greater than or equal to T melt When the particle type is changed from solid particle to fluid particle, the flow calculation begins. When the temperature is less than T melt , the particle type is solid particle, fixed and motionless.

[0035] 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:

[0036] (4.1) Based on the particle’s current velocity and coordinate set, obtain the Laplace operator of the velocity;

[0037] (4.2) Collect the viscosity coefficient of each particle based on the temperature of each particle;

[0038] (4.3) According to the viscosity coefficient, Laplace operator and preset shear force parameters of each particle, the first temporary velocity and first temporary coordinate generated by the viscosity term are obtained;

[0039] (4.4) Obtain the surface tension of each particle based on the preset potential-based model;

[0040] (4.5) obtaining a second temporary velocity and a second temporary coordinate of each particle based on the surface tension, the preset time node, the preset gravity, the first temporary velocity, and the first temporary coordinate;

[0041] (4.6) Calculate the pressure of each particle based on its second temporary velocity and second temporary coordinate;

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

[0043] Further,

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

[0045]

[0046] represents the velocity vector of particle i at the kth time step;

[0047] represents the velocity vector of particle j at the kth time step;

[0048] represents the position vector of particle j;

[0049] represents the position vector of particle i;

[0050] λ represents the average value of the square of the distance between particles;

[0051] w represents the weight function;

[0052] represents the distance r between particles j and i; r e is the effective radius, take r e =2.1l0;

[0053] n 0 represents the standard particle number density at the initial moment;

[0054] The viscosity coefficient μ of each particle i ;

[0055]

[0056] A, B, C are constants and have different values ​​for different quartz materials. w (i) represents the temperature of the particle.

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

[0058]

[0059] in, is the coordinate component of the tangential shear force of the spherical head model;

[0060] represents the first temporary velocity vector of particle i;

[0061] Represents the first temporary coordinate vector of particle i;

[0062] represents the first temporary velocity vector of particle j;

[0063] represents the first temporary coordinate vector of particle j;

[0064] represents the velocity vector of particle i at the kth time step;

[0065] represents the velocity vector of particle j at the kth time step;

[0066] represents the coordinate vector of particle i at the kth time step;

[0067] Δt represents the time step;

[0068] μ represents the dynamic viscosity coefficient;

[0069] ρ represents the density of the particles;

[0070] d represents the dimension, for three-dimensional problems, d=3;

[0071] represents the average value of the square of the distance between particles i and j;

[0072] represents the position vector of particle j;

[0073] represents the position vector of particle i;

[0074] represents the distance r between particles j and i; r e is the effective radius, take r e =2.1l0;

[0075] n 0 represents the standard particle number density at the initial moment;

[0076] Λ i represents the set of all particles around particle i excluding ghost particles;

[0077] w represents the weight function;

[0078] <n′> i represents the particle number density after considering virtual particles;

[0079] <n> i represents the particle number density;

[0080] l0 represents the spatial step length, that is, the distance between particles at the initial moment;

[0081] The shear force is calculated using engineering calculation methods, as shown below:

[0082]

[0083] in, is the velocity at the outer edge of the boundary layer

[0084] ψ represents the elicitation factor;

[0085] q or represents the cold wall heat flux; h r represents the total enthalpy of the air flow;

[0086] p r represents the Prandtl number;

[0087] θ represents the angle between the tangent direction of the spherical head model and the incoming flow direction.

[0088] Furthermore, the surface tension of each particle is obtained in step (4.4) according to a preset potential-based model, specifically:

[0089] The acceleration due to surface tension is calculated as follows:

[0090]

[0091] Where A represents the set of i;

[0092] B represents the set of particles surrounding particle i;

[0093] C represents surface tension;

[0094] σ is the Stefan-Boltzmann constant;

[0095] r ij represents the distance between particles i and j;

[0096] r e represents the effective radius of the force between particles, which is taken as 3.1l0 here;

[0097] r min Indicates the minimum distance between particles at the initial moment;

[0098]

[0099] Among them, P(r ij ) represents the potential energy between particles i and j;

[0100] E st =∑ ij P(r ij );

[0101] Among them, E st represents the sum of the potential energies of all particles around particle i;

[0102]

[0103] in, represents the resultant force of particles around particle i on it;

[0104] w ij represents the weight function;

[0105]

[0106] in, represents the particle acceleration obtained based on the potential based model;

[0107] Indicates the mass of the particle.

[0108] 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:

[0109] Second temporary speed Calculate using the following method:

[0110]

[0111] in, is the acceleration of surface tension, is the acceleration due to gravity.

[0112] Second temporary coordinates Calculate using the following method:

[0113]

[0114] Furthermore, the pressure of each particle is calculated according to the second temporary velocity and the second temporary coordinate of each particle in (4.6), specifically:

[0115] According to the second temporary coordinate, the temporary value of the particle number density is calculated The calculation method is as follows:

[0116]

[0117] The pressure is calculated using the following Poisson equation for pressure:

[0118]

[0119] in, represents the pressure of particle j at step k+1;

[0120] represents the pressure of particle i at step k+1;

[0121] γ represents the empirical coefficient;

[0122] represents the particle number density considering virtual particles;

[0123] P free Indicates the external gas pressure;

[0124] n 0 represents the standard particle number density at the initial moment;

[0125] Represents the gradient.

[0126] Furthermore, (4.7) calculates the velocity and coordinates of each particle based on the pressure of each particle, specifically:

[0127] After obtaining the pressure value of the k+1 time step, the velocity and coordinates of the k+1 time step are obtained using the following method:

[0128]

[0129]

[0130] in, represents the velocity vector of particle i at the k+1 time step;

[0131] Represents the second temporary velocity vector of particle i;

[0132] ρ 0 represents the particle density;

[0133] represents the gradient of particle pressure;

[0134] in, Represents the position vector of particle i at the k+1th step.

[0135] Furthermore, the surface morphology of the axisymmetric hemispherical model is output according to the temperature, velocity and coordinates of each particle, specifically:

[0136] Output the type of each particle according to the temperature of each particle;

[0137] Output the position information of each particle according to its coordinates and velocity;

[0138] According to the particle type, flow state and the position information, the surface morphology of the axisymmetric hemispherical model is output.

[0139] 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 0.1n 0 , to determine whether the position information of each particle is within the calculation domain. When the particle number density is less than 0.1n 0 When the particle flows out of the boundary, it is considered that the particle flows out of the boundary. After 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; the ghost particles are represented as ghost particles, which simulate the fluid particles leaving the calculation domain after being subjected to shear force.

[0140] The beneficial effects of the present invention are:

[0141] (1) The present invention creates axisymmetric spherical head models of different particles within a preset computational domain, applies ablation to the particles over a certain period of time, and causes some particles to escape from the computational domain, thereby simulating the ablation state of the protective layer during flight and determining the particle evolution process within the computational domain. Furthermore, since some particles escaping the computational domain do not require computation, it is not necessary to perform computation on all particles, thus avoiding large-scale computations and improving computational efficiency and accuracy.

[0142] (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 scenario, thereby improving the prediction accuracy of the material ablation behavior. BRIEF DESCRIPTION OF THE DRAWINGS

[0143] Figure 1 Flow chart of the method of the present invention;

[0144] Figure 2 Schematic diagram of the axisymmetric particle model before ablation;

[0145] Figure 3 Schematic diagram of the axisymmetric particle model after ablation. DETAILED DESCRIPTION

[0146] When an aircraft is in flight, the thermal protection layer at different locations on the aircraft surface is affected by factors such as temperature, pressure, and speed, resulting in different degrees of ablation at different locations on the quartz composite thermal protection layer. The one-dimensional ablation prediction of traditional algorithms or the calculation of the ablation setback at the stagnation point position of the spherical head model cannot calculate the surface tension, overload, asymmetric ablation, etc. of the molten liquid layer. Existing technologies make it difficult to simulate the ablation setback at different locations in the molten state based on calculations, which affects the prediction accuracy of the quartz composite thermal protection layer.

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

[0148] In the embodiments of the present application, the fully implicit particle method MPT is obtained by improving 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 the two-dimensional or three-dimensional ablation morphology of the surface, as well as the ablation retreat at different positions.

[0149] Specifically, within the computational domain, different types of particles are used to create a particle model, and a potential-based model is used to impose an interaction force between particles so that the distance between particles can be maintained 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. Figure 2 shown.

[0150] Specifically, using aerodynamic 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. The particle type changes with the change in 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 the velocity and coordinate of each particle is calculated to obtain the first temporary velocity and first temporary coordinate. According to the potential-based model, the surface tension of each particle is obtained. The second temporary velocity and temporary coordinate of each particle are calculated using surface tension and gravity. 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, second temporary coordinate and pressure.

[0151] 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.

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

[0153] like Figure 1 As shown, the present invention proposes an ablation simulation method for an axisymmetric model of a quartz composite material under aerodynamic heating conditions, and the specific steps are as follows:

[0154] Step 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.

[0155] 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.

[0156] In an embodiment of the present application, an axisymmetric hemispherical head model is created in a computational domain by using different types of particles, including 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 heating, 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 as virtual particles to avoid collisions between these particles.

[0157] Step 2: Based on the preset aerodynamic heating conditions, the thermal conductivity, density, specific heat capacity and initial temperature of each particle are collected;

[0158] Aerodynamic heating is applied to the heating surface of an axisymmetric hemispherical model within the computational domain, and a fixed temperature is set. As the aerodynamic heating process progresses, the temperature of the solid particles gradually rises. The thermal conductivity, density, specific heat capacity, and initial temperature of each solid particle are collected over a period of time.

[0159] For example, aerodynamic heating can be used to simulate the effect of temperature on the protective layer during flight. Aerodynamic heating conditions include the heat flow, airflow, and pressure to which the aircraft is subjected within a certain period of time. When an aircraft is flying at high speed, it has a large relative velocity with the airflow. When gas molecules collide with the aircraft, their speed will stagnate, and the kinetic energy of the gas will be converted into thermal energy, generating temperatures 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.

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

[0161] Calculate the temperature of each particle, specifically:

[0162]

[0163] 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; Δt represents the time step; k(i) represents the thermal conductivity of the i-th particle; ρ(i) represents the density of the i-th particle; c p (i) represents the specific heat capacity of the i-th particle; d represents the dimension; are the position vectors of the jth and ith particles respectively; then represents the distance r between particles j and i; r e is the effective radius, take r e =2.1l0;n 0 represents the standard particle number density at the initial moment; w represents the weight function; represents the position vector of particle j; represents the position vector of particle i; q n (i) represents the net heat flux entering the surface; l0 represents the spatial step size, that is, the distance between particles at the initial moment;

[0164] <n′> i represents the particle number density after considering virtual particles;

[0165] <n> i represents the particle number density,

[0166] <n′> i =max( <n> i ,n 0 ), for surface particles,<n′> i > <n> i , internal particles<n′> i = <n> i , the particle number density is used to achieve aerodynamic heating only on the surface;

[0167] represents the average value of the square of the distance between particles.

[0168] The net heat flux q entering the particle surface n (i) Calculate using the following method:

[0169]

[0170] Among them, q or is the cold wall heat flux; h r is the total enthalpy of the air flow; h w is the wall enthalpy; ε is the radiation coefficient; σ is the Stefan-Boltzmann constant; ψ is the ejection factor; T w represents the temperature of the particle; T envir Indicates the ambient temperature, v -δ is the loss velocity, v w is the evaporation rate, ΔH L is the heat loss effect, is the evaporation heat effect, is the coefficient of Lees distribution of heat flux on the Lees spherical surface, and θ is the spherical center angle;

[0171] is the coefficient of Lees distribution of heat flux on the Lees spherical surface, θ is the spherical center angle, The calculation is done using the following formula:

[0172]

[0173] Where γ represents the specific heat ratio of the gas, M ∞ is the incoming flow Mach number;

[0174] The ejection factor ψ is calculated as follows:

[0175]

[0176] 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; P is the gas mass flow rate generated by carbon sublimation; v is the vapor pressure; Pe is the external gas pressure; M av is the ratio of the molecular weight of air to the molecular weight of SiO2; for laminar and turbulent flow, S f Equal to 0.62 and 0.2 respectively.

[0177] 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.

[0178] Set the solid melting temperature T melt is the temperature threshold. When the particle temperature is greater than or equal to T melt When the particle type is changed from solid particle to fluid particle, the flow calculation begins. When the temperature is less than T melt , the particle type is solid particle, fixed and motionless.

[0179] Step 4: Calculate 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;

[0180] 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.

[0181] It should be noted that in practical scenarios, the velocity and coordinates of each particle in the computational domain need to be output every cycle. Therefore, for each cycle, the current velocity in this application is the velocity of the particle calculated in the previous cycle. For the first cycle, the current velocity is a preset initial value, typically 0.

[0182] 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, drawing on the idea of ​​virtual particles in the MPS method, a shear force is applied to the upper surface of the axisymmetric model, and the first temporary velocity and first temporary coordinate of each particle are calculated based on the viscosity term.

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

[0184] Specifically, the velocity and coordinates of each particle in the calculation domain are calculated 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, including the following steps:

[0185] (4.1) Based on the particle’s current velocity and coordinate set, obtain the Laplace operator of the velocity;

[0186]

[0187] represents the velocity vector of particle i at the kth time step;

[0188] represents the velocity vector of particle j at the kth time step;

[0189] represents the position vector of particle j;

[0190] represents the position vector of particle i;

[0191] λ represents the average value of the square of the distance between particles;

[0192] w represents the weight function;

[0193] represents the distance r between particles j and i; r e is the effective radius, take r e =2.1l0;

[0194] n 0 represents the standard particle number density at the initial moment;

[0195] (4.2) According to the temperature of each particle, collect the viscosity coefficient μ of each particle i ;

[0196] Drawing on the idea of ​​virtual particles in the MPS method, a method is proposed to impose a net heat flow boundary on the surface that takes into account the influence of surface wall enthalpy correction, radiation heat dissipation and gas induced factors. Heat transfer, phase change and temperature-varying viscosity coefficient 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 using the temperature-varying viscosity coefficient expression of silicon-based materials.

[0197]

[0198] A, B, C are constants and have different values ​​for different quartz materials. w (i) represents the temperature of the particle.

[0199] (4.3) According to the viscosity coefficient, Laplace operator and preset shear force parameters of each particle, the first temporary velocity and first temporary coordinate generated by the viscosity term are obtained;

[0200] The changes in particle velocity and position caused by the viscosity term and shear force are calculated. Drawing on the idea of ​​virtual particles in the MPS method, 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.

[0201]

[0202] in, is the coordinate component of the tangential shear force of the spherical head model;

[0203] Represents the first temporary velocity vector of particle i;

[0204] Represents the first temporary coordinate vector of particle i;

[0205] represents the first temporary velocity vector of particle j;

[0206] represents the first temporary coordinate vector of particle j;

[0207] represents the velocity vector of particle i at the kth time step;

[0208] represents the velocity vector of particle j at the kth time step;

[0209] represents the coordinate vector of particle i at the kth time step;

[0210] Δt represents the time step;

[0211] μ represents the dynamic viscosity coefficient;

[0212] ρ represents the density of the particles;

[0213] d represents the dimension, for three-dimensional problems, d=3;

[0214] represents the average value of the square of the distance between particles i and j;

[0215] represents the position vector of particle j;

[0216] represents the position vector of particle i;

[0217] represents the distance r between particles j and i; r e is the effective radius, take r e =2.1l0;

[0218] n 0 represents the standard particle number density at the initial moment;

[0219] Λ i represents the set of all particles around particle i excluding ghost particles;

[0220] w represents the weight function;

[0221] Mn′> i represents the particle number density after considering virtual particles;

[0222] <n> i represents the particle number density;

[0223] l0 represents the spatial step length, that is, the distance between particles at the initial moment;

[0224] The shear force is calculated using engineering calculation methods, as shown below:

[0225]

[0226] in, is the velocity at the outer edge of the boundary layer

[0227] ψ represents the elicitation factor;

[0228] q or represents the cold wall heat flux; h r represents the total enthalpy of the air flow;

[0229] p r represents the Prandtl number;

[0230] θ represents the angle between the tangent direction of the spherical head model and the incoming flow direction.

[0231] (4.4) Obtain the surface tension of each particle based on the preset potential-based model;

[0232] 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.

[0233] Surface tension, also known as contraction force, refers to the fact that the forces acting on molecules at an interface differ from those acting on molecules within the bulk of the liquid. While the net force acting on a liquid molecule within a liquid is zero due to the forces acting on it from the surrounding liquid molecules, this is not the case for liquid molecules on the surface. Because the attraction of the gaseous molecules above it is less than that of the liquid molecules within, the net force acting on the molecule is non-zero and points perpendicularly to the interior of the liquid, resulting in the liquid surface having a tendency to shrink. There are two main ways to describe surface tension: for grid methods, the CSF surface tension model is used. For the particle method used in this application, which is a gridless method, a potential-based model can be used.

[0234] 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.

[0235] The acceleration due to surface tension is calculated as follows:

[0236]

[0237] Where A represents the set of i;

[0238] B represents the set of particles surrounding particle i;

[0239] C represents surface tension;

[0240] σ is the Stefan-Boltzmann constant;

[0241] r ij represents the distance between particles i and j;

[0242] r e represents the effective radius of the force between particles, which is taken as 3.1l0 here;

[0243] r min Indicates the minimum distance between particles at the initial moment;

[0244]

[0245] Among them, P(r ij ) represents the potential energy between particles i and j;

[0246] E st =∑ ij P(r ij );

[0247] Among them, E st represents the sum of the potential energies of all particles around particle i;

[0248]

[0249] in, represents the resultant force of particles around particle i on it;

[0250] w ij represents the weight function;

[0251]

[0252] in, represents the particle acceleration obtained based on the potential based model;

[0253] Indicates the mass of the particle.

[0254] (4.5) obtaining a second temporary velocity and a second temporary coordinate of each particle based on the surface tension, the preset time node, the preset gravity, the first temporary velocity, and the first temporary coordinate;

[0255] Second temporary speed Calculate using the following method:

[0256]

[0257] in, is the acceleration of surface tension, is the acceleration due to gravity.

[0258] Second temporary coordinates Calculate using the following method:

[0259]

[0260] (4.6) Calculate the pressure of each particle based on its second temporary velocity and second temporary coordinate; calculate the temporary value of the particle number density based on the second temporary coordinate. The calculation method is as follows:

[0261]

[0262] The pressure is calculated using the following Poisson equation for pressure:

[0263]

[0264] in, represents the pressure of particle j at step k+1;

[0265] represents the pressure of particle i at step k+1;

[0266] γ represents the empirical coefficient;

[0267] represents the particle number density considering virtual particles;

[0268] P free Indicates the external gas pressure;

[0269] n 0 represents the standard particle number density at the initial moment;

[0270] Represents the gradient.

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

[0272] Specifically:

[0273] After obtaining the pressure value of the k+1 time step, the velocity and coordinates of the k+1 time step are obtained using the following method:

[0274]

[0275] in, represents the velocity vector of particle i at the k+1 time step;

[0276] Represents the second temporary velocity vector of particle i;

[0277] ρ 0 represents the particle density;

[0278] represents the gradient of particle pressure;

[0279] in, Represents the position vector of particle i at the k+1th step.

[0280] 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.

[0281] include:

[0282] Output the type of each particle according to the temperature of each particle;

[0283] Output the position information of each particle according to its coordinates and velocity;

[0284] According to the particle type, flow state and the position information, the surface morphology of the axisymmetric hemispherical model is output.

[0285] 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 0.1n 0 , to determine whether the position information of each particle is within the calculation domain. When the particle number density is less than 0.1n 0 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.

[0286] 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 allows you to view the calculation results at each time, and you can also create animations to observe the changes during the calculation process.

[0287] Parts of the present invention that are not described in detail belong to common knowledge among those skilled in the art.< / n> < / n> < / n> < / n> < / n> < / n> < / n> < / n> < / n>

Claims

1. A method for simulating the ablation morphology of an axisymmetric model of quartz material, characterized in that: include: (1) Using solid particles and dummy particles, an axisymmetric hemispherical model is constructed to form a computational domain; and the coordinates of each particle in the axisymmetric hemispherical model are collected 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 computational 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) outputting the surface morphology of the axisymmetric hemispherical model according to the temperature, velocity, and coordinates of each particle, and completing the simulation of the ablation morphology of the axisymmetric model of the quartz material; The solid particles represent solid particles, simulating the state of the axisymmetric hemispherical model before heating; the dummy particles represent virtual particles, used to ensure that aerodynamic heating conditions are applied only on the hemispherical surface; 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.

2. 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 an axisymmetric hemispherical model within the computational domain, and a fixed temperature is set. As the aerodynamic heating process progresses, the temperature of the solid particles gradually rises. The thermal conductivity, density, specific heat capacity, and initial temperature of each solid particle are collected over a period of time.

3. The method for simulating ablation morphology of an axisymmetric model of quartz material according to claim 2, 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; Δt represents the time step; k(i) represents the thermal conductivity of the i-th particle; ρ(i) represents the density of the i-th particle; c p (i) represents the specific heat capacity of the i-th particle; d represents the dimension; are the position vectors of the jth and ith particles respectively; then represents the distance r between particles j and i; r e is the effective radius, take r e =2.1l0;n 0 represents the standard particle number density at the initial moment; w represents the weight function; represents the position vector of particle j; represents the position vector of particle i; q n (i) represents the net heat flux entering the surface; l0 represents the spatial step size, that is, the distance between particles at the initial moment; <n′> i represents the particle number density after considering virtual particles; <n> i represents the particle number density, < / n> <n′> i =max( <n> i ,n 0 ), for surface particles,<n′> i > <n> i , internal particles<n′> i = <n> i , the particle number density is used to achieve aerodynamic heating only on the surface;< / n> < / n> < / n> represents the average value of the square of the distance between particles.

4. The method for simulating ablation morphology of an axisymmetric model of quartz material according to claim 3, characterized in that: The net heat flux q entering the particle surface n (i) Calculate using the following method: Among them, q or is the cold wall heat flux; h r is the total enthalpy of the air flow; h w is the wall enthalpy; ε is the radiation coefficient; σ is the Stefan-Boltzmann constant; ψ is the ejection factor; T w represents the temperature of the particle; T envir Indicates the ambient temperature, v -δ is the loss velocity, v w is the evaporation rate, ΔH L is the heat loss effect, is the evaporation heat effect, is the coefficient of Lees distribution of heat flux on the Lees spherical surface, and θ is the spherical center angle; is the coefficient of Lees distribution of heat flux on the Lees spherical surface, θ is the spherical center angle, The calculation is done using the following formula: Where γ represents the specific heat ratio of the gas, M ∞ is the incoming flow Mach number; The ejection factor ψ is calculated as follows: in, is the mass flow rate of the ejected gas generated by carbon ablation and pyrolysis of the original material, is the mass flow rate of pyrolysis gas; The mass flow rate of carbon monoxide gas produced by carbon oxidation; P is the gas mass flow rate generated by carbon sublimation; v is the vapor pressure; Pe is the external gas pressure; M av is the ratio of the molecular weight of air to the molecular weight of SiO2; for laminar and turbulent flow, S f Equal to 0.62 and 0.2 respectively.

5. The method for simulating ablation morphology of an axisymmetric model of quartz material according to claim 3, 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; Set the solid melting temperature T melt is the temperature threshold. When the particle temperature is greater than or equal to T melt When the particle type is changed from solid particle to fluid particle, the flow calculation begins. When the temperature is less than T melt , the particle type is solid particle, fixed and motionless.

6. The method for simulating ablation morphology of an axisymmetric model of quartz material according to claim 3, characterized in that: 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 Laplace operator of the velocity; (4.2) Collect the viscosity coefficient of each particle based on the temperature of each particle; (4.3) According to the viscosity coefficient, Laplace operator and preset shear force parameters of each particle, the first temporary velocity and first temporary coordinate generated by the viscosity term are obtained; (4.4) Obtain the surface tension of each particle based on the preset potential-based model; (4.5) obtaining a second temporary velocity and a second temporary coordinate of each particle based on 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 its second temporary velocity and second temporary coordinate; (4.7) Based on the pressure of each particle, the velocity and coordinates of each particle are calculated.

7. The method for simulating ablation morphology of an axisymmetric model of quartz material according to claim 6, 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; w represents the weight function; represents the distance r between particles j and i; r e is the effective radius, take r e =2.1l0; n 0 represents the standard particle number density at the initial moment; The viscosity coefficient μ of each particle i ; A, B, C are constants and have different values ​​for different quartz materials. w (i) represents the temperature of the particle.

8. The method for simulating ablation morphology of an axisymmetric model of quartz material according to claim 7, characterized in that: In step (4.3), the first temporary velocity and first temporary coordinate generated by the viscosity term are obtained based on the viscosity coefficient, Laplace operator and preset shear force parameters of each particle. 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; Δt represents the time step; μ represents the dynamic viscosity coefficient; ρ represents the density of the particles; d represents the 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 r between particles j and i; r e is the effective radius, take r e =2.1l0; n 0 represents the standard particle number density at the initial moment; Λ i represents the set of all particles around particle i excluding ghost particles; w represents the weight function; <n′> i represents the particle number density after considering virtual particles; <n> i represents the particle number density;< / n> l0 represents the spatial step length, that is, the distance between particles at the initial moment; The shear force is calculated using engineering calculation methods, as shown below: in, is the velocity at the outer edge of the boundary layer ψ represents the elicitation factor; q or represents the cold wall heat flux; h r represents the total enthalpy of the air flow; p r represents the Prandtl number; θ represents the angle between the tangent direction of the spherical head model and the incoming flow direction.

9. The method for simulating ablation morphology of an axisymmetric model of quartz material according to claim 8, 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 surrounding particle i; C represents surface tension; σ is the Stefan-Boltzmann constant; r ij represents the distance between particles i and j; r e represents the effective radius of the force between particles, which is taken as 3.1l0 here; r min Indicates the minimum distance between particles at the initial moment; Among them, P(r ij ) represents the potential energy between particles i and j; E st =∑ ij P(r ij ); Among them, E st represents the sum of the potential energies of all particles around particle i; in, represents the resultant force of particles around particle i on it; w ij represents the weight function; in, represents the particle acceleration obtained based on the potential based model; Indicates the mass 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.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:

11. The method for simulating ablation morphology of an axisymmetric model of quartz material according to claim 10, characterized in that: The pressure of each particle is calculated according to the second temporary velocity and the second temporary coordinate of each particle in (4.6), 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 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 considering virtual particles; P free Indicates the external gas pressure; n 0 represents the standard particle number density at the initial moment; Represents the gradient.

12. The method for simulating ablation morphology of an axisymmetric model of quartz material according to claim 11, characterized in that: The velocity and coordinates of each particle are calculated based on the pressure of each particle in (4.7), specifically: After obtaining the pressure value of the k+1 time step, the velocity and coordinates of 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; ρ 0 represents the particle density; represents the gradient of particle pressure; in, Represents the position vector of particle i at the k+1th step.

13. The method for simulating ablation morphology of an axisymmetric model of quartz material according to claim 12, 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 its coordinates and velocity; According to the particle type, flow state and the position information, the surface morphology of the axisymmetric hemispherical model is output.

14. A method for simulating ablation morphology of an axisymmetric model of quartz material according to any one of claim 13, 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 0.1n 0 , to determine whether the position information of each particle is within the calculation domain. When the particle number density is less than 0.1n 0 When the particle flows out of the boundary, it is considered that the particle flows out of the boundary. After 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; the ghost particles are represented as ghost particles, which simulate the fluid particles leaving the calculation domain after being subjected to shear force.

Images