CN111859529A - A Multi-Grid Disturbance Domain Update Acceleration Method for Numerical Simulation of Flow Around Aircraft - Google Patents

Abstract

The invention discloses a multi-grid disturbance domain update acceleration method for numerical simulation of flow around an aircraft. Aiming at the problem that a large number of invalid calculations exist in the existing multi-grid method when simulating the flow around an aircraft, by adjusting the fine grid and coarse grid Convective and viscous dynamic computational domains are adopted in the above, so that the flow governing equations are solved in the convective dynamic computational domain containing only the unconverged disturbed elements in each iteration step, and the viscous effect is only contained in the viscous effect dominant element. The viscous dynamic calculation domain of , realizes the solution idea of only solving the unconverged disturbed element and only considering the viscous effect in the local area. By effectively reducing the calculation amount of a single iteration step and the number of iteration steps required to reach a convergent state, the invention can significantly improve the calculation efficiency of the multi-grid method for simulating the flow around the aircraft.

Description

A Multi-Grid Disturbance Domain Update Acceleration Method for Numerical Simulation of Flow Around Aircraft

technical field

The invention relates to the technical field of computational fluid dynamics, in particular to a multi-grid disturbance domain update acceleration method for numerical simulation of flow around an aircraft.

Background technique

With the rapid development of computer technology, computational fluid dynamics numerical simulation has become one of the important technical supports for the aerodynamic design of contemporary aircraft. Although numerical simulation can accurately predict the static and dynamic aerodynamic characteristics of aircraft, for the aerodynamic design of aircraft, there is still a problem that the computational efficiency of numerical simulation cannot meet the increasing demand for practical applications. Therefore, it is of great significance and engineering application value to propose an efficient method for numerical simulation of flow around aircraft to improve the efficiency of aircraft design iteration and shorten the development cycle.

Multigrid method is one of the mainstream acceleration techniques for numerical simulation of flow around aircraft. The method can be used not only for the iterative solution of the steady surrounding flow field in the prediction of the static aerodynamic characteristics of the aircraft, but also for the inner iterative solution in the iterative propulsion of the unsteady surrounding flow field in the prediction of the dynamic aerodynamic characteristics of the aircraft. The multi-grid method is an accelerated convergence technique that cyclically solves the flow control equations between coarse and fine grids. The acceleration principle has two aspects: First, the coarse grid can advance the solution with a larger iterative step size, so that it can be solved. Drive the fine mesh to reach the convergence state faster; secondly, the cyclic solution between the coarse and fine meshes helps to attenuate the high-frequency error that hinders the convergence of the fine mesh, so that the overall error of the fine mesh can be attenuated faster .

At present, the multi-grid method still has a drawback in the numerical simulation of the flow around the aircraft, that is, each iteration step must update and solve all the grid elements in the preset calculation domain of the aircraft flow field grid. This global update solution method ignores the characteristics of the evolution of the flow field around the aircraft during the solution process, resulting in different degrees of invalid calculation. For a subsonic aircraft, the mathematical properties of the flow are elliptical. Therefore, in order to avoid the unphysical response caused by the finite boundary of the numerical simulation, the preset computational domain of the aircraft flow field grid needs to be in the range of dozens of times the reference length of the aircraft. Do a simulation. However, as the flow field around the aircraft gradually approaches the steady state, the region that does not satisfy the flow control equations also shrinks accordingly, and the preset computational domain of the aircraft flow field grid will also tend to be redundant. For supersonic vehicles, the mathematical properties of the flow are hyperbolic, that is, any point in the flow field is only affected by the flow in its upstream dependent domain. Therefore, in the iterative process, the region located downstream of the flow field around the aircraft will not converge before its upstream flow field, and the upstream flow that has reached a steady state will not be affected by the update of the unconverged downstream flow. It can be seen that in the entire iterative solution process, the global update and solution of all grid elements in the preset calculation domain of the aircraft flow field grid will generate a large number of invalid calculations, which will significantly reduce the computational efficiency of the numerical simulation of the flow around the aircraft. Therefore, if the invalid calculation in the existing method can be eliminated, the computational efficiency of the numerical simulation of the flow around the aircraft will be significantly improved, thereby effectively solving the problem that the computational efficiency of numerical simulation in the aerodynamic design of the aircraft cannot meet the increasing practical application requirements.

SUMMARY OF THE INVENTION

In view of this, the present invention provides a multi-grid disturbance domain update acceleration method for numerical simulation of flow around an aircraft, so as to solve the problem of a large number of ineffective calculations when simulating flow around an aircraft in the existing multi-grid method.

The invention provides a multi-grid disturbance domain update acceleration method for numerical simulation of flow around an aircraft, which is characterized in that it includes the following steps:

S1: Read in data, including the fine mesh of the aircraft flow field, the preset computational domain, boundary conditions and calculation settings;

S2: Generate a coarse grid with a specified number of layers by thinning the fine grid;

S3: initialize the flow field in the preset computational domain on the coarse grid and the fine grid according to the incoming flow condition or according to a given flow field;

S4: According to the flow field initialization method, establish the convective dynamic calculation domain and the viscous dynamic calculation domain;

S5: Divide the data to be stored for all layers of grid solution into two types; the first type of data is the inherent information of the unit, and the static data structure is used to store the information of all grid units in the preset computational domain, including the grid Coordinates and flow field variables; the second type of data is the information related to the solution update, and the dynamic data structure is used to store only the information of the grid cells in the convective dynamic calculation domain, including the update amount of the conservation quantity and the local iteration step size;

S6: According to the boundary condition type, assign value to the conservation quantity of the boundary virtual grid of the current layer grid;

S7: Divide the residual term of the flow control equation into two types: inviscid term and viscous term; wherein, the inviscid term includes the forcing function of convective flux and coarse grid, and the inviscid term is in the convective dynamic computational domain Calculated in ; the viscous term includes viscous flux and turbulence model equation source term, and the viscous term is calculated in the viscous dynamic computational domain;

S8: in the convective dynamic calculation domain, solve the update amount of the conserved quantity, and update the flow field variable;

S9: determine whether the current grid layer is a fine grid; if yes, go to step S10; if not, go to step S21;

S10: Determine whether the fine mesh reaches the convergence condition; if not, execute step S11; if yes, execute step S23;

S11: Determine whether all the boundary elements of the convective dynamic calculation domain have been subjected to inviscid disturbance one by one; if so, go to step S12; Judging; after traversing all boundary elements of the convective dynamic computational domain, perform step S13;

S12: Measure the propagation direction of the inviscid perturbation, and include the unit located in the propagation direction in the adjacent units of the boundary unit into the convective dynamic calculation domain, and return to step S11 to perform the next boundary unit of the convective dynamic calculation domain. judge;

S13: Judging, one by one, whether all boundary elements of the convective dynamic calculation domain satisfy the four conditions of converged, located in compressible flow, located at the most upstream and no longer affected by other units in the convective dynamic calculation domain; , then go to step S14; if not, go back to step S13 to judge the next boundary unit of the convective dynamic computational domain; after traversing all the boundary units of the convective dynamic computational domain, go to step S16;

S14: Remove the boundary element from the convective dynamic calculation domain, and determine whether the boundary element exists in the viscous dynamic calculation domain; if so, execute step S15; The next boundary element of the convective dynamic computational domain is judged as described above;

S15: remove the boundary element from the viscous dynamic calculation domain, and return to step S13 to judge the next boundary element of the convective dynamic calculation domain;

S16: Determine whether all boundary elements of the viscous dynamic calculation domain are dominated by the viscous effect one by one; if so, go to step S17; if not, return to step S16 to judge the next boundary element of the viscous dynamic calculation domain ; After traversing all boundary elements of the viscous dynamic computational domain, execute step S18;

S17: Include the unit located in the convective dynamic calculation domain in the adjacent units of the boundary unit into the viscous dynamic calculation domain, and return to step S16 to judge the next boundary element of the viscous dynamic calculation domain;

S18: According to the convective dynamic calculation domain and the viscous dynamic calculation domain of the fine grid, update the convective dynamic calculation domain and the viscous dynamic calculation domain of the coarse grid of all layers;

