A numerical simulation method for supersonic shock wave turbulence interference

Abstract

The invention belongs to the technical field of computational fluid dynamics, and discloses a numerical simulation method for supersonic shock wave turbulence interference, comprising the following steps: S100 initializing flow physical quantities in the Reynolds-averaged Navier-Stokes equations and turbulence characteristic quantities in the #imgabs0#SST turbulence model equation; S200 solving the RANS equation to obtain the #imgabs2# at the #imgabs1#th moment; calculating the scale-reduced dimensionless adverse pressure gradient quantity #imgabs3# according to the correlation value obtained in S200; calculating the adjustment function #imgabs5# according to the #imgabs4# value obtained in S300 to obtain an improved turbulent kinetic energy generation term #imgabs6#; solving the improved #imgabs8#SST turbulence model equation according to the #imgabs7# obtained in S400 to obtain #imgabs9#; substituting the #imgabs10# obtained in S500 into S200, and solving it iteratively until the RANS equation residual is Res n Less than a certain value e or reach the maximum number of steps N = N max When , the iterative solution stops. According to the average value of the latest flow physical quantity output by S600 #imgabs11#, the numerical simulation ends.

Description

Technical Field

[0001] The invention belongs to the technical field of computational fluid dynamics, and in particular relates to a numerical simulation method for supersonic shock wave turbulence interference. Background Art

[0002] When supersonic flow passes through the obvious deformation of the aircraft surface (such as control surfaces, inlet compression surfaces, etc.), complex shock wave interference will occur, usually accompanied by separation and reattachment of boundary layer flow, significant enhancement of aerodynamic heat, etc. Therefore, how to accurately and quickly predict complex shock wave interference is of great significance to the overall aerodynamic evaluation of the aircraft.

[0003] The Reynolds-Averaged Navier-Stokes (RANS) method is a commonly used method in the industrial sector during the aerodynamic design of aircraft, and is usually used for aerodynamic evaluation of aircraft. Compared with direct numerical simulation or large eddy simulation, it has the advantages of fewer computational grids and higher computational efficiency. One of the key technologies of the RANS method is the turbulence model, among which the two-equation model (such as the Shear Stress Transport turbulence model, or SST model) is one of the most widely used models in aerospace engineering. This model has good performance and application in typical flows such as attached boundary layers, shear layers, and small separations, but it cannot yet predict supersonic complex shock wave turbulence with strong interference flows.

[0004] Existing computational experience shows that the standard SST model over-predicts the separation zone induced by supersonic strong shock wave interference, while the existing improved SST model has limited improvement in the prediction of aerodynamic forces of different forms of shock wave interference flows, especially in three-dimensional flows, and its universality needs to be enhanced. Summary of the invention

[0005] The object of the present invention is to provide a numerical simulation method for supersonic shock wave turbulence interference in order to improve the inaccurate separation zone, surface pressure and friction resistance near the shock wave interference, thereby improving the overall aerodynamic prediction of the aircraft.

[0006] The technical solution adopted by the present invention is as follows: a numerical simulation method for supersonic shock wave turbulence interference, the numerical simulation method comprising the following steps: S100: initializing the average value of flow physical quantities in the Reynolds-averaged Navier-Stokes RANS equation, and Turbulence characteristic quantities in the SST turbulence model; S200: Solve the RANS equation to obtain the n The relevant value at time ; S300: Calculate the dimensionless adverse pressure gradient of the scale reduction according to the correlation value obtained in S200 ; S400: According to the result obtained in S300 Value, calculate the adjustment function , the improved turbulent kinetic energy generation term is obtained ; S500: According to the result obtained in S400 , solve the improved The SST turbulence model equation is obtained ; S600: According to the result obtained in S500 Substitute into step S200 and iterate to solve until the RANS equation residual Res n Less than a certain value ε or reach the maximum number of steps N = N max When , the iterative solution stops; S700: according to the result of S600, output the latest N The average value of the flow physical quantity at the moment ends the numerical simulation.

