A numerical simulation method for supersonic shock wave turbulence interference

Abstract

The present invention belongs to the technical field of computational fluid dynamics, and discloses a numerical simulation method for supersonic shock-turbulence interaction, including the following steps: S100 Initialize the flow physical quantities in the Reynolds-averaged Navier-Stokes equations and the turbulent characteristic quantities in the #imgabs0# SST turbulence model equations; S200 Solve the RANS equations to obtain #imgabs2# at the #imgabs1#-th moment; Calculate the scale-normalized dimensionless inverse pressure gradient quantity #imgabs3# according to the relevant values obtained in S200; Calculate the adjustment function #imgabs5# according to the #imgabs4# value obtained in S300 to obtain the improved turbulent kinetic energy production term #imgabs6#; Solve the improved #imgabs8# SST turbulence model equations according to the #imgabs7# obtained in S400 to obtain #imgabs9#; Substitute the #imgabs10# obtained in S500 into S200 and perform iterative solution until the residual of the RANS equations Res n is less than a certain value ε or reaches the maximum number of steps N = N max at which point, stop the iterative solution. Output the average value #imgabs11# of the latest flow physical quantities according to S600 to end the numerical simulation.

Description

Technical Field

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

[0002] When supersonic flow passes through obvious deformed parts of the aircraft surface (such as control surfaces, inlet compression surfaces, etc.), complex shock interference phenomena will occur, usually accompanied by phenomena such as separation and reattachment of boundary layer flow, and significant enhancement of aerodynamic heating. Therefore, how to accurately and quickly predict complex shock interference is of great significance for 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 process of aircraft, and is usually used for the aerodynamic evaluation of aircraft. Compared with methods such as 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, and the two-equation model (such as the Shear Stress Transport turbulence model, i.e., the 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 well predict supersonic complex shock-turbulence strong interference flows.

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

[0005] The object of the present invention is to provide a numerical simulation method for supersonic shock-turbulence interference in view of the above problems, hoping to improve the problems of inaccurate separation zones, surface pressures, and frictional drag near shock interference, and further improve the overall aerodynamic force prediction of the aircraft.

[0006] The technical solution adopted by the present invention is as follows: A numerical simulation method for supersonic shock-turbulence interference, the numerical simulation method includes the following steps: S100: Initialize the average values of the flow physical quantities in the Reynolds-averaged Navier-Stokes RANS equation, and the turbulence characteristic quantities in the SST turbulence model; S200: Solve the RANS equation to obtain the relevant values at the n th moment ; S300: Calculate the scale-normalized dimensionless adverse pressure gradient according to the relevant values obtained in S200 ; S400: Calculate the adjustment function according to the value obtained in S300 , an improved turbulent kinetic energy generation term is obtained ; S500: According to what is obtained in S400 , solve the improved SST turbulence model equation to obtain ; S600: According to what is obtained in S500 Substitute it into step S200 and perform iterative solution until the RANS equation residual Res n is less than a certain value ε or the maximum number of steps is reached N = N max ; then stop the continued iterative solution; S700: According to the result of S600, output the average value of the flow physical quantity at the latest N moment, and end the numerical simulation.

[0007] Furthermore, in the step S100, it includes: initializing the relevant quantities in the RANS equation, including: average density , average velocity , average pressure , and the turbulent kinetic energy , specific dissipation rate and turbulent eddy viscosity coefficient in the SST turbulence model; where, "―" represents the time-averaged quantity, and "~" represents the mass-averaged quantity.

[0008] Furthermore, in the step S200, it includes:

[0009] (1)

[0010] where the molecular viscous stress tensor and the Reynolds stress tensor are:

[0011] (2)

[0012] The heat flux density is:

[0013] (3)

[0014] where, "―" represents the time-averaged quantity, "~" represents the mass-averaged quantity, is the average density, is the average velocity, 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, is the time, or is the physical space coordinate, , is the tensor index. When taking 1, 2, and 3 respectively, it represents , , , is the molecular viscous stress tensor, is the Reynolds stress tensor; the constant coefficients , , .

[0015] Furthermore, the step S300 includes:

[0016] Calculate the scale-normalized dimensionless inverse pressure gradient , and the specific formula is as follows:

[0017] (4)

[0018] where, is the modulus of the pressure gradient, is the modulus of the velocity, is the average density, is the molecular dynamic viscosity coefficient.

[0019] Furthermore, the step S400 calculates the adjustment function according to the value obtained in the step S300:

[0020] (5)

[0021] (6)

[0022] where, and are constants, , , is the scale-normalized dimensionless inverse pressure gradient, is expressed as the spatial position coordinate;

[0023] For the turbulent kinetic energy in the right end of the equation, correct the turbulent kinetic energy production term , and its expression form is:

[0024] (7)

[0025] Constrain the production term as follows:

[0026] (8)

[0027] wherein is the destruction term in the turbulent kinetic energy equation, and the expression is , the adjustment function and the transition function are both functions of the spatial position and are in the boundary layer flow with zero pressure gradient; is the Reynolds stress tensor, is the mean velocity, is the physical space coordinate; when the adjustment function is = 1, it is restored to the standard turbulent kinetic energy generation term;

[0028] The adjustment function is introduced, and the average values of the latest flow physical quantities are solved in the improved SST equation set, and the accuracy of the SST model prediction is improved according to the average values.

[0029] Furthermore, the steps S500 to S700 specifically include the following: substituting the turbulent kinetic energy generation term obtained in step S400 into the SST turbulence model equation for iterative solution: (9)

[0030] (10)

[0031] wherein, 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, is the inner layer mode and the outer layer mode branch switching function, is the physical space coordinate;

[0032] The constant coefficient is the corresponding constant in the inner layer mode and the corresponding constant in the outer layer mode are mixed according to the following relationship:

[0033] (11)

[0034] wherein, the constant coefficient , take any one;

[0035] It should be noted that the specific values of the constant coefficients are the same as those in the original SST model (Menter F R. Two - equation eddy - viscosity turbulence models for engineering applications, AIAA Journal, 1994, 32: 1598 - 1605).

[0036] The above formulas (9) and (10) are solved by using the implicit LUSGS method for iteration to obtain , and the superscript represents the physical quantity at the time step;

[0037] Calculate the turbulent eddy - viscosity coefficient , and its specific form is as follows:

[0038] (12)

[0039] where is the Bradshaw constant, taking , is the average density, is the average density at the n - th time step, is the strain - rate tensor modulus , is the turbulent eddy - viscosity coefficient in the mixing function.

[0040] Furthermore, substitute the turbulent eddy - viscosity coefficient into formula (1) for iterative calculation until the time step, at which time the RANS equation iteration satisfies the preset convergence condition or reaches the maximum set number of calculation steps , then stop the calculation.

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

[0042] Furthermore, in the equation, MUSCL format is used for spatial discretization, combined with the minmod limiter to capture shock waves, the flux format is selected as AUSMPW +, and LU - SGS method is used for time discretization.

[0043] It should be noted that the above spatial and time discretization methods are existing technologies.

[0044] Furthermore, the constants and Calibrate by directly simulating the data of strong shock interference induced by compressed corners, where The constant is used to adjust the increased amplitude of the turbulent kinetic energy generation term, and the constant is used to adjust the smoothness of the transition function.

[0045] In summary, due to the adoption of the above technical solutions, the beneficial effects of the present invention are as follows:

[0046] 1. By introducing a regulation function into the generation term in the turbulent kinetic energy equation, which is calibrated based on the DNS data of the compressed corner flow, to improve the prediction accuracy of the turbulent kinetic energy;

[0047] 2. The improved SST model can improve the prediction accuracy of the aerodynamic force of different forms of shock-turbulence interference flows, especially enhancing the universality of the model in three-dimensional configurations. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] Figure 1 is the flowchart of the method of the present invention;