S19: According to the convective dynamic calculation domain, reallocate the storage space of the second type of data in all coarse grids and fine grids;

S20: Interpolate the conserved amount to a coarser grid, return to step S6, and enter the calculation of the coarser grid;

S21: Determine whether a coarser grid needs to be calculated; if yes, go back to step S20; if not, go to step S22;

S22: Interpolate the correction amount of the conservation amount back to a finer grid, and determine whether the finer grid is a fine grid; if so, return to step S6, and enter the calculation of the next iteration step; if not, return to step S6, enter Calculation of finer grids;

S23: output the result.

In a possible implementation manner, in the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention, in step S3, according to the incoming flow conditions or according to a given flow field, the coarse grid and the fine grid are Initialize the flow field in the preset computational domain on the grid, including:

Initializing according to the incoming flow condition, assigning the conservation quantities of all grid cells in the preset computational domain on the coarse grid and the fine grid as the incoming flow value;

According to the initialization of the given flow field, the conservation quantities of all grid cells in the preset calculation domain on the fine grid are assigned as the given flow field values, and the values of all grid cells in the preset calculation domain on the coarse grid The conserved quantities are obtained by interpolating the conserved quantities of all grid cells on the fine grid.

In a possible implementation manner, in the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention, step S4, according to the flow field initialization method, establishes a convective dynamic calculation domain and a viscous dynamic calculation domain, which specifically includes:

According to the initialization of the inflow conditions, for the convective dynamic calculation domain, the coarsest grid layer takes the 1-layer element next to the wall as the initial element, and the fine grid takes the L-layer adjacent element of the wall as the initial element, and L satisfies:

L=2 ^l (1)

Among them, l represents the number of coarse grid layers, and l is a positive integer, which is given by the calculation settings read in step S1; the initial elements of the coarser grids are generated according to the initial elements of the fine grid; for the viscous dynamic computational domain, all The grid layer takes the 1-layer unit next to the wall as the initial unit;

According to the initialization of the given flow field, for the convective dynamic computational domain, the initial cell of the fine grid is the disturbed cell whose flow characteristics in the given flow field are inconsistent with the incoming flow conditions, satisfying:

||WW _∞ ||/||ΔW ⁽¹⁾ || _max ＞ε _a,c (2)

Among them, W represents the conserved quantity; W _∞ represents the conserved quantity of the incoming flow condition; ||ΔW ⁽¹⁾ || _max represents the opposite flow

is the maximum value of the update amount of the conservation quantity of all units in the viscous computational domain in the first iteration step; ε _a,c represent the new threshold for convection, which is 10 ^-4 ≤ε _a,c ≤10 ^-6 ; for the viscous dynamic computational domain, fine The initial element of the grid is the viscous effect dominant element in the fine-grid convection dynamic computational domain, which satisfies:

表示第1迭代步的粘性效应衡量参数；I,J,K分别表示网格方向；

表示对流通量Jacobian矩阵沿i方向的谱半径，

表示粘性通量Jacobian矩阵沿i方向的谱半径；ε_a,v表示粘性新增阈值，取10^-3≤ε_a,v≤10^-2；所有粗网格层的对流动态计算域依据细网格的对流动态计算域生成，所有粗网格层的粘性动态计算域依据细网格的粘性动态计算域生成。Among them, Ψ represents the viscous effect measurement parameter, which is the ratio of viscous disturbance to mass flow without viscous disturbance;

Represents the viscous effect measurement parameter of the first iteration step; I, J, K represent the grid direction respectively;

represents the spectral radius of the Jacobian matrix of convective flux along the i direction,

Represents the spectral radius of the viscous flux Jacobian matrix along the i direction; ε _a,v represents the new viscosity threshold, which is 10 ^-3 ≤ε _a,v ≤10 ^-2 ; the convection dynamic computational domain of all coarse mesh layers is based on the fine mesh The convective dynamic computational domain of the grid is generated, and the viscous dynamic computational domain of all coarse grid layers is generated according to the viscous dynamic computational domain of the fine grid.

In a possible implementation manner, in the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention, in step S7, the residual term of the flow control equation is divided into two types: inviscid term and viscous term; The inviscid term includes the convective flux and the forcing function of the coarse grid, and the inviscid term is calculated in the convective dynamic computational domain; the viscous term includes the viscous flux and the turbulence model equation source term, and the viscous term is in the viscous flow. Computation in the dynamic computational domain, including:

The flow control equation is expressed as:

为无粘项，

为粘性项；在对流动态计算域中求解残差项的无粘项，在粘性动态计算域中求解残差项的粘性项。where t is the time, F _c is the convective flux, F _v is the viscous flux, Q _F is the forcing function of the coarse mesh, Q _T is the source term of the turbulence model equation, |Ω| is the volume of the mesh element, N _f represents the number of element faces, and ΔS represents the area of the element face; the right-hand term of equation (5) is collectively referred to as the residual term of the flow control equation,

for non-sticky terms,

is the viscous term; the inviscid term of the residual term is solved in the convective dynamic computational domain, and the viscous term of the residual term is solved in the viscous dynamic computational domain.

In a possible implementation manner, in the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention, step S11 is to judge whether all boundary elements of the convective dynamic calculation domain have been subjected to inviscid disturbance one by one, specifically including :

Judge whether all boundary elements of the convective dynamic computational domain satisfy one by one:

||ΔW||/||ΔW ⁽¹⁾ || _max ＞ε _a,c (6)

Among them, ||ΔW|| represents the modulo value of the update quantity of the conservation quantity of the boundary element of the convective dynamic computational domain in the current iteration step;

In step S12, the propagation direction of the inviscid disturbance is measured, and the unit located in the propagation direction among the adjacent units of the boundary unit is included in the convective dynamic calculation domain, which specifically includes:

The inviscid disturbance propagates at the speed of sound relative to the flow, and along the propagation direction of the inviscid disturbance, the propagation velocity of the inviscid disturbance is positive; let q represent the unit direction vector, then the inviscid disturbance propagates along the q direction as:

u·q+a＞0 (7)

Among them, u represents the flow velocity vector, a represents the speed of sound; q takes the unit vector of the grid center of the boundary element of the convective dynamic calculation domain pointing to a certain grid point of the boundary element, if q satisfies the formula (7), it will be combined with the boundary of the convective dynamic calculation domain. Cells that share this grid are included in the convective dynamic computational domain.

In a possible implementation manner, in the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention, in step S13, it is judged one by one for all the boundary elements of the convective dynamic calculation domain whether it satisfies the conditions of converged, compressible In the flow, at the most upstream and no longer affected by the other units in the convective dynamic calculation domain are four conditions; if so, go to step S14; The unit makes a judgment; it specifically includes the following steps:

S131: Determine whether all boundary elements of the convective dynamic computational domain satisfy one by one:

||ΔW||/||ΔW ⁽¹⁾ || _max <ε _d (8)

Wherein, ε _d represents the deletion threshold, and ε _d ≤ 10 ⁻⁷ ; if yes, then go to step S132; if not, return to step S131 to judge the next boundary element of the convective dynamic computational domain;

S132: judge whether the Mach number of the boundary unit is greater than 0.3; if so, execute step S133; if not, return to step S131 to judge the next boundary unit of the convective dynamic computational domain;

S133: Determine whether the boundary unit satisfies:

Wherein, q takes the unit vector of the cell center of the boundary cell of the convective dynamic calculation domain pointing to the cell center of the boundary cell immediately adjacent to the cell center, and θ _d represents the tolerance angle of the upstream cell, which is 5° _≤θd≤10 °; if so, execute step S134; If not, return to step S131 to judge the next boundary element of the convective dynamic computational domain;

S134: For supersonic inviscid flow, no further judgment is required, and step S14 is directly executed;

For subsonic flow and viscous flow, determine whether the boundary element satisfies:

||ΔW||+||Δ(ΔW)||<ε _d (10)

Among them, ||Δ(ΔW)|| represents the modulo value of the update quantity of the boundary element conservation quantity affected by the adjacent elements in the convective dynamic calculation domain of the boundary element; Δt represents the iteration step size, C _CFL represents the CFL number of the time-advance format, ΔR ⁱ Represents the influence of the adjacent elements in the convective dynamic computational domain on the residual term of the boundary element along the i direction;

