阅读说明：本技术 一种固体火箭发动机燃面退移计算方法 (Solid rocket engine combustion surface retreating calculation method ) 是由 魏然 石德磊 李露 孙林 鲍福廷 刘旸 惠卫华 于 2021-07-27 设计创作，主要内容包括：本发明提供了一种固体火箭发动机燃面退移计算方法,包括：构建三角面片集合；构建最短距离场函数构建笛卡尔网格上的燃速分布场r(x)；修改场函数在每个节点坐标x处的数值；基于新的场函数生成用三角面片表达的等值面,将该等值面中的所有三角面片的集合记为I；进行几何布尔运算；计算固体火箭发动机在当前时刻的总燃气生成率；计算固体火箭发动机燃烧室内当前时刻的压强p；燃烧表面面积变化曲线、燃气生成率曲线,以及压强p随时间变化的曲线,即为燃面退移计算结果。该方法整个计算流程不涉及数值差分运算,节省了算力,消除了在棱线、角点附近区域进行数值差分运算所可能引入的计算误差。(The invention provides a method for calculating the combustion surface retreating of a solid rocket engine, which comprises the following steps: constructing a triangular patch set; constructing shortest distance field functions Constructing a burning rate distribution field r (x) on a Cartesian grid; modifying field functions The value at each node coordinate x; based on new field functions Generating representations by triangular patches Iso-surface, all triangles in the iso-surfaceThe collection of the patches is marked as I; performing geometric Boolean operation; calculating the total gas generation rate of the solid rocket engine at the current moment; calculating the pressure p at the current moment in the combustion chamber of the solid rocket engine; the combustion surface area change curve, the gas generation rate curve and the pressure p curve changing along with time are combustion surface retreating calculation results. The whole calculation process of the method does not relate to numerical difference operation, so that the calculation force is saved, and the calculation errors possibly introduced by the numerical difference operation in the areas near the edge lines and the corner points are eliminated.)

1. A method for calculating the combustion surface recession of a solid rocket engine is characterized by comprising the following steps:

step 1, dispersing a propellant charge geometric model of a solid rocket engine into triangular patch models, and classifying all triangular patches into two triangular patch sets; the set of triangular patches starting to burn at the moment of ignition is marked as B, and all triangular patches in B form a thin-sheet-shaped curved surface; the set of all triangular patches in the model is marked as G, and all triangular patches in G form a closed geometric body;

step 2, dividing a Cartesian grid into three-dimensional envelope boxes where the input solid rocket engine is located, and recording the node coordinate of each node of the grid as x;

step 3, for the triangular patch set B, on each node of the Cartesian grid, enabling the triangular patch set B to be in a triangular shapeCalculating the shortest distance from the node coordinate x to any triangular patch in the triangular patch set B by using an AABB quick distance field algorithm to form a shortest distance field function defined on all nodes of a Cartesian grid

Step 4, calculating the burning rate r of the propellant at each node coordinate x based on a solid propellant burning rate model to form a burning rate distribution field r (x) on a Cartesian grid;

wherein r (x) ap ^ n, a and n are respectively the burning rate coefficient and the pressure index of the propellant at the node coordinate x and are obtained by physical parameters of the propellant; p is the pressure in the combustion chamber of the solid rocket engine at the current moment;

step 5, on each node, performing the calculationModifying field functionsThe value at each node coordinate x;

step 6, using Marching cubes isosurface extraction algorithm based on new field functionGenerating representations by triangular patchesThe isosurface records the set of all triangular patches in the isosurface as I;

step 7, performing geometric Boolean operation, cutting off the part of the isosurface I outside the closed geometric body G, and updating the I set to only contain the triangular patch in the cut isosurface; at the moment, the area of the sheared I is the area of the combustion surface of the propellant at the current moment;