[0007] Furthermore, the step S100 includes: initializing the relevant quantities in the RANS equation, including: average density , average speed , average pressure ,as well as Turbulent Kinetic Energy in SST Turbulence Model , specific dissipation rate and the turbulent eddy viscosity coefficient ; Among them, “-” represents the time average quantity, and “~” represents the mass average quantity.

[0008] Furthermore, the step S200 includes: (1) The molecular viscosity stress tensor and Reynolds stress tensor are: (2) Heat flux for: (3) Among them, “―” represents the time average, and “~” represents the mass average. is the average density, is the average speed, is the average pressure, is the total energy per unit mass, is the average temperature, is the molecular dynamic viscosity coefficient, is the turbulent eddy viscosity coefficient, is the turbulent kinetic energy, is the Kronecker tensor, For time, or is the physical space coordinate, , is a tensor index, and when it is 1, 2, or 3, it represents , , , is the molecular viscous stress tensor, is the Reynolds stress tensor; constant coefficient , , .

[0009] Furthermore, the step S300 includes: Computation of the dimensionless adverse pressure gradient with reduced scale , the specific formulas are as follows: (4) in, The norm of the pressure gradient, is the modulus of velocity, is the average density, is the molecular dynamic viscosity coefficient.

[0010] Further, the step S400 is based on the Value, calculate the adjustment function : (5) (6) in, and is a constant, , , is the dimensionless adverse pressure gradient after scaling, Expressed as spatial position coordinates; Turbulent kinetic energy The turbulent kinetic energy generation term on the right side of the equation The correction is expressed as: (7) The generated items are constrained as follows: (8) in is the destruction term in the turbulent kinetic energy equation, and its expression is , adjustment function and transition function They are all functions of spatial position and flow in the boundary layer with zero pressure gradient; is the Reynolds stress tensor, is the average speed, is the physical space coordinate; when the adjustment function is = 1, it returns to the standard turbulent kinetic energy generation term.

[0011] Introducing the adjustment function , in the improved The SST equations are solved to obtain the average value of the latest flow physical quantities, and the accuracy of SST model prediction is improved based on the average value.

[0012] Furthermore, the steps S500 to S700 specifically include the following steps: the turbulent kinetic energy generation term obtained in step S400 Substitution The SST turbulence model equations are solved iteratively: (9) (10) in, is the average density, is the molecular dynamic viscosity coefficient, is the turbulent kinetic energy, is the specific dissipation rate, is the turbulent eddy viscosity coefficient, Inner mode and the outer mode Branch switching function, is the physical space coordinate; Constant coefficient Inner mode The corresponding constant in With outer mode The corresponding constant in Mixed according to the following relationship: (11) Among them, the constant coefficient ,Pick Any one; It should be noted that the constant coefficient The specific value of is consistent with that in the original SST model (Menter F R. Two-equationeddy-viscosity turbulence models for engineering applications, AIAA Journal, 1994, 32: 1598-1605).

[0013] The implicit LUSGS method is used to iterate the above formulas (9) and (10) to obtain , superscript Expressed as Physical quantity of the time step; Calculate the turbulent eddy viscosity coefficient , the specific form is as follows: (12) in is the Bradshaw constant, take , is the average density, is the average density of n time steps, is the strain rate tensor Model , is the turbulent eddy viscosity coefficient The mixing function in .

[0014] Furthermore, the turbulent eddy viscosity coefficient Substitute into formula (1) and perform iterative calculation until the Time step, at which the RANS equation iteration meets the preset convergence condition or reaches the maximum set number of calculation steps , the calculation stops.

[0015] Furthermore, according to the improved equation, the average value of the latest flow physical quantity is output , end the calculation.

[0016] Furthermore, in the equation, the spatial discretization adopts the MUSCL format, which is used with the minmod limiter to capture shock waves, the flux format is AUSMPW+, and the time discretization adopts the LU-SGS method.

[0017] It should be noted that the above spatial and temporal discretization methods are prior art.

[0018] Furthermore, the constant and The calibration is performed by directly simulating the data of the strong interference problem of shock waves induced by compression corners, where The constant is used to adjust the amplitude of the increase in turbulent kinetic energy generation. Used to adjust the ease of the transition function.