For subsonic inviscid elements, the influence of the adjacent elements in the convective dynamic computational domain on the residual term of the boundary element along the ⁱ direction is expressed as:

Among them, ΔF _c represents the change of convective flux; the subscript i+1 represents the adjacent unit of the boundary element of the convective dynamic computing domain along the positive i direction, and the subscript i-1 represents the boundary element of the convective dynamic computing domain along the negative i direction ΔS _i+1/2 represents the area of the boundary element of the convective dynamic computational domain along the positive i direction; ΔS _i-1/2 represents the area of the boundary element of the convective dynamic computational domain along the negative i direction;

For a viscous element, the influence ΔR ⁱ of the adjacent element in the convective dynamic computational domain on the residual term of the boundary element along the i direction is expressed as:

If yes, go to step S14; if no, go back to step S131 to judge the next boundary unit of the convective dynamic computational domain.

In a possible implementation manner, in the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention, step S16 is to judge whether all boundary elements of the viscous dynamic calculation domain are dominated by the viscous effect one by one, which specifically includes:

Judge whether all boundary elements of the viscous dynamic computational domain satisfy one by one:

Among them, Φ represents the scaling factor of the newly added viscosity threshold, ||ΔW|| _max represents the maximum value of the update amount of the conservation quantity of all elements in the convective dynamic computational domain in the current iteration step, (||ΔW|| _max ) _min represents the th 1 Iteration step to the minimum value of the maximum value of the single-step conservation amount update amount in the current iteration step.

In the multi-grid disturbance domain update acceleration method for the numerical simulation of flow around an aircraft provided by the present invention, the residual estimation only calculates the inviscid term of the residual in the In the dynamic calculation domain, the viscous term of the residual is calculated, and the time integration is only performed in the convective dynamic calculation domain of each layer grid. Therefore, the present invention adopts two types of dynamic calculation domains, convective and viscous, in each layer grid. When the grid method simulates the flow around the aircraft, the calculation amount of each layer of grids significantly improves the calculation efficiency of the multi-grid method for the numerical simulation of the flow around the aircraft, thus effectively solving the problem that the calculation efficiency of numerical simulation is difficult to meet the increasing practical application of aircraft aerodynamic design. question of needs. The data to be stored is divided into two categories: the inherent information of all units and the information related to solving and updating; wherein, for the first type of data, the information of all units in the preset computational domain is stored in a static data structure; for the second type of data data, only the information of grid cells in the convective dynamic computing domain of each grid layer is stored in a dynamic data structure, and the storage space of the second type of data is adjusted according to the convective dynamic computing domain of each grid layer in time at each iteration step. Therefore, , compared with the existing multi-grid method, the invention can effectively reduce the storage space requirement of the numerical simulation of the flow around the aircraft. For the fine grid of the aircraft flow field, only the disturbed elements are added to the convective dynamic calculation domain according to the disturbance propagation characteristics; by weighing the influence of the viscous effect, only the grid elements dominated by the viscous effect are added to the viscous dynamic calculation of the fine grid In the domain; using the solution convergence feature, without affecting the convergence rate, the converged elements are removed from the two types of dynamic calculation domains of convection and viscosity; Convective and viscous two types of dynamic computational domains of the grid; therefore, the present invention can always keep all grid layers in the flow field of the aircraft to solve only the unconverged disturbed element, and only consider the viscosity in the viscous effect dominant element, thus The ineffective calculation of the existing multi-grid method in the numerical simulation of the flow around the aircraft can be effectively avoided, and the problem that the numerical simulation calculation efficiency is difficult to meet the increasing practical application requirements of the aerodynamic design of the aircraft can be solved.

Description of drawings

1 is a flowchart of a multi-grid disturbance domain update acceleration method provided in Embodiment 1 of the present invention;

2 is a schematic diagram of generating a coarse grid from a multi-grid according to a fine grid in Embodiment 1 of the present invention;

3 is a schematic diagram of a multi-grid cycle in Embodiment 1 of the present invention;

Fig. 4 utilizes the multi-grid disturbance domain update acceleration method provided by the embodiment of the present invention to solve the flow field of the transonic NACA0012 airfoil inviscid flow and the evolution process diagram of the convective dynamic calculation domain;

Fig. 5a is the grid quantity change curve of the convective dynamic computational domain for solving the inviscid flow around the transonic NACA0012 airfoil by utilizing the multi-grid disturbance domain update acceleration method provided by the embodiment 1 of the present invention;

Fig. 5b is the convergence curve of solving the inviscid flow around the transonic NACA0012 airfoil using the multi-grid disturbance domain update acceleration method provided by the embodiment 1 of the present invention;

6 is a diagram showing the evolution process of the flow field and two types of dynamic calculation domains of convection and viscosity for solving the flow field of the transonic RAE2822 airfoil turbulent flow around the airfoil using the multi-grid disturbance domain update acceleration method provided in Embodiment 1 of the present invention;

Fig. 7a is the grid quantity change curve of two types of dynamic computational domains, convective and viscous, for solving the turbulent flow around the transonic RAE2822 airfoil by using the multi-grid disturbance domain update acceleration method provided by Embodiment 1 of the present invention;

Fig. 7b is a convergence curve for solving the turbulent flow around a transonic RAE2822 airfoil by using the multi-grid disturbance domain update acceleration method provided in Embodiment 1 of the present invention.

Detailed ways

The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are merely illustrative and not intended to limit the present invention.

The solution process of simulating the flow around an aircraft by using the multi-grid method has the following characteristics: (1) the disturbance of the flow field starts at the point where the governing equations cannot be satisfied; (2) the disturbance carrying discontinuous information will gradually propagate to the direction at a finite speed with time. (3) In compressible flow, the convergence of the upstream flow field will not be later than that of the downstream flow field; (4) Regardless of separation or not, the viscous effect only dominates in a limited region effect. In order to make full use of the above-mentioned four characteristics and avoid invalid calculation in the existing multi-grid method, the present invention provides a multi-grid disturbance domain update acceleration method for numerical simulation of flow around an aircraft. The non-converged convergent element is solved, and the viscosity is only considered in the viscous effect-dominant element, which effectively reduces the calculation amount of a single iteration step and the number of iteration steps required to reach a convergent state, thereby significantly improving the calculation efficiency of the multi-grid method for simulating the flow around the aircraft .

The invention provides a multi-grid disturbance domain update acceleration method for numerical simulation of flow around an aircraft, comprising the following steps:

S1: Read in data, including the fine mesh of the aircraft flow field, the preset computational domain, boundary conditions and calculation settings;

S2: Generate a coarse grid with a specified number of layers through a sparse fine grid;

S3: Initialize the flow field in the preset computational domain on the coarse grid and fine grid according to the incoming flow conditions or according to the given flow field;

S4: According to the flow field initialization method, establish the convective dynamic calculation domain and the viscous dynamic calculation domain;

S5: Divide the data that needs to be stored for solving all layers of grids into two categories; the first type of data is the inherent information of the unit, and the static data structure is used to store the information of all grid units in the preset computational domain, including grid coordinates and Flow field variables; the second type of data is the information related to the solution update. The dynamic data structure is used to store only the information of the grid cells in the convective dynamic calculation domain, including the update amount of the conservation quantity and the local iteration step size;

S6: According to the boundary condition type, assign value to the conservation quantity of the boundary virtual grid of the current layer grid;

S7: The residual term of the flow control equation is divided into two types: inviscid term and viscous term; among them, the inviscid term includes the forcing function of convective flux and coarse grid, and the inviscid term is calculated in the convective dynamic calculation domain; the viscous term terms include viscous flux and turbulence model equation source terms, and the viscous term is calculated in the viscous dynamic computational domain;

S8: In the convective dynamic calculation domain, solve the update quantity of the conserved quantity and update the flow field variables;

S9: determine whether the current grid layer is a fine grid; if yes, go to step S10; if not, go to step S21;

