Nonlinear numerical iterative solution method for full-speed domain numerical simulation of complex aircraft

Abstract

The present invention discloses a nonlinear numerical iterative solution method for full-speed domain numerical simulation of complex aircraft, including: generating anisotropic grids in the calculation domain based on metric tensors, constructing a numerical discrete format and RANS turbulence model for compressible flow, and using the Gauss-Seidel method to iteratively calculate the nonlinear budget sub-group of the RANS partial differential equations, thereby reducing the high-frequency errors that may occur in the anisotropic grid during the iteration process. In order to quickly access grid nodes and ensure the accuracy of node interpolation solutions, a node template enhancement strategy and an efficient coloring traversal algorithm are implemented; and a generalized conjugate residual method is used to perform static residual evaluation of forward and backward iterative solutions, thereby greatly improving the computational efficiency and robustness. The present invention is good at processing anisotropic mixed grids, which can not only ensure the iterative robustness of numerical solutions on large aspect ratio boundary layer grids, but also ensure the convergence efficiency of numerical solutions on large-size coarse grids, thereby ensuring the efficiency and stability of the entire nonlinear equation solving iterative process.

Description

Technical Field

[0001] The present invention relates to the field of computational aerodynamics, and in particular to improving the efficiency and computational robustness of numerical simulation solutions for aerospace vehicles, and specifically to a nonlinear numerical iterative solution method for full-speed domain numerical simulation of complex aircraft. Background Art

[0002] Low-cost research and development and modification design of aircraft require efficient computational aerodynamic analysis and design tools. On the one hand, computational aerodynamic tools use high-performance parallel computing platforms and numerical simulation methods to simulate the aerodynamic shape performance of aircraft, which can greatly reduce the number of wind tunnel tests and save test costs; on the other hand, computational aerodynamic numerical analysis tools assisted by parametric design platforms can shorten the aerodynamic design cycle of aircraft and significantly improve the design efficiency of aircraft. Computational aerodynamic analysis and design tools need to solve the Navier-Stokes control equations of flow motion, which is an unsteady, strongly nonlinear and highly coupled set of partial differential equations. Theoretical solutions are very difficult and are usually solved numerically. However, the computational efficiency and robustness of numerical solutions are closely related to the complexity of the aerodynamic flow field of actual problems.

[0003] The surface of subsonic, supersonic and hypersonic aircraft has very complex flow structures, and flow phenomena such as turbulent transition, boundary layer separation and shock / expansion waves often occur. In order to make a detailed description and accurate analysis of these flow field structures, in addition to high-efficiency and high-fidelity numerical solution technology, special processing of the computational grid is also required. Anisotropic hexahedral or prismatic grids are generated in the near-wall area of ​​the aircraft, and isotropic tetrahedral grids are generated in the area far from the wall. It is precisely because of the mixed application of anisotropic and isotropic grids that the convergence efficiency and iterative robustness of the numerical solution process are severely restricted.

[0004] The present invention proposes a highly efficient and robust nonlinear numerical iterative method, and couples it with anisotropic grid generation technology and high-performance parallel technology to form a new numerical simulation solver for complex aircraft aerodynamics. It is used to improve the processing speed of computational aerodynamic analysis and design tools for anisotropic grids, and to improve the robustness, computational efficiency and accuracy of Reynolds average numerical solutions. This technology can be effectively applied to the aerodynamic design of complex aircraft. Summary of the invention

[0005] The present invention proposes a nonlinear numerical iterative solution method for full-speed domain numerical simulation of complex aircraft, with the aim of improving the processing speed of computational aerodynamic analysis and design tools for anisotropic grids in complex aircraft, and improving the robustness, computational efficiency and accuracy of the numerical solution.

[0006] The nonlinear numerical iterative solution method for full-speed domain numerical simulation of complex aircraft includes:

[0007] Firstly, anisotropic mesh is generated for the computational domain of complex aircraft based on the metric tensor constructed based on the flow field gradient characteristics.

[0008] Construct a RANS turbulence model and discretize each term of the Reynolds averaged RANS equation;

[0009] Based on preprocessing operator iteration, node template enhancement and generalized conjugate residual evaluation control, the nonlinear operator of the Reynolds-averaged RANS equation is iteratively solved using the Gauss–Seidel method.

[0010] In the above technical solution, further, the anisotropic grid is generated based on the metric tensor constructed based on the flow field gradient characteristics, specifically including:

