A high dynamic flow field simulation method for large-scale component maneuvers during rocket stage separation

Abstract

The present invention relates to a high-dynamic flow field simulation method for large-scale component maneuvers during separation between rocket stages, including: generating component surface grids, assigning motion type labels to grid cells in each region, and setting the opening time corresponding to each virtual boundary; completing simulation control parameter settings and starting flow simulation; traversing the opening time corresponding to each virtual boundary, and determining the flux on the virtual boundary at the current calculation time step; calculating the physical quantities in the unit at the next time step according to the flux on the virtual boundary at the current calculation time step, integrating the aerodynamic force and torque acting on each component, and calculating the center of mass translation position and rotation around its main inertia axis of each component at the next time step; determining the node positions of grid cells in various motion type regions in the calculation domain, and completing grid motion and grid deformation. The method solves the problems of large calculation amount, insufficient robustness, and switching of boundary condition types over time in the grid deformation dynamic grid.

Description

Technical Field

[0001] The invention belongs to the technical field of computational fluid dynamics and relates to a high dynamic flow field simulation method for maneuvering large-scale components between rocket stages. Background Art

[0002] In the field of aerospace, unsteady flow processes involving large-scale movement of components are an important type of fluid-motion coupling problem, which is widely present in processes such as interstage separation and mount separation. The main method for dealing with such boundary motion problems with a large range of motion is the overlapping moving grid technology. This method first divides a background grid covering the entire calculation domain, and divides an independent grid around each moving component that follows the movement of the component. The evolution of the flow field is achieved by interpolating the overlapping areas between the grids. The second method is to use a method that combines grid deformation with local reconstruction. During the separation simulation process, a certain range around the moving boundary can be set as a grid deformation zone. The grid in this area deforms with the movement of the boundary. When severely distorted grid units appear, the grid in the deformation zone is regenerated, that is, local grid reconstruction is performed. In this process, the flow parameters are mapped from the old grid through interpolation operations. The overlapping dynamic mesh method and the hybrid deformation-reconstruction method have their own advantages: the former can adapt to the unsteady flow simulation of moving parts in various scenarios, but the solution process requires the introduction of grid overlapping area partitioning, hole digging, point finding, interpolation, and parallel overlapping grid implementation due to parallel computing, which brings a large amount of calculation and complex grid processing, and brings great challenges to program development and optimization. At the same time, the overlapping dynamic mesh method is more sensitive to the grid unit scale, and often requires the unit scale of the grid overlapping area to be consistent with each other, which brings high requirements to the pre-processing process. The hybrid dynamic mesh technology of deformation and reconstruction is easier to implement in the pre-processing process and the development of the calculation program because it is always carried out within a set of grids during the simulation process and has the ability of automatic grid reconstruction. However, an important bottleneck of this method is that if the grid deformation process is difficult to maintain a good grid unit quality, the calculation process will converge slowly, resulting in frequent grid reconstruction, making it difficult for the calculation efficiency and calculation accuracy to meet the design requirements.

[0003] In recent years, with the proposal and application of new launch vehicle interstage separation technology, more moving parts with large scale differences in different dimensions have appeared, such as shell fragments, nozzles, etc., which are usually thin-walled structures. If the overlapping dynamic mesh method is used, the grid unit scales of the component grid block and the background grid block are quite different. Therefore, more manual intervention is required in the shell edge area and gap area during the pre-processing process, which increases the difficulty of simulation work and prolongs the cycle. Therefore, it is necessary to establish a mesh deformation method that can maintain a high mesh unit quality, so that the deformation-reconstruction hybrid dynamic mesh method can adapt to the large-scale motion problem of thin-walled structural components. Summary of the invention

[0004] The technical problem solved by the present invention is: in view of the difficult problems existing in the grid deformation dynamic grid, such as large amount of calculation in the grid deformation process, insufficient robustness, and switching of boundary condition types over time, the present invention proposes a high dynamic flow field simulation method for the maneuvering of large-scale components separated between rocket stages, realizes efficient modeling of the complex scenario of the maneuvering process of large-scale components separated between rocket stages, and improves the simulation efficiency and robustness.

[0005] The technical solution of the present invention is: a high dynamic flow field simulation method for large-scale component maneuvers during separation between rocket stages, comprising the following steps:

[0006] S1. Generate the surface mesh of the component, assign the motion type label to the mesh cells in each area, and set the opening time corresponding to each virtual boundary;

[0007] S2. Complete the simulation control parameter settings, including the boundary condition parameters of the far field and the engine jet inlet, the flux solution method, the unit surface physical quantity reconstruction method, realize the physical quantity initialization in the unit, and start the flow simulation;

[0008] S3, traversing the opening time corresponding to each virtual boundary, determining the flux on the virtual boundary at the current calculation time step, including calculating the flux on the virtual boundary and other general unit boundaries;

[0009] S4. Calculate the physical quantities in the unit at the next time step according to the flux on the virtual boundary at the current calculation time step, integrate the aerodynamic force and torque acting on each component, and calculate the center of mass translation position and rotation around its main inertia axis of each component at the next time step;

[0010] S5. Determine the node positions of the grid cells in various motion types in the computational domain according to the center of mass translation position of each component in the next time step and the rotation around its principal inertia axis, complete the grid movement and grid deformation, and repeat steps S3 to S5 to complete the high dynamic flow field simulation of component maneuvers.

[0011] Furthermore, the component surface mesh is generated, the motion type label is assigned to each area mesh unit, and the opening time corresponding to each virtual boundary is set, including:

[0012] s11, generate component surface mesh;

[0013] s12, establish virtual boundaries at the gaps between components and the holes on the surface of components;

[0014] s13, using the front advancing method to generate the component surface mesh boundary layer mesh;

[0015] s14, using the Delaynay method, generate grid cells in the area outside the boundary layer of the component surface in the computational domain by region: if the area outside the boundary layer of the component surface grid contains virtual boundaries, the virtual boundaries are used as the double-sided boundaries of the grid generation, and the entire area is divided into two areas on the left and right sides, and grid cells are generated for each area;

[0016] s15, assign motion type labels to each regional grid cell, including the following types:

[0017] Static area: the units in this area remain stationary during the calculation process, and the motion type labels are batch assigned 0;

[0018] Rigid body motion area, in which the unit will follow a certain component to perform rigid body motion, and the motion type label is consistent with the component number;

[0019] Local deformation area: The elements in this area can be deformed during calculation, and the motion type label is assigned a value of -1 in batches;

[0020] s16, record the number information of the adjacent areas for each virtual boundary: i The number of a virtual boundary is represented by a binary number (ZoneLeft_i, ZoneRight_i), which represents the number of the grid area on the left and right sides of the virtual boundary respectively;

[0021] s17, setting the corresponding opening time for each virtual boundary.

[0022] Furthermore, in step s14, the area number is named in the form of N×10+1, where N is a positive integer.

[0023] Furthermore, the step of setting the corresponding opening time for each virtual boundary in step s17 includes:

[0024] According to the process of separation between rocket stages:

[0025] Set the opening time of the virtual boundary corresponding to the hole on the surface of the interstage separation body component as the opening time of the interstage flame exhaust window;

[0026] Set the opening time of the virtual boundary formed by the gap between the interstage separation parts as the cutting time of the pyrotechnics;

[0027] The opening time of the virtual boundary formed by the gap between the rocket upper stage component, the lower stage component and the separation body component is set as the free separation time;

[0028] Sort the opening time on each virtual boundary and record the opening time of the last virtual boundary, which is recorded as .

[0029] Furthermore, the flux on the virtual boundary at the current calculation time step is determined, including:

[0030] Using logical variables Record whether the virtual boundary opening state needs to be changed. If the current calculation step physical time Greater than the opening time of the last virtual boundary , then there is no need to update the virtual boundary opening state. False; otherwise, execute the following process:

[0031] Using variables Record the result of the predicted opening time of a virtual boundary calculated in the previous time step, using the variable Record the predicted opening time of the next virtual boundary, and traverse each numbered i Virtual border opening time :

[0032] if Less than the current calculation step physical time , then there is no need to consider i The virtual boundary, number i The value of is increased by 1, and the opening time of the next virtual boundary is traversed ; Otherwise, execute the judgment: if Less than the predicted opening time of the next virtual boundary , then the opening time of the next virtual boundary will be predicted Reassign to ;