step 8, multiplying the area of each triangular surface patch t forming the I by the combustion speed at the position, and then multiplying by the density of a propellant, namely the gas generation rate of the position of the current time t, summing the gas generation rate data corresponding to all the triangular surface patches to obtain the total gas generation rate of the solid rocket engine at the current time, and recording the combustion surface area and all the gas generation rate data;

step 9, calculating the pressure p at the current moment in the combustion chamber of the solid rocket engine based on the numerical model of the inner ballistic curve according to the combustion surface area and the gas generation rate data;

step 10, regarding the sheared I as the surface of combustion at the beginning, replacing the set B with the set I, calculating the shortest distance from the node to the set B of the triangular patches by using an AABB quick distance field algorithm on each node of the Cartesian grid to form a new shortest distance field function defined on all nodes of the Cartesian gridTraverse all mesh nodes, check for correspondencesWhether less than zero; if so, the method will be usedIs modified toFinally, using the modifiedReplace original

Step 11, based on the new one obtained in step 10And the new one obtained in step 9p, circularly executing the steps 4-10; until I does not contain any triangular patch; at the moment, the propellant grains are completely burnt out;

the combustion surface area change curve, the gas generation rate curve and the pressure p curve changing along with time are combustion surface retreating calculation results.

2. The method for calculating the combustion surface recession of the solid rocket engine according to claim 1, wherein in step 2, a three-dimensional envelope box where the input solid rocket engine is located can be further divided into a conformal structure or a non-structural grid, and in this case, in step 6, other isosurface extraction algorithms of a grid type are selected.

3. The method for calculating the combustion surface recession of a solid-rocket engine according to claim 1, wherein the step 7 is specifically:

traversing each triangle t in the I, and judging whether the t is positioned inside the G or not;

if t is inside G, no operation is performed;

if t is outside the closed geometry represented by G, then t is deleted from I;

if t intersects with triangles in G, t is divided into a part inside G and a part outside G along the intersecting line, the part inside G is divided into triangles, and the triangles are added into I, and t is deleted from I.

Technical Field

The invention belongs to the field of simulation calculation of a solid rocket power system, and particularly relates to a method for calculating the combustion surface recession of a solid rocket engine.

Background

The combustion surface retreating of the solid rocket engine refers to the process that during the operation of the solid rocket engine, the part of the solid propellant charge exposed to high-temperature fuel gas is combusted layer by layer and is converted from solid to gas. In order to calculate the thrust-time curve of the solid rocket engine, the time-varying process of the combustion area and the gas generation rate of the solid propellant charge during the combustion surface recession process must be calculated, and the calculation task is called the combustion surface recession calculation of the solid rocket engine.

The shape of the solid propellant grains is quite complex and only part of the surface is exposed to the combustion gases at the start of combustion; during the layer-by-layer combustion of the solid propellant charge, there are also differences in the combustion speed (continuous or abrupt changes) throughout the space, which factors ultimately make the process extremely complicated in its shape change and accompanied by changes in the topology. Common calculation methods for the combustion surface retreating of the solid rocket engine include an analytical analysis method, a CAD calculation method, a moving grid method, a shortest distance field method, a level set method and the like.

Wherein, analytical analysis method: manual analysis is carried out manually to carry out combustion surface retreating and moving expansion calculation, and for solid propellant grains with high geometrical shape complexity and combustion speed distribution, correct results are difficult to obtain through manual analysis.

CAD calculation method: the combustion surface retreating and expanding calculation is carried out by depending on the offset calculation function provided by commercial CAD software, and when the topological structure of the solid propellant grain changes along with the layer-by-layer combustion, calculation errors or failures occur at a probability, and the condition of non-uniform distribution of the combustion speed cannot be calculated.

The dynamic grid method comprises the following steps: combustion surface regression and expansion are carried out by a moving mesh method provided by commercial finite element software. For solid propellant grains with high geometrical shape complexity and the topological structure of the solid propellant grains changing along with layer-by-layer combustion, the calculation has high probability of failure or inaccurate result generation.

Shortest distance field method: similar to the level set method, but cannot be calculated for the case of non-uniform distribution of combustion speed.