[0049] Figure 2 is the comparison diagram of the skin friction coefficient of the flat plate with a 13° incident oblique shock of the present invention;

[0050] Figure 3 is the comparison diagram of the midline pressure of the double strut surface of the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0051] The present invention will be described in detail below with reference to the accompanying drawings.

[0052] In order to make the objectives, technical solutions and advantages of the present invention clearer, the present invention will be further described in detail below with reference to 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 used to limit the present invention.

[0053] As Figure 1 shown, by introducing a regulation function k into the generation term in the equation to improve the prediction accuracy of the SST model.

[0054] The numerical simulation method includes the following steps:

[0055] S100: Initialize the average values of the flow physical quantities in the Reynolds-averaged Navier-Stokes (RANS) equations, and the turbulent characteristic quantities in the SST turbulence model; Initialize the relevant quantities in the RANS equations, including: average density average velocity , average pressure , and the turbulent kinetic energy in the SST turbulence model , specific dissipation rate and turbulent eddy viscosity coefficient ; where, "―" represents the time-averaged quantity, and "~" represents the mass-averaged quantity.

[0056] S200: Solve the RANS equations to obtain the relevant values at the n moment; (1)

[0057] where the molecular viscous stress tensor and the Reynolds stress tensor are:

[0058] (2)

[0059] Heat flux density is:

[0060] (3)

