A solid propellant burning surface shifting method based on the material point method

Abstract

This invention provides a solid propellant burning surface migration method based on the material point method. The method includes the following steps: constructing implicit functions characterizing the burned region, the burning interface, and the unburned region; generating a .k text file containing grid node information and boundary division information using a mesh generation tool; importing the .k file and creating a grid data instance; creating a grid node scalar data instance and discretizing the solid propellant by mesh generation; setting the iteration time step, initial conditions, boundary conditions, and burning surface migration velocity; calculating the grid node scalar gradient at the current time step; updating the grid node scalar data; calculating the amount of solid propellant consumed during combustion within time dt; updating the internal ballistic parameters; and performing time-progression updates based on the grid scalar gradient, time step size, and combustion rate. This invention does not exhibit numerical discontinuities in simulating the large deformation of the propellant combustion process and accurately and reliably simulates the propellant burning surface-wall thickness relationship.

Description

Technical Field

[0001] This invention relates to the field of solid propellant technology, and more specifically to a solid propellant burning surface shifting method based on the material point method. Background Technology

[0002] Rocket engines play a crucial role in propulsion systems, with solid rocket engines (SROs) being widely used in various fields such as rockets and missiles, aerospace, meteorological observation, and disaster mitigation due to their advantages. Because the propellant charge in a solid rocket engine is fixed in advance, its internal trajectory is difficult to adjust during operation. Therefore, accurately predicting the internal trajectory of a solid rocket engine is of great significance. Since the internal trajectory largely depends on the propellant design, the propellant charge is a core aspect of solid rocket engine design. Finding a universal method for calculating the burnout surface displacement of three-dimensional propellants has significant engineering application value.

[0003] Due to the high cost and inherent risks of propellant manufacturing, numerical simulation methods for propellant combustion have become crucial for calculating propellant burnout due to their speed, safety, and low cost. However, traditional finite element mesh-based methods for propellant combustion suffer from mesh distortion, making it difficult to obtain accurate and reliable simulations of the burnout-wall thickness relationship when the propellant topology varies significantly. Summary of the Invention

[0004] In view of this, the present invention provides an algorithm for solid propellant burning surface migration based on the material point method. This method overcomes the mesh distortion defects of the traditional finite element method, and does not exhibit numerical discontinuity problems in simulating the large deformation of the propellant combustion process, thus accurately and reliably simulating the propellant burning surface-wall thickness relationship.

[0005] To achieve the above objectives, the technical solution of the present invention includes the following steps:

[0006] S1 constructs an implicit function characterizing the burned region, the combustion interface, and the unburned region.

[0007] S2 uses a mesh generation tool to mesh the solid propellant, generating a .k text file containing mesh node information and boundary delineation information; the .k file is then imported and a mesh data instance is created.

[0008] S3 creates scalar data instances of grid nodes, resulting in the discretized grid.

[0009] S4 is configured with the following parameters: iteration time step size, initial conditions, boundary conditions, and burning surface movement speed.

[0010] S5 calculates the scalar gradient of the grid nodes at the current time step and updates the scalar data of the grid nodes.

[0011] S6 calculates the amount of solid propellant consumed during combustion within time dt and updates the internal ballistic parameters.

[0012] S7 performs time-advancing updates based on the grid scalar gradient, time step, and burn rate.

[0013] Furthermore, implicit functions scalar Values ​​are assigned based on the physical processes described by the transport equations.

[0014] scalar It has no actual physical meaning; it is only used to characterize the current state of a spatial node, that is, to represent three states: being in a burned region, on a burning interface, or in an unburned region. Therefore, implicit functions can be constructed. As shown in the following formula:

[0015]

[0016] In the formula, Ω - Indicates the inside of an enclosed space. Ω represents the boundary of a closed space. + Indicates the outside of an enclosed space. This represents the coordinates of a point in the spatial computational domain.

[0017] Furthermore, in S2, the solid propellant is meshed using a mesh generation tool, specifically as follows:

[0018] Based on the geometric model of the solid propellant charge, discrete mesh nodes are generated by dividing the geometric space of the charge using modeling tools.

[0019] Furthermore, in S4, initial and boundary conditions are set, including:

[0020] Based on the initial spatial distribution of the burning surface, the scalars on the mesh nodes are initialized according to the following formula, and the corresponding boundary conditions are set:

[0021]

[0022] In the formula, Ω - Indicates the inside of an enclosed space. Ω represents the boundary of a closed space. + Indicates the outside of an enclosed space. This represents the coordinates of a point in the spatial computational domain.

