CN111651831A

CN111651831A - Partition disturbance domain updating calculation method for constant-viscosity compressible streaming of aircraft

Info

Publication number: CN111651831A
Application number: CN202010284225.XA
Authority: CN
Inventors: 蒋崇文; 胡姝瑶; 高振勋; 李椿萱
Original assignee: Beihang University
Current assignee: Beihang University
Priority date: 2020-04-13
Filing date: 2020-04-13
Publication date: 2020-09-11
Anticipated expiration: 2040-04-13
Also published as: CN111651831B

Abstract

The invention discloses a partitioned disturbance domain updating calculation method for a constant-viscosity compressible streaming of an aircraft, and belongs to the field of computational fluid mechanics. The method comprises the steps of reading in a file; initializing a flow field; establishing an initial convection dynamic domain and an initial viscosity dynamic domain; allocating a storage space; assigning a value to the conservation quantity of the boundary virtual grid according to the type of the boundary condition; residual error estimation; time integration; increasing the convection dynamic domain; narrowing the convective, viscous dynamic domain; reallocating storage space; increasing the viscous dynamic domain; repeating the steps until a convergence condition is reached; and outputting the result. According to the method, the convection dynamic domain and the viscosity dynamic domain which can be adaptively increased/reduced are adopted, only disturbed units which are not solved and are not converged are solved, the viscosity effect is considered only in partial disturbed units, the calculation amount and the total iteration step number of a single iteration step can be effectively reduced, and therefore the convergence rate of the constant viscosity compressible streaming of the solving aircraft is improved.

Description

Partition disturbance domain updating calculation method for constant-viscosity compressible streaming of aircraft

Technical Field

The invention belongs to the field of computational fluid mechanics, and particularly relates to a partitioned disturbance domain updating and calculating method suitable for aircraft steady viscous compressible streaming numerical simulation.

Background

Although the solution of the viscous flow control equation N-S equation can be realized by greatly improving the hardness performance of the computer, the numerical simulation of the viscous streaming of the real aircraft shape by adopting the grid with enough resolution is still time-consuming, and the practical application requirements of the disciplines such as fluid mechanics, aircraft design and the like are difficult to meet. The method is characterized in that a time propulsion method is generally adopted in the current computational fluid dynamics to simulate the constant viscous compressible streaming of an aircraft; the method can simultaneously solve subsonic velocity and supersonic velocity flow, and is suitable for large Mach number and Reynolds number range. However, the viscous flow control equation has a large calculation amount, a large amount of viscous flow simulation grids, a slow iterative solution by a time-marching method, and the like, so that it takes a long time to directly solve the complete viscous flow control equation N-S equation by the method, and it is difficult to meet the actual application requirements of the disciplines such as fluid mechanics, aircraft design, and the like. Therefore, an efficient numerical simulation technology for the steady viscous streaming of the aircraft has been a research hotspot of computational fluid mechanics.

Besides the acceleration technology which is universal with the inviscid flow simulation, the acceleration technology principle which is special for the existing viscous flow numerical simulation mainly comprises two main categories of simplifying a mathematical model and calculating in a partitioning mode. The first method is to reduce the amount of computation by replacing the N-S equation with a more solvable simplified equation as the governing equation for the flow. Flow with negligible viscous effects can be governed by Euler equations, all-potential equations, small perturbation equations, etc. The control equation of the laminar flow of the appendage can be simplified into a boundary layer equation, a viscous shock wave layer equation, a parabolic N-S equation and the like. Turbulent flow may introduce RANS, large vortex simulation, etc. mode treatment. The second method is based on the boundary layer theory, divides the non-adhesive and adhesive areas, and realizes the non-adhesive subarea calculation considering the adhesive effect only in the near-wall area and the wake area by combining the control equations, thereby further reducing the calculation amount of the adhesive flow numerical simulation on the basis of the first accelerating method.

Most of the existing inviscid-sticky partition calculation has the following defects: 1) the applicability is limited: the simplified N-S equation has limited applicability, and the simplified form of the N-S equation is established on the premise of more assumptions, for example, the full potential equation is only established for non-swirl and isentropic non-viscous flow, and the boundary layer equation and the parabolic N-S equation are only suitable for appendage viscous flow; 2) the programming is complex to realize: different field variables are solved by different numerical methods for different control equations, and the information exchange among different regions needs special treatment to eliminate the jump among different mathematical models; 3) partitioning difficulty: most of partitioning technologies divide a calculation domain into regular regions according to experience and physical concepts, and the method is not only difficult to define reasonable regions, but also is not suitable for solving dynamic problems with constantly changing ranges; 4) the calculation efficiency is low: the conventional numerical simulation method is used for globally updating all grid units on a static preset calculation domain, and flow field evolution characteristics in the solving process are not considered, so that a large amount of invalid calculations are caused.

To address the above-mentioned shortcomings, the literature, "Hu s., Jiang c., Gao z., Lee c.zonaldistry Region Update Method for step compatible Viscous flow [ J ]. Computer Physics Communications,2019,244: 97-116" proposes a Viscous flow numerical simulation acceleration technique without Viscous-Viscous dynamic partition solution. However, the method for determining the viscosity domain proposed in this document cannot take into account the influence of the wall surface condition and the flow state, and the key parameters for controlling the viscosity domain range cannot be adaptively adjusted. For various aerospace aircrafts with different flight profiles, the key parameters of the aerospace aircrafts need to be adjusted greatly according to experience, and calculation fails once the key parameters are selected improperly.

Disclosure of Invention

Technical problem to be solved by the invention and technical scheme adopted by the technical problem

a. Aiming at the problems of limited applicability and complex programming realization of most inviscid-inviscid partition calculation, the invention controls inviscid, laminar flow and turbulent flow by Euler and N-S, RANS equations respectively. The advantage of these three types of equation combinations is: first, these three types of equations may be applied to numerical simulations of inviscid, laminar and turbulent flow, respectively, without additional constraints; and secondly, the three equations have the same field variable and can be solved by adopting a unified numerical method, so that the information conversion among the regions can be avoided, and the realization difficulty is greatly reduced.

b. Aiming at the problem that the partitioning of most of inviscid-sticky partition calculation is difficult, the invention dynamically and adaptively partitions inviscid domains participating in calculation and sticky domains needing to consider the sticky effect in the flow field grid of the aircraft, and the inviscid domains and the sticky domains can be adaptively updated in each iteration step.

c. Aiming at the problem that the determination method of the viscous domain can not consider the wall surface condition, the flow state influence and the problem that the control key parameter of the viscous domain range can not be adaptively adjusted, the adaptive adjustment of the viscous domain is realized by means of standardizing viscous effect measurement parameters, introducing a convergence rate factor and the like in the determination condition of the viscous domain, the robustness and the adaptability of the algorithm can be further improved while the viscous domain is more consistent with the physical law, and the viscous domain can be simulated by the same parameter for the flow-around problems of different aircrafts.