Level set method: during calculation, the surface which is burning and the surface which is not burning need to be distinguished and processed in each step, and the calculation speed is slow. When sharp corners exist in the solid propellant grain, in the area of the external angle bisector of the edge line and the angular point, because the gradient estimation based on the numerical difference theory has a systematic error, the combustion surface calculation in the area is inaccurate, so that the shape of the combustion surface in the area is distorted. When the burning rate distribution is not uniform, each step of calculation needs to be performed with multiple reinitialization operations, and the calculation speed is low. When the burning rate distribution has a sudden change, the numerical difference near the discontinuous point is continuously changed instead of the sudden change, so that the numerical dissipation phenomenon and the floating effect are serious, and the calculation accuracy is lowered.

Therefore, in order to solve the problem of combustion surface recession calculation of the solid rocket engine with a complex shape and complex combustion speed distribution, the invention provides a combustion surface recession calculation method of the solid rocket engine.

Disclosure of Invention

In order to overcome the defects of the prior art, the invention provides a method for calculating the combustion surface recession of a solid rocket engine.

In order to achieve the above purpose, the invention provides the following technical scheme:

a method for calculating the combustion surface recession of a solid rocket engine comprises the following steps:

step 1, dispersing a propellant charge geometric model of a solid rocket engine into triangular patch models, and classifying all triangular patches into two triangular patch sets; the set of triangular patches starting to burn at the moment of ignition is marked as B, and all triangular patches in B form a thin-sheet-shaped curved surface; the set of all triangular patches in the model is marked as G, and all triangular patches in G form a closed geometric body;

step 2, dividing a Cartesian grid into three-dimensional envelope boxes where the input solid rocket engine is located, and recording the node coordinate of each node of the grid as x;

and 3, for the triangular patch set B, calculating the shortest distance from a node coordinate x to any triangular patch in the triangular patch set B by using an AABB (Axis-aligned bounding boxes) fast distance field algorithm on each node of the Cartesian grid to form a shortest distance field function defined on all nodes of the Cartesian grid

Step 4, calculating the burning rate r of the propellant at each node coordinate x based on a solid propellant burning rate model to form a burning rate distribution field r (x) on a Cartesian grid;

wherein r (x) ap ^ n, a and n are respectively the burning rate coefficient and the pressure index of the propellant at the node coordinate x and are obtained by physical parameters of the propellant; p is the pressure in the combustion chamber of the solid rocket engine at the current moment;

step 5, on each node, performing the calculationModifying field functionsThe value at each node coordinate x;

step 6, Marching cubes and the like are usedValue-plane extraction algorithm based on new field functionsGenerating representations by triangular patchesThe isosurface records the set of all triangular patches in the isosurface as I;

step 7, performing geometric Boolean operation, cutting off the part of the isosurface I outside the closed geometric body G, and updating the I set to only contain the triangular patch in the cut isosurface; at the moment, the area of the sheared I is the area of the combustion surface of the propellant at the current moment;

step 8, multiplying the area of each triangular surface patch t forming the I by the combustion speed at the position, and then multiplying by the density of a propellant, namely the gas generation rate of the position of the current time t, summing the gas generation rate data corresponding to all the triangular surface patches to obtain the total gas generation rate of the solid rocket engine at the current time, and recording the combustion surface area and all the gas generation rate data;

step 9, calculating the pressure p at the current moment in the combustion chamber of the solid rocket engine based on an inner ballistic curve numerical model (such as a quasi-steady-state or transient zero-dimensional inner ballistic model) according to the combustion surface area and the gas generation rate data;

step 10, regarding the sheared I as the surface of combustion at the beginning, replacing the set B with the set I, calculating the shortest distance from the node to the set B of the triangular patches by using an AABB quick distance field algorithm on each node of the Cartesian grid to form a new shortest distance field function defined on all nodes of the Cartesian gridTraverse all mesh nodes, check for correspondencesWhether less than zero; if so, the method will be usedIs modified toFinally, using the modifiedReplace original