[0023] Furthermore, in S5, when calculating the scalar gradient of the grid nodes at the current time step, the following criteria should be followed based on the actual combustion characteristics of solid propellants:

[0024] The burning surface of solid propellants can only move from the burned area to the unburned area;

[0025] A burned area cannot be turned back into an unburned area;

[0026] The gradient of the burned region, i.e., the grid node with a scalar value of -1, is always zero, expressed as: In the formula, Represents the gradient value. For the i-th grid node in t n scalar of time step Let be the gradient of the implicit function of grid node p at the nth iteration step;

[0027] When calculating the scalar gradient at grid nodes with values ​​in (0,1], to ensure that the propellant burning surface can influence the unburned region and to avoid spurious convergence where the node scalar remains constant, central difference or upwind schemes should be used for gradient calculation. Taking a one-dimensional grid node as an example, when using the central difference scheme... When using the first-order upwind pattern

[0028]

[0029] In the formula, For the i-th grid node in t n scalar of time step For the (i+1)th grid node in t n+1 scalar of time step Let Δx be the gradient of the implicit function of grid node p at the nth iteration step, and Δx be the gradient of the implicit function of grid node i+1 at step t. n+1 The coordinates of the time step and the i-th grid node at time step t n The coordinate difference of the time step, where u is the burning rate of the solid propellant.

[0030] Furthermore, in S7, when performing time-progression updates based on the grid scalar gradient, time step, and burn rate, the following guidelines should be followed:

[0031] Since it is impossible for a burned area grid node to revert to an unburned state, the scalar value of a grid node with a scalar value of -1 remains unchanged.

[0032] The grid nodes on the burning surface will inevitably enter the burned state after the burning surface moves, and their scalar update value will directly change from 0 to -1, and the gradient value will be -1.

[0033] In the solid propellant space, the grid nodes can only be in three states: burned, burning, and unburned, i.e., the three values ​​of -1, 0, and 1 specified in the above formula. However, in the numerical calculation process, the scalar value of the grid node near the burning surface in the unburned area affected by the movement of the burning surface can be in the intermediate state of (0,1).

[0034] Without considering erosion combustion and assuming uniform pressure distribution within the combustion chamber, we assume that a n The combustion interface is equal everywhere.

[0035] The time step size for time iteration should satisfy the CFL condition to avoid iteration divergence.

[0036] For the case where the scalar of the neighboring grid node is (0,1], it indicates that the burning surface has not yet moved to the neighboring node, and the gradient of the node here is 0, thus avoiding false movement of the burning surface.

[0037] For mesh nodes on the boundary, if they are in a burned area or on a burning surface, their state can be updated directly. If they are in an unburned area, virtual nodes can be extended outward along the normal, and their value will still be 1.

[0038] Beneficial effects:

[0039] 1. This invention provides a solid propellant burning surface migration method based on the material point method. The shape function of this method is established based on the background mesh, and a finite element shape function is introduced. Its computational cost is much smaller than other meshless methods. This method overcomes the mesh distortion defect of the traditional finite element method and does not have numerical discontinuity problems in simulating the large deformation of the propellant combustion process. It can accurately and reliably simulate the propellant burning surface-wall thickness relationship.

[0040] 2. The present invention provides a solid propellant burning surface shifting method based on the material point method. The shape function provided by this method has interpolation characteristics, so the boundary conditions of the propellant are easier to apply. The mass points in the material point method are fixedly connected with the background grid nodes, eliminating the need for a large amount of time to search, thus greatly improving the computational efficiency.

[0041] 3. The present invention provides a solid propellant burning surface shifting method based on the material point method. The critical time step of the material point method is determined by the size of the background grid cells. Compared with the general meshless method, which is determined by the minimum spacing between particles, the time step of the material point method can remain stable because the grid size is constant.

[0042] 4. This invention provides a solid propellant burning surface migration method based on the material point method. Due to the single-valued mapping between particles and background grid nodes, the material point method prevents penetration between particles even without any contact algorithm. Therefore, it has significant advantages over existing algorithms in solid propellant burning surface migration. Attached Figure Description

[0043] Figure 1 This is a flowchart of the solid propellant burning surface shifting algorithm based on the material point method in an embodiment of the present invention. Detailed Implementation

[0044] The present invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0045] This invention is a solid propellant burning surface shift algorithm based on the material point method, the specific process of which is as follows: Figure 1 As shown, it includes the following steps:

[0046] Step 1: Construct implicit functions representing the burned region, the combustion interface, and the unburned region.