[0011] All metric tensors in the flow field constitute the metric space, which is:

[0012]

[0013] Where M is a symmetric positive definite matrix, called the metric tensor matrix, and x=[x, y, z] is the coordinate of the grid point in Euler space;

[0014] The metric tensor matrix M is a metric that describes the shape and size characteristics of the grid unit. Its mathematical expression is:

[0015]

[0016] The change information of the flow field gradient characteristics is expressed by the Hessian matrix To describe, (i=1, 2, 3) and They are The eigenvalues ​​and eigenvectors of the matrix; using this property of the metric tensor matrix, a conventional unstructured grid is first generated in the metric space, then converted to the Euclidean space, and finally an anisotropic grid is generated.

[0017] Furthermore, each term of the Reynolds-averaged RANS equation is discretized. Specifically, the convection term is discretized using Roe flux difference splitting or van Leer flux vector splitting, and a second-order upwind format is used; the diffusion term is discretized using a second-order central format, so that the accuracy of the entire average flow calculation reaches second order.

[0018] Furthermore, the Reynolds averaged RANS discretization equation is a nonlinear system of equations, and its general form can be expressed as:

[0019] R ( Q ) = 0 Q = [ r , u , v , w ,T ]

[0020] Where Q is the solution flux about density, velocity and temperature, ρ, u, v, w, T are density, velocity in x, y and z directions, and temperature respectively, and R is the nonlinear expression of the conserved flux Q, called the nonlinear operator.

[0021] Furthermore, in order to iteratively calculate the solution flux, the nonlinear operator R(Q) is expressed in a preprocessed form as:

[0022]

[0023] Among them, Q n and Q n+1 They represent the numerical solution fluxes of the nth and n+1th iteration calculations respectively; represents the Jacobian matrix form of the nonlinear operator with respect to the flux Q; V is the control volume of the grid cell; is the pseudo time step of the numerical calculation;

[0024] Furthermore, in the process of iterating the nonlinear preprocessing operator, the Gauss–Seidel method is used for iterative solution. In this process, the coloring algorithm is used to determine the order of grid cell access during iteration, and the specific process is as follows: at the beginning, all grid cells are assigned color 0; then each grid cell is traversed in turn, and an attempt is made to assign color 1; but if any cell adjacent to the current grid cell has been assigned color n=1, the attempt to assign fails, and the color index is increased by 1, that is, n=n+1, and the attempt is repeated until the color is successfully assigned to the current cell; when iterating the nonlinear preprocessing operator, the grid cell with color n=1 is first visited, and then the grid cell with color n=2 is visited, and so on; in this way, the stability of the Gauss–Seidel iteration is improved.

[0025] Furthermore, the nonlinear operator represents the discretization of the Reynolds average equation on the grid unit, wherein the physical variables are stored at the center of each grid unit, and the physical solution on the grid node n is obtained by weighted interpolation of the physical solution of the center k of the grid unit connected to the node, and all the grid units connected to the node constitute a node template; if the node template cannot meet any of the three distance threshold conditions, the node template needs to be expanded for enhancement, and the operation process is to add each grid unit adjacent to the grid unit in the node template to the template to become an extended template, and evaluate again whether the extended template meets one of the three conditions until it is met; the three distance threshold conditions are:

[0026] ① ;② ③

[0027] in, is the vector from node n to the center k of any unit in the template, represents the vector size, n corresponds to the grid node, k corresponds to the center of the template unit, , , Respectively Take the vectors corresponding to the minimum, maximum, and intermediate values. ε represents the minimum distance threshold, which is usually taken as the minimum distance from the center of the grid unit to the viscous surface.

[0028] Furthermore, generalized conjugate residual evaluation control is used in the iterative process to evaluate the static residual after each iteration, and then the solutions of all units are updated. The convergence of the static residual is measured by the residual norm before and after the iteration. If the ratio of the current iteration step to the initial residual root mean square norm is less than the set target residual 0.1, the nonlinear operator iteration is completed and terminated; if at the end of the iteration, the root mean square norm of the static residual increases by more than 2 times, the nonlinear operator iteration fails, and the under-relaxation parameters and pseudo time step are adaptively adjusted.

[0029] Furthermore, before each iteration, the physical solution needs to be checked for regularity. When a non-physical solution appears, it needs to be normalized or the solution of the previous iteration step needs to be retained. After the regularization check is completed, the next round of iteration step can be advanced.