The invention provides a partitioned disturbance domain updating and calculating method for a constant-viscosity compressible streaming of an aircraft, which comprises the following steps of:

s1: reading in a file;

the read-in file comprises a preset calculation domain, grid coordinates and boundary conditions of the aircraft flow field grid, calculation setting and the like;

s2: initializing an aircraft flow-around flow field according to an incoming flow condition and a given flow field;

s3: establishing an initial convection dynamic domain and an initial viscosity dynamic domain of an aircraft flow field grid;

aiming at a mode of initialization according to the incoming flow condition, an initial convection dynamic domain and an initial viscosity dynamic domain are both established according to a wall boundary;

aiming at the initialization mode according to the given flow field, the initial convection dynamic domain comprises all units with flow characteristics inconsistent with the incoming flow in the given flow field, and the initial viscous dynamic domain comprises a viscous effect dominant unit in the initial convection dynamic domain;

s4: allocating a storage space;

the data required to be stored for solving are divided into two types of data: the first type data is inherent information of a unit, comprises grid coordinates, flow field variables and the like, and is stored in a static data structure in a unit in a preset calculation domain; the second kind of data is information related to solving and updating, and comprises a conservative quantity updating quantity, a local time step length, a solving last step conservative quantity updating quantity and the like, and only the second kind of data of units in the convection dynamic domain is stored by adopting a dynamic data structure;

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

s6: residual error estimation;

dividing a residual error term of the flow control equation into an inviscid term and a viscosity term, solving the inviscid term in a convection dynamic domain, and solving the viscosity term in a viscosity dynamic domain;

s7: time integration;

in the convection dynamic domain, updating the conservation quantity by adopting a time advance format;

s8: increasing the convection dynamic domain by adopting a newly added unit operation, and specifically realizing the increase of the convection dynamic domain through the following 2 sub-steps:

s81: judging whether the convection dynamic domain boundary unit is subjected to inviscid disturbance;

s82: if the convection dynamic domain boundary unit is not disturbed, measuring the propagation direction of the disturbance, and bringing the influenced adjacent unit into the convection dynamic domain;

s9: reducing the convection dynamic domain and the viscosity dynamic domain by adopting a deleting unit operation, and specifically realizing the reduction of the convection dynamic domain through the following 4 judgment conditions:

judgment 1: judging whether the convective dynamic domain boundary unit is converged;

and (3) judging: judging whether the convection dynamic domain boundary unit is positioned in the compressible flow;

and 3, judgment: judging whether the convection dynamic domain boundary unit is positioned at the most upstream;

and 4, judgment: judging whether the boundary unit of the convection dynamic domain is not influenced by other units in the convection dynamic domain any more;

when the convection dynamic domain boundary unit simultaneously meets the 4 judgment conditions, removing the convection dynamic domain boundary unit from the convection dynamic domain; when the convection dynamic domain boundary unit exists in the sticky dynamic domain at the same time, removing the convection dynamic domain boundary unit from the sticky dynamic domain;

s10: reallocating storage space of the second type of data;

s11: increasing the viscous dynamic domain by adopting a newly added unit operation, and specifically realizing the increase of the viscous dynamic domain through the following 2 sub-steps:

s111: judging whether the boundary unit of the viscous dynamic domain is dominated by a viscous effect;

s112: if the viscous dynamic domain boundary unit is dominated by the viscous effect, all the adjacent units of the viscous dynamic domain boundary unit which are positioned in the convection dynamic domain are taken into the viscous dynamic domain;

s12: repeating steps S5-S11 until a convergence condition is reached;

the convergence condition comprises that the solution reaches the specified maximum iteration step number, or the residual error or the updating amount is smaller than a given convergence threshold value;

s13: and outputting the result.

Further, in step S2, when initializing according to the inflow conditions, assigning the conservative values of all the units in the aircraft flow field grid as inflow values; when initializing according to a given flow field, assigning the conservation quantities of all units in the aircraft flow field grid as given flow field values.

Further, the step S3 specifically includes the following steps:

1) for the way of initialization according to incoming flow conditions:

the initial convection dynamic domain comprises a plurality of layers of adjacent units of the wall surface, and the initial viscous dynamic domain comprises 1 layer of units adjacent to the wall surface;

2) for the way of initialization according to a given flow field:

flow characteristics of cells in the initial convective dynamic domain satisfy

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

Wherein W represents a conservative amount; subscript ∞ indicates incoming flow conditions; | Δ W⁽¹⁾||_maxThe maximum value of all unit conservation-constant updating quantities in the step 1 is solved;_a,cthe newly added threshold value of convection is taken as 10^-4≤_a,c≤10^-6；

The initial viscous dynamic domain should contain a viscous effect dominant unit in the initial convective dynamic domain, the effect of the viscous effect is determined according to a standardized viscous effect metric, a viscous effect metric parameter Ψ is defined to represent the ratio of viscous disturbance to non-viscous disturbance mass flow, and when the disturbance speed is represented by a grid scale weighted value of the disturbance speed, the viscous effect metric parameter Ψ is represented as the ratio of the spectral radii of the convective flux and the viscous flux Jacobian matrix, i.e., the ratio of the convective flux to the viscous flux Jacobian matrix

In the formula, I, J and K represent grid directions;

respectively representing the spectral radii of the convection flux and viscous flux Jacobian matrixes along the direction I, I is I, J and K;

measuring parameters by using viscosity effect of adjacent wall surface units in the 1 st iteration step

Performing standardization to standardize viscosity effect measurement parameter

Is shown as

In the formula, subscripts max and min represent the maximum value and the minimum value, respectively;

the normalized viscous-effect metric for the viscous-effect-dominant cell is at a larger magnitude, so the cells in the initial viscous dynamic domain satisfy

In the formula (I), the compound is shown in the specification,_a,vindicating the newly added threshold value of viscosity, take 10^-3≤_a,v≤10^-2。

Further, the step S6 specifically includes the following steps:

the flow control equation in semi-discrete form is expressed as:

wherein W represents a conservative amount; t represents time; F. q represents flux and source terms, respectively; omega, N_fAnd Delta S respectively represents the volume of the grid unit, the number of surfaces and the area of the unit surface; subscripts c, v, T, and axisy denote convection, viscosity, average magnitude of turbulent pulsating field, and axial symmetry, respectively; the right-hand term of the equation in the formula (5) is collectively called the residual of the control equation, and F in the residual is_c、F_c,TAnd Q_axisy,cCollectively referred to as non-sticky terms, F_v、F_v,T、Q_axisy,vAnd Q_TCollectively referred to as sticky terms, and solved in the convection dynamics domain and the sticky dynamics domain, respectively.