[0047] Due to the complex burning surface shape of solid propellants, and the abrupt changes in topology leading to discontinuities in the burning surface, complex topological changes occur during burning surface migration. Therefore, this embodiment uses a material point method based on the Eulerian description. The expression for the propellant burning surface migration material point equation is as follows:

[0048]

[0049] Implicit functions are expressed as follows

[0050]

[0051] In the formula, It is an implicit function characterizing a closed entity, a propellant grain, in three-dimensional space, Ω. - Indicates the inside of an enclosed space. Ω represents the boundary of a closed space. + Indicates the outside of an enclosed space. This represents the coordinates of a point in the spatial computational domain.

[0052] During propellant combustion, the parallel-layer combustion law applies. Therefore, the propellant combustion boundary layer can be approximated as the boundary points moving at a constant velocity along the normal direction. At this point, we have: Further results can be obtained In the formula, This indicates the burning rate of the solid propellant, where u represents the burning rate value of the solid propellant. Represents a unit vector. Represents the gradient value. For the scalar gradient of the grid nodes, This represents a scalar value representing a grid node within time t.

[0053] Step 2: Generate a .k text file containing grid node information and boundary division information using a grid generation tool; import the .k file and create a grid data instance.

[0054] Mesh generation software was used to create a mesh for the propellant grain model, and boundary conditions were set. To ensure accurate calculation results and prevent discontinuities, the mesh was generally set to an unstructured mesh. Preferably, to ensure accurate calculation results, the mesh distortion rate was controlled within 10%. In this embodiment, due to the size limitations of the protruding toothed propellant portion, in order to preserve and restore the shape characteristics of the geometric model as much as possible in space, the maximum mesh size was 1 mm, resulting in a total of 20 million volume meshes and 3.3 million mesh nodes.

[0055] Import the .k file using the MeshData handle class and create a mesh data instance;

[0056] Pair By performing first-order precision time discretization, we can obtain the following equation:

[0057] in For the i-th grid node in t n scalar of time step For the i-th grid node in t n+1 scalar of time step For the i-th grid node in t n The scalar gradient at time step u n For solid propellant in t n The burning rate of the time step.

[0058] The MeshData handle class, which processes .k files, primarily reads the corresponding fields from the .k file containing pistol grain mesh information and boundary delimitation information, and generates a node coordinate matrix, an element node number matrix, and a mesh node adjacency information matrix that meet the computational requirements. The fun_node_nb_id() function, which calculates the mesh node adjacency information matrix, is a preprocessing procedure for subsequently reconstructing and calculating the scalar gradients of the mesh nodes.

[0059] Step 3: Create a grid node scalar data instance using the NodeData handle class, divide the grid to discretize the solid propellant, and obtain the discretized grid result.

[0060] This class is used to store handles for scalar and scalar gradient data of mesh nodes and their calculation methods. It sets variables: NodeResults (matrix storing scalar data of mesh nodes); Temp_GXN (matrix storing temporary scalar gradient data between mesh nodes in the X direction); Temp_GYN (matrix storing temporary scalar gradient data between mesh nodes in the Y direction); Temp_GZN (matrix storing temporary scalar gradient data between mesh nodes in the Z direction); NodeGX (matrix storing scalar data of mesh nodes in the X direction); NodeGY (matrix storing scalar data of mesh nodes in the Y direction); NodeGZ (matrix storing scalar data of mesh nodes in the Z direction); h_mesh_data (handle to the mesh data); tlist (time point for computation iteration); N_step (total iteration time steps); i_step (current time step); time (current time); dtime (current time step size); and u_pp (current propellant burning rate).

[0061] The MeshData class, which handles grid node data, primarily processes and stores scalar values ​​on grid nodes during time iterations. It then reconstructs the scalar gradients of unstructured grid nodes based on these node scalar values. The basic principle of its algorithm is to calculate and analyze the scalar gradients between the current grid node and its adjacent nodes in the same cell, and then reconstruct the scalar gradients on the current grid node using the least squares method.

[0062] Step 4: Set the following parameters: iteration time step size, initial conditions, boundary conditions, and burning surface movement speed.

[0063] Set the iteration time step size: Set a function named fun_set_timestep(obj,tlist) to solve for the time step, based on the initially set total combustion time and time step. Preferably, in this embodiment, the combustion time is set to 0.03s and the time step size is 0.001s, and this time step is placed into the time step matrix.

[0064] Initial conditions are set: The node data is initialized using the `fun_nodedata_init(obj,phi)` function. Based on the time step, the corresponding temporary scalar gradient data between storage grid nodes and the zero matrices of the X, Y, and Z directions of the storage grid node scalar data are set for use in subsequent data updates.