[0033] The predicted opening time of the next virtual boundary and the result calculated at the previous time step Compare them. If the two are equal, there is no need to change the opening state of a virtual boundary. is false; if the two are not equal, you need to reset the open state on the virtual boundary. is true;

[0034] Calculate the flux on all element faces in the computational domain If a cell face is also turned on on the virtual boundary, it is treated as an internal cell face, otherwise it is treated as a boundary cell face.

[0035] Furthermore, when True to reset the open state on the virtual boundary as follows:

[0036] Traverse each number in order i Virtual border opening time :if Less than the current calculation step time , then continue to traverse numbered i The area numbers on the left and right sides of the virtual boundary (ZoneLeft_i, ZoneRight_i); otherwise, enter the next loop, number i The value of is increased by 1;

[0037] Traverse all element faces in the computational domain: Consider traversing the number j The zone numbers on the left and right sides of the unit are (ZoneLeft_j, ZoneRight_j). To reduce repeated logical judgments on the left and right sides, only check whether ZoneLeft_i×ZoneRight_i is equal to ZoneLeft_j× ZoneRight_j. If they are equal, the number is j The unit is set to on.

[0038] Furthermore, the node positions of the grid units in various motion types in the computational domain are determined according to the center of mass translation position and the rotation around its principal inertia axis of each component in the next time step, including:

[0039] For the grid cells in the stationary type region, their node positions and the previous time step remain unchanged;

[0040] For following i The area where the component performs rigid body motion, its node coordinates Satisfies the following equation:

[0041]

[0042] in, For the i The centroid of the components, For the i The angular velocity of each component about the three principal axes of inertia;

[0043] For the local deformation area, the vertex spring model is used to realize the mesh deformation movement. For the vertex spring, the resultant force on the mesh point is expressed as:

[0044]

[0045] in, To connect nodes i , j The stiffness coefficient of the spring, Is a node The position vector of is the total number of unit nodes in the computational domain, Is with the node The number of connected nodes; for all nodes, the above formula is expressed as a linear system. Assuming that the initial force on the internal nodes is zero, the force on the nodes is expressed as the following system:

[0046]

[0047] in Representative i Line k The matrix coefficients of the columns are calculated as follows:

[0048] .

[0049] Furthermore, the stubbornness coefficient The calculation method is as follows:

[0050]

[0051] in Representative Node i and nodes j The distance is the boundary reinforcement factor, represents the torsion effect, Is a constant value.

[0052] Furthermore, the boundary reinforcement factor The calculation method is as follows:

[0053]

[0054] in, For Node i and nodes j The shortest distance between the center of the connected unit edge and each solid wall boundary, is an artificial empirical distance parameter, which is taken as twice the average size of the component surface mesh unit.

[0055] Furthermore, the torsion effect The calculation method is as follows:

[0056]

[0057] in, Represents a connected node i and nodes j The number of all triangular cells to which the cell edges belong, To connect nodes i and nodes j Among all the triangular elements to which the element edge belongs, m The minimum interior angle of a triangular element.

[0058] The beneficial effects of the present invention compared with the prior art are:

[0059] (1) The present invention adopts mesh motion deformation unit classification to realize a batch-based, customizable dynamic mesh unit type marking method. Under the premise of meeting the simulation accuracy, the area of ​​the unit that needs to be meshed is limited to a limited range, thereby reducing the computational overhead of the mesh deformation processing process. The area that needs to be meshed in the simulation of the maneuvering process of large-scale components separated between rocket stages is reduced to a range that is acceptable to engineering applications.

[0060] (2) The present invention adopts a virtual boundary type switching method to realize an automatic switching method suitable for the boundary types of gaps between separated components between rocket stages and holes in the components themselves. Without redefining the grid boundary type and without significantly increasing the amount of grid boundary type processing calculations, the complex timing of the separation process between rocket stages can be modeled more efficiently, which can reduce the dependence on human participation of simulation designers.

[0061] (3) The present invention adopts a boundary effect correction based on a quadratic smooth function of the wall distance and a torsion correction based on the inner angle of the triangular unit, thereby realizing a more robust strengthening method for the spring model coefficient near the boundary, which can effectively adapt to the complex grid deformation process of the large-scale component maneuvering process of the rocket stage separation, and improve the simulation robustness. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] Figure 1 This is a flow chart of a high dynamic flow field simulation method for large-scale component maneuvers during separation between rocket stages of the present invention. DETAILED DESCRIPTION