[0030] The nonlinear numerical iterative solution method for full-speed domain numerical simulation of complex aircraft proposed in the present invention generates anisotropic grids in the calculation domain based on the metric tensor, constructs the numerical discrete format and RANS turbulence model of compressible flow, and uses the line implicit method to iteratively calculate the nonlinear budget sub-group of the RANS partial differential equations, thereby reducing the high-frequency errors that may occur in the anisotropic grid during the iteration process. The node template enhancement strategy and the efficient coloring traversal algorithm are used to quickly access the grid nodes and ensure the accuracy of the node interpolation solution. At the same time, the generalized conjugate residual method is used to perform static residual evaluation of the forward and backward iterative solutions, so as to adjust the calculation control parameters in time, thereby greatly improving the calculation efficiency and robustness. The present invention can be deployed on a parallel computing cluster based on the Linux kernel, and can be applied to the efficient numerical simulation and fast iterative analysis of the aerodynamic shape of complex aircraft. The main beneficial effects are:

[0031] 1) This method is suitable for characterizing complex flow phenomena such as turbulent transition, boundary layer separation, and shock / expansion waves on the surface of subsonic, supersonic and hypersonic vehicles, especially with high computational accuracy in describing anisotropic flow field structures.

[0032] 2) This method is good at handling anisotropic grids, which can not only ensure the iterative robustness of the numerical solution on the large aspect ratio boundary layer grid, but also ensure the convergence efficiency of the numerical solution on the large size coarse grid, thereby ensuring the efficiency and stability of the entire iterative process of solving the nonlinear equations.

[0033] 3) The present invention improves the robustness, computational efficiency and computational accuracy of the nonlinear iterative solution process of the Reynolds mean flow equation and the turbulence equation through three strategies: preprocessing operator iteration strategy, enhanced node template and generalized conjugate residual evaluation control. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] Figure 1 : is a flowchart of iterative calculation of flow field using the nonlinear numerical iterative solution method of the present invention;

[0035] Figure 2 : is a schematic diagram of a geometric calculation model of Example 1 of the present invention;

[0036] Figure 3 : is a graph showing the variation of the residual error of the mean flow with the number of iteration steps in Example 1 of the present invention;

[0037] Figure 4 : is a visualization cloud diagram of the flow field result on the wing surface of a supersonic aircraft according to Example 2 of the present invention, which is the Mach number distribution on the wing surface;

[0038] Figure 5 : The convergence results of the iterative residual (a) and the maximum viscosity coefficient (b) of Example 2 of the present invention are compared with those of the commercial software;

[0039] Figure 6 : This is Example 3 of the present invention regarding the calculation grid and flow field result visualization of the external flow field of the XBL-100 light sport aircraft.

[0040] Figure 7 : This is Example 4 of the present invention regarding the distribution of the spaceplane calculation grid (left) and the Mach number isosurface (right). DETAILED DESCRIPTION

[0041] The technical solution of the present invention is further described in detail below in conjunction with the accompanying drawings and specific examples.

[0042] The present invention provides a nonlinear numerical iterative solution method for full-speed domain numerical simulation of complex aircraft, which is mainly used for accurate numerical simulation of subsonic or supersonic complex flow phenomena, and can significantly improve the convergence efficiency and computational robustness of nonlinear flow fields. According to a specific embodiment of the present invention, its main implementation process includes the following (see Figure 1 ):

[0043] Step 1: Generate anisotropic grids for the external flow field of the aircraft based on the metric tensor. The metric tensor reflects the density information of the flow field domain grid and constitutes the calculation domain of the numerical simulation.

[0044] Step 2: Perform Reynolds averaging on the NS equations, construct a turbulence model to close the RANS control equations, and construct a numerical discrete format to discretize the equation group.

[0045] Step 3: Iterate and solve the nonlinear preprocessing operator on the computational domain grid, and use the Gauss–Seidel iterative method to reduce the high-frequency error of the anisotropic grid. Use the coloring algorithm to traverse and access the grid cells of the entire computational domain to improve the stability of the Gauss–Seidel method iteration.

[0046] Step 4: Interpolation template construction of grid node solution. The physical solution on the grid node is obtained by weighted interpolation of the solution on the center of the grid unit connected to the node, so all the units connected to the node constitute a template. Evaluate whether the node template meets the three distance threshold conditions, and use the enhancement strategy to expand the node template that does not meet any of the conditions.