S10: Determine whether the fine mesh reaches the convergence condition; if not, execute step S11; if yes, execute step S23;

S11: Judge whether all boundary elements of the convective dynamic calculation domain have been subjected to inviscidity disturbance one by one; if so, go to step S12; if not, return to step S11 to judge the next boundary element of the convective dynamic calculation domain; After all boundary elements of the domain are dynamically calculated, step S13 is performed;

S12: Measure the propagation direction of the inviscid disturbance, and include the unit located in the propagation direction in the adjacent units of the boundary unit into the convective dynamic calculation domain, and return to step S11 to judge the next boundary unit of the convective dynamic calculation domain;

S13: Judging, one by one, whether all the boundary elements in the convective dynamic computational domain satisfy the four conditions: converged, located in the compressible flow, located in the most upstream, and no longer affected by other elements in the convective dynamic computational domain; if so, go to step S13 S14; if not, return to step S13 to judge the next boundary element of the convective dynamic computational domain; after traversing all boundary elements of the convective dynamic computational domain, execute step S16;

S14: remove the boundary element from the convective dynamic calculation domain, and determine whether the boundary element exists in the viscous dynamic calculation domain; if so, go to step S15; The next boundary unit is judged;

S15: remove the boundary element from the viscous dynamic calculation domain, and return to step S13 to judge the next boundary element of the convective dynamic calculation domain;

S16: Determine whether all the boundary elements of the viscous dynamic calculation domain are dominated by the viscous effect one by one; if so, go to step S17; if not, return to step S16 to judge the next boundary element of the viscous dynamic calculation domain; After calculating all boundary elements of the domain, perform step S18;

S17: Include the unit located in the convective dynamic calculation domain in the adjacent units of the boundary unit into the viscous dynamic calculation domain, and return to step S16 to judge the next boundary element of the viscous dynamic calculation domain;

S18: According to the convective dynamic calculation domain and the viscous dynamic calculation domain of the fine grid, update the convective dynamic calculation domain and the viscous dynamic calculation domain of the coarse grid of all layers;

S19: Reallocate the storage space of the second type of data in all coarse grids and fine grids according to the convective dynamic computational domain;

S20: Interpolate the conserved amount to a coarser grid, return to step S6, and enter the calculation of the coarser grid;

S21: Determine whether a coarser grid needs to be calculated; if yes, go back to step S20; if not, go to step S22;

S22: Interpolate the correction amount of the conservation amount back to a finer grid, and determine whether the finer grid is a fine grid; if so, return to step S6, and enter the calculation of the next iteration step; if not, return to step S6, enter Calculation of finer grids;

S23: output the result.

The specific implementation of the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention will be described in detail below through a specific embodiment.

Embodiment 1: As shown in FIG. 1 , the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention specifically includes the following steps:

The first step: data reading.

The read-in data includes a fine mesh of the aircraft flow field, a preset computational domain, boundary conditions, and computational settings.

Step 2: Generate a coarse mesh.

Generates a coarse mesh with a specified number of layers by sparse fine meshing. For structured grids, the commonly used method for generating coarse grids is to sparse the fine grids twice along all grid directions. Figure 2 illustrates the method for generating coarse grids by taking one dimension as an example. In Figure 2, h represents the fine grid, 2h and 4h represent the first layer and the second layer coarse grid, respectively. That is, each cell of the first layer coarse grid contains 2 fine grid cells, and each cell in the second layer coarse grid contains 2 first layer coarse grid cells.

The third step: flow field initialization.

The flow field initialization is divided into two ways: initialization according to incoming flow conditions and initialization according to a given flow field. When initialized according to inflow conditions, the conserved quantities of all grid cells in the preset computational domain on coarse and fine grids are assigned as inflow values. When initialized according to the given flow field, the conservation quantities of all grid cells in the preset computational domain on the fine grid are assigned to the given flow field value, and the conservation quantities of all grid cells in the pre-set computational domain on the coarse grid Interpolated from the conserved quantities of all grid cells on a fine grid. Taking one dimension as an example, the conservation of coarse grid cells satisfies:

Among them, W _2h represents the conserved quantity of the coarse mesh cell, W _h represents the conservation quantity of the fine mesh cell; |Ω| _2h represents the volume of the coarse mesh cell, and |Ω| _h represents the fine mesh The volume of grid cells; 2I-1 represents the 2I-1th cell along the grid direction, and 2I represents the 2Ith cell along the grid direction.

The fourth step: establish the convective dynamic computational domain and the viscous dynamic computational domain.

(1) When initialized according to the inflow conditions, both the convective dynamic calculation domain and the viscous dynamic calculation domain are established according to the wall boundary. The disturbance carrying flow field change information is generated from the surface of the object, and the units dominated by viscosity are also close to the surface of the object. Therefore, the convective dynamic computational domain takes several layers of adjacent elements on the wall as the initial units. In order to avoid the influence of the growth rate of the convective dynamic computational domain on the convergence rate of the multi-grid method, the first layer of elements next to the wall is taken as the initial element for the coarsest grid, and for the fine grid under the structural grid framework, the convective The dynamic calculation domain needs to take the L-layer adjacent element of the wall as the initial element, and L satisfies:

L=2 ^l (2)

Among them, l represents the number of coarse grid layers, and l is a positive integer, which is given by the calculation settings read in step S1; the initial elements of the coarser grids are generated according to the initial elements of the fine grid; for the viscous dynamic computational domain, all The grid layer is taken as the initial element of the 1-layer element next to the wall surface, so as to avoid the invalid calculation caused by the rapid growth of the viscous dynamic computational domain.

(2) When initialized according to a given flow field, the convective dynamic computational domain and the viscous dynamic computational domain are established according to the flow characteristics of the given flow field. For the convective dynamic computational domain, the initial elements of the fine grid are the disturbed elements in the given flow field, and the flow characteristics of these elements do not match the incoming flow conditions and satisfy:

||WW _∞ ||/||ΔW ⁽¹⁾ || _max ＞ε _a,c (3)

Among them, W _∞ represents the conserved quantity of the incoming flow condition; ||ΔW ⁽¹⁾ || _max represents all units in the convective dynamic computational domain

The maximum value of the update quantity of the conservation quantity of the element in the first iteration step; ε _a,c represents the new threshold of convection, which is 10 ^-4 ≤ε _a,c ≤10 ^-6 .

For the viscous dynamic computational domain, the initial element of the fine grid is the viscous effect dominant element in the fine grid convective dynamic computational domain. The viscous effect measure parameter Ψ is defined so that it represents the ratio of the mass flow rate between viscous disturbances and inviscid disturbances. The viscous effect measurement parameter of the viscous effect-dominant element should be in a relatively large order, so the elements in the viscous dynamic computational domain should satisfy:

表示第1迭代步的粘性效应衡量参数；I,J,K分别表示网格方向；

表示对流通量Jacobian矩阵沿i方向的谱半径，

表示粘性通量Jacobian矩阵沿i方向的谱半径；ε_a,v表示粘性新增阈值，取10^-3≤ε_a,v≤10^-2。Among them, Ψ represents the viscous effect measurement parameter, which is the ratio of viscous disturbance to mass flow without viscous disturbance;

Represents the viscous effect measurement parameter of the first iteration step; I, J, K represent the grid direction respectively;

represents the spectral radius of the Jacobian matrix of convective flux along the i direction,

Represents the spectral radius of the viscous flux Jacobian matrix along the i direction; ε _a,v represents the new threshold of viscosity, which is 10 ^-3 ≤ε _a,v ≤10 ^-2 .

The convective dynamic computational domains of all coarse grid layers are generated according to the convective dynamic computational domains of fine grids, and the viscous dynamic computational domains of all coarse grid layers are generated from the viscous dynamic computational domains of fine grids.

Step 5: Allocate storage space.

The data needed to be stored in all layers of grid solution are divided into two categories; the first type of data is the inherent information of the unit, and the static data structure is used to store the information of all grid units in the preset computational domain, including grid coordinates and flow field. Variables; the second type of data is the information related to the solution update. The dynamic data structure is used to store only the information of the grid cells in the convective dynamic calculation domain, including the update amount of the conserved quantity and the local iterative step size.