[0063] In order to solve the problems of large computational complexity, insufficient robustness and time-dependent switching of boundary condition types in the grid deformation dynamic grid, the present invention proposes a batch-based and customized unit type marking method, a virtual boundary type switching method and a new comprehensive correction method for boundary correction and torsion correction, thereby achieving efficient modeling of the complex scenario of the maneuvering process of large-scale components separated between rocket stages and improving the simulation efficiency and robustness.

[0064] like Figure 1 As shown, a high dynamic flow field simulation method for large-scale component maneuvers between rocket stages includes the following steps:

[0065] Step 1: Generate the surface mesh of the component, assign the motion type label to the mesh cells in each area, and set the opening time corresponding to each virtual boundary;

[0066] Step 2: Complete the simulation control parameter settings, including the boundary condition parameters of the far field and the engine jet inlet, the flux solution method, the unit surface physical quantity reconstruction method, realize the initialization of the physical quantity in the unit, and start the flow simulation;

[0067] Step 3: traverse the opening time corresponding to each virtual boundary, and determine the flux on the virtual boundary at the current calculation time step, including calculating the flux on the virtual boundary and other general unit boundaries;

[0068] Step 4: Calculate the physical quantities in the unit at the next time step according to the flux on the virtual boundary at the current calculation time step, integrate the aerodynamic forces and moments acting on each component, and calculate the center of mass translation position and rotation around its principal inertia axis of each component at the next time step;

[0069] Step 5. Determine the node positions of the grid units in various motion types in the computational domain according to the center of mass translation position of each component in the next time step and the rotation around its principal inertia axis, complete the grid movement and grid deformation, and repeat steps 3 to 5 to complete the high dynamic flow field simulation of component maneuvers.

[0070] The present invention will be further described below in conjunction with the embodiments.

[0071] Example 1

[0072] This embodiment proposes a high dynamic flow field simulation method for large-scale component maneuvers between rocket stages, and the steps are as follows:

[0073] Step S1, generate a component surface mesh, assign a unit motion type label to each area mesh unit, and set the opening time corresponding to each virtual boundary.

[0074] Step S1 specifically includes the following sub-steps:

[0075] Sub-step s11, using mesh pre-processing software to generate component surface mesh.

[0076] Sub-step s12, establishing virtual boundaries at gaps between components and holes on the surfaces of components.

[0077] The purpose of this step is to facilitate the creation of variable type boundary conditions on the created virtual boundary in the subsequent steps, and to use the virtual boundary to distinguish different dynamic mesh types.

[0078] Sub-step s13, using the front advancement method to generate the component surface mesh boundary layer mesh.

[0079] Due to the use of the front advancement method, the generated mesh unit types are prisms, hexagons and pyramids, and do not include tetrahedral units. Specifically, in the process of generating the boundary layer mesh, if the surface of the component is adjacent to the virtual boundary, the adjacent boundary is used as the matching boundary for generating the boundary layer mesh, that is, the unit height of the boundary layer mesh and the adjacent unit height magnification factor are applied to the virtual boundary to form a mesh unit that matches the boundary layer mesh on the component surface. In addition, in order to avoid interference between the boundary layer mesh areas of adjacent components during the generation process, 1 to 2 layers of isotropic tetrahedral unit buffers need to be reserved.

[0080] Sub-step s14, generating grid cells by region. Specifically, using the Delaynay method, grid cells in the region outside the component surface boundary layer in the computational domain are generated by region, with tetrahedron and pyramid cells being used preferentially.

[0081] Specifically, if the area outside the mesh boundary layer of the component surface mesh contains the virtual boundary in sub-step s13, the virtual boundary is used as the double-side boundary of the mesh generation, and the entire area is divided into two areas on the left and right sides and mesh units are generated respectively. In particular, the area number is named in the manner of N×10+1, that is, named in the manner of 11, 21, 31, ...

[0082] Sub-step s15, assigning unit motion type labels to the grid units in each region. According to the motion and deformation modes of the grid units during the calculation process, the units in the flow field calculation domain during the component motion simulation are divided into the following types:

[0083] Static area: the units in this area remain stationary during the calculation process, and the unit motion type labels in this area are batch assigned to 0;

[0084] Rigid body motion area, the unit in this area will follow a certain part to perform rigid body motion. The unit motion type label in this area is consistent with the part number. For example, the unit motion type label of the area that follows part No. 1 to perform rigid body motion is assigned to 1 in batches;

[0085] Local deformation area: The elements in this area can be deformed during the calculation process, and the element motion type labels in this area are batch assigned a value of -1.

[0086] Sub-step s16, recording the numbering information of the adjacent areas for each virtual boundary.

[0087] Specifically, no. i The number of a virtual boundary is represented by a binary number (ZoneLeft_i, ZoneRight_i), which represents the numbers of the grid zones on the left and right sides of the virtual boundary, for example, (11, 21) represents the virtual boundary formed by zone 11 and zone 21.

[0088] Sub-step s17, setting the corresponding opening time for each virtual boundary. Specifically, according to the process of separation between rocket stages:

[0089] Set the opening time of the virtual boundary corresponding to the hole on the surface of the interstage separation body component as the opening time of the interstage flame exhaust window;

[0090] Set the opening time of the virtual boundary formed by the gap between the interstage separation parts as the cutting time of the pyrotechnics;

[0091] The opening time of the virtual boundary formed by the gap between the rocket upper stage component, the lower stage component and the separation body component is set as the free separation time;

[0092] Sort the opening time on each virtual boundary and record the opening time of the last virtual boundary, which is recorded as .

[0093] Step S2, complete the simulation control parameter setting, including the far field and engine jet inlet boundary condition parameters, flux solution method, unit surface physical quantity reconstruction method, realize the physical quantity initialization in the unit, and start the flow simulation.

[0094] Step S3, traverse the opening time corresponding to each virtual boundary, and determine the flux on the virtual boundary at the current calculation time step, including calculating the flux on the virtual boundary and other general unit boundaries.

[0095] As the calculation process proceeds, the multiple virtual boundaries set will be gradually opened at several key time points. In order to avoid confirming the open state of the unit interface on all virtual boundaries at each time step, the present invention uses a logical judgment process to confirm whether the open state of the current virtual boundary needs to be changed.

[0096] Step S3 specifically includes the following sub-steps:

[0097] Sub-step s31, confirming whether the opening state of the virtual boundary needs to be changed according to the physical time of the current calculation step. Specifically,

[0098] Use a logical variable Records whether the virtual boundary opening state needs to be changed.

[0099] If the current calculation step physical time Greater than the opening time of the last virtual boundary , then there is no need to update the virtual boundary opening state. False; otherwise, execute the following process:

[0100] Use a variable The result of calculating the predicted opening time of a virtual boundary in the previous time step is recorded. In this embodiment, the initial value is set to -1.0s.

[0101] Use a variable The predicted opening time of a virtual boundary is recorded, which is 1.0 s in this embodiment.

[0102] Traverse each number in order i Virtual border opening time :

[0103] if Less than the current calculation step physical time , then there is no need to consider i The virtual boundary, number i The value of is increased by 1, and the opening time of the next virtual boundary is traversed Otherwise, continue with the following process:

[0104] if Less than the predicted opening time of the next virtual boundary , then the opening time of the next virtual boundary will be predicted Reassign to ;

[0105] The predicted opening time of the next virtual boundary and the result calculated at the previous time step Compare them. If the two are equal, there is no need to change the opening state of a virtual boundary. is false; if the two are not equal, you need to reset the open state on the virtual boundary. is true.

[0106] Sub-step s32, if the logical variable If true, you need to reset the virtual boundary to enable state, including:

[0107] Traverse each number in order i Virtual border opening time :

[0108] if Less than the current calculation step time , then continue to traverse numbered i The area numbers on the left and right sides of the virtual boundary (ZoneLeft_i, ZoneRight_i); otherwise, enter the next loop, number i The value of is increased by 1;

[0109] Traverse all element faces in the computational domain:

[0110] Consider the traversal number j The zone numbers of the left and right sides of the unit (ZoneLeft_j, ZoneRight_j);