[0061] where, "―" represents the time-averaged quantity, "~" represents the mass-averaged quantity, is the average density, is the average velocity, 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, is the time, or is the physical space coordinate, , is the tensor index, and when taking 1, 2, 3 respectively, it represents , , , is the molecular viscous stress tensor, is the Reynolds stress tensor;

[0062] Constant coefficient , , .

[0063] S300: Calculate the scale-normalized dimensionless inverse pressure gradient ;

[0064] Calculation of the Dimensionless Inverse Pressure Gradient with Scale Regularization , and the specific formula is as follows:

[0065] (4)

[0066] Wherein, is the modulus of the pressure gradient, is the modulus of the velocity, is the average density, is the molecular dynamic viscosity coefficient.

[0067] S400: Calculate the adjustment function according to the value obtained in S300 to obtain the improved turbulent kinetic energy generation term ;

[0068] Calculate the adjustment function according to the value obtained in step S300:

[0069] (5)

[0070] (6)

[0071] Wherein, and are constants, , , is the dimensionless inverse pressure gradient with scale regularization, is expressed as the spatial position coordinate;

[0072] Correct the turbulent kinetic energy generation term in the right - hand side of the turbulent kinetic energy equation, and its expression form is:

[0073] (7)

[0074] Constrain the generation term as follows:

[0075] (8)

[0076] Wherein is the destruction term in the turbulent kinetic energy equation, and the expression is , the adjustment function and the transition function are both functions of the spatial position and in the boundary - layer flow with zero pressure gradient; is the Reynolds stress tensor, is the average velocity, is the physical space coordinate; when the adjustment function is When = 1, it is restored to the standard turbulent kinetic energy generation term;

[0077] S500: According to what is obtained in S400 , solve the improved SST turbulence model equation to obtain ;

[0078] S600: Substitute what is obtained in S500 into step S200 and perform iterative solution until the RANS equation residual Res n is less than a certain value ε or reaches the maximum number of steps N = N max at which point, stop the continued iterative solution;

[0079] S700: According to the result of S600, output the average value of the flow physical quantity at the latest N moment to end the numerical simulation.

[0080] Substitute the turbulent kinetic energy generation term obtained in step S400 into the SST turbulence model equation for iterative solution:

[0081] (9)

[0082] (10)

[0083] where, 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, is the inner layer mode and the outer layer mode branch switching function, is the physical space coordinate;

[0084] The constant coefficient is the corresponding constant in the inner layer mode and the corresponding constant in the outer layer mode are mixed according to the following relationship:

[0085] (11)

[0086] where, the constant coefficient , take any one of them;

[0087] The above formulas (9) and (10) are solved iteratively using the implicit LUSGS method to obtain , and the superscript represents the physical quantity at the time step;

[0088] Calculate the turbulent eddy viscosity coefficient , and its specific form is as follows:

[0089] (12)

[0090] where is the Bradshaw constant, taking , is the average density, is the average density at the nth time step, is the strain rate tensor of the modulus , is the turbulent eddy viscosity coefficient in the mixing function.

[0091] Substitute the turbulent eddy viscosity coefficient into formula (1) for iterative calculation until the time step, at which time the RANS equation iteration satisfies the preset convergence condition or reaches the maximum set number of calculation steps , then stop the calculation.

[0092] Output the average value of the flow physical quantity at the latest N moment , and end the calculation.

[0093] The method introduces a regulation function for the generation term in the turbulent kinetic energy equation, and solves for the improved turbulent eddy viscosity coefficient in the improved SST equation set to improve the prediction accuracy of the SST model.

[0094] Example 1

[0095] As Figure 2 shown, this example is an incident oblique shock flat plate configuration, used to verify the numerical prediction ability of the modified SST turbulence model for the flow field with a strong adverse pressure gradient. The oncoming flow conditions are a Mach number of 2.89, an oncoming flow Reynolds number , and the incident shock flat plate angle is 13°.

[0096] Figure 2The comparison between the distribution of the friction drag coefficient of the object surface in the 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 than the standard SST turbulence model, which is close to the experimental results. That is, the SST turbulence model has improved the accuracy of predicting two-dimensional shock wave turbulence interference.

[0097] Example 2

[0098] If Figure 3 As shown in the figure, 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 .

[0099] Figure 3 The pressure distribution at the midline of the double-bracket configuration 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.

[0100] The above 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 , is the average density, is the average speed, is the average pressure; 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 ; The step S400 is based on the data obtained in 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; S500: Based on the results from S400 , solve the improved The SST turbulence model equation is obtained , is the turbulent kinetic energy, is the specific dissipation rate, is the turbulent eddy viscosity coefficient, Expressed as time; 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 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 .

6. 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.

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

8. 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.

9. The numerical simulation method for supersonic shock wave turbulence interference according to claim 1, 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