Further, the step S8 specifically includes the following steps:

s81, judging whether the convection dynamic domain boundary unit is disturbed without adhesion;

and if the conservative updating amount of the disturbed unit is in a non-negligible magnitude, the convection dynamic domain boundary unit which has been disturbed without adhesion satisfies the following conditions:

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

in the formula, | Δ W | | | represents a module value of a conservation quantity update quantity of the convection dynamic domain boundary unit for solving the current iteration step; | Δ W⁽¹⁾||_maxThe maximum value of all unit conservation-constant updating quantities in the step 1 is solved;_a,cthe newly added threshold value of convection is taken as 10^-4≤_a,c≤10^-6；

S82, if the convection dynamic domain boundary unit is disturbed without adhesion, measuring the propagation direction of the disturbance without adhesion, and bringing the influenced adjacent unit into the convection dynamic domain;

the inviscid disturbance is transmitted at the sound velocity relative to the flow, and along any direction, if the transmission speed of the disturbance is positive, the adjacent units of the convection dynamic domain boundary unit along the direction are also influenced by the disturbance;

let the unit vector q be a unit vector pointing in any direction from the center of the grid of the convection dynamic domain boundary cell, the propagation of the disturbance along the direction of the unit vector q is expressed as:

u·q+a＞0 (7)

wherein u represents a flow velocity vector; a represents the speed of sound and a represents,

the unit vector q is one of the following a-c:

a. the unit vector q is the unit external normal direction of the unit surface of the convection dynamic domain boundary unit;

b. the unit vector q is a unit vector pointing to the grid center of the boundary cell of the convection dynamic domain from the grid center of the boundary cell of the convection dynamic domain to the grid center of the adjacent cell;

c. the unit vector q is a unit vector pointing to the grid point of the grid center of the convection dynamic domain boundary cell.

Further, the 4 determinations in step S9 are specifically as follows:

judgment 1: when the conservative update amount of the convection dynamic domain boundary cell is in a small order, i.e., when the following equation (8) is satisfied,

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

in the formula, | Δ W | | | represents a module value for solving the conservation-constant update quantity of the convection dynamic domain boundary unit of the current iteration step; | Δ W⁽¹⁾||_maxThe maximum value of all unit conservation-constant updating quantities in the 1 st iteration step is solved;_dindicating a deletion threshold, take_d≤10^-7，

Judging that the convective dynamic domain boundary unit has converged;

and (3) judging: according to the compressible flow definition, when the Mach number of the convection dynamic domain boundary unit is larger than 0.3, judging that the convection dynamic domain boundary unit is positioned in the compressible flow;

and 3, judgment: let the unit vector q represent the unit vector pointing to the adjacent cell center of the boundary cell center of the convection dynamic domain, then the adjacent cell satisfies the requirement that the downstream cell of the boundary cell of the convection dynamic domain is the downstream cell

In the formula, u representsA flow velocity vector; theta_dRepresenting the tolerance angle of the upstream unit, and taking theta more than or equal to 5 degrees_d≤10°，

When all the adjacent units in the convection dynamic domain boundary unit in the convection dynamic domain satisfy the formula (9), the convection dynamic domain boundary unit is considered to be positioned at the most upstream;

and 4, judgment: in supersonic non-viscous flow, as the mathematical property of the flow control equation is hyperbolic, that is, any point in the flow field is not influenced by the downstream flow, the supersonic non-viscous unit meeting judgment 3 also meets judgment 4;

in subsonic flow and viscous flow, after considering the influence of the adjacent unit in the convection dynamic domain on the conservation quantity updating quantity of the boundary unit of the convection dynamic domain, the conservation quantity updating quantity of the boundary unit of the convection dynamic domain is still in a smaller magnitude, and the boundary unit of the convection dynamic domain is considered to be not influenced by other units any more, namely, the requirement of meeting the influence of other units is met

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

Where Δ (Δ W) represents the influence of the immediate neighbor cells in the flow dynamics domain on the conservative update amount of the flow dynamics domain boundary cells, and is expressed as:

in the formula, Δ t represents a time step; c_CFLCFL number representing a time advance format; | Ω | represents the volume of the grid cell; i, J, K represent grid directions; Δ RⁱRepresenting the influence of the adjacent units in the flow dynamic domain on the residual error of the flow dynamic domain boundary unit along the direction I, wherein I is I, J and K;

for subsonic inviscid cells, the effect Δ R of the residuals of the neighboring cells in the convective dynamic domain on the convective dynamic domain boundary cells along the grid directionⁱIs shown as

In the formula, Δ W represents a conservative update amount; Δ F_cIndicating convective flux changesThe amount, i.e. the difference between the convection flux of the current step and the previous step; Δ S represents the area of the cell face; subscript i + -1 represents the immediate vicinity of the convective dynamic domain boundary cell along the positive and negative i directions; subscript i + -1/2 denotes the cell faces of the convective dynamic domain boundary cell in the positive and negative i-directions;

representing the spectral radius of the convection flux Jacobian matrix along the i direction;

for viscous units

In the formula (I), the compound is shown in the specification,

representing the spectral radius of the viscous flux Jacobian matrix along the i-direction.

Further, the specific determination process of the sub-step S111 in the step S11 is as follows:

introducing a convergence rate factor to scale the newly added viscosity threshold on the basis of the standardized viscosity effect measurement parameters, and defining a convergence rate factor phi:

in the formula, | Δ W | non-woven phosphor_maxThe maximum value of all unit conservation-constant updating quantities of the current iteration step is solved; (| Δ W | | luminous flux)_max)_minRepresents the minimum value of the maximum value of the single-step conservation quantity updating quantity in the 1 st to the current iteration steps, when the solution is converged, i.e., | | delta W | luminance_max＝(||ΔW||_max)_minPhi is 1.0; when the solution divergence or convergence rate is oscillating, i.e., | | Δ W | | non-woven cells_max>(||ΔW||_max)_min，Φ＜1.0；

The viscous dynamic domain boundary unit is mainly satisfied by the viscous effect

When the solution converges, i.e. | | Δ W | | non-woven calculation_max＝(||ΔW||_max)_minThe new threshold of actual viscosity is_a,v(ii) a When the solution divergence or convergence rate is oscillating, i.e., | | Δ W | | non-woven cells_max>(||ΔW||_max)_minThe actual viscosity increase threshold used is less than_a,vTherefore, the problem of calculation failure caused by insufficient viscous dynamic domain can be avoided by increasing the viscous dynamic domain.

The invention has the advantages of

(1) The partition disturbance domain updating calculation method provided by the invention can obviously improve the convergence rate.