[0111] To reduce repeated logical judgments on the left and right sides, it is only necessary to judge whether ZoneLeft_i×ZoneRight_i is equal to ZoneLeft_j×ZoneRight_j. If they are equal, the number is j The unit is set to the on state. Here, the method of using N×10+1 as the area number in sub-step s16 is used.

[0112] Sub-step s33, calculate the flux on all element faces in the computational domain Specifically, if a cell face is also in an open state on a virtual boundary, the cell face is processed as an internal cell face, otherwise it is processed as a boundary cell face.

[0113] Step S4, calculate the physical quantity in the unit at the next time step according to the flux on the virtual boundary at the current calculation time step, integrate the aerodynamic force and torque acting on each component, and calculate the center of mass translation position and rotation around its main inertia axis of each component at the next time step.

[0114] Step S4 specifically includes the following sub-steps:

[0115] Step S41, perform a time integration based on the flux on the unit surface to obtain the physical quantities on the unit in the calculation domain at the next time step, including pressure ,temperature ,density , the speed in three directions .

[0116] Step S42, using the updated physical quantities on the units in the calculation domain, spatially integrate the unit on the component surface to obtain the aerodynamic force and aerodynamic moment on each component. Specifically, it includes is the aerodynamic force in three directions acting on the center of mass of the ith component, and To act on i The moment about the principal axis of inertia of a component.

[0117] Step S43, based on the aerodynamic force acting on each component, calculate the mass center translation velocity, position, and rotation attitude and rotation velocity around the principal inertia axis of each component in the next time step.

[0118] Specifically, the center of mass translation is given by Newton's laws of motion:

[0119]

[0120] in, For the i The quality of each component, For the i The centroid of the components, To act on i Aerodynamic forces in three directions on the center of mass of a component.

[0121] The rotation of a rigid body about its principal axes of inertia is given by the following Euler equations:

[0122]

[0123] in, To act on i The moments of the three principal inertia axes of each component are: For the i The angular velocity of the components around the three principal axes of inertia, For the i The moment of inertia of a component about its three principal axes of inertia.

[0124] Step S5, determine the node positions of the grid units in various motion type areas in the calculation domain according to the center of mass translation position of each component in the next time step and the rotation around its main inertia axis, complete the grid movement and grid deformation, repeat steps S3~S5, and complete the high dynamic flow field simulation of component maneuvers.

[0125] Step S51: For the grid cells in the static type area, the node positions thereof remain unchanged from the previous time step.

[0126] Step S52, for following i The area where the component performs rigid body motion, its node coordinates Satisfies the following equation:

[0127]

[0128] Step S53: For the local deformation area, the vertex spring model is used to make it more suitable for realizing the mesh deformation movement, and the stiffness coefficient is no longer taken as a constant. For the vertex spring, the resultant force on the mesh point can be expressed as:

[0129]

[0130] in, To connect nodes i , j The stiffness coefficient of the spring, Is a node The position vector of is the total number of unit nodes in the computational domain, Is with the node The number of connected nodes. For all nodes, the above formula can be expressed as a linear system. Considering a linear spring system, without loss of generality, assuming that the initial force on the internal nodes is zero, the node force can be expressed as the following system:

[0131]

[0132] in Representative i Line k The matrix coefficients of the columns are calculated as follows:

[0133]

[0134] Since the linear equation is diagonally dominant, the Jacobi iteration or Seidel iteration format converges. In this embodiment, generally 3 to 4 iterations are required to achieve satisfactory accuracy.

[0135] In order to improve the applicability of the grid deformation method to the large-scale maneuverability problem of components, the present invention adopts a stiffness coefficient calculation model that integrates two correction methods: boundary correction and torsion correction. It can be specifically expressed as:

[0136]

[0137] in Representative Node i and nodes j distance. is the boundary reinforcement factor, represents the torsion effect, .

[0138] In order to adapt to the three-dimensional multi-component movement, the present invention proposes a new boundary reinforcement factor and torsion effect calculation formula based on the above calculation model. Specifically, the boundary reinforcement factor Based on the distance from the vertex to its nearest wall, it is calculated according to the following formula with quadratic smoothness:

[0139]

[0140] in, For Node i and nodes j The shortest distance between the center of the connected unit edge and each solid wall boundary, is an artificial empirical distance parameter, and in this embodiment is taken as 2.0 times the average size of the component surface grid unit.