[0047] Step 5: Use the generalized conjugate residual method to evaluate the static residual after each iteration, and then update the solution of all units. The convergence of the static residual is measured by the residual norm before and after the iteration. If the residual meets the control condition, the preconditioning operator iteration is completed and the nonlinear solution is updated; otherwise, the preconditioning operator iteration fails, and the under-relaxation parameters and pseudo time step are adaptively adjusted, and the nonlinear preconditioning operator iteration is repeated.

[0048] Step 6: When the numerical calculation reaches the set number of iterations and flow time, the flow field outputs the results.

[0049] The specific process is as follows:

[0050] (1) Generation of anisotropic grids. The biggest feature of anisotropic grids is that the grid size control has obvious directionality. In areas where the flow field outside the entire aircraft changes greatly, such as near the boundary layer or the shock wave position, a small-sized encrypted grid is generated to improve the accuracy of capturing the gradient flow field, while a large-sized uniform grid is used in the external field area far away from the aircraft to save computing resources. In the present invention, a metric tensor is used to describe the directionality of the anisotropic grid. The metric tensor at any grid point is represented by a symmetric positive definite matrix, which contains the grid size information of the point in all directions. If the grid size information in each direction is the same, it becomes an isotropic grid. In essence, all metric tensors in the flow field constitute a metric space, so all size calculations for generating anisotropic grids are performed in the metric space.

[0051] Take the two-dimensional ellipse model as an example to illustrate the measurement of anisotropic size. Define any point P(x,y) in Euler space, then the equation of the circle can be expressed as: The three control variables that define the ellipse are the major semi-axis h1, the minor semi-axis h2 and the angle θ, so the ellipse equation can be expressed as:

[0052] x T [ cos i − sin i sin i cos i ] [ 1 / h 1 0 0 1 / h 2 ] [ cos i sin i − sin i cos i ] x = 1

[0053] make

[0054] M = [ cos i − sin i sin i cos i ] [ 1 / h 1 0 0 1 / h 2 ] [ cos i sin i − sin i cos i ]

[0055] but

[0056]

[0057] Where M is a symmetric positive definite matrix, called the metric tensor matrix. It can be seen that the circle equation in the Euler space is transformed into an ellipse after being modified by the metric tensor matrix, so that the size in all directions along the origin O changes from uniform to anisotropic. Define the metric distance between any two points A and B in the metric space The expression is:

[0058]

[0059] in, is the vector from A to B, Yes The metric tensor matrix at . The metric tensor matrix M is a metric that describes the shape and size characteristics of the mesh unit and is the basis for creating high-quality anisotropic meshes. Its mathematical expression is:

[0060]

[0061] The change information of the flow field gradient characteristics is expressed by the Hessian matrix To describe, and They are The eigenvalues ​​and eigenvectors of a matrix, is a flow field variable, usually k can represent pressure, velocity or temperature, x represents the grid point coordinates, and the subscripts i, j=1,2,3 represent the three coordinate directions. It is by utilizing this characteristic of the metric tensor matrix that a conventional unstructured grid is first generated in the metric space, then converted to the Euclidean space, and finally an anisotropic grid is generated.

[0062] (2) Construction of compressible discrete format and RANS turbulence model. For the numerical simulation of unsteady complex aircraft external flow field, the Reynolds-averaged RANS method is used to numerically solve the basic Navier-Stokes control equations of fluid motion. The Reynolds-averaged RANS method requires an ensemble average of the flow field turbulence to obtain the Reynolds-averaged RANS equation. Since the Reynolds stress term caused by the pulsating quantity will inevitably be generated during the ensemble averaging process, the Boussinesq hypothesis is introduced, which assumes that the turbulent Reynolds stress is proportional to the strain, and the proportionality coefficient is expressed by the turbulent viscosity coefficient.

[0063] The calculation of the turbulent viscosity coefficient is constructed by a mathematical model. In the numerical solution method of the present invention, the turbulence model adopts a one-equation SA model or a two-equation SST k-ω equation, which is particularly suitable for characterizing the calculation of the external flow field of complex aircraft.