[0019] In summary, due to the adoption of the above technical solution, the beneficial effects of the present invention are: 1. Through the turbulent kinetic energy The generated term in the equation introduces the adjustment function , which is calibrated based on DNS data of compression corner flow to improve the prediction accuracy of turbulent kinetic energy; 2. Through improved The SST model can improve the prediction accuracy of different forms of shock wave turbulence interfering with the aerodynamic force of the flow, especially the universality of the model is enhanced in three-dimensional configurations. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] Figure 1 is a flow chart of the method of the present invention; Figure 2 This is a comparison diagram of the friction resistance coefficient of the flat plate with 13° incident oblique shock wave of the present invention; Figure 3 For the present invention Comparison of midline pressure on the double support plate surface. DETAILED DESCRIPTION

[0021] The present invention will be described in detail below in conjunction with the accompanying drawings.

[0022] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention.

[0023] like Figure 1 As shown, through k The generated term in the equation introduces the adjustment function , in order to improve the prediction accuracy of the SST model.

[0024] The numerical simulation method comprises the following steps: S100: Initializes the average values ​​of flow quantities in the Reynolds-averaged Navier-Stokes RANS equations, and Turbulence characteristic quantities in the SST turbulence model; initialization of related quantities in the RANS equation, including: average density , average speed , average pressure ,as well as Turbulent Kinetic Energy in SST Turbulence Model , specific dissipation rate and the turbulent eddy viscosity coefficient ; Among them, “-” represents the time average quantity, and “~” represents the mass average quantity.

[0025] S200: Solve the RANS equation and get n The relevant value at the moment; (1) The molecular viscosity stress tensor and Reynolds stress tensor are: (2) Heat flux for: (3) Among them, “―” represents the time average, and “~” represents the mass average. is the average density, is the average speed, is the average pressure, is the total energy per unit mass, is the average temperature, is the molecular dynamic viscosity coefficient, is the turbulent eddy viscosity coefficient, is the turbulent kinetic energy, is the Kronecker tensor, For time, or is the physical space coordinate, , is a tensor index, and when it is 1, 2, or 3, it represents , , , is the molecular viscous stress tensor, is the Reynolds stress tensor; Constant coefficient , , .

[0026] S300: Calculate the dimensionless adverse pressure gradient of the scale reduction according to the correlation value obtained in S200 ; Computation of the dimensionless adverse pressure gradient with reduced scale , the specific formulas are as follows: (4) in, The norm of the pressure gradient, is the modulus of velocity, is the average density, is the molecular dynamic viscosity coefficient.

[0027] S400: Based on the results obtained in S300 Value, calculate the adjustment function , the improved turbulent kinetic energy generation term is obtained ; According to the step S300 Value, calculate the adjustment function : (5) (6) in, and is a constant, , , is the dimensionless adverse pressure gradient after scaling, Expressed as spatial position coordinates; Turbulent kinetic energy The turbulent kinetic energy generation term on the right side of the equation The correction is expressed as: (7) The generated items are constrained as follows: (8) in is the destruction term in the turbulent kinetic energy equation, and its expression is , adjustment function and transition function They are all functions of spatial position and flow in the boundary layer with zero pressure gradient; is the Reynolds stress tensor, is the average speed, is the physical space coordinate; when the adjustment function is = 1, it returns to the standard turbulent kinetic energy generation term.

[0028] S500: Based on the results from S400 , solve the improved The SST turbulence model equation is obtained ; S600: Based on what was learned from S500 Substitute into step S200 and iterate to solve until the RANS equation residual Res n Less than a certain value ε or reaches the maximum number of steps N = N max When , the iterative solution is stopped; S700: Output the latest data according to the results of S600 N The average value of the flow physical quantity at the moment ends the numerical simulation.