[0065] Setting boundary conditions: The function fun_set_boundary(obj,phi) is used to obtain boundary condition data information. It uses a nested loop algorithm to obtain the boundary node data values ​​and assign them to the scalar value phi.

[0066] Setting the burning surface advance velocity: Based on the internal ballistic prediction, an initial burning velocity is set for the burning surface advance to facilitate subsequent iterations. Preferably, in this embodiment, the initial burning velocity is set to 10 mm / s. The selection of the initial burning velocity value does not affect the final result, but it will affect the intermediate calculation process. An unreasonable initial burning velocity will increase the number of iterations during the calculation process, thereby increasing the calculation time.

[0067] Step 5: Calculate the scalar gradient of the grid nodes at the current time step and update the scalar data of the grid nodes;

[0068] Set up the function fun_nodedata_gradient(obj), which is a method of the NodeData handle class, to calculate the scalar gradient of a node.

[0069] (1) Write a nested loop to traverse all nodes on the grid and calculate the scalar gradient between the node and its neighboring nodes: if the number of the neighboring node is less than the number of the current node, it means that the scalar gradient between the two points has been calculated before and there is no need to calculate the scalar gradient again.

[0070] (2) Set up a matrix to store the scalar value of the current node at the current time step: If the scalar values ​​of the two nodes are both less than or equal to 0, or both are greater than 0, it means that the two nodes are in the burned region or in the unburned region, there is no burning surface shift between the two nodes, and the scalar gradient is zero; if the scalar values ​​of the two nodes have opposite signs, it means that one node is in the burned region and the other is in the unburned region, there is burning surface shift between the two nodes, and the scalar gradient is not zero;

[0071] (3) Calculate the distance vector between the current node and its neighboring nodes and the scalar gradient between the current node and its neighboring nodes;

[0072] (4) Count the non-zero elements in each row of Temp_GXN, Temp_GYN, and Temp_GZN, sum them, and then average them.

[0073] Update grid node scalar data: Write the function `fun_nodedata_update(obj)` to update grid node scalar data:

[0074] Start updating the scalar values ​​on the grid nodes over time. If the scalar value on the grid node is less than or equal to -1, it means that the node is already in the burned zone and does not need to be updated. The next time step will directly set it to -1. For "dead zone" nodes, assign the minimum scalar value of the grid node to -1 and update to the next time step.

[0075] Step 6: Calculate the amount of solid propellant consumed during time dt; update the internal ballistic parameters;

[0076] The stored node gradients are obtained and used to calculate the updated internal ballistic parameters for each time step in the solid rocket motor's internal trajectory.

[0077] Step 7: Perform time-advanced updates based on the grid scalar gradient, time step, and burn rate;

[0078] The grid scalar gradient, time step, and burn thickness per unit time (i.e., burn rate) obtained through iteration are updated in the inner trajectory to obtain the above parameter values ​​at each time step. After obtaining the data, the trend is analyzed and summarized.

[0079] The combustion surface shift algorithm based on material point mass (MPM) simulated the combustion process of the propellant in the end-pack form. The inner surface, outer surface and toothed side of the propellant in the circular hole burned simultaneously. The inner surface of the circular hole continuously pushed outward, the outer surface continuously pushed inward, and the toothed side continuously pushed towards the middle. The thickness of the propellant continuously decreased until the inner surface intersected with the outer surface, and the entire propellant was completely burned.

[0080] Compared with traditional methods for propellant combustion calculation based on finite element meshes, the advantages of this invention are as follows: Based on Eulerian-described material points, this invention can obtain more accurate combustion surface state judgment and accurate normal vectors even when the propellant topology changes significantly. The discrete mesh remains constant, and the combustion surface migration process is implicitly contained within the propagation of the implicit scalar transport equation, naturally handling complex situations such as combustion surface intersections and topological changes. This solves the problem of strong abrupt changes and discontinuities in combustion results in previous combustion surface migration calculations.

[0081] In the description of this specification, the references to terms such as "one embodiment," "some embodiments," "illustrative embodiment," "example," "specific example," or "some examples," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the invention. In this specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples.

[0082] Although embodiments of the invention have been shown and described, those skilled in the art will understand that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the claims and their equivalents.