[0064] For the discretization of each term of the Reynolds averaged RANS equation, the convection term uses Roe flux difference splitting or van Leer flux vector splitting, and is discretized using the second-order accuracy upwind format; the diffusion term is discretized using the second-order center format, so that the accuracy of the entire average flow calculation reaches the second order. Compressible gas satisfies the ideal gas state equation, the gas dynamic viscosity is affected by temperature, and its viscosity coefficient satisfies the Sutherland formula.

[0065] (3) Nonlinear iterative solution control strategy. The nonlinear operator R(Q) represents the discretized form of the Reynolds average equation on the grid cells, where the physical variables are stored at the center of each grid cell, and the physical solution or solution gradient at the grid node is obtained by interpolating the physical solution at the cell center, and the physical solution on the boundary is obtained by extrapolating the solution inside the computational domain. The preprocessing form of the nonlinear operator R(Q) is:

[0066]

[0067] Here, Q n and Q n+1 They represent the numerical solution fluxes of the nth and n+1th iteration calculations respectively; represents the Jacobian matrix form of the nonlinear operator with respect to the flux Q; V is the control volume of the grid cell; is the pseudo time step of the numerical computation.

[0068] In the iterative solution process of the nonlinear preprocessing operator, high-frequency errors are the main source of residual errors in the discretized equations. Therefore, the Gauss–Seidel method is used for iterative solution. Based on the block iteration strategy, the preprocessing operator solutions of all grid cells consistent with the line direction are updated, which alleviates the defects of the Gauss–Seidel method in processing large aspect ratio anisotropic grids and effectively reduces high-frequency errors in all directions.

[0069] The efficiency and stability of the preprocessing operator iteration depends on the order in which the grid cells are visited during iteration. The present invention proposes a coloring algorithm, and the specific process is as follows: at the beginning, all cells are assigned color 0; then each cell is traversed in turn, and color 1 is first tried; if any adjacent cell has been assigned color n=1, the attempt fails. The color index is increased by 1, that is, n=n+1, and the attempt is repeated until the color is successfully assigned to the cell; when the preprocessing operator is iterated, the color n=1 cell is first visited, then the color n=2 cell is visited, and so on. The coloring algorithm ensures that at the end of the coloring process, no two adjacent cells have the same color, which improves the stability of the Gauss–Seidel method iteration.

[0070] (4) Enhancement of node templates. The physical solution on a grid node is obtained by weighted interpolation of the solutions on the centers of the grid cells connected to the node, so all cells connected to the node constitute a template. However, if the node template cannot meet any of the following three conditions, the node template needs to be enhanced.

[0071] ① ;② ③

[0072] here, is the vector from node n to the center k of any unit in the template, represents the vector size, n corresponds to the grid node, k corresponds to the center of the template unit, , , Respectively Take the vectors corresponding to the minimum, maximum, and intermediate values. ε represents the minimum distance threshold, which is usually taken as the minimum distance from the center of the grid unit to the sticky surface. Node template enhancement is to expand the existing template. The operation process is to add the cells adjacent to each template unit to the template to become an extended template, and evaluate again whether the extended template meets the three conditions until the above conditions are met.

[0073] (5) Generalized conjugate residual method to control solution iteration. For the system of equations processed by the preconditioner, the generalized conjugate residual method is used for iterative calculation. The preconditioner generates a new search direction for the generalized conjugate residual. The static residual after each iteration is evaluated, and then the solution of all units is updated. The convergence of the static residual is measured by the residual norm before and after the iteration. If the ratio of the current iteration step to the initial residual root mean square norm is less than the set target residual 0.1, the preconditioner iteration is completed and terminated; if at the end of the iteration, the root mean square norm of the static residual increases by more than 2 times, the preconditioner iteration fails, and the under-relaxation parameter and pseudo time step are adaptively adjusted.

[0074] The maximum number of iterations allowed in the preprocessing operator and the maximum number of search directions allowed in the generalized conjugate residual method are specified. The residual reduction target for the preprocessing operator, generalized conjugate residual, and nonlinear solution update is specified to find the optimal under-relaxation parameter and CFL. In order to reduce the computational cost of the static residual, the static residual evaluation can be performed only at the 10th iteration, the 20th iteration, and every 15 iterations thereafter, which can improve the computational efficiency.

[0075] (6) Regularity check. Before each iteration, the physical solution needs to be checked for regularity, mainly including flow field density, temperature, and internal energy. When a non-physical solution appears, it needs to be normalized or the solution of the previous iteration step needs to be retained. After the regularization check is completed, the next iteration step can be carried out.