In the conventional global updating method, the most time-consuming residual estimation and time integration are performed on all grid cells in a static preset calculation domain, and the total solution time is about 99%. In contrast, in the partitioned disturbance domain update calculation method of the present invention, the residual estimation of step S6 calculates only the non-sticky term of the residual in the convection dynamic domain, and calculates only the sticky term of the residual in the sticky dynamic domain; the time integration of step S7 is also performed only in the convective dynamic domain. The convection dynamic domain only contains disturbed units which are not converged in the preset calculation domain, and the viscous dynamic domain only contains viscous effect dominant regions in the convection dynamic domain. Therefore, the partitioned disturbance domain updating calculation method can effectively improve the calculation efficiency of the constant viscous compressible streaming numerical simulation of the aircraft.

(2) The partition disturbance domain updating calculation method provided by the invention can save storage space.

In the conventional global update method, all data of all grid cells in a static preset calculation domain are stored in a static data structure. In contrast, the partition perturbation domain updating calculation method of the invention divides the data to be stored into two types of intrinsic information of units and information related to solving and updating; for the first type of data, storing information of units in a preset calculation domain in a static data structure; for the second type of data, only information for cells in the stream dynamic domain is stored in a dynamic data structure. In step S10, the storage space of the second type data is adjusted according to the convection dynamic domain in each iteration. Therefore, the partitioned disturbance domain updating calculation method can effectively reduce the memory requirement of the aircraft constant viscosity compressible streaming numerical simulation.

(3) The partition disturbance domain updating calculation method provided by the invention can adaptively adjust the viscous domain, distinguish the influence of the flow state and the wall surface condition and determine the viscous effect leading area which accords with the physics.

In the standardized viscous effect measurement parameters defined by the formula (3), the viscous effect measurement parameters defined by the formula (2) are standardized by taking the viscous effect measurement parameters of the adjacent wall surface units in the step 1, so that the same parameter value in the partitioned disturbance domain updating calculation method has good applicability to flow problems of large Mach number and Reynolds number, and the problem that algorithm parameters need to be adjusted according to experience for different aircraft working conditions is effectively avoided.

Taking the normalized viscosity effect measurement parameter defined by the formula (3)

The condition of laminar flow is standardized,

the turbulent flow condition is standardized, so that the viscous effect dominant region under the turbulent flow condition determined by the partition disturbance domain updating calculation method is thicker than the viscous effect dominant region under the laminar flow condition under the same condition, and the physical law that the turbulent flow boundary layer grows faster than the laminar flow boundary layer along the flow direction under the same condition is met.

In step S11, the increase of the sticky dynamic domain takes into account the convergence rate factor defined by equation (14), so that the partition perturbation domain update calculation method can resolve the influence of the wall surface characteristics, and can implement adaptive adjustment of the sticky domain to avoid the calculation failure problem caused by insufficient sticky domain.

(4) The partition disturbance domain updating calculation method provided by the invention can effectively avoid invalid calculation.

The newly added unit operation of the partition disturbance domain updating calculation method of the invention is that in steps S8 and S11, only the disturbed unit and the viscous effect leading unit are respectively added into the convection dynamic domain and the viscous dynamic domain according to the disturbance propagation characteristics. In step S9, the delete unit operation may utilize the solution convergence feature to remove the converged units from the solution update in time without affecting the convergence rate. Therefore, the partitioned disturbance domain updating calculation method can always solve only the disturbed units with unconverged solution and consider the viscosity only in the viscous effect leading unit, so that the invalid calculation in the constant viscous compressible flow surrounding solution of the traditional aircraft is effectively avoided.

Drawings

FIG. 1 is a flow chart of a method for calculating a zonal disturbance domain update of a constant viscous compressible flow around an aircraft according to the present invention;

fig. 2 is a schematic diagram of a dynamic domain evolution process of a partitioned perturbation domain update calculation method according to embodiment 1 of the present invention;

fig. 3 is an evolution process diagram of a flow field and a dynamic domain for solving a supersonic velocity plate problem by the partitioned disturbance domain updating calculation method in embodiment 2 of the present invention;

fig. 4 is an evolution process diagram of a flow and dynamic domain for solving a subsonic speed NACA0012 airfoil laminar flow streaming problem by the partition disturbance domain update calculation method in embodiment 3 of the present invention;

fig. 5a) and b) are dynamic domain grid quantity change curves and convergence curves for solving the problems of supersonic cone laminar flow and turbulent flow streaming by the partition disturbance domain updating calculation method in embodiment 4 of the present invention, respectively.

Detailed Description

The flow field evolution in the aircraft steady viscous compressible flow-around numerical simulation process has the following characteristics:

1) the disturbance of the flow field begins from the place that the constant control equation can not meet;

2) the disturbance carrying the intermittent information is gradually transmitted to the surrounding flow field at a limited speed along with time;

3) the original state of the undisturbed area is maintained

4) In compressible flow, the convergence of the upstream flow field is not later than that of the downstream flow field;

5) whether or not detached, the viscous effect only dominates over a limited area.

In order to fully utilize the 5-point characteristics and avoid invalid calculation in numerical simulation of the constant viscous compressible streaming of the traditional aircraft, the invention provides a partitioned disturbance domain updating calculation method for the constant viscous compressible streaming of the aircraft. According to the method, the convection dynamic domain and the viscosity dynamic domain which can be adaptively increased/reduced are adopted, only disturbed units which are not solved and are not converged are solved, the viscosity effect is considered only in part of disturbed units, the calculation amount and the total iteration step number of a single iteration step can be effectively reduced, and therefore the solving speed of the steady viscous compressible streaming of the aircraft is improved.

As shown in fig. 1, the method for updating and calculating the partitioned disturbance domain of the aircraft steady-state viscous compressible streaming specifically includes the following steps:

s1: and (6) reading in a file.

The read-in files comprise preset calculation domains, grid coordinates and boundary conditions of the aircraft flow field grid, calculation settings and the like.

S2: and initializing a flow field.

The initialization of the flow field is divided into two modes of initialization according to the inflow condition and initialization of the given flow field. When initializing according to an incoming flow condition, assigning the conservative quantities of all units in the aircraft flow field grid as incoming flow values; when initializing according to a given flow field, assigning the conservation quantities of all units in the aircraft flow field grid as given flow field values.

S3: and establishing an initial convection dynamic domain and an initial viscosity dynamic domain of the aircraft flow field grid.

a. When initialized according to incoming flow conditions, perturbations carrying flow field variation information are generated from the surface of the object, with the unit dominated by viscosity being in close proximity to the surface of the object. Therefore, a plurality of layers of adjacent units, such as 5-10 layers, of the wall surface are taken in the initial convection dynamic domain; the initial viscous dynamic domain takes 1 layer of cells next to the wall.