Step 11, based on the new one obtained in step 10And new p obtained in step 9, and circularly executing the steps 4-10; until I does not contain any triangular patch; at the moment, the propellant grains are completely burnt out;

the combustion surface area change curve, the gas generation rate curve and the pressure p curve changing along with time are combustion surface retreating calculation results.

Preferably, in the step 2, a three-dimensional envelope box in which the input solid rocket engine is located may be further divided into a conformal structure or a non-structural grid, and in this case, the step 6 selects another isosurface extraction algorithm of a grid type.

Preferably, the step 7 specifically includes:

traversing each triangle t in the I, and judging whether the t is positioned inside the G or not;

if t is inside G, no operation is performed;

if t is outside the closed geometry represented by G, then t is deleted from I;

if t intersects with triangles in G (the intersection decision process can be accelerated by using AABB (Axis-aligned bounding boxes)), then t is divided into a part inside G and a part outside G along the intersection line, the part inside G is divided into triangles (depending on the specific shape), and the triangles are added into I, and t is deleted from I.

The method for calculating the combustion surface recession of the solid rocket engine has the following beneficial effects:

(1) the whole calculation flow does not relate to numerical difference operation, so that the calculation force is saved, and the numerical dissipation phenomenon and the floating effect of the numerical difference operation in the combustion surface shape distortion and the vicinity of the combustion speed discontinuity point brought by the edge line and the angular point external angle bisector area are avoided.

(2) Using the method of regenerating the distance field at each step, rather than performing a re-initialization, avoids the computational power consumption of re-initialization (when the burn rate distribution is more uneven, regenerating the distance field using the AABB fast distance field algorithm is much faster than solving the re-initialization equation).

(3) Ignoring the form factor of the non-burning surface of the solid propellant charge when calculating the change in geometry; after the implicit calculation of the change of the geometric shape is finished, the part which does not exist in the Marching cubes result is cut off by a geometric Boolean calculation method, and the calculation power consumption caused by distinguishing a burning surface from a non-burning surface is avoided.

Drawings

In order to more clearly illustrate the embodiments of the present invention and the design thereof, the drawings required for the embodiments will be briefly described below. The drawings in the following description are only some embodiments of the invention and it will be clear to a person skilled in the art that other drawings can be derived from them without inventive effort.

FIG. 1 is a classification method of a triangular patch model of propellant charge of a solid rocket engine;

FIG. 2 is a Cartesian grid divided at the location of a three-dimensional envelope box in which the solid rocket motor charges are located;

FIG. 3 is an iso-surface extracted from the shortest distance field function value based on the Marching cubes iso-surface extraction algorithm;

figure 4 is an iso-surface with the portion of the propellant grain outside having been cut away.

Detailed Description

In order that those skilled in the art will better understand the technical solutions of the present invention and can practice the same, the present invention will be described in detail with reference to the accompanying drawings and specific examples. The following examples are only for illustrating the technical solutions of the present invention more clearly, and the protection scope of the present invention is not limited thereby.

Example 1

The invention provides a method for calculating the combustion surface retreating of a solid rocket engine, which comprises the following steps:

step 1, dispersing a propellant charge geometric model of the solid rocket engine into a triangular patch model, and classifying all triangular patches into two triangular patch sets. Specifically, in this embodiment, the method for classifying the triangular patch specifically includes: the set of triangular patches that start to burn at the moment of ignition is denoted as B, and all triangular patches in B form a thin-sheet-shaped curved surface. The set of all triangle patches in the model is denoted as G (all triangle patches here include all triangle patches, B is a subset of G), and all triangle patches in G constitute a closed geometry. As shown in fig. 1, the bold lines and the dashed lines together make up the geometric model of the solid rocket engine propellant charge (set G), but the set "surface initially burned" contains only the dashed line (set B).