[0076] (7) Parallel computing platform. The carrier of the entire numerical iterative solution process is a parallel computing cluster based on the Linux kernel. MPI realizes data communication between computing nodes and ensures efficient iteration of the solution on the grid nodes of the entire computing domain.

[0077] Implementation Cases:

[0078] The calculation parameters of the following case are set as follows: the SA model is used for the turbulence model; the mean flow convection term uses the second-order accuracy Roe FDS format, and the diffusion term uses the second-order accuracy central format; the convection term of the SA model uses the first-order accuracy upwind format, and the diffusion term uses the second-order accuracy central format; the unit gradient of the source term of the SA model uses the surface weighted surface gradient; these gradients are modified using line implicit line mapping; the maximum number of GS iterations of the preprocessing operator is 500; the linear static residual reduction target of the preprocessing operator is 0.5; the linear residual reduction target is 0.92; the nonlinear residual reduction target is 0.96; the failure condition of the preprocessing operator is that the linear residual increases by more than 2 times, or the linear residual at the end of the iteration exceeds 0.95 times the maximum residual recorded during the iteration.

[0079] Example 1 is a subsonic flow problem around a convex structure in a two-dimensional channel. The free flow Mach number Ma is 0.2, the angle of attack is 0°, and its geometric model is as follows: Figure 2 The entire computational domain uses four grid numbers from coarse to fine: 3649, 56833, 226305, and 903169 for comparative analysis.

[0080] Figure 3 This is the convergence of the mean flow residual monitoring after 500,000 nonlinear iterations, which illustrates the iterative convergence effect of the method of the present invention. It can be seen that from coarse grid to fine grid, the grid size gradually increases, and the rate of decrease of the calculation residual gradually slows down, among which the flow field residual corresponding to the coarse grid decreases the fastest. It can be seen that the nonlinear numerical iterative solution method proposed in the present invention has a very high iterative convergence efficiency for coarse grids, and as the number of grids increases, the calculation convergence speed will decrease.

[0081] Example 2: Flow around a supersonic aircraft wing. The angle of attack α=3°, the free stream Mach number Ma is 1.9, and there is a shock wave near the upper wing suction surface. The average chord length of the wing is 0.64m, the wingspan is 1.12m, and the root-tip ratio is 0.56. Figure 4 and Figure 5 The results of iterative calculation of the solution method of the present invention are shown. Figure 4 The capturing effect of the method of the present invention on supersonic flow and shock wave is described. Figure 5 This reflects that the iterative convergence effect of the method of the present invention is better; it can be seen that a local supersonic area appears at the leading edge of the upper wing, and after the shock wave, the flow velocity gradually becomes a subsonic flow. From the residual convergence, it can be seen that after about 1600s, the residual of the average flow variable tends to a stable state, and the iterative solution reaches convergence, and its calculation convergence efficiency is significantly greater than the calculation result of FLUENT commercial software.

[0082] Example 3: Calculation of the external flow field of an XBL-100 light sport aircraft. The flight altitude is 3000m, the free stream Mach number Ma is 0.15, and the angle of attack is 10°. Figure 6 The anisotropic mesh and pressure field distribution visualization results of the aircraft surface are shown, where (a): the actual shape; (b): the mesh; (c): the pressure coefficient distribution of the entire aircraft surface; (d): the local pressure coefficient distribution (the light and dark on the aircraft surface in (c) and (d) represent different pressures, respectively).

[0083] Example 4: Numerical simulation of the reentry process of a hypersonic spaceplane. The flight altitude is 70,000 m, and the free stream Mach number is 21. Figure 7The anisotropic grid and Mach number flow field distribution on the aircraft surface and its vicinity are displayed. It can be seen that the anisotropic adaptive grid technology proposed in the present invention based on the construction of the flow field characteristic Hessian matrix can very effectively capture the hypersonic shock wave system structure and its interaction, greatly reducing the number of calculation grids and accelerating the calculation progress. At the same time, it also effectively overcomes the calculation difficulties of the low-density area on the leeward side and wake area of ​​the high-altitude high-speed aircraft, and accelerates the calculation convergence speed of the over-expansion area. This implementation example successfully demonstrated the numerical simulation results of the nonlinear numerical iterative solution method proposed in the present invention for the re-entry process of a hypersonic aerospace vehicle.