b. When initialized according to a given initial flow field, two types of dynamic domains, convection and viscous, are established according to the flow characteristics of the given initial flow field.

For the initial convective dynamic domain, its cells should be disturbed cells in a given initial flow field, whose flow characteristics should not be consistent with the incoming flow conditions. Therefore, the flow characteristics of the cells in the initial convective dynamic domain should be satisfied

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

Wherein W represents a conservative amount; subscript ∞ indicates incoming flow conditions;_a,cthe newly added threshold value of convection is taken as 10^-4≤_a,c≤10^-6(ii) a In order to avoid the influence of the grid, the numerical format and other factors on the initial convection dynamic domain range, the maximum value of the all-unit conservative updating quantity in the step 1 is solved by the formula (1) | Δ W⁽¹⁾||_maxNormalization is performed.

The initial viscous dynamic domain should contain a viscous effect dominant cell in the initial convective dynamic domain, the effect of the viscous effect being determined from a normalized viscous effect metric defining a viscous effect metric Ψ representing the ratio of viscous disturbance to non-viscous disturbance mass flow. When the disturbance velocity is expressed as a grid scale weighted value of the disturbance velocity, the viscous effect metric Ψ is expressed as the ratio of the spectral radii of the convective flux and the viscous flux Jacobian matrix, i.e., the ratio

In the formula, I, J and K represent grid directions;

in order to reduce the influence of factors such as flow characteristics, numerical formats, grids and the like on the viscous dynamic domain, parameters are measured by the viscous effect of the adjacent wall surface units in the step 1 of solving

Equation (2) is normalized. Considering that turbulent boundary layer will be more laminar than laminar boundary layer under the same conditionThe boundary layer grows faster in the flow direction, taking into account the influence of the flow state on the viscous dynamic domain

The maximum value of (c) is normalized to the laminar flow situation,

is normalized to the turbulence, the viscosity effect metric is normalized

Can be expressed as

In the formula, subscripts max and min represent the maximum value and the minimum value, respectively.

The normalized viscous-effect metric for the viscous-effect-dominated cell should be at a large magnitude, so the cells in the initial viscous dynamic domain should satisfy

S4: allocating a storage space;

the data required to be stored for solving are divided into two types of data: the first type is the inherent information of the unit, and the static data structure is adopted to store the information of the unit in the preset calculation domain, such as grid coordinates, flow field variables and the like; the second type is information related to solving and updating, and only information of units in the flow dynamic domain is stored by adopting a dynamic data structure, such as the conservative updating quantity, the local time step length, the conservative quantity of the last iteration step to be solved, and the like.

For the first type of data, the partition disturbance domain updating calculation method adopts a static array to store the information of the units in the preset calculation domain. For example, let I, J, K denote the grid direction, I_max、J_max、K_maxEach representing the number of cells in each grid direction. For a block I_max×J_max×K_maxIf there are 10 quantities to be recorded per cell, a real four-dimensional array A can be used₁(10,I_max,J_max,K_max) And storing.

For the second type of data, the partition disturbance domain updating calculation method adopts a dynamic data structure, and only stores the information of the units in the flow dynamic domain. For example, let I, J, K denote the grid direction, and cells with the same J, K designation are collectively referred to as a row. A special data type is defined for storing information for a row in the dynamic domain. The special data types include: a. a single linked list-records the I label range of the line unit in the current step in the dynamic domain updating; b. an integral two-dimensional array, wherein the I label range of the unit in the previous step is stored when the dynamic domain is updated, and the I label range of the unit in the current step is stored after the dynamic domain is updated; c. a real two-dimensional array, storing the second type of data for the row, corresponding to the range of the I-label of the integer two-dimensional array. For a block I_max×J_max×K_maxIf there are 10 quantities to record per cell, a special data type two-dimensional array A is used₂(J_max,K_max)，A₂The medium real number type two-dimensional array is B (10, I)₁:I₂)(I₁、I₂Respectively represent the minimum and maximum values of the I indices of the units in the row (J, K), 1 ≦ I₁,I₂≤I_max)。

S5: and assigning a value to the conservation quantity of the boundary virtual grid according to the type of the boundary condition.

S6: and estimating a residual error.

In order to effectively avoid invalid calculation in the conventional global updating method and reduce the calculation amount of a single iteration step, the invention divides the residual error of the control equation into a non-sticky item and a sticky item, and only solves the non-sticky item of the residual error in the convection dynamic domain and the sticky item of the residual error in the sticky dynamic domain.

The flow control equation in semi-discrete form can be expressed as:

wherein t represents time; F. q represents flux and source terms, respectively; omega, N_fAnd Delta S respectively represents the volume of the grid unit, the number of surfaces and the area of the unit surface; subscripts c, v, T, and axisy denote convection, viscosity, average magnitude of turbulent pulsating field, and axial symmetry, respectively; the right-hand term of the equation in the formula (5) is collectively called the residual of the control equation, and F in the residual is_c、F_c,TAnd Q_axisy,cCollectively referred to as non-sticky terms, F_v、F_v,T、Q_axisy,vAnd Q_TCollectively called sticky terms, and respectively solved in a convection dynamic domain and a sticky dynamic domain, and the adopted calculation method is consistent with a conventional global updating method.

S7: and (4) time integration.

In the convection dynamic domain, the left term of equation (5) is calculated by using a time-marching format, and the adopted time integration method is completely consistent with the conventional global updating method.

S8: the convection dynamic domain is increased.

The unit (I, J, K) represents the boundary unit of the convection dynamic domain, and the newly added unit operation in the partition disturbance domain updating calculation method realizes the increase of the convection dynamic domain through the following 2 sub-steps:

s81, judging whether the unit (I, J, K) is disturbed without adhesion;

the conservative updating amount of the disturbed unit is in a non-negligible order, and the disturbed unit satisfies

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

In the formula, | Δ W | | | represents a module value for solving the conservation-constant update quantity of the convection dynamic domain boundary unit of the current iteration step; | Δ W⁽¹⁾||_maxThe maximum value of all unit conservation-constant updating quantities in the step 1 is solved;_a,cthe newly added threshold value of convection is taken as 10^-4≤_a,c≤10^-6；

S82, if the unit (I, J, K) is not disturbed, measuring the propagation direction of the disturbance, and including the affected adjacent unit of the unit (I, J, K) in the convection dynamic domain.

The inviscid disturbance propagates at the speed of sound relative to the flow, and in either direction, if the propagation speed of the disturbance is positive, the immediately adjacent cells in that direction of cell (I, J, K) will also be affected by the disturbance. Let q be the unit vector where the cell center of the cell (I, J, K) points in either direction, the propagation of the perturbation along the direction of the unit vector q can be expressed as:

u·q+a＞0 (7)

wherein u represents a flow velocity vector; a represents the speed of sound.

Alternatives to the unit vector q include:

a. the unit vector q is the unit out-of-unit normal of the cell plane of the cell (I, J, K);

taking the I forward direction as an example, if the outer normal direction n of the interface of the cell (I, J, K) and the cell (I +1, J, K)_I+1/2If equation (7) is satisfied, the unit (I +1, J, K) is incorporated into the convection dynamics domain.

b. The unit vector q is the unit vector in which the cell center of the cell (I, J, K) points to its immediate cell center;

taking the I forward direction as an example, if the unit vector of the cell (I, J, K) pointing to the cell (I +1, J, K) satisfies equation (7), the cell (I +1, J, K) is included in the convective dynamic domain.

c. The unit vector q is a unit vector in which the cell center of the unit (I, J, K) points to the cell point;

if the unit vector of cell (I, J, K) whose center points to a certain lattice point satisfies equation (7), the immediate neighboring cell sharing that lattice point with cell (I, J, K) is included in the convective dynamic domain.

S9: reducing the convection, viscous dynamic domain.

In order to reduce the calculation amount of a single iteration step as much as possible on the premise of not influencing the convergence rate, for the boundary units (I, J, K) of the convection dynamic domain, the deletion unit operation in the partition disturbance domain updating calculation method of the invention realizes the reduction of the convection dynamic domain through the judgment conditions of the following 4 sub-steps:

s91, judging whether the unit (I, J, K) is converged;

the conservative update amount of the converged unit should be in a small order, i.e., satisfy

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

In the formula, | Δ W | | | represents a module value for solving the conservation-constant update quantity of the convection dynamic domain boundary unit of the current iteration step; | Δ W⁽¹⁾||_maxThe maximum value of all unit conservation-constant updating quantities in the step 1 is solved;_dindicating a deletion threshold, take_d≤10^-7。

S92: determining whether the unit (I, J, K) is in a compressible flow;

the mach number of the cell in compressible flow should be greater than 0.3, according to compressible flow definition.

S93: determining whether the unit (I, J, K) is located most upstream;

let q denote the unit vector in which the cell center of a cell (I, J, K) points to its immediate cell center, then the immediate cell should satisfy that the immediate cell is the downstream cell of the cell (I, J, K)

Wherein u represents a flow velocity vector; 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; a unit (I, J, K) is considered to be most upstream if all immediately adjacent units in the convective dynamic domain satisfy equation (9).

S94: judging whether the unit (I, J, K) is no longer influenced by other units in the convection dynamic domain;

in the supersonic non-viscous flow, the mathematical property of the flow control equation is hyperbolic, that is, any point in the flow field is not affected by the downstream flow, so the supersonic non-viscous unit meeting the step S93 naturally meets the step S94.

In subsonic flow and viscous flow, if the influence of the adjacent units on the conservative updating quantities of the units (I, J, K) is considered, and the conservative updating quantities of the units (I, J, K) still meet the convergence condition, the units (I, J, K) are considered to be not influenced by other units any more, namely meet the requirement of meeting the requirement

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

Where Δ (Δ W) represents the effect of the immediate neighbor cells in the flow dynamics domain on the amount of cell (I, J, K) conservation updates, which can be expressed as:

in the formula, Δ t represents a local time step; c_CFLCFL number representing a time advance format; | Ω | represents the volume of the grid cell; i, J, K represent grid directions; Δ RⁱRepresenting the influence of the adjacent units in the flow dynamic domain on the residual error of the flow dynamic domain boundary unit along the direction I, wherein I is I, J and K;

for subsonic inviscid cells, the effect Δ R of the neighboring cells on the residual of convection dynamic domain boundary cells along the grid directionⁱCan be expressed as

In the formula, Δ W represents a conservative update amount; Δ F_cRepresenting the convection flux variation, namely the difference between the current step and the previous step; subscript i + -1 represents the immediate vicinity of the convective dynamic domain boundary cell along the positive and negative i directions; subscript i + -1/2 denotes the cell faces of the convective dynamic domain boundary cell in the positive and negative i-directions;

representing the spectral radius of the convective flux Jacobian matrix along the i direction.

For viscous units

In the formula (I), the compound is shown in the specification,

If the unit (I, J, K) satisfies the above 4 judgment conditions at the same time, the delete unit operation will remove it from the convection dynamic domain. If the unit (I, J, K) is present in the sticky dynamic domain at the same time, then the delete unit operation will also remove it from the sticky dynamic domain in order to ensure that the sticky dynamic domain is always included in the convective dynamic domain.

S10: the second type of data storage space is reallocated.

And for the line (J, K) in the stream dynamic domain, if the I label range of the line (J, K) is changed, the storage space of the line for storing the second type data is redistributed according to the updated I label range.

S11: increasing the viscous dynamic domain.

Let the cell (I, J, K) represent a boundary cell of the sticky dynamic domain, and the newly added cell operation implements the augmentation of the sticky dynamic domain by 2 sub-steps:

s111, judging whether the unit (I, J, K) is dominated by a viscosity effect;

if the normalized sticking effect measure of a cell (I, J, K) is of non-negligible magnitude, the cell must take into account the sticking effect. In order to effectively distinguish the influence of wall surface characteristics and realize the self-adaptive adjustment of a viscous domain to avoid the problem of calculation failure caused by insufficient viscous domain, a convergence rate factor is introduced to scale the newly added viscous threshold on the basis of standardized viscous effect measurement parameters, so that the self-adaptive capacity of the partition disturbance domain updating calculation method is improved through the change of the convergence rate in the solution.

Defining a convergence speed factor phi

In the formula, | Δ W | non-woven phosphor_maxThe maximum value of all unit conservation-constant updating quantities of the current iteration step is solved; (| Δ W | | luminous flux)_max)_minRepresents the minimum value of the maximum value of the single-step conservation constant updating quantity in the 1 st to the current iteration steps, when the solution is converged, i.e., | | delta W | counting_max＝(||ΔW||_max)_minPhi is 1.0; when the solution divergence or convergence rate is oscillating, i.e., | | Δ W | | non-woven cells_max>(||ΔW||_max)_minPhi is less than 1.0. Taking into account the fact that the values are formatted in different waysThe same problem may have obvious convergence rate difference, and in order to avoid the viscous dynamic domain from being too sensitive to the numerical format and the convergence process, the formula (14) adds | | Δ W⁽¹⁾||_maxThe term is measured as the logarithm of the maximum of the conservative update quantity.

The units (I, J, K) are subject to a viscous effect which is to be satisfied