For the first type of data, a static data structure is used to store the information of all grid cells in the preset computational domain. For example, let I, J, K denote the grid direction, and I _max , J _max , and K _max respectively denote the number of elements in each grid direction. For a grid of I _max ×J _max ×K _max , if each element has If 10 quantities need to be recorded, a real four-dimensional array A ₁ (10, I _max , J _max , K _max ) can be used for storage.

For the second type of data, the dynamic data structure is used to store only the information of grid cells in the convective dynamic computational domain. For example, let I, J, K denote the grid direction, the cells with the same J, K labels are collectively called a row, and a special data type is defined to store the information of a row in the two types of dynamic calculation domains of convection and viscosity. This special data type Including: a. a single-item list - record the I label range of the row unit in the current iteration step in the update of the two types of dynamic calculation fields; b. an integer two-dimensional array - store the last value in the update of the two types of dynamic calculation fields The range of the I label of the row unit in the iteration step, and the range of the I label of the row unit of the current iteration step is stored after the two types of dynamic calculation domains are updated; c. A real number two-dimensional array - to store the unit located in the convective dynamic calculation domain of the row The second type of data corresponds to the range of the I label of the integer two-dimensional array. For a grid of I _max ×J _max ×K _max , if each unit has 10 quantities to be recorded, a special data type two-dimensional array A ₂ (J _max ,K _max ) is used, and the real number type two in A ₂ is used. The dimensional array is B(10, I ₁ :I ₂ ), I ₁ and I ₂ represent the minimum and maximum values of the I labels of the cells in the row (J, K), respectively, 1<I ₁ , I ₂ <I _max .

In existing multigrid methods, all data for all grid cells in all grid layers are stored in static data structures. In contrast, the multi-grid disturbance domain update acceleration method provided by the present invention divides the data to be stored into two categories: inherent information of all units and information related to solution update; wherein, for the first type of data, the The static data structure stores the information of all grid cells in the preset computing domain; for the second type of data, the dynamic data structure only stores the information of the grid cells in the convective dynamic computing domain of each grid layer, and the fifteenth step will In each iteration step, the storage space of the second type of data is adjusted according to the convective dynamic computational domain of each grid layer in time. Therefore, compared with the existing multi-grid method, the present invention can effectively reduce the storage space requirement.

Step 6: Boundary condition processing.

According to the boundary condition type, assign the value of the conservation quantity of the boundary virtual grid of the current layer grid.

Step 7: Residual estimation.

The flow control equation is expressed as:

为无粘项，

为粘性项；在对流动态计算域中求解残差项的无粘项，在粘性动态计算域中求解残差项的粘性项。where t is the time, F _c is the convective flux, F _v is the viscous flux, Q _F is the forcing function of the coarse grid, Q _T is the source term of the turbulence model equation, N _f is the number of element faces, and ΔS is the is the area of the element face; the right-hand term of equation (6) is collectively referred to as the residual term of the flow control equation,

for non-sticky terms,

is the viscous term; the inviscid term of the residual term is solved in the convective dynamic computational domain, and the viscous term of the residual term is solved in the viscous dynamic computational domain.

Step 8: Time integration.

In the convective dynamic computational domain, the left-hand term of Equation (6) is calculated using the time advance scheme; the time integration method used is consistent with the existing method.

Among the existing multi-grid methods, the two most time-consuming steps of residual estimation and time integration need to be performed on all grid elements of each grid in the static preset computational domain, accounting for about 99% of the total solution time. In contrast, in the multi-grid disturbance domain update acceleration method provided by the present invention, the residual estimation in the seventh step only calculates the inviscid term of the residual in the convective dynamic calculation domain of each layer of grids, and only calculates the inviscid term of the residual in each layer of grids. The viscous term of the residual is calculated in the viscous dynamic calculation domain of ; the time integration of the eighth step is also only performed in the convective dynamic calculation domain of each layer of grids. Therefore, the present invention can simultaneously reduce the calculation amount of each layer of grid, and significantly improve the calculation efficiency of the multi-grid method.

Step 9: Determine whether the current grid layer is a fine grid; if so, perform Step 10; if not, perform Step 17.

Step 10: judge whether the fine mesh meets the convergence condition; if not, execute step 11; if yes, execute step 19.

Step 11: Increase the convective dynamic computational domain.

Judging whether all the boundary elements of the convective dynamic calculation domain have been subjected to inviscid perturbation one by one, that is, to judge whether all the boundary elements of the convective dynamic calculation domain satisfy:

||ΔW||/||ΔW ⁽¹⁾ || _max ＞ε _a,c (7)

Among them, ||ΔW|| represents the modulo value of the update amount of the conservation quantity of the boundary element of the convective dynamic calculation domain in the current iteration step; if the boundary element of the convective dynamic calculation domain satisfies Equation (7), the propagation direction of the inviscid disturbance is measured. , and the unit located in the propagation direction in the adjacent units of the boundary unit is included in the convective dynamic calculation domain, and returns to the eleventh step to judge the next boundary unit of the convective dynamic calculation domain; if the boundary element of the convective dynamic calculation domain is not If formula (7) is satisfied, return to step 11 to judge the next boundary element of the convective dynamic computational domain; after traversing all boundary elements of the convective dynamic computational domain, perform the twelfth step.

The inviscid disturbance propagates at the speed of sound relative to the flow, and along the propagation direction of the inviscid disturbance, the propagation velocity of the inviscid disturbance is positive; let q represent the unit direction vector, then the inviscid disturbance propagates along the q direction as:

u·q+a＞0 (8)

Among them, u is the flow velocity vector, a is the speed of sound; q is the unit vector of the grid center of the boundary element of the convective dynamic calculation domain pointing to a certain grid point of the boundary element. Cells that share this grid are included in the convective dynamic computational domain.

The twelfth step: reduce the convective dynamic computational domain and the viscous dynamic computational domain.

In order to reduce the calculation amount of a single iteration step as much as possible without affecting the convergence rate, the boundary elements of the convective dynamic computational domain can be reduced. The following conditions are judged for all boundary elements of the convective dynamic computational domain one by one to reduce the convective dynamic computational domain:

(1) Judge whether all the boundary elements of the convective dynamic calculation domain have converged one by one, that is, judge whether all the boundary elements of the convective dynamic calculation domain satisfy:

||ΔW||/||ΔW ⁽¹⁾ || _max <ε _d (9)

Wherein, ε _d represents the deletion threshold, and takes ε _d ≤10 ^-7 ; if the boundary element satisfies the formula (9), continue to judge (2); if not, judge the next boundary element of the convective dynamic computational domain;

(2) Determine whether the boundary element is in a compressible flow, that is, determine whether the Mach number of the boundary element is greater than 0.3; if so, continue to judge (3); The next boundary unit is judged;

(3) Judging whether the boundary element is located at the most upstream, that is, judging whether the boundary element satisfies:

Among them, q is taken as the unit vector of the cell center of the boundary element of the convective dynamic calculation domain pointing to the cell center of the adjacent cell, and θd represents the tolerance angle of the upstream element, which is 5° _≤θd≤10 °; if the boundary element satisfies the _formula ( 10), then continue to judge (4); if not, return to judge (1) to judge the next boundary element of the convective dynamic computational domain;

(4) Determine whether the boundary element is no longer affected by other elements in the convective dynamic calculation domain. For supersonic inviscid flow, no judgment is required and it is automatically satisfied; for subsonic flow and viscous flow, determine whether the boundary element satisfies:

||ΔW||+||Δ(ΔW)||<ε _d (11)

Among them, ||Δ(ΔW)|| represents the modulo value of the update quantity of the boundary element conservation quantity affected by the adjacent elements in the convective dynamic calculation domain of the boundary element; Δt represents the iteration step size, C _CFL represents the CFL number of the time-advance format, ΔR ⁱ Represents the influence of the adjacent elements in the convective dynamic computational domain on the residual term of the boundary element along the i direction; for subsonic inviscid elements, the influence of the adjacent elements in the convective dynamic computational domain on the residual term of the boundary element along the i direction ΔR ⁱ represents for:

Among them, ΔF _c represents the change of convective flux; the subscript i+1 represents the adjacent unit of the boundary element of the convective dynamic computing domain along the positive i direction, and the subscript i-1 represents the boundary element of the convective dynamic computing domain along the negative i direction ΔS _i+1/2 represents the area of the boundary element of the convective dynamic computational domain along the positive i direction; ΔS _i-1/2 represents the area of the boundary element of the convective dynamic computational domain along the negative i direction; For a viscous element, the influence ΔR ⁱ of the adjacent element in the convective dynamic computational domain on the residual term of the boundary element along the i direction is expressed as:

If the boundary element satisfies equations (11) and (12), then the boundary element has been converged, located in the compressible flow, located in the most upstream and no longer affected by other elements in the convective dynamic computational domain. If the conditions are met, the boundary element is removed from the convective dynamic calculation domain. If the boundary element exists in the viscous dynamic calculation domain at the same time, the boundary element is removed from the viscous dynamic calculation domain and returns to step 12. Judge the next boundary element of the convective dynamic computational domain; if the boundary element does not satisfy equations (11) and (12), return to the twelfth step to judge the next boundary element of the convective dynamic computational domain. After traversing all boundary elements of the convective dynamic computational domain, perform the thirteenth step.

The thirteenth step: increase the viscous dynamic computational domain.

Determine whether all boundary elements of the viscous dynamic computational domain are dominated by the viscous effect one by one, that is, to judge whether all the boundary elements of the viscous dynamic computational domain satisfy:

Among them, Φ represents the scaling factor of the newly added viscosity threshold, ||ΔW|| _max represents the maximum value of the update amount of the conservation quantity of all elements in the convective dynamic computational domain in the current iteration step, (||ΔW|| _max ) _min represents the th 1 Iterative step to the minimum value of the maximum value of the update value of the single-step conservation quantity in the current iteration step; if the boundary element of the viscous dynamic computational domain satisfies the equations (15) and (16), the adjacent elements of the boundary element are located in the convection. The elements in the dynamic calculation domain are included in the viscous dynamic calculation domain, and return to step 13 to judge the next boundary element of the viscous dynamic calculation domain; if the boundary elements of the viscous dynamic calculation domain do not satisfy equations (15) and (16) , then return to the thirteenth step to judge the next boundary element of the viscous dynamic computational domain; after traversing all the boundary elements of the viscous dynamic computational domain, execute the fourteenth step. The viscous term of the flow governing equation is elliptical. Therefore, all the adjacent elements of the viscous perturbed element can be added to the viscous dynamic computational domain without enlarging the convective dynamic computational domain.

Step 14: Update the coarse grid dynamic computational domain.

According to the convective dynamic computational domain of the fine grid, the convective dynamic computational domain of the coarse grid of all layers is updated, and the viscous dynamic computational domain of the coarse grid of all layers is updated according to the viscous dynamic computational domain of the fine grid. For the structured grid, the corresponding relationship between the coarse grid and the fine grid of the dynamic computational domain is shown in Figure 2. If a fine grid cell is included in the dynamic computational domain of the fine grid, it corresponds to the coarse grid of the cell. Grid cells are also added to the dynamic computational domain of the coarse grid.

In the eleventh step of the present invention, for fine grids, only the disturbed units are added to the convective dynamic calculation domain according to the disturbance propagation characteristics; in the thirteenth step, only the viscous effect dominates The grid elements are added to the viscous dynamic computational domain of the fine grid; in the twelfth step, using the solution convergence feature, without affecting the convergence rate, the converged elements are transferred from the convection and viscous dynamic computational domains in time. In the fourteenth step, according to the dynamic computational domain of the fine grid, the two types of dynamic computational domains of convection and viscosity of all coarse grids are updated. Therefore, the present invention can always keep only the unconverged disturbed element to be solved in all grid layers, and only consider the viscosity in the viscous effect dominant element, thereby effectively avoiding the invalid calculation in the existing multi-grid method.

Step 15: Reallocate the storage space of the second type of data.

According to the convective dynamic computational domain, reallocate the storage space for the second type of data in all coarse and fine grids.

Step 16: Interpolate the conserved quantity to a coarser grid, return to step 6, and enter the calculation of a coarser grid.

The interpolation of the conserved quantities from the fine grid to the coarse grid is the same as the interpolation in the third step for the given flow field initialization.

Step 17: Determine whether a coarser grid needs to be calculated; if a coarser grid needs to be calculated, go back to step 16; if there is no need to calculate a coarser grid, execute step 18;

The commonly used circulation mode of multi-grid is shown in Figure 3. The left subgraph represents a "V" type cycle, and the right subgraph represents a "W" type cycle. In Figure 3, h represents a fine grid, and 2h represents the first layer of coarse mesh. grid, 4h represents layer 2 coarse grid. Taking the "V"-shaped loop as an example, if the current calculation layer is a fine grid, the seventeenth step is to select a coarser grid to be calculated; if the current calculation layer is the coarsest grid, then the seventeenth step is to select no need. Calculates a coarser grid.

Step 18: Interpolate the conservation amount correction back to a finer grid, and determine whether the finer grid is a fine grid; if the finer grid is a fine grid, return to step 6 and enter the next iteration step If the finer mesh is not a fine mesh, go back to step 6 and enter into the calculation of the finer mesh.

Step 19: Output the results.

The inviscid flow around the NACA0012 airfoil with Mach number 0.755 and the turbulent flow around the transonic RAE2822 airfoil with Mach number 0.729 are simulated as follows.

(1) Simulate the inviscid flow around the NACA0012 airfoil with a Mach number of 0.755.

FIG. 4 illustrates the evolution process of the convective dynamic computational domain of the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention. In Fig. 4, η _c represents the ratio of the grid size of the convective dynamic computational domain to the preset computational domain, and N _max represents the total number of iteration steps required to reach a convergent state. The flow field is initialized with flow conditions and disturbances are generated from the wall boundary. At the beginning of the solution, the convective dynamic computational domain gradually expands from the wall to the entire preset computational domain. In the middle and late stages of the solution, the convective dynamic computational domain begins to gradually shrink from the far field to the wall and from the upstream to the downstream as the solution converges, and finally only the near-wall region and its downstream region remain. It can be seen from the change curve of the grid quantity in the convective dynamic calculation domain of the multi-grid disturbance domain update acceleration method shown in FIG. 5a that the multi-grid disturbance domain update acceleration method provided by the present invention is suitable for multi-grids with different grid layers. , the evolution process of the convective dynamic computational domain is basically the same under the condition of different grid layers.

Fig. 5b compares the convergence curve of the multi-grid disturbance domain update acceleration method provided by the present invention and the existing global update multi-grid method. The comparison results show that, for the multi-grid methods with different layers, compared with the existing methods, the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention can reduce the calculation amount of a single iteration step and the total number of iteration steps. Significantly improve the efficiency of numerical simulation. Among them, when the present invention adopts 1 layer of coarse grid, compared with the existing multi-grid method, 52.1% of the calculation time can be saved; when the present invention adopts 2 layers of coarse grid, compared with the existing multi-grid method 45.6% of computing time can be saved.

(2) Simulate the turbulent flow around the transonic RAE2822 airfoil with Mach number 0.729.

FIG. 6 shows the dynamic calculation domain evolution process of the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention. The flow field is initialized from the flow conditions, therefore, the convective dynamic computational domain and the viscous dynamic computational domain are still established based on the wall surface. The growth rate of the convective dynamic computational domain is much greater than that of the viscous dynamic computational domain. It can be seen from the situation of 1% N _max in Figure 6 that when the convective dynamic computational domain has covered the entire preset computational domain, the viscous dynamic computational domain still remains near the wall, including only the near-wall area and its downstream area. 7a, it can be seen from the grid quantity change curves of the two types of dynamic calculation domains, convective and viscous, that the multi-grid disturbance domain update acceleration method provided by the present invention is suitable for viscous flow, convective dynamic Both the computational domain and the viscous dynamic computational domain follow the evolution law of first increasing and then shrinking.