[0141] Torsion Effect Considering the amplification effect in the large torsion process, the calculation is performed according to the following formula:

[0142]

[0143] in, Represents a connected node i and nodes j The number of all triangular cells to which the cell edges belong, To connect nodes i and nodes j Among all the triangular elements to which the element edge belongs, m The minimum internal angle of a triangular element. Compared with the original method, only the connection nodes are considered i and nodes j The method adopted by the present invention further constrains the torsion degree of the triangular unit and improves the robustness of the mesh deformation during the maneuver simulation of a large range of components.

[0144] Although the present invention has been disclosed as above in the form of a preferred embodiment, it is not intended to limit the present invention. Any person skilled in the art may make possible changes and modifications to the technical solution of the present invention by using the methods and technical contents disclosed above without departing from the spirit and scope of the present invention. Therefore, any simple modifications, equivalent changes and modifications made to the above embodiments based on the technical essence of the present invention without departing from the content of the technical solution of the present invention shall fall within the protection scope of the technical solution of the present invention.

[0145] The contents not described in detail in the specification of the present invention belong to the common knowledge of those skilled in the art.

Claims

1. A high dynamic flow field simulation method for large-scale component maneuvers during rocket stage separation, characterized in that: The following steps are involved: S1. Generate the surface mesh of the component, assign the motion type label to the mesh cells in each area, and set the opening time corresponding to each virtual boundary; The method for setting the opening time of each virtual boundary is as follows: According to the process of separation between rocket stages: Set the opening time of the virtual boundary corresponding to the hole on the surface of the interstage separation body component as the opening time of the interstage flame exhaust window; Set the opening time of the virtual boundary formed by the gap between the interstage separation parts as the cutting time of the pyrotechnics; The opening time of the virtual boundary formed by the gap between the rocket upper stage component, the lower stage component and the separation body component is set as the free separation time; Sort the opening time on each virtual boundary and record the opening time of the last virtual boundary, which is recorded as ; S2. Complete the simulation control parameter settings, including the boundary condition parameters of the far field and the engine jet inlet, the flux solution method, the unit surface physical quantity reconstruction method, realize the physical quantity initialization in the unit, and start the flow simulation; S3, traverse the opening time corresponding to each virtual boundary, determine the flux on the virtual boundary at the current calculation time step, including calculating the flux on the virtual boundary and other general unit boundaries; specifically: Using logical variables Record whether the virtual boundary opening state needs to be changed. If the current calculation step physical time Greater than the opening time of the last virtual boundary , then there is no need to update the virtual boundary opening state. False; otherwise, execute the following process: Using variables Record the result of the predicted opening time of a virtual boundary calculated in the previous time step, using the variable Record the predicted opening time of the next virtual boundary, and traverse each numbered i Virtual border opening time : if Less than the current calculation step physical time , then there is no need to consider i The virtual boundary, number i The value of is increased by 1, and the opening time of the next virtual boundary is traversed ; Otherwise, execute the judgment: if Less than the predicted opening time of the next virtual boundary , then the opening time of the next virtual boundary will be predicted Reassign to ; The predicted opening time of the next virtual boundary and the result calculated at the previous time step Compare them. If the two are equal, there is no need to change the opening state of a virtual boundary. is false; if the two are not equal, you need to reset the open state on the virtual boundary. is true; Calculate the flux on all element faces in the computational domain , if a cell face is also in the open state on the virtual boundary, then the cell face is treated as an internal cell face, otherwise it is treated as a boundary cell face; S4. Calculate the physical quantities in the unit at the next time step according to the flux on the virtual boundary at the current calculation time step, integrate the aerodynamic force and torque acting on each component, and calculate the center of mass translation position and rotation around its main inertia axis of each component at the next time step; S5. Determine the node positions of the grid cells in various motion types in the computational domain according to the center of mass translation position of each component in the next time step and the rotation around its principal inertia axis, complete the grid movement and grid deformation, and repeat steps S3 to S5 to complete the high dynamic flow field simulation of component maneuvers.