[0084] The above embodiments are only used to illustrate the design ideas and features of the present invention, and their purpose is to enable those skilled in the art to understand the content of the present invention and implement it accordingly. The protection scope of the present invention is not limited to the above embodiments. Therefore, any equivalent changes or modifications made based on the principles and design ideas disclosed by the present invention are within the protection scope of the present invention.

Claims

1. A nonlinear numerical iterative solution method for full-speed domain numerical simulation of complex aircraft, characterized by: include: Firstly, an anisotropic grid is generated for the flow field outside the aircraft based on the metric tensor constructed based on the flow field gradient characteristics; Specifically include: All metric tensors in the flow field constitute the metric space, which is: x T Mx=1 Where M is a symmetric positive definite matrix, called the metric tensor matrix, and x = [x, y, z] is the coordinate of the grid point in Euler space; The metric tensor matrix M is a metric that describes the shape and size characteristics of the grid unit. Its mathematical expression is: M=V T diag{|λ1|,|λ2|,|λ3|}V The change information of the flow field gradient characteristics is described by the Hessian matrix H, λ i (i=1, 2, 3) and V are the eigenvalues ​​and eigenvectors of the H matrix respectively; using this property of the metric tensor matrix, a conventional unstructured grid is first generated in the metric space, then converted to the Euclidean space, and finally an anisotropic grid is generated; A RANS turbulence model is constructed, and each item of the Reynolds average RANS equation is discretized; the discretization of each item of the Reynolds average RANS equation is specifically: the convection term is discretized using Roe flux difference splitting or van Leer flux vector splitting, and a second-order accuracy upwind format is used; the diffusion term is discretized using a second-order center format, so that the accuracy of the entire average flow calculation reaches the second order; In the computational domain, based on the iteration of preprocessing operators, node template enhancement and generalized conjugate residual evaluation control, the Gauss–Seidel method is used to iteratively solve the nonlinear operator of the Reynolds-averaged RANS equation; in the iterative process of the nonlinear preprocessing operator, the coloring algorithm is used to determine the order of grid cell access during iteration, and the specific process is as follows: at the beginning, all grid cells are assigned color 0; then each grid cell is traversed in turn, and an attempt is made to assign color 1; but if any cell adjacent to the current grid cell has been assigned color n=1, the assignment attempt fails, and the color index is increased by 1, that is, n=n+1, and the attempt is repeated until the color is successfully assigned to the current cell; when iterating the nonlinear preprocessing operator, the grid cell with color n=1 is first visited, and then the grid cell with color n=2 is visited, and so on; in this way, the stability of the Gauss–Seidel method iteration is improved.

2. The nonlinear numerical iterative solution method for full-speed domain numerical simulation of complex aircraft according to claim 1 is characterized in that: The nonlinear operator represents the discretization of the Reynolds average equation on the grid unit, where the physical variables are stored at the center of each grid unit, and the physical solution on the grid node n is obtained by weighted interpolation of the physical solution of the center k of the grid unit connected to the node, and all the grid units connected to the node constitute a node template; if the node template cannot meet any of the three distance threshold conditions, the node template needs to be expanded for enhancement, and the operation process is to add each grid unit adjacent to the grid unit in the node template to the template to become an extended template, and evaluate again whether the extended template meets one of the three conditions until it is met; the three distance threshold conditions are: ①||r nk ||>ε; Among them, r nk is the vector from node n to the center k of any unit in the template, || || represents the vector size, n corresponds to the grid node, k corresponds to the center of the template unit, r n,min 、r n,max 、r n,1 / 2 Respectively represent ||r nk ||Take the vectors corresponding to the minimum value, maximum value, and middle value, and ε represents the minimum distance threshold.

3. The nonlinear numerical iterative solution method for full-speed domain numerical simulation of complex aircraft according to claim 1 is characterized in that: During the iteration process, generalized conjugate residual evaluation control is used to evaluate the static residual after each iteration, and then the solutions of all units are updated. The convergence of the static residual is measured by the residual norm before and after the iteration. If the ratio of the current iteration step to the initial residual root mean square norm is less than the set target residual 0.1, the nonlinear operator iteration is completed and terminated; if at the end of the iteration, the root mean square norm of the static residual increases by more than 2 times, the nonlinear operator iteration fails, and the under-relaxation parameters and pseudo time step are adaptively adjusted.

Images