[0083] In summary, the above are merely preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for advancing the burning surface of solid propellants based on the material point method, characterized in that, Includes the following steps: S1 constructs an implicit function characterizing the burned region, the combustion interface, and the unburned region. ; S2 uses a mesh generation tool to mesh the solid propellant, generating a .k text file containing mesh node information and boundary delineation information; the .k file is then imported and a mesh data instance is created. S3 creates scalar data instances of grid nodes to obtain the discretized result of the grid; S4 is configured with the following parameters: iteration time step size, initial conditions, boundary conditions, and burning surface movement speed; S5 calculates the scalar gradient of the grid nodes at the current time step and updates the scalar data of the grid nodes; S6 calculates the amount of solid propellant consumed during combustion within time dt and updates the internal ballistic parameters; S7 performs time-advancing updates based on the grid scalar gradient, time step, and burn rate; In S5, when calculating the scalar gradient of the grid nodes at the current time step, the following criteria should be followed based on the actual combustion characteristics of solid propellants: The burning surface of solid propellants can only move from the burned area to the unburned area; A burned area cannot be turned back into an unburned area; The gradient of the burned region, i.e., the grid node with a scalar value of -1, is always zero, expressed as: In the formula, Represents the gradient value. For the i-th grid node in scalar of time step For the first Grid nodes at each iteration step The gradient of the implicit function; When calculating the scalar gradient at grid nodes with values ​​in (0,1], to ensure that the propellant burning surface can influence the unburned region and to avoid spurious convergence where the node scalar remains constant, central difference or upwind schemes should be used for gradient calculation. Taking a one-dimensional grid node as an example, when using the central difference scheme... The first-order upwind pattern is sometimes used. In the formula, For the i-th grid node in scalar of time step For the (i+1)th grid node in scalar of time step For the first Grid nodes at each iteration step The gradient of the implicit function, For the (i+1)th grid node in The coordinates of the time step and the i-th grid node are in The coordinate difference of the time step, This represents the burning rate of the solid propellant.

2. The solid propellant burning surface shifting method based on the material point method as described in claim 1, characterized in that, The implicit function scalar Values ​​are assigned based on the physical processes described by the transport equations; scalar It has no actual physical meaning; it is only used to characterize the current state of a spatial node, that is, to represent three states: being in a burned region, on a burning interface, or in an unburned region. Therefore, implicit functions can be constructed. As shown in the following formula: In the formula, Indicates the inside of an enclosed space. Indicates the boundary of a closed space. Indicates the outside of an enclosed space. This represents the coordinates of a point in the spatial computational domain.

3. The solid propellant burning surface shifting method based on the material point method as described in claim 1, characterized in that, In step S2, the solid propellant is meshed using a mesh generation tool, specifically as follows: Based on the geometric model of the solid propellant charge, discrete mesh nodes are generated by dividing the geometric space of the charge using modeling tools.

4. The solid propellant burning surface shifting method based on the material point method as described in claim 1, characterized in that, In step S4, initial conditions and boundary conditions are set, including: Based on the initial spatial distribution of the burning surface, the scalars on the mesh nodes are initialized according to the following formula, and the corresponding boundary conditions are set: In the formula, Indicates the inside of an enclosed space. Indicates the boundary of a closed space. Indicates the outside of an enclosed space. This represents the coordinates of a point in the spatial computational domain.

5. The solid propellant burning surface shifting method based on the material point method as described in claim 1, characterized in that, In S7, when performing time-progression updates based on the grid scalar gradient, time step, and burn rate, the following criteria should be followed: Since it is impossible for a burned area grid node to revert to an unburned state, the scalar value of a grid node with a scalar value of -1 remains unchanged. The grid nodes on the burning surface will inevitably enter the burned state after the burning surface moves, and their scalar update value will directly change from 0 to -1, and the gradient value will be -1. In the solid propellant space, the grid nodes can only be in three states: burned, burning, and unburned, i.e., the three values ​​of -1, 0, and 1 specified in the above formula. However, in the numerical calculation process, the scalar value of the grid node near the burning surface in the unburned area affected by the movement of the burning surface can be in the intermediate state of (0,1). Without considering erosive combustion and assuming uniform pressure distribution within the combustion chamber, it is assumed that... The same value is found everywhere at the combustion interface; The time step size for time iteration should satisfy the CFL condition to avoid iteration divergence; For the case where the scalar of the neighboring grid node is (0,1], it indicates that the burning surface has not yet moved to the neighboring node, and the gradient of the node here is 0, thus avoiding false movement of the burning surface; For mesh nodes on the boundary, if they are in a burned area or on a burning surface, their state can be updated directly. If they are in an unburned area, virtual nodes can be extended outward along the normal, and their value will still be 1.

Images