Wide-speed-domain boundary layer transition prediction method using algebraic intermittent factor

Abstract

The invention discloses a wide-speed-domain boundary layer transition prediction method using an algebraic intermittent factor, and the method comprises the steps: constructing an intermittent factor function in an algebraic form through calculating a transition indication variable and a transition criterion which can be solved locally in a flow field based on a universal Reynolds average Navier-Stokes computational fluid mechanics method; and a generation item and a dissipation item of the turbulent energy transport equation of the turbulence model are further corrected through an intermittent factor function, so that accurate prediction of transition from the subsonic velocity to the hypersonic velocity boundary layer is realized.

Description

Technical Field

[0001] The present invention relates to the field of hypersonic boundary layer transition prediction, and in particular to a wide speed range boundary layer transition prediction method. Background Art

[0002] In recent years, with the continuous upgrading of demands for high-speed transportation and round-trip travel between the earth and the sky, aircraft have gradually developed in the direction of large airspace, wide speed range, and reusability. In such aircraft, boundary layer transition is a common flow phenomenon. Boundary layer transition characterizes the evolution of flow from laminar state to turbulent state, which significantly affects the friction resistance of the aircraft wall, heat exchange and flow separation, and is of great significance to the design of the aerodynamic shape, thermal protection system and propulsion system of the aircraft. Therefore, in the design of wide-speed hypersonic aircraft, boundary layer transition prediction is a key technical issue.

[0003] Transition prediction techniques usually include experimental research, direct numerical simulation, large eddy simulation, flow stability analysis, transition model, etc. Among the above methods, the transition model has the advantages of easy use, good robustness and low computing resource requirements, and is a more commonly used transition prediction method in current engineering practice. However, the flow scenarios faced by wide-speed hypersonic vehicles are very complex, with many transition mechanisms and influencing factors. Compared with a single subsonic / transonic scenario, the traditional transition model based on multiple transport equations faces the defects of poor computational robustness, low prediction accuracy, and high difficulty in model implantation when applied to wide-speed domain flows, which hinders the transition prediction and aerodynamic design of wide-speed hypersonic vehicles.

[0004] Therefore, it is necessary to propose a new technical solution to solve the above technical problems. Summary of the invention

[0005] In view of the above problems, the present invention provides a wide speed domain boundary layer transition prediction method using algebraic intermittency factors, the purpose of which is to accurately predict the subsonic to hypersonic boundary layer transition.

[0006] In order to achieve the above object, the technical solutions that can be adopted in the present invention are as follows:

[0007] A method for predicting boundary layer transition in a wide velocity domain using an algebraic intermittent factor comprises the following steps:

[0008] (1) Using the average flow information, the vortex Reynolds number Re is calculated. v ;

[0009] (2) Through the vortex Reynolds number Re v , calculate the local momentum thickness Reynolds number Re used to characterize the transition flow field θL ;

[0010] (3) According to the input free flow turbulence Tu∞ and the Mach number M at the outer edge of the boundary layer calculated in the previous step e Calculate the transition criterion Re θt The value of

[0011] (4) According to the local momentum thickness Reynolds number Re of the flow field θL and transition criterion Re θt , construct the algebraic intermittent factor γ;

[0012] (5) Use the intermittent factor γ to modify the turbulence model, including the generation term P of the transport equation for the turbulent kinetic energy k k , dissipation term D k , the mixed function F1 is used for correction:

[0013] P k =γP k,SST

[0014] D k ={f lim min[max(γ,0.1),1]+(1-fl lim )}D k,SST

[0015] Where P k,SST and D k,SST They represent the generation term and dissipation term of the turbulent kinetic energy k transport equation in the original k-ωSST turbulence model respectively; the function fl lim To ensure the reasonable dissipation of turbulent kinetic energy in free flow;

[0016]

[0017] Among them, F 1SST represents the mixing function in the original k-ωSST turbulence model, F3 represents the limiter function, and R y Represents the distance d based on the wall w and the Reynolds number of the turbulent kinetic energy k;