And 2, as shown in fig. 2, dividing a cartesian grid into three-dimensional envelope boxes in which the input solid rocket engine is located, and recording the node coordinates of each node of the grid as x.

And 3, for the triangular patch set B, calculating the shortest distance from a node coordinate x to any triangular patch in the triangular patch set B by using an AABB (Axis-aligned bounding boxes) fast distance field algorithm on each node of the Cartesian grid to form a shortest distance field function defined on all nodes of the Cartesian grid

And 4, calculating the burning rate r of the propellant at each node coordinate x based on the solid propellant burning rate model without considering whether the node belongs to a burning surface or a non-burning surface to form a burning rate distribution field r (x) on a Cartesian grid.

Wherein, r (x) ap ^ n, a and n are respectively the burning rate coefficient and the pressure index of the propellant at the node coordinate x, and are obtained by physical parameters of the propellant. And p is the pressure in the combustion chamber of the solid rocket engine at the current moment.

Step 5, on each node, performing the calculationModifying field functionsThe value at each node coordinate x.

Step 6, using Marching cubes isosurface extraction algorithm based on new field functionGenerating representations by triangular patchesThe iso-surface is represented by the set of all triangular patches in the iso-surface as I, I is composed of a large number of triangular patches, as shown in fig. 3.

And 7, performing geometric Boolean operation, cutting off the part of the isosurface I outside the closed geometric body G, and updating the I set to only contain the triangular patch in the cut isosurface. As shown in fig. 4, the iso-surface from which the portion outside the propellant grain has been cut off is also the surface of combustion at the present time, and at this time, the area of the cut I is the area of the surface of combustion of the propellant at the present time, which specifically includes the following steps:

go through each triangle t in I, and determine if t is inside G.

If t is inside G then nothing is done.

If t is outside the closed geometry represented by G, t is deleted from I.

If t intersects with triangles in G (the intersection decision process can be accelerated by using AABB (Axis-aligned bounding boxes)), then t is divided into a part inside G and a part outside G along the intersection line, the part inside G is divided into triangles (depending on the specific shape), and the triangles are added into I, and t is deleted from I.

And 8, multiplying the area of each triangular surface patch t forming the I by the combustion speed at the position, and then multiplying by the density of a propellant, namely the gas generation rate of the position of the current moment t, summing the gas generation rate data corresponding to all the triangular surface patches to obtain the total gas generation rate of the solid rocket engine at the current moment, and recording the combustion surface area and all the gas generation rate data.

And 9, calculating the pressure p at the current moment in the combustion chamber of the solid rocket engine based on an internal ballistic curve numerical model (such as a quasi-steady-state or transient zero-dimensional internal ballistic model) according to the combustion surface area and the gas generation rate data.

Step 10, regarding the sheared I as the surface of combustion at the beginning, replacing the set B with the set I, calculating the shortest distance from the node to the set B of the triangular patches by using an AABB quick distance field algorithm on each node of the Cartesian grid to form a new shortest distance field function defined on all nodes of the Cartesian gridTraverse all mesh nodes, check for correspondencesWhether less than zero. If so, the method will be usedIs modified toFinally, using the modifiedReplace original

Step 11, based on the new one obtained in step 10And new p obtained in step 9, and the steps 4-10 are executed circularly. Until I does not contain any triangular patches. At this point, the propellant grains are completely burned.

The combustion surface area change curve, the gas generation rate curve and the pressure p curve changing along with time are combustion surface retreating calculation results.

The whole calculation flow of the method does not relate to numerical difference operation, so that the calculation force is saved, and the numerical dissipation phenomenon and the floating effect of the numerical difference operation in the outer angle bisector area caused by the shape distortion of the combustion surface and the vicinity of the discontinuous point of the combustion speed are avoided.

The above-mentioned embodiments are only preferred embodiments of the present invention, and the scope of the present invention is not limited thereto, and any simple modifications or equivalent substitutions of the technical solutions that can be obviously obtained by those skilled in the art within the technical scope of the present invention are included in the scope of the present invention.