When the solution converges, i.e. | | Δ W | | non-woven calculation_max＝(||ΔW||_max)_minThe new threshold of actual viscosity is_a,v(ii) a When the solution divergence or convergence rate is oscillating, i.e., | | Δ W | | non-woven cells_max>(||ΔW||_max)_minThe actual viscosity newly-increased threshold value adopted is smaller than_a,vTherefore, the problem of calculation failure caused by insufficient viscous dynamic domain can be avoided by increasing the viscous dynamic domain.

S112: if the unit (I, J, K) is dominated by the viscous effect, all its immediate neighbors in the convective dynamic domain are included in the viscous dynamic domain.

Because the viscous term of the control equation is elliptical, the operation of the newly added unit adds all the immediately adjacent units of the units which are subjected to viscous disturbance into the viscous dynamic domain on the premise of not increasing the convection dynamic domain.

S12: steps S5 to S11 are repeated until a convergence condition is reached.

S13: and outputting the result.

The invention is further described below with reference to the accompanying drawings and examples, it being understood that the examples described below are intended to facilitate the understanding of the invention, and are not intended to limit it in any way.

Example 1

Simulating the problem of a laminar flow with a Mach number of 2.5 flowing through a wedge angle of 5 degrees, and FIG. 2 shows the evolution process of a convection dynamic domain and a viscous dynamic domain of a partition disturbance domain updating calculation method under two operations of adding and deleting units, wherein η_c、η_vRespectively representing convection, viscosity dynamic domain and preset calculation domainRatio of the grid quantities of (A), N_maxRepresenting the number of iteration steps required to solve for convergence. In this embodiment, the flow field is initialized with incoming flow conditions and disturbances are generated from wall boundaries. In step 1 of the solution, the initial convective dynamic domain contains 10 layers of adjacent cells of the wall and the initial viscous dynamic domain contains the immediate cells of the wall. In the early stage of the solution, the two types of dynamic domains are gradually expanded by adding unit operation. At step 330, as shown in fig. 2b), the convection and viscous dynamic domains reach the maximum range at the same time, the maximum convection dynamic domain covers the oblique shock wave and the entire wave-back region thereof, and the maximum viscous dynamic domain is elongated and covers only the thin region of the near-wall region. Subsequently, as shown in fig. 2c) and d), the two types of dynamic domains shrink simultaneously from upstream. Since supersonic inviscid flow only needs to consider the influence of its upstream flow field, the first removed is the unit in the convective dynamic domain that does not consider the viscous effect. At 48% N_maxThen, as shown in fig. 2e), only the viscous dominant cell remains in the convection dynamic domain, i.e. the convection and viscous dynamic domains have completely coincided. At the end of the solution, only the cells near the wall and with lower speed remain in the dynamic domain, as shown in fig. 2 f).

The dynamic domain grid quantity variation curve shown in fig. 2g) shows that, during the whole solving process, both dynamic domains are increased and then reduced, η_cAnd η_vThe peak values of (a) are only 0.657 and 0.519, that is, only 66.2% of cells in the preset calculation domain participate in the solution at the same time, and only 51.9% of cells take the viscosity effect into account at the most. The convergence curve of the prior art global updating method compared with the partitioned perturbation domain updating calculation method of the invention in fig. 2 proves that, except for the initial stage of the solution, the convergence rate of the invention in each iteration step is consistent with that of the global updating method, and the contraction of the calculation domain does not cause adverse effect on the solution convergence.

Example 2

For the case of two wall surfaces of an insulating wall and a cold wall, the Mach number is 3.5, and the Reynolds number is 5.83 × 10⁶The problem of laminar flow through a flat plate is simulated, and fig. 3 shows that under two wall conditions, the partition disturbance domain updating calculation method of the invention is_a,c＝10^-5、_a,v＝10^-2、_d＝10^-7In both wall cases, the convective dynamic domain and the viscous dynamic domain are both established according to the wall boundary and gradually expand from the lower boundary of the pre-set computational domain to the upper boundary in the early stage of the solution, as shown in line 1 of FIG. 3, line 2 of FIG. 3 compares the maximum convective dynamic domains for both wall cases, since the thickness of the boundary layer is smaller and the intensity of the shock wave of the generated leading edge is smaller for the cold wall case than for the thermal insulation wall case, η for the cold wall case_cThe peak is slightly smaller. Line 3 of fig. 3 compares the maximum viscous dynamic range for both wall cases. The maximum viscous dynamic range for both wall cases covers only a very thin area near the wall. The newly added unit operation defines a thicker viscous dynamic domain for the case of an insulating wall with a thicker boundary layer than for the case of a cold wall, the viscous dynamic domain having a maximum thickness of only 65% of the case of an insulating wall in the case of a cold wall. In both wall cases, the maximum thickness of the viscous dynamic domain is about 1.55 times its nominal thermal boundary layer thickness. Therefore, the partition disturbance domain updating calculation method has the capability of adaptively distinguishing the wall influence, and can define the viscous effect dominant region conforming to the physical law.

In particular, the zone disturbance domain update calculation method of the present invention achieves a time savings of 50.1% and 44.8% and an iteration step reduction of 5.51% and 3.66% for the adiabatic wall and cold wall cases, respectively, as compared to the global update method. Therefore, by reducing the calculation amount of a single iteration step and the total iteration step number, the partition perturbation domain updating calculation method can remarkably improve the numerical simulation efficiency of the steady viscous compressible flow.

Example 3

For NACA0012 airfoil, the incoming flow Mach number is 0.5 and the incoming flow Reynolds number is 10⁴And the simulation is carried out under the attack angle of 0 degree, and figure 4 shows that the partition disturbance domain updating calculation method adopts two time advancing formats of 3-order R-K explicit and LU-SGS implicit_a,c＝10^-5、_a,v＝10^-2、_d＝10^-7Dynamic domain evolution process of time. During the pre-solution period, the convection and viscous dynamics domains grow together from the wall boundaries as the disturbance propagates, as is the case for the 3 rd order R-K format at step 173.In the subsonic case, the convective dynamic domain will grow to the entire predefined computational domain, while the maximum viscous dynamic domain still contains only a small portion of the predefined computational domain, η_vThe peak versus apparent and implicit formats are 0.692 and 0.673, respectively. The lower the reynolds number, the thicker the boundary layer. Thus, the viscous dynamic domain of example 3 covers a thicker near-wall region than example 2. As time progresses, the convective dynamic domain gradually shrinks towards the airfoil, removing the converging elements located in the far field from the solution update, as shown in fig. 4, column 3. Eventually, only the low speed cells near the walls and wake remain in the dynamic domain, as shown in the enlarged view of FIG. 4.