[0018] (6) Using a general CFD method in any computational fluid dynamics solver, the turbulence model equations modified by the above steps are solved to obtain the flow velocity, density, temperature, force coefficient, and wall heat flux data.

[0019] Furthermore, in step (1), Re v The calculation formula is:

[0020]

[0021] Where ρ is the density, d w is the wall distance, Ω is the modulus of vorticity, and μ is the molecular viscosity.

[0022] Furthermore, in step (2), the local momentum thickness Reynolds number Re θL as follows:

[0023]

[0024] Among them, F ratio Re represents the maximum value of the vortex Reynolds number Re v,max The momentum thickness Reynolds number Re obtained by integration θ The specific formula is as follows:

[0025]

[0026] Among them, c1, c2 and c3 are formula coefficients, M e represents the Mach number at the outer edge of the boundary layer, T e represents the temperature at the outer edge of the boundary layer, T w represents the wall temperature; for isothermal walls, T w Determined by input parameters; for an adiabatic wall, T w Take adiabatic wall temperature γ g is the specific heat ratio, and Pr is the laminar Prandtl number.

[0027] Furthermore, in step (2), M e and T e The following isentropic relationship is used for solution:

[0028]

[0029]

[0030] Among them, p ∞ is the free flow pressure, ρ ∞ represents the density of the free flow, U ∞ is the speed of the free stream, a ∞ represents the sound velocity of the free flow; p represents the local pressure of each grid point in the flow field calculation, ρ e represents the density at the outer edge of the boundary layer, U e represents the velocity at the outer edge of the boundary layer and a e represents the speed of sound at the outer edge of the boundary layer, and R represents the gas constant.

[0031] Furthermore, in step (3), the transition criterion Re θt The calculation formula is:

[0032]

[0033] Re θt,high =62.12(Tu ∞ +0.032)-0.745 M e

[0034]

[0035] Among them, Re θt,low represents the free stream Mach number M ∞ Transition criterion value under ≤1.2, Re θt,hig represents the free stream Mach number M ∞ >Transition criterion value under 1.2.

[0036] Furthermore, in step (4), the algebraic intermittent factor γ is constructed as follows:

[0037] γ=γ1+γ2

[0038]

[0039] γ2=min(10Φ1Φ2γ1,3)

[0040]

[0041] It can be seen from the above formula that when the intermittent factor γ = 0, it means that the flow is laminar, and when 0 < γ < 1, it means that the flow is in the transition stage; the intermittent factor γ consists of the first part γ1 and the second part γ2; Φ1 and Φ2 are two limiter functions used to shield the second part γ2 of the intermittent factor in the fully developed turbulent zone and the viscous bottom layer of the boundary layer; F turb is based on the turbulent viscosity μ t and the dissipation function of the molecular viscosity μ.

[0042] Compared with the prior art, the present invention has the following beneficial effects:

[0043] It is possible to predict the transition through local information of the flow field without adding additional transport equations or non-local variable solving strategies, thus achieving accurate prediction of subsonic to hypersonic boundary layer transition and aerodynamic heating.

[0044] The design method provided by the present invention can be stored on a storage medium as a computer program, and includes the following technical solutions: an electronic device, comprising: one or more processors; and a storage device for storing one or more programs, when the one or more programs are executed by the one or more processors, the one or more processors implement the above-mentioned prediction method.

[0045] as well as:

[0046] A computer readable medium stores a computer program, which implements the above prediction method when executed by a processor. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] Figure 1 Schematic diagram of the flow chart of the novel algebraic transition model proposed in the present invention.

[0048] Figure 2 Schematic diagram of the computational grid and boundary conditions for the plate example.

[0049] Figure 3 The prediction results of the new algebraic transition model proposed in this invention and the traditional two-equation transport transition model in the subsonic plate transition example are compared.

[0050] Figure 4 The prediction results of the new algebraic transition model proposed in the present invention and the traditional two-equation transport transition model in the supersonic flat-plate transition example are compared.