Fig. 7b compares the convergence curves of the multigrid disturbance domain update acceleration method provided by the present invention and the existing global update multigrid method. The comparison results show that the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention is the same as the existing multi-grid method, and a convergent solution can be obtained. Benefiting from the reduction of two types of dynamic computational domains of convection and viscosity and considering only local viscosity, the above-mentioned multi-grid disturbance domain update acceleration method provided by the present invention obtains a significant increase in computational efficiency and can save 40.7% of computational time.

In the multi-grid disturbance domain update acceleration method for the numerical simulation of flow around an aircraft provided by the present invention, the residual estimation only calculates the inviscid term of the residual in the In the dynamic calculation domain, the viscous term of the residual is calculated, and the time integration is only performed in the convective dynamic calculation domain of each layer grid. Therefore, the present invention adopts two types of dynamic calculation domains, convective and viscous, in each layer grid. When the grid method simulates the flow around the aircraft, the calculation amount of each layer of grids significantly improves the calculation efficiency of the multi-grid method for the numerical simulation of the flow around the aircraft, thus effectively solving the problem that the calculation efficiency of numerical simulation is difficult to meet the increasing practical application of aircraft aerodynamic design. question of needs. The data to be stored is divided into two categories: the inherent information of all units and the information related to solving and updating; wherein, for the first type of data, the information of all units in the preset computational domain is stored in a static data structure; for the second type of data data, only the information of grid cells in the convective dynamic computing domain of each grid layer is stored in a dynamic data structure, and the storage space of the second type of data is adjusted according to the convective dynamic computing domain of each grid layer in time at each iteration step. Therefore, , compared with the existing multi-grid method, the invention can effectively reduce the storage space requirement of the numerical simulation of the flow around the aircraft. For the fine grid of the aircraft flow field, only the disturbed elements are added to the convective dynamic calculation domain according to the disturbance propagation characteristics; by weighing the influence of the viscous effect, only the grid elements dominated by the viscous effect are added to the viscous dynamic calculation of the fine grid In the domain; using the solution convergence feature, without affecting the convergence rate, the converged elements are removed from the two types of dynamic calculation domains of convection and viscosity; Convective and viscous two types of dynamic computational domains of the grid; therefore, the present invention can always keep all grid layers in the flow field of the aircraft to solve only the unconverged disturbed element, and only consider the viscosity in the viscous effect dominant element, thus The ineffective calculation of the existing multi-grid method in the numerical simulation of the flow around the aircraft can be effectively avoided, and the problem that the numerical simulation calculation efficiency is difficult to meet the increasing practical application requirements of the aerodynamic design of the aircraft can be solved.

It will be apparent to those skilled in the art that various modifications and variations can be made in the present invention without departing from the spirit and scope of the invention. Thus, provided that these modifications and variations of the present invention fall within the scope of the claims of the present invention and their equivalents, the present invention is also intended to include these modifications and variations.

Claims (7)

1. A multiple grid disturbance domain updating and accelerating method for aircraft bypass flow numerical simulation is characterized by comprising the following steps:

s1: reading data, including fine grids of an aircraft flow field, a preset calculation domain, boundary conditions and calculation settings;

s2: generating a coarse grid with a specified number of layers by thinning the fine grid;

s3: initializing flow fields in preset calculation domains on the coarse grid and the fine grid according to an incoming flow condition or a given flow field;

s4: establishing a convection dynamic calculation domain and a viscosity dynamic calculation domain according to a flow field initialization mode;

s5: dividing data to be stored for solving all layers of grids into two types; the first kind of data is inherent information of a unit, and a static data structure is adopted to store information of all grid units in the preset calculation domain, wherein the information comprises grid coordinates and flow field variables; the second kind of data is information related to solving and updating, and only the information of grid cells in the convection dynamic calculation domain is stored by adopting a dynamic data structure, wherein the information comprises a conservative updating quantity and a local iteration step length;

S6: assigning a value to the conservation quantity of the boundary virtual grid of the current layer grid according to the type of the boundary condition;

s7: dividing residual terms of a flow control equation into two types of non-sticky terms and sticky terms; wherein the inviscid terms comprise forcing functions of the convection flux and the coarse grid, the inviscid terms being computed in the convection dynamics computation domain; the viscosity term comprises a viscosity flux and turbulence model equation source term, the viscosity term being calculated in a viscosity dynamics computation domain;

s8: in the convection dynamic calculation domain, solving a conservative updating quantity and updating a flow field variable;

s9: judging whether the current grid layer is a fine grid or not; if yes, go to step S10; if not, go to step S21;

s10: judging whether the fine grid reaches a convergence condition; if not, go to step S11; if yes, go to step S23;

s11: judging whether all boundary units of the convection dynamic calculation domain are subjected to inviscid disturbance one by one; if yes, go to step S12; if not, returning to the step S11, and judging the next boundary unit of the convection dynamic calculation domain; after traversing all the boundary cells of the convection dynamic calculation domain, executing step S13;

s12: measuring the propagation direction of the inviscid disturbance, bringing the cells positioned in the propagation direction in the adjacent cells of the boundary cell into the convection dynamic calculation domain, returning to the step S11, and judging the next boundary cell in the convection dynamic calculation domain;

S13: judging whether all boundary units of the convection dynamic calculation domain meet four conditions of converged, in compressible flow, at the most upstream and no longer influenced by other units in the convection dynamic calculation domain; if yes, go to step S14; if not, returning to the step S13, and judging the next boundary unit of the convection dynamic calculation domain; after traversing all the boundary cells of the convection dynamic calculation domain, executing step S16;

s14: removing the boundary unit from the convection dynamic calculation domain, and judging whether the boundary unit exists in the viscosity dynamic calculation domain; if yes, go to step S15; if not, returning to the step S13, and judging the next boundary unit of the convection dynamic calculation domain;

s15: removing the boundary cell from the viscous dynamic calculation domain, returning to step S13, and determining the next boundary cell in the convective dynamic calculation domain;

s16: judging whether all boundary units of the viscosity dynamic calculation domain are dominated by viscosity effect one by one; if yes, go to step S17; if not, returning to the step S16, and judging the next boundary unit of the viscosity dynamic calculation domain; after traversing all boundary cells of the sticky dynamic computation domain, executing step S18;

S17: bringing the cell located in the convection dynamic calculation domain in the immediate cell of the boundary cell into the viscous dynamic calculation domain, returning to step S16, and determining the next boundary cell of the viscous dynamic calculation domain;

s18: updating the convection dynamic calculation domain and the viscosity dynamic calculation domain of all layers of the coarse grids according to the convection dynamic calculation domain and the viscosity dynamic calculation domain of the fine grids;

s19: reallocating storage spaces of second-class data in all the coarse grids and the fine grids according to the convection dynamic calculation domain;

s20: interpolating the conservative quantity to a coarser grid, returning to the step S6, and entering the calculation of the coarser grid;

s21: judging whether a coarser grid needs to be calculated; if yes, return to step S20; if not, go to step S22;

s22: interpolating the conservative quantity correction back to the finer grid, and judging whether the finer grid is a fine grid or not; if yes, returning to the step S6, and entering the calculation of the next iteration step; if not, returning to the step S6 to enter the calculation of the finer grid;

s23: and outputting the result.

2. The method for accelerating updating of multiple grid disturbance domains according to claim 1, wherein step S3, initializing flow fields in preset computation domains on the coarse grid and the fine grid according to an incoming flow condition or according to a given flow field, specifically comprises:

Initializing according to incoming flow conditions, and assigning the conservative quantities of all grid units in a preset calculation domain on the coarse grid and the fine grid as incoming flow values;

and initializing according to a given flow field, assigning the conservation quantities of all grid units in a preset calculation domain on the fine grid to be given flow field values, and obtaining the conservation quantities of all grid units in the preset calculation domain on the coarse grid according to the conservation quantity interpolation of all grid units on the fine grid.

3. The method for accelerating updating of multiple grid disturbance domains according to claim 1, wherein step S4 is to establish a convection dynamic calculation domain and a viscosity dynamic calculation domain according to a flow field initialization method, and specifically includes:

initializing according to incoming flow conditions, regarding a convection dynamic calculation domain, taking 1 layer of cells adjacent to a wall surface as initial cells by the coarsest layer of grids, taking L layers of adjacent cells of the wall surface as initial cells by the fine grids, wherein L satisfies the following conditions:

L＝2^l(1)

wherein l represents the number of coarse grid layers, is a positive integer, and is given by the calculation setting read in step S1; generating initial units of the coarse layer grids according to the initial units of the fine grids; for the viscous dynamic calculation domain, all grid layers take 1 layer of units adjacent to the wall surface as initial units;

According to initialization of a given flow field, for a convection dynamic calculation domain, an initial unit of a fine grid is an interfered unit with flow characteristics in the given flow field inconsistent with incoming flow conditions, and the following conditions are met:

||W-W_∞||/||ΔW⁽¹⁾||_max＞_a,c(2)

wherein W represents a conservative amount; w_∞A conservation quantity representing an incoming flow condition; | Δ W⁽¹⁾||_maxRepresenting the maximum value of the conservative updating quantity of all units in the convection dynamic calculation domain in the 1 st iteration step;_a,cthe newly added threshold value of convection is taken as 10^-4≤_a,c≤10^-6(ii) a For the viscous dynamic calculation domain, the initial unit of the fine grid is a viscous effect dominant unit in the fine grid convection dynamic calculation domain, and the following conditions are met:

wherein Ψ represents a viscous effect measurement parameter as a ratio of viscous disturbance to non-viscous disturbance mass flow;

a viscosity effect measurement parameter representing the 1 st iteration; i, J and K respectively represent grid directions;

represents the spectral radius of the convection flux Jacobian matrix along the i direction,

representing the spectral radius of the viscous flux Jacobian matrix along the i direction;_a,vindicating the newly added threshold value of viscosity, take 10^-3≤_a,v≤10^-2(ii) a And the convection dynamic calculation domains of all the coarse grid layers are generated according to the convection dynamic calculation domain of the fine grid, and the viscosity dynamic calculation domains of all the coarse grid layers are generated according to the viscosity dynamic calculation domain of the fine grid.

4. The multiple mesh perturbation domain updating acceleration method according to claim 3 characterized in that step S7, the residual terms of the flow control equation are divided into two types, non-sticky terms and sticky terms; wherein the inviscid terms comprise forcing functions of the convection flux and the coarse grid, the inviscid terms being computed in the convection dynamics computation domain; the viscosity term comprises a viscosity flux and a turbulence model equation source term, and is calculated in a viscosity dynamic calculation domain, and specifically comprises the following steps:

The flow control equation is expressed as:

wherein t represents time, F_cDenotes the convective flux, F_vDenotes viscous flux, Q_FForcing function, Q, representing a coarse grid_TRepresenting the source term of the turbulence model equation, | Ω | representing the volume of the grid cell, N_fThe number of the unit surfaces is shown, and Delta S represents the area of the unit surfaces; the right-hand term of equation (5) is collectively referred to as the residual term of the flow control equation,

is a non-sticky item, and is provided with a sticky part,

is a viscosity term; solving the non-sticky terms of the residual terms in the flow dynamics computation domain and solving the sticky terms of the residual terms in the sticky dynamics computation domain.

5. The method as claimed in claim 4, wherein the step S11 of determining whether all boundary cells of the convective dynamic computation domain have been disturbed without sticking one by one includes:

judging whether all boundary units of the convection dynamic calculation domain meet the following requirements one by one:

||ΔW||/||ΔW⁽¹⁾||_max＞_a,c(6)

wherein, | Δ W | | | represents a modulus of a constant update quantity of a boundary unit of the convection dynamic calculation domain in the current iteration step;

step S12, measuring the propagation direction without viscous disturbance, and incorporating the cell located in the propagation direction in the cells immediately adjacent to the boundary cell into the convection dynamic calculation domain, which specifically includes:

The inviscid disturbance is transmitted at the sound velocity relative to the flow, and the transmission speed of the inviscid disturbance is positive along the transmission direction of the inviscid disturbance; let q denote the unit direction vector, then the propagation of the inviscid perturbation along the q direction is expressed as:

u·q+a＞0 (7)

wherein u represents a flow velocity vector, and a represents a sound velocity; and q is a unit vector of a grid point of the boundary unit of the convection dynamic calculation domain pointing to a certain grid point of the boundary unit, and if q satisfies the formula (7), an adjacent unit sharing the grid point with the boundary unit of the convection dynamic calculation domain is brought into the convection dynamic calculation domain.

6. The multiple-grid-disturbance-domain updating acceleration method according to claim 5, wherein step S13 is implemented to determine whether four conditions of converged, in compressible flow, in most upstream, and no longer influenced by other units in the convective dynamic computation domain are simultaneously satisfied for all boundary cells of the convective dynamic computation domain one by one; if yes, go to step S14; if not, returning to the step S13, and judging the next boundary unit of the convection dynamic calculation domain; the method specifically comprises the following steps:

s131: judging whether all boundary units of the convection dynamic calculation domain meet the following requirements one by one:

||ΔW||/||ΔW⁽¹⁾||_max＜_d(8)

wherein,_dindicating a deletion threshold, take _d≤10^-7(ii) a If yes, go to step S132; if not, returning to the step S131, and judging a next boundary unit of the convection dynamic calculation domain;

s132: judging whether the Mach number of the boundary unit is larger than 0.3; if yes, go to step S133; if not, returning to the step S131, and judging a next boundary unit of the convection dynamic calculation domain;

s133: judging whether the boundary unit satisfies the following conditions:

wherein, q is a unit vector of grid center of the convection dynamic calculation domain boundary unit pointing to the adjacent cell center of the boundary unit, theta_dRepresenting the tolerance angle of the upstream unit, and taking theta more than or equal to 5 degrees_dLess than or equal to 10 degrees; if yes, go to step S134; if not, returning to the step S131, and judging a next boundary unit of the convection dynamic calculation domain;

s134: for supersonic non-viscous flow, step S14 is directly executed without further judgment;

for subsonic and viscous flows, it is determined whether the boundary cell satisfies:

||ΔW||+||Δ(ΔW)||＜_d(10)

wherein, | Δ (Δ W) | | represents a module value of the boundary unit conservative update quantity affected by the adjacent unit in the convection dynamic calculation domain of the boundary unit; Δ t denotes the iteration step, C_CFLNumber of CFLs, Δ R, representing a time advance formatⁱRepresenting the influence of the residual error items of the adjacent units on the boundary units along the i direction in the flow dynamic calculation domain;

For subsonic inviscid cells, the effect Δ R in the i direction of the residual terms of the adjacent cells to the boundary cell in the convection dynamics computation domainⁱExpressed as:

wherein, Δ F_cRepresenting the amount of change in convective flux; the subscript i +1 represents the immediate vicinity of the boundary cell of the convective dynamics computation domain along the positive i direction, the subscript i-1 represents the immediate vicinity of the boundary cell of the convective dynamics computation domain along the negative i direction, Δ S_i+1/2The area of a boundary unit of the convection dynamic calculation domain along a positive i direction unit surface is represented; delta S_i-_1/2The area of a boundary unit of the convection dynamic calculation domain along a negative i direction unit surface is represented;

for sticky cells, the effect Δ R in the i direction of the residual term of the immediate cell on the boundary cell in the flow dynamics computation domainⁱExpressed as:

if yes, go to step S14; if not, the process returns to step S131 to determine the next boundary unit of the convection dynamic calculation domain.

7. The method for accelerating updating of multiple grid disturbance domains according to claim 6, wherein step S16, determining whether all boundary cells of the sticky dynamic computation domain are dominated by a sticky effect one by one, specifically comprises:

judging whether all boundary units of the viscosity dynamic calculation domain meet the following conditions one by one:

Wherein phi represents the scaling factor of the newly added viscosity threshold, | | Δ W | | magnetism_maxRepresenting the maximum value of the conservative update quantity of all units in the convection dynamic calculation domain in the current iteration step (| | Δ W | | count_max)_minRepresents the minimum value of the maximum value of the conservation constant updating quantity of the single step from the 1 st iteration step to the current iteration step.

Images