[0029] The turbulent kinetic energy generation term obtained in step S400 Substitution The SST turbulence model equations are solved iteratively: (9) (10) in, is the average density, is the molecular dynamic viscosity coefficient, is the turbulent kinetic energy, is the specific dissipation rate, is the turbulent eddy viscosity coefficient, Inner mode and the outer mode Branch switching function, is the physical space coordinate; Constant coefficient Inner mode The corresponding constant in With outer mode The corresponding constant in Mixed according to the following relationship: (11) Among them, the constant coefficient ,Pick Any one; The implicit LUSGS method is used to iterate the above formulas (9) and (10) to obtain , superscript Expressed as Physical quantity of the time step; Calculate the turbulent eddy viscosity coefficient , the specific form is as follows: (12) in is the Bradshaw constant, take , is the average density, is the average density of n time steps, is the strain rate tensor Model , is the turbulent eddy viscosity coefficient The mixing function in .

[0030] The turbulent eddy viscosity coefficient Substitute into formula (1) and perform iterative calculation until the Time step, at which the RANS equation iteration meets the preset convergence condition or reaches the maximum set number of calculation steps , the calculation stops.

[0031] Output latest N The average value of the flow physical quantity at any moment , end the calculation.

[0032] The method is to The generated term in the equation introduces the adjustment function , in the improved The SST equations are solved to obtain the improved turbulent eddy viscosity coefficient , in order to improve the prediction accuracy of the SST model.

[0033] Example 1 like Figure 2As shown, this embodiment is an incident oblique shock wave flat plate configuration, which is used to verify the numerical prediction ability of the modified SST turbulence model for the flow field with strong adverse pressure gradient. The incoming flow condition is Mach number 2.89 and the incoming flow Reynolds number , the incident shock wave plate angle is 13°.

[0034] Figure 2 A comparison between the distribution of the surface friction drag coefficient in a 13° incident oblique shock wave / flat plate and the experimental results is given. It can be seen that the separation zone predicted by the standard SST turbulence model is too large, which is quite different from the experimental results. The modified SST turbulence model proposed in this embodiment predicts a smaller separation zone size than the standard SST turbulence model, which is close to the experimental results. This improves the prediction accuracy of the SST turbulence model for two-dimensional shock wave turbulence interference.

[0035] Example 2 like Figure 3 As shown, this embodiment is a three-dimensional configuration of double leading edge support plates, which is used to verify the effectiveness and wide applicability of the modified SST turbulence model in three-dimensional complex shock wave interference. The configuration of this example is Symmetrical double leading edge support plate, the incoming flow condition is Mach number 3.92, the incoming flow Reynolds number .

[0036] Figure 3 The pressure distribution at the midline position of the double-bracket configuration surface is given. It can be seen that the results calculated by the modified SST turbulence model are closer to the experimental results at the separation point. Secondly, the numerical prediction value of the low-pressure area after the cross shock wave is also closer to the experimental value.

[0037] The above descriptions are only two typical embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A numerical simulation method for supersonic shock wave turbulence interference, characterized in that: The numerical simulation method comprises the following steps: S100: Initializes the average values ​​of flow quantities in the Reynolds-averaged Navier-Stokes RANS equations, and Turbulence characteristic quantities in the SST turbulence model; S200: Solve the RANS equation and get n The relevant value at time ; S300: Calculate the dimensionless adverse pressure gradient of the scale reduction according to the correlation value obtained in S200 ; S400: Based on the results obtained in S300 Value, calculate the adjustment function , the improved turbulent kinetic energy generation term is obtained ; S5 00: According to S400 , solve the improved The SST turbulence model equation is obtained ; S600: Based on what was learned from S500 Substitute into step S200 and iterate to solve until the RANS equation residual Res n Less than a certain value ε or reaches the maximum number of steps When , the iterative solution is stopped; S700: Output the latest data according to the results of S600 N The average value of the flow physical quantity at the moment ends the numerical simulation.

2. The numerical simulation method for supersonic shock wave turbulence interference according to claim 1 is characterized in that: The step S100 includes: initializing the relevant quantities in the RANS equation, including: average density , average speed , average pressure ,as well as Turbulent Kinetic Energy in SST Turbulence Model , specific dissipation rate and the turbulent eddy viscosity coefficient ; Among them, "-" represents the time average quantity, and "~" represents the mass average quantity.