[0051] Figure 5 The prediction results of the new algebraic transition model proposed in this invention and the traditional two-equation transport transition model in the hypersonic flat-plate transition example are compared. DETAILED DESCRIPTION

[0052] In order to make the purpose, design process, technical method and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the accompanying drawings.

[0053] The present invention discloses a method for predicting boundary layer transition in a wide speed range using an algebraic intermittent factor, comprising the following steps:

[0054] (1) Using the average flow information, the vortex Reynolds number Re is calculated. v as follows:

[0055]

[0056] Where ρ is the density, d w is the wall distance, Ω is the modulus of vorticity, and μ is the molecular viscosity.

[0057] (2) Further through the vortex Reynolds number Re v , calculate the local momentum thickness Reynolds number Re used to characterize the transition flow field θL as follows:

[0058]

[0059] Among them, F ratio Re represents the maximum value of the vortex Reynolds number Re v,max The momentum thickness Reynolds number Re obtained by integration θ The specific formula is as follows:

[0060]

[0061] Among them, c1, c2 and c3 are formula coefficients, M e represents the Mach number at the outer edge of the boundary layer, T e represents the temperature at the outer edge of the boundary layer, T w represents the wall temperature. For an isothermal wall, T w Determined by input parameters; for an adiabatic wall, T w Take adiabatic wall temperature γ g is the specific heat ratio, and Pr is the laminar Prandtl number. For ideal air, γ g It is usually set to 1.4, and Pr is usually set to 0.72. In addition, < e and T e They represent the Mach number and temperature at the outer edge of the boundary layer respectively, and are solved using the following isentropic relationship:

[0062]

[0063] Among them, p ∞ is the free flow pressure, ρ ∞ represents the density of the free flow, U ∞ is the speed of the free stream, a ∞ represents the speed of sound of the free flow. These four quantities are given by the boundary conditions of the calculation input. p represents the local pressure of each grid point in the flow field calculation, ρ e represents the density at the outer edge of the boundary layer, U e represents the velocity at the outer edge of the boundary layer and a e represents the speed of sound at the outer edge of the boundary layer. These four variables are obtained in CFD (Computational Fluid Dynamics) calculations. R represents the gas constant, which is given by the fluid properties and is usually taken as 287 for ideal air.

[0064] (3) According to the input free flow turbulence Tu ∞ and the Mach number M at the outer edge of the boundary layer calculated in the previous step e Calculate the transition criterion Re θt Value

[0065]

[0066] Re θt,high =62.12(Tu ∞ +0.032) -0.745 M e

[0067]

[0068] Among them, Re θt,low Indicates low-speed flow (i.e. free-flow Mach number M ∞ ≤1.2) under the transition criterion value, Reθt,high Indicates high-speed flow (i.e. free-flow Mach number M ∞ >1.2) under the transition criterion value.

[0069] (4) According to the local momentum thickness Reynolds number Re calculated in step (2) θL and the transition criterion Re calculated in step (3) θt , construct the algebraic intermittent factor γ in the following form:

[0070] γ=γ1+γ2

[0071]

[0072] γ2=min(10Φ1Φ2γ1,3)

[0073]

[0074] After this step, the intermittent factor γ used to judge the transition can be obtained. When γ = 0, it means that the flow is laminar, when 0 < γ < 1, it means that the flow is in the transition stage, and when γ = 1, it means that the flow develops into full turbulence. The intermittent factor γ consists of the first part γ1 and the second part γ2. Φ1 and Φ2 are two limiter functions used to shield the second part γ2 of the intermittent factor in the fully developed turbulent zone and the viscous bottom layer of the boundary layer; F turb is based on the turbulent viscosity μ t The dissipation function of the molecular viscosity μ can make the transition criterion value decrease rapidly after the turbulence is generated, thereby ensuring the accurate generation of the intermittent factor.