Compared with a global updating method, the partitioned perturbation domain updating method provided by the invention respectively obtains time saving amounts of 42.8% and 41.7% for the explicit and implicit conditions. Wherein the time consumption for residual estimation and time integration is reduced by 43% and 39% in both formats. Therefore, the same value of the parameter in the partitioned disturbance domain updating calculation method has good applicability to the flow problems of large Mach number and Reynolds number, and the problem that the algorithm parameter needs to be adjusted according to experience aiming at different problems can be effectively avoided.

Example 4

For Mach number 6 and Reynolds number 1.05 × 10⁷The problem that laminar flow and turbulent flow pass through a 5-degree half cone angle cone is simulated, and the updating calculation method of the subarea disturbance domain is respectively shown in figures 5a) to b)_a,c＝10^-5、_a,v＝10^-2、_d＝10^-7The dynamic domain grid quantity change curve and the convergence curve are obtained, the area near the cone top, which is the area difficult to converge in the whole flow field due to the factors of dense grids, large flow gradient and the like, and the convergence condition of the whole flow field can be immediately met when the area converges_cThe N curves both exhibit a gradual followed by a sharp decline η comparing the two flow conditions_vIt is known that the maximum viscous dynamic range determined by the operation of the newly added unit in the turbulent flow situation is obviously thicker than that in the laminar flow situation, which is consistent with the physical law that the turbulent boundary layer grows faster in the flow direction than the laminar boundary layer. Thus, the present invention partition perturbation domain updateThe calculation method has the capability of adaptively distinguishing the influence of the flow state on the viscous effect dominant region.

Compared with the global updating method, the partition disturbance domain updating calculation method provided by the invention respectively obtains 32.2% and 25.5% of time saving amount, 4.93% and 5.76% of maximum memory requirement saving amount and 11.7% and 9.5% of total iteration step saving amount for laminar flow and turbulent flow conditions. Therefore, for the problem of supersonic incoming flow, the partition disturbance domain updating calculation method provided by the invention can be used for remarkably improving the numerical simulation efficiency of the constant viscosity compressible flow and reducing the storage requirement.

It will be apparent to those skilled in the art that various modifications and improvements can be made to the embodiments of the present invention without departing from the inventive concept thereof, and these modifications and improvements are intended to be within the scope of the invention.

Claims

1. A partitioned disturbance domain updating calculation method for a constant-viscosity compressible streaming of an aircraft is characterized by comprising the following steps:

s1: reading in a file;

s4: allocating a storage space;

the data required to be stored for solving are divided into two types of data: the first type data is inherent information of the unit, comprises grid coordinates and flow field variables, and is stored in a static data structure in a unit in a preset calculation domain; the second kind of data is information related to solving and updating, and comprises a conservative quantity updating quantity, a local time step length and a conservative quantity updating quantity in the last step of solving, and only the second kind of data of the units in the dynamic domain of the flow is stored by adopting a dynamic data structure;

s6: residual error estimation;

s7: time integration;

s10: reallocating storage space of the second type of data;

s12: repeating steps S5-S11 until a convergence condition is reached;

s13: and outputting the result.

2. The method according to claim 1, characterized in that in step S2, when initializing according to inflow conditions, the conservative of all cells in the aircraft flow field grid is assigned as an inflow value; when initializing according to a given flow field, assigning the conservation quantities of all units in the aircraft flow field grid as given flow field values.

3. The method according to claim 1, wherein step S3 is implemented as follows:

1) for the way of initialization according to incoming flow conditions:

2) for the way of initialization according to a given flow field:

flow characteristics of cells in the initial convective dynamic domain satisfy

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

In the formula, W represents conservationAn amount; subscript ∞ indicates incoming flow conditions; | Δ W⁽¹⁾||_maxThe maximum value of all unit conservation-constant updating quantities in the 1 st iteration step is solved;_a,cthe newly added threshold value of convection is taken as 10^-4≤_a,c≤10^-6；

The initial viscous dynamic domain comprises a viscous effect leading unit in the initial convection dynamic domain, the effect of the viscous effect is determined according to a standardized viscous effect parameter, a viscous effect parameter psi is defined and represents the ratio of viscous disturbance to non-viscous disturbance mass flow, when the disturbance speed is represented by a grid scale weighted value of the disturbance speed, the viscous effect parameter psi is represented as the ratio of the convection flux and the spectral radius of a viscous flux Jacobian matrix, namely

In the formula, I, J and K represent grid directions;

Is shown as

4. The method according to claim 1, wherein step S6 is implemented as follows:

the flow control equation in semi-discrete form is expressed as:

5. The method according to claim 1, wherein step S8 is implemented as follows:

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

in the formula, | Δ W | | | represents a module value of a conservation quantity update quantity of the convection dynamic domain boundary unit for solving the current iteration step; | Δ W⁽¹⁾||_maxThe maximum value of all unit conservation-constant updating quantities in the 1 st iteration step is solved;_a,cthe newly added threshold value of convection is taken as 10^-4≤_a,c≤10^-6；

u·q+a＞0 (7)

the unit vector q is one of the following a-c:

6. The method according to claim 1, wherein the 4 determinations in step S9 are specifically as follows:

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

Judging that the convective dynamic domain boundary unit has converged;

Wherein u represents a flow velocity vector; theta_dRepresenting the tolerance angle of the upstream unit, and taking theta more than or equal to 5 degrees_d≤10°；

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

in the formula (I), the compound is shown in the specification,Δ t represents the local time step; c_CFLCFL number representing a time advance format; | Ω | represents the volume of the grid cell; i, J, K represent grid directions; Δ RⁱRepresenting the influence of the adjacent units in the flow dynamic domain on the residual error of the flow dynamic domain boundary unit along the direction I, wherein I is I, J and K;

In the formula, Δ W represents a conservative update amount; Δ F_cRepresenting the amount of change of the convection flux, namely the difference of the convection flux of the current iteration step and the previous iteration step; Δ S represents the area of the cell face; subscript i + -1 represents the immediate vicinity of the convective dynamic domain boundary cell along the positive and negative i directions; subscript i + -1/2 denotes the cell faces of the convective dynamic domain boundary cell in the positive and negative i-directions;

for viscous units

In the formula (I), the compound is shown in the specification,

7. The method according to claim 3, wherein the specific determination procedure of the sub-step S111 in the step S11 is as follows:

When the solution converges, i.e. | | Δ W | | non-woven calculation_max＝(||ΔW||_max)_minThe new threshold of actual viscosity is_a,v(ii) a When the solution divergence or convergence rate is oscillating, i.e., | | Δ W | | non-woven cells_max>(||ΔW||_max)_minThe actual viscosity increase threshold used is less than_a,vTherefore, the problem of calculation failure caused by insufficient viscous dynamic domain is avoided by increasing the viscous dynamic domain.