3. The numerical simulation method for supersonic shock wave turbulence interference according to claim 1 is characterized in that: The step S200 includes: (1) The molecular viscosity stress tensor and Reynolds stress tensor are: (2) Heat flux for: (3) Among them, "―" represents the time average, "~" represents the mass average, is the average density, is the average speed, is the average pressure, is the total energy per unit mass, is the average temperature, is the molecular dynamic viscosity coefficient, is the turbulent eddy viscosity coefficient, is the turbulent kinetic energy, is the Kronecker tensor, For time, or is the physical space coordinate, , is a tensor index, and when it is 1, 2, or 3, it represents , , , is the molecular viscous stress tensor, is the Reynolds stress tensor; Constant coefficient , , .

4. The numerical simulation method for supersonic shock wave turbulence interference according to claim 1 is characterized in that: The step S300 includes: Computation of the dimensionless adverse pressure gradient with reduced scale , the specific formulas are as follows: (4) in, The norm of the pressure gradient, is the modulus of velocity, is the average density, is the molecular dynamic viscosity coefficient.

5. The numerical simulation method for supersonic shock wave turbulence interference according to claim 1 is characterized in that: The step S400 is based on the Value, calculate the adjustment function : (5) (6) in, and is a constant, , , is the dimensionless adverse pressure gradient after scaling, Expressed as spatial position coordinates; Turbulent kinetic energy The turbulent kinetic energy generation term on the right side of the equation The correction is expressed as: (7) The generated items are constrained as follows: (8) in is the destruction term in the turbulent kinetic energy equation, and its expression is , adjustment function and transition function They are all functions of spatial position and flow in the boundary layer with zero pressure gradient; is the Reynolds stress tensor, is the average speed, is the physical space coordinate; when the adjustment function is = 1, it returns to the standard turbulent kinetic energy generation term.

6. The numerical simulation method for supersonic shock wave turbulence interference according to claim 1 is characterized in that: The steps S500 to S700 specifically include the following steps: Substitution The SST turbulence model equations are solved iteratively: (9) (10) in, is the average density, is the molecular dynamic viscosity coefficient, is the turbulent kinetic energy, is the specific dissipation rate, is the turbulent eddy viscosity coefficient, Inner mode and the outer mode Branch switching function, is the physical space coordinate; Constant coefficient Inner mode The corresponding constant in With outer mode The corresponding constant in Mixed according to the following relationship: (11) Among them, the constant coefficient Pick Any one; The implicit LUSGS method is used to iterate the above formulas (9) and (10) to obtain , superscript Expressed as The physical quantity of the time step; Calculate the turbulent eddy viscosity coefficient , the specific form is as follows: (12) in is the Bradshaw constant, take , is the average density, is the average density of n time steps, is the strain rate tensor Model , is the turbulent eddy viscosity coefficient The mixing function in .

7. The numerical simulation method for supersonic shock wave turbulence interference according to claim 3 is characterized in that: The turbulent eddy viscosity coefficient Substitute into formula (1) and perform iterative calculation until the Time step, at which the RANS equation iteration meets the preset convergence condition or reaches the maximum set number of calculation steps , the calculation stops.

8. The numerical simulation method for supersonic shock wave turbulence interference according to claim 7 is characterized in that: Output latest N The average value of the flow physical quantity at any moment , end the calculation.

9. The numerical simulation method for supersonic shock wave turbulence interference according to claim 3 is characterized in that: The MUSCL format is used for spatial discretization, and the minmod limiter is used to capture shock waves. The flux format is AUSMPW+, and the LU-SGS method is used for time discretization.

10. The numerical simulation method for supersonic shock wave turbulence interference according to claim 5, characterized in that: constant and The calibration is performed by directly simulating the data of the strong interference problem of shock waves induced by compression corners, where The constant is used to adjust the amplitude of the increase in turbulent kinetic energy generation. Used to adjust the ease of the transition function.

Images