2. A high dynamic flow field simulation method for large-scale component maneuvers between rocket stages according to claim 1, characterized in that: Generate the surface mesh of the component, assign motion type labels to the mesh cells in each area, and set the opening time corresponding to each virtual boundary, including: s11, generate component surface mesh; s12, establish virtual boundaries at the gaps between components and the holes on the surface of components; s13, using the front advancing method to generate the component surface mesh boundary layer mesh; s14, using the Delaynay method, generate grid cells in the area outside the boundary layer of the component surface in the computational domain by region: if the area outside the boundary layer of the component surface grid contains virtual boundaries, the virtual boundaries are used as the double-sided boundaries of the grid generation, and the entire area is divided into two areas on the left and right sides, and grid cells are generated for each area; s15, assign motion type labels to each regional grid cell, including the following types: Static area: the units in this area remain stationary during the calculation process, and the motion type labels are batch assigned 0; Rigid body motion area, in which the unit will follow a certain component to perform rigid body motion, and the motion type label is consistent with the component number; Local deformation area: The elements in this area can be deformed during calculation, and the motion type label is assigned a value of -1 in batches; s16, record the number information of the adjacent areas for each virtual boundary: i The number of a virtual boundary is represented by a binary number (ZoneLeft_i, ZoneRight_i), which represents the number of the grid area on the left and right sides of the virtual boundary respectively; s17, setting the corresponding opening time for each virtual boundary.

3. A high dynamic flow field simulation method for large-scale component maneuvers between rocket stages according to claim 2, characterized in that: In step s14, the area number is named in the form of N*10+1, where N is a positive integer.

4. A high dynamic flow field simulation method for large-scale component maneuvers between rocket stages according to claim 3, characterized in that: when True to reset the open state on the virtual boundary as follows: Traverse each number in order i Virtual border opening time :if Less than the current calculation step time , then continue to traverse numbered i The area numbers on the left and right sides of the virtual boundary (ZoneLeft_i, ZoneRight_i); otherwise, enter the next loop, number i The value of is increased by 1; Traverse all element faces in the computational domain: Consider traversing the number j The zone numbers of the left and right sides of the unit are (ZoneLeft_j, ZoneRight_j). To reduce repeated logical judgments on the left and right sides, only check whether ZoneLeft_i*ZoneRight_i and ZoneLeft_j* ZoneRight_j are equal. If they are equal, the number is j The unit is set to on.

5. A high dynamic flow field simulation method for large-scale component maneuvers between rocket stages according to claim 4, characterized in that: The node positions of the grid cells in various motion types in the computational domain are determined based on the translation position of the center of mass of each component in the next time step and the rotation around its principal inertia axis, including: For the grid cells in the stationary type region, their node positions and the previous time step remain unchanged; For following i The area where the component performs rigid body motion, its node coordinates Satisfies the following equation: in, For the i The centroid of the components, For the i The angular velocity of each component about the three principal axes of inertia; For the local deformation area, the vertex spring model is used to realize the mesh deformation movement. For the vertex spring, the resultant force on the mesh point is expressed as: in, To connect nodes i , j The stiffness coefficient of the spring, Is a node The position vector of is the total number of unit nodes in the computational domain, Is with the node The number of connected nodes; for all nodes, the above formula is expressed as a linear system. Assuming that the initial force on the internal nodes is zero, the force on the nodes is expressed as the following system: in Representative i Line k The matrix coefficients of the columns are calculated as follows: 。 6. A high dynamic flow field simulation method for large-scale component maneuvers between rocket stages according to claim 5, characterized in that: Stubbornness coefficient The calculation method is as follows: in Representative Node i and nodes j The distance is the boundary reinforcement factor, represents the torsion effect, Is a constant value.

7. A high dynamic flow field simulation method for large-scale component maneuvers between rocket stages according to claim 6, characterized in that: Boundary Enhancement Factor The calculation method is as follows: in, For Node i and nodes j The shortest distance between the center of the connected unit edge and each solid wall boundary, is an artificial empirical distance parameter, which is taken as twice the average size of the component surface mesh unit.

8. A high dynamic flow field simulation method for large-scale component maneuvers between rocket stages according to claim 6, characterized in that: Torsion Effect The calculation method is as follows: in, Represents a connected node i and nodes j The number of all triangular cells to which the cell edges belong, To connect nodes i and nodes j Among all the triangular elements to which the element edge belongs, m The minimum interior angle of a triangular element.

Images