[0075] (5) In order to obtain the effect of transition on aerodynamic force, aerodynamic heat and other factors, the intermittent factor γ can be further used to correct the turbulence model after step (4). Taking the k-ωSST turbulence model as an example, the generation term P of the transport equation of turbulent kinetic energy k is k and the dissipation term D k To make corrections:

[0076] P k =γP k,SST

[0077] D k ={f lim min[max(γ,0.1),1]+(1-fl lim )}D k,SST

[0078] Where P k,SST and D k,SST They represent the generation and dissipation terms of the turbulent kinetic energy k transport equation in the original k-ωSST turbulence model. Function fl lim Used to ensure the reasonable dissipation of turbulent kinetic energy in the free flow:

[0079]

[0080] Where S is the modulus of the strain rate. Finally, the mixing function F1 in the original k-ωSST turbulence model is modified as follows:

[0081] F1=max(F 1SS ,F3),

[0082] where F 1SST represents the mixing function in the original k-ωSST turbulence model, F3 represents the limiter function, R y Represents the distance d based on the wall w and the Reynolds number of the turbulent kinetic energy k.

[0083] (6) By using a general computational fluid dynamics (CFD) method in any computational fluid dynamics (CFD) solver to solve the turbulence model equations modified by the above steps, the flow velocity, density, temperature, force coefficient, and wall heat flux information can be obtained.

[0084] The following is an example of a subsonic to hypersonic flat plate calculation to verify the specific effect of the present invention. The calculation conditions refer to the experimental conditions in the public literature, as shown in Table 1. The calculation grid and boundary conditions are as follows: Figure 2 The number of flow direction × normal grids in the upstream and downstream of the inclined wedge is 21 × 151 and 301 × 151 respectively, and the dimensionless height y of the first layer grid is + is 1.

[0085]

[0086] The transition model proposed in the present invention can be implanted in any computational fluid dynamics solver. Here, the Fenglei software open sourced by the China Aerodynamics Research and Development Center is taken as an example to illustrate the specific implementation process and calculation results of the present invention. First, a new subroutine is established in the structural grid Reynolds-averaged Navier-Stokes equation (RANS) solver of the Fenglei software to program and implement the technical solution steps (1)-(4), so as to calculate the algebraic intermittent factor function. The flow field variable information required in the calculation process can be directly obtained from the structural grid RANS solver of the Fenglei software. Then the calculated algebraic intermittent factor function is transmitted to the k-ωSST turbulence model module of the structural grid RANS solver of the Fenglei software, and the technical solution step (5) is programmed and implemented to correct the k-ωSST turbulence model. After the above programming is completed, the structural grid RANS solver of the Fenglei software is called to read the grid to start the CFD flow field calculation, that is, the technical solution step (6). When the calculation is completed, the Fenglei software can automatically output the flow field calculation results.

[0087] Figures 3 to 5 The flat plate wall friction coefficient and Stanton number distribution calculated by Fenglei software using the new algebraic transition model proposed in this invention are given, and the γ-Re distribution in Fenglei software is also given. θ The calculation results of the two-equation transport transition model are compared with the experimental results in the open literature. It can be seen that the new algebraic transition model proposed in this invention is similar to the traditional γ-Re under subsonic conditions. θ The two-equation transport transition model obtained similar transition prediction results; under supersonic and hypersonic conditions, the method proposed in the present invention obtained results that were more consistent with the experimental values, indicating the beneficial effects of the present invention.

[0088] In addition, there are many specific implementation methods and approaches of the present invention, and the above is only a preferred implementation of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be regarded as the protection scope of the present invention.

Claims

1. A method for predicting boundary layer transition in a wide velocity domain using an algebraic intermittent factor, characterized in that: The following steps are included: (1) Using the average flow information, the vortex Reynolds number Re is calculated. v ; (2) Through the vortex Reynolds number Re v , calculate the local momentum thickness Reynolds number Re used to characterize the transition flow field θL ; (3) According to the input free flow turbulence Tu ∞ and the Mach number M at the outer edge of the boundary layer calculated in the previous step e Calculate the transition criterion Re θt The value of (4) According to the local momentum thickness Reynolds number Re of the flow field θL and transition criterion Re θt , construct the algebraic intermittent factor γ; (5) Use the intermittent factor γ to modify the turbulence model, including the generation term P of the transport equation for the turbulent kinetic energy k k , dissipation term d k , the mixed function F1 is used for correction: P k =γP k,SST D k ={f lim min[max(γ,0.1),1]+(1-f lim )}D k,SST Where P k,SST and D k,SST They represent the generation term and dissipation term of the turbulent kinetic energy k transport equation in the original k-ωSST turbulence model respectively; the function lim To ensure the reasonable dissipation of turbulent kinetic energy in free flow; where F 1SST represents the mixing function in the original k-ωSST turbulence model, F3 represents the limiter function, and R y Represents the distance d based on the wall w and the Reynolds number of the turbulent kinetic energy k; (6) Using a general CFD method in any computational fluid dynamics solver, the turbulence model equations modified by the above steps are solved to obtain the flow velocity, density, temperature, force coefficient, and wall heat flux data.

2. The design method according to claim 1, characterized in that: In step (1), Re v The calculation formula is: Where ρ is the density, d w is the wall distance, Ω is the modulus of vorticity, and μ is the molecular viscosity.

3. The design method according to claim 2, characterized in that: In step (2), the local momentum thickness Reynolds number Re θL as follows: Among them, F ratio Re represents the maximum value of the vortex Reynolds number Re v,max The momentum thickness Reynolds number Re obtained by integration θ The specific formula is as follows: Among them, c1, c2 and c3 are formula coefficients, M e represents the Mach number at the outer edge of the boundary layer, T e represents the temperature at the outer edge of the boundary layer, T w represents the wall temperature; for isothermal walls, T w Determined by input parameters; for an adiabatic wall, T w Take adiabatic wall temperature γ g is the specific heat ratio, and Pr is the laminar Prandtl number.

4. The design method according to claim 3, characterized in that: In step (2), M e and T e The following isentropic relationship is used for solution: Among them, p ∞ is the free flow pressure, ρ ∞ represents the density of the free flow, U ∞ is the speed of the free stream, a ∞ represents the sound velocity of the free flow; p represents the local pressure of each grid point in the flow field calculation, ρ e represents the density at the outer edge of the boundary layer, U e represents the velocity at the outer edge of the boundary layer and a e represents the speed of sound at the outer edge of the boundary layer, and R represents the gas constant.

5. The design method according to claim 4, characterized in that: In step (3), the transition criterion Re θt The calculation formula is: Re θt,high =62.12(You ∞ +0.032) -0.745 M e Among them, Re θt,low represents the free stream Mach number M ∞ Transition criterion value under ≤1.2, Re θt,high represents the free stream Mach number M ∞ >Transition criterion value under 1.

2.

6. The design method according to claim 5, characterized in that: In step (4), the algebraic intermittent factor γ is constructed as follows: γ=γ1+…2 γ2=min(10Φ1Φ2γ1,3) It can be seen from the above formula that when the intermittent factor γ = 0, it means that the flow is laminar, and when 0 < γ < 1, it means that the flow is in the transition stage; the intermittent factor γ consists of the first part γ1 and the second part γ2; Φ1 and Φ2 are two limiter functions used to shield the second part γ2 of the intermittent factor in the fully developed turbulent zone and the viscous bottom layer of the boundary layer; F turb is based on the turbulent viscosity μ t and the dissipation function of the molecular viscosity μ.

7. An electronic device comprising: one or more processors; And the storage device is used to store one or more programs, when the one or more programs are executed by the one or more processors, the one or more processors implement the design method described in any one of claims 1 to 6.

8. A computer readable medium having a computer program stored thereon, characterized in that: When the program is executed by a processor, the design method according to any one of claims 1 to 6 is implemented.

Images