Method for adaptively controlling computational fluid mechanics flux format dissipation

Abstract

The invention discloses a method for adaptively controlling and calculating fluid mechanics flux format dissipation, which comprises the following steps of: when fluid mechanics of an aircraft to be subjected to high Reynolds number flow simulation is solved, solving a flux by utilizing a discrete control equation set in each iteration step; obtaining a value of a physical quantity on a grid unit in a flow field of a next iteration step, and determining left and right temperature values, left and right pressure values and left and right Mach numbers of each surface based on the value of the physical quantity; calculating left and right pressure gradients of each surface according to the pressure values of the left and right grid units of each surface; calculating a compression factor on each surface, and determining a discontinuous detection factor of each grid unit based on the compression factor; and introducing a discontinuous detection threshold, converting the discontinuous detection factor of the grid unit into a discontinuous detection factor of each surface on the grid unit when the circulation amount of the discrete control equation set is calculated by adopting a Roe format, and determining whether the flow field is interrupted or not according to the discontinuous detection factor of each surface and the discontinuous detection threshold.

Description

Technical Field

[0001] The present invention relates to the technical field of numerical simulation in computational fluid dynamics, and specifically relates to a method for adaptively controlling the dissipation of a computational fluid dynamics flux format. Background Art

[0002] In computational fluid dynamics, the dissipation characteristics of numerical formats play a crucial role in the accuracy and stability of simulation results. An ideal numerical format should be able to maintain low dissipation in smooth regions to capture small-scale flow structures, while increasing dissipation in discontinuous regions (such as shock waves) to ensure numerical stability. However, traditional numerical formats often struggle to meet both of these requirements simultaneously. For example, high-precision central formats perform well in smooth regions but are prone to non-physical oscillations in discontinuous regions; while upwind formats with larger dissipation can stably capture shock waves but introduce excessive dissipation in smooth regions, reducing simulation accuracy. Considering that most numerical flux formats can be regarded as a combination of a central format and a dissipation term, researchers can flexibly adjust the combination of the central format and the dissipation term according to the specific needs of the problem. For example, reducing the dissipation term and increasing the weight of the central format in smooth flow regions to capture small-scale flow structures; while increasing the weight of the dissipation term in discontinuous flow regions to suppress numerical oscillations and ensure computational stability. By reasonably dynamically adjusting the magnitude and distribution of the dissipation term, we can design high-precision and high-stability numerical formats suitable for simulating complex flow problems such as high Mach number and high Reynolds number flows.

[0003] To dynamically adjust the magnitude of the dissipation term, the method of this application proposes an adaptive dissipation mechanism based on a discontinuity detection function. By introducing an adaptive dissipation control parameter into the dissipation term of the numerical format, the numerical format only increases dissipation in discontinuous regions, thereby minimizing numerical dissipation in smooth regions while ensuring numerical stability. This method can effectively solve the limitations of traditional numerical formats in complex flow simulations and provide a more accurate and efficient numerical simulation tool for complex problems such as high Reynolds number flows and the interaction between shock waves and turbulence.

[0004] The adaptive dissipation adjustment method can make the numerical format dynamically adjust its dissipation characteristics according to local flow characteristics by introducing a regulation factor related to the local flow field. This method requires the ability to accurately detect discontinuities in the flow and only increase the dissipation of the numerical format in these regions. Therefore, developing effective discontinuity detectors and adaptive dissipation formats is crucial for improving the accuracy and reliability of numerical simulations. However, for complex geometries and flow conditions, such as high Reynolds number boundary layer flows, existing discontinuity detection methods are still not perfect, which may cause the numerical format to increase dissipation in smooth regions, further reducing numerical accuracy. Summary of the Invention

[0005] The object of the present invention is to provide a method for adaptively controlling the dissipation of the flux format in computational fluid dynamics to overcome the problems existing in the prior art.

[0006] To achieve the above object, the present invention adopts the following technical solutions:

[0007] A method for adaptively controlling the dissipation of the flux format in computational fluid dynamics, comprising:

[0008] When performing hydrodynamic solution for an aircraft to be simulated with high Reynolds number flow, by solving the convective flux using the discrete control equations at each iteration step, obtaining the values of physical quantities on the grid cells in the flow field at the next iteration step, and based on the values of the physical quantities, determining the left and right values of each face included in each grid cell, determining the left and right temperature values and left and right pressure values of each face based on the left and right values of each face, and judging the positive and negative of the left and right temperature values and left and right pressure values of each face to ensure they are positive values, and on this basis, determining the left and right Mach numbers of each face;

[0009] Calculating the left and right pressure gradients of each face according to the pressure values of the left and right grid cells of each face; wherein, the left and right pressure gradients of the face are equal to the pressure gradients at the cell centers of the left and right two grid cells adjacent to the face;

[0010] Calculating the compression factor on each face based on the left and right pressure gradients, left and right pressure values, and left and right Mach numbers of each face, and determining the discontinuity detection factor of each grid cell based on the compression factor;

[0011] Introducing a discontinuity detection threshold, when calculating the convective flux of the discrete control equations using the Roe format, converting the discontinuity detection factor of the grid cell into the discontinuity detection factor of each face on the grid cell, and determining whether there is an interruption in the flow field through the discontinuity detection factor of each face and the discontinuity detection threshold.

[0012] Further, after performing the (n - 1)-th step iterative solution of the convective flux using the discrete control equations, the values of physical quantities on the grid cells in the flow field at the n-th iteration step can be obtained where the parameters ρ, u, v, w, E represent density, velocities in the x, y, and z directions, and energy respectively, the parameter superscript n represents the n-th iteration step, and the subscript i represents the i-th grid cell;

[0013] Calculating the left and right values of each face included in each grid cell based on the values of the physical quantities where the parameter superscripts L and R represent the left value calculated using the grid cell on the left side of the face and the right value calculated using the grid cell on the right side of the face respectively; the subscripts i, j represent the i-th and j-th grid cells adjacent to the left and right of the currently calculated face, and the double subscript i, j represents that the parameter is the physical quantity on the face between grid cell i and grid cell j;

[0014] Among them, the left and right temperatures and the left and right pressure values The formulas used in the calculation are respectively:

[0015]

[0016] In the above formula, It means that the parameter in the brackets is the physical quantity in the left value of the face, It means that the parameter in the brackets is the physical quantity in the right value of the face; γ is the specific heat ratio of the gas, U represents the local velocity,

[0017] Furthermore, judging the positive and negative of the left and right temperature values and the left and right pressure values of each face to ensure that they are positive values, including:

[0018] Detect whether there is a negative value in the left and right pressure values of the face. If there is, the pressure takes the preset value ε. If not, it remains unchanged; for the left and right temperature values of the face, if the left and right temperature values are negative, take the absolute value, otherwise remain unchanged.

[0019] Furthermore, the calculation process of the left and right Mach numbers of each face is:

[0020] First, it is necessary to calculate the left and right sound speed values of the face according to the left and right temperatures of the face where R is the gas constant, and calculate the left and right Mach numbers of the face respectively represent the left and right local velocities of the face; then according to the normal vector n of the face i,j Calculate the left and right normal Mach numbers of the face and the left and right tangential Mach numbers

[0021] Furthermore, calculating the left and right pressure gradients of each face according to the left and right pressure values of the grid cells of each face, including:

[0022] For the finite volume method of the second-order format, the calculation method of the pressure gradient at the cell center is the same as the calculation method of the conserved quantity gradient. Based on the pressure value on the i-th grid cell at the n-th iteration step Use the Green-Gauss method or the least squares method to calculate the pressure gradient at the cell center of the i-th grid cell

[0023] Furthermore, calculate the compression factor α on the face at the n-th iteration step ij :

[0024]

[0025] In the above formula, α ijrepresents the compression factor on the currently calculated surface, which is the surface between the i-th and j-th grid cells adjacent left and right; p i , p j represent the pressure values on the i-th and j-th grid cells respectively, represent the pressure gradients at the cell centers of the i-th and j-th grid cells respectively, Ω i , Ω j represent the volumes of the i-th and j-th grid cells respectively, ||·|| represents the modulus operation, and ∈ is a preset pressure threshold.

[0026] Furthermore, calculate the discontinuity detection factor α for each grid cell i :

[0027]

[0028] In the above formula, α i represents the discontinuity detection factor of the i-th grid cell, N f represents the number of faces of the i-th grid cell.

[0029] Furthermore, the discontinuity detection factor for each face adopts the maximum value of the discontinuity detection factors of the two adjacent left and right grid cells;

[0030] When α f < α thr , it is considered that there is a discontinuity in the flow field, and at this time the Roe format dissipation term remains unchanged; when α f > α thr , it is considered that the local flow field is smooth, and at this time the Roe format dissipation term is multiplied by a coefficient β, usually taking a value of 0 < β < 1; where α f is the discontinuity detection factor of the f-th face, and α thr is the discontinuity detection threshold.

[0031] A terminal device includes a processor, a memory, and a computer program stored in the memory; when the processor executes the computer program, the method for adaptively controlling the dissipation of the computational fluid dynamics flux format is implemented.

[0032] A computer-readable storage medium stores a computer program; when the computer program is executed by a processor, the method for adaptively controlling the dissipation of the computational fluid dynamics flux format is implemented.

[0033] Compared with the prior art, the present invention has the following technical features:

[0034] 1. The implementation of this method is relatively simple. Only a sub-function for solving the discontinuity detection function value needs to be added when calculating the numerical flux each time, and the modification to the original program is relatively small;

[0035] 2. This method achieves the goal of dynamically adjusting the numerical format dissipation according to local flow field changes. The additional empirical parameters have less dependence on specific cases, reducing the program's reliance on the user's experience and improving the robustness and computational efficiency of the numerical simulation process. Description of the Drawings

[0036] Figure 1 is the adaptive dissipation numerical format process of the present invention;

[0037] Figure 2 is the contour map of the pressure coefficient distribution on the symmetry plane of the test case adopted by the present invention, where (a) is the overall figure and (b) is the partial enlarged figure;

[0038] Figure 3 is the contour map of the discontinuity detection factor distribution on the symmetry plane of the test case adopted by the present invention;

[0039] Figure 4 is the contour map of the vorticity distribution at the same moment obtained by the original Roe format;

[0040] Figure 5 is the contour map of the vorticity distribution at the same moment obtained by the adaptive dissipation Roe format of the present invention application. Detailed Implementation Manner

[0041] Aiming at the accuracy problems of shock wave capturing and turbulent small-scale structure simulation in high Reynolds number flows of aircraft, a method for implementing an adaptive numerical format for dynamically adjusting numerical dissipation is proposed. First, a three-dimensional unstructured computational grid is constructed based on the geometric parameters of the aircraft, and working condition parameters such as the incoming flow Mach number and angle of attack are set to initialize the physical quantity field of the flow field. During the iterative solution of the Navier-Stokes equations, the dissipation adaption is achieved through the following steps: 1) Real-time detect the abnormal values of negative pressure and negative temperature at the grid interface and perform physical corrections; 2) Use the least squares method or the Green-Gauss method to calculate the pressure gradients of the left and right cells at the interface, and construct a compression factor in combination with the normal / tangential Mach number to quantify the flow discontinuity intensity; 3) Design a cell discontinuity detection factor, fuse the pressure gradient and Mach number characteristics to enhance the sensitivity to weak shock waves; 4) Dynamically adjust the coefficient of the numerical format dissipation term based on the detection factor, maintaining strong dissipation in the shock wave region to ensure stability and reducing numerical dissipation in the smooth flow region to improve the resolution. Finally, the shock wave structure, surface pressure coefficient distribution and aerodynamic parameters of the aircraft's flow field are output, solving the problem that it is difficult to balance accuracy and stability in complex flows with traditional numerical formats.

[0042] See the attached Figure 1 , a method for adaptively controlling the dissipation of computational fluid dynamics flux formats, comprising the following steps:

[0043] Step 1: For the aircraft to be simulated for high Reynolds number flow, a three-dimensional computational domain is constructed based on the three-dimensional model of the aircraft. Commercial grid generation software such as Pointwise or ICEM CFD is used to mesh the computational domain. The types of grid cells used in the meshing include tetrahedrons, hexahedrons, or hybrid grids. Then, the Navier-Stokes equations are discretized to obtain the discrete control equations. The discrete control equations perform spatial discretization of the mass, momentum, and energy conservation equations by the finite volume method, which is used to calculate the evolution process of the flow field physical quantities (including density, velocities in the x, y, and z directions, and energy) over time, and the updated flow field physical quantities are output at each iteration step.

[0044] After performing the (n - 1)-th step iterative solution of the flux using the discrete control equations, the values of the physical quantities on the grid cells in the flow field at the n-th iteration step can be obtained. Among them, the parameters ρ, u, v, w, and E represent density, velocities in the x, y, and z directions, and energy respectively. The superscript n of the parameter represents the n-th iteration step, and the subscript i represents the i-th grid cell.

[0045] Based on the values of the physical quantities on the grid cells Calculate the left and right values of each face contained in each grid cell. Among them, the superscripts L and R of the parameter represent the left value calculated using the grid cell on the left side of the face and the right value calculated using the grid cell on the right side of the face respectively; the subscripts i and j represent the i-th and j-th grid cells adjacent to the left and right of the currently calculated face. The single subscript i represents that the parameter is the physical quantity of grid cell i, and the double subscript i, j represents that the parameter is the physical quantity on the face between grid cell i and grid cell j. The same applies hereinafter. Calculate the left and right temperature values of the face based on the left and right values of each face of the grid cell. And the left and right pressure values Then, detect whether there are negative values in the left and right temperature values and the left and right pressure values of the face; if the left and right temperature values are negative, take their absolute values; if the left and right pressure values are negative, replace them with a preset value; then, calculate the left and right Mach numbers of the face using the left and right temperature values of the face. Among them, the left and right temperature values, the left and right pressure values, and the left and right Mach numbers are the temperature, pressure values, and Mach numbers calculated using the grid cells on the left and right sides of the face respectively.

[0046] Among them, the left and right temperatures And the left and right pressure values The formulas used in the calculation are respectively:

[0047]

[0048]

[0049] In the above formula, It indicates that the parameter in the brackets is a physical quantity in the left value of the face, indicating that the parameter in the brackets is a physical quantity in the right value of the face; γ is the specific heat ratio of the gas, generally taken as 1.4 for air, and U represents the local velocity,

[0050] In one embodiment of the present invention, it is detected whether there are negative values in the left and right pressure values of the face. If there are, the pressure is taken as 10 -16 , and if not, it remains unchanged, that is For the left and right temperature values of the face, if there are negative values in the left and right temperature values, the absolute value is taken otherwise it remains unchanged. Among them When the superscript is L, it represents the left pressure value, and when the superscript is R, it represents the right pressure value, and the meaning of is the same.

[0051] Next, it is necessary to calculate the left and right Mach numbers of the face through the left and right temperature values

[0052] First, it is necessary to calculate the left and right sound speed values of the face according to the temperatures on the left and right of the face where R is the gas constant, and calculate the left and right Mach numbers of the face respectively represent the left and right local velocities of the face; then according to the normal vector n of the face i,j calculate the left and right normal Mach numbers of the face and the left and right tangential Mach numbers

[0053] It should be particularly noted that for the second-order finite volume method, the left and right values of the face can be reconstructed according to the gradients at the cell centers of the left grid cell and the right grid cell.

[0054] It should be noted that n is a positive integer greater than or equal to 1. Before performing iterative solution, in addition to setting the Roe format entropy correction coefficient for the parameters related to the adaptive Roe format dissipation, the following parameters also need to be set in the input file of the flow field solver: the current case discontinuity detection threshold α thr and the smooth region dissipation coefficient β, etc.; according to the inventor's experience, for the vast majority of aircraft aerodynamic simulation problems, the above parameters can be set as follows in sequence: α thr = 0.5, β = 0.1.

[0055] In addition, when the iteration step n = 1, at this time, the smooth region dissipation coefficient β of each grid cell in the flow field is 1.0. Generally speaking, in order to ensure the stability of numerical calculation during initial calculation, the user can enable the adaptive dissipation function after a certain number of iteration steps (the user can decide according to the working conditions of the case and the grid quality).

[0056] Step 2: Calculate the left and right pressure gradients of each face based on the pressure values of the left and right grid cells of each face.

[0057] After the processing in Step 1, there are no negative temperature and pressure values in all grid cells in the flow field; for the second-order finite volume method, the left and right pressure gradients of a face are equal to the pressure gradients at the cell centers of the left and right two grid cells adjacent to this face; the calculation method of the pressure gradient at the cell center is the same as that of other conserved quantity gradients, based on the pressure value on the i-th grid cell at the n-th iteration step The pressure gradient at the cell center of the i-th grid cell can be calculated using the Green-Gauss method or the least squares method

[0058] In the embodiments of the present invention, since the pressure gradient quantity is not required in the numerical solution process of the flow field, the general solver does not solve the pressure gradient. Therefore, it is necessary to explain the calculation of the pressure gradient here. The solver uses the second-order finite volume method to solve the flow control equations, calculates the gradients of the solution variables through the Green-Gauss method and the least squares method, and generally uses a limiter in order to maintain the stability of the numerical solution near the shock wave. However, the present invention does not use a limiter in the process of solving the pressure gradient. The reason is that the limiter will make the gradient of the physical quantity become smaller or even zero near the discontinuity during use, which obviously does not conform to the original intention of using the pressure gradient to detect the discontinuity. Therefore, the calculation of the pressure gradient involved in the present invention can directly use the gradient calculation module in the solver, and a similar method can also be directly used for other solvers, and finally the pressure gradient at the cell center of each unit i at the n-th iteration step is obtained In addition, for the second-order finite volume method, the pressure gradient at the cell center is the same as the value at the face center. Therefore, it can be directly considered that the left and right pressure gradient values of the face are the gradient values at the cell centers of the left and right cells.

[0059] Step 3: When calculating the convective flux at the n-th iteration step of the discrete control equation set, calculate the compression factor α on the face at the n-th iteration step ij :

[0060]

[0061]

[0062] In the above formula, α ij represents the compression factor on the face being currently calculated, and this face is the face between the i-th and j-th grid cells adjacent left and right; p i , p j represent the pressure values on the i-th and j-th grid cells respectively, represent the pressure gradients at the cell centers of the i-th and j-th grid cells respectively, and Ω i , Ωj represent the volumes of the \(i\)-th and \(j\)-th grid cells respectively, \(\|\cdot\|\) represents the modulus operation, and \(\epsilon\) is a preset pressure threshold. In this embodiment, \(\epsilon = 10\). -6 .

[0063] The compression factor \(\alpha\) in the above formula ij fully considers the strong discontinuities in the flow field, including the sharp changes in physical quantities brought about by shock waves and strong shear layers, while the detection of weak shock waves is further enhanced by introducing the function .

[0064] The value of the compression factor for each face is obtained through the above formula, and then the discontinuity detection factor \(\alpha\) of each grid cell is calculated i :

[0065]

[0066] In the above formula, \(\alpha\) i represents the discontinuity detection factor of the \(i\)-th grid cell, and \(N\) f represents the number of faces of the \(i\)-th grid cell. It should be noted here that the value range of the compression factor \(\alpha\) on the face is \(0 < \alpha\) ij < 1, and the value of the discontinuity detection factor \(\alpha\) of the grid cell obtained through the above formula ij is necessarily less than 1. i

[0067] Step 4, through the discontinuity detection factor of each grid cell, it can be judged whether there is a discontinuity in the flow field. This requires introducing a discontinuity detection threshold \(\alpha\) thr , which is used as the criterion for judging whether there is a discontinuity in the flow field. The method of this application takes the Roe format as an example to elaborate in detail the implementation process of the adaptive dissipation method, but its application is not limited to the Roe format and can be extended to other numerical formats; when calculating the flux of the \(n\)-th iteration step of the discrete control equations using the Roe format, operations are performed on the faces of the grid cells. Therefore, it is necessary to convert the discontinuity detection factor of the grid cell into the discontinuity detection factor of each face on the grid cell:

[0068] And the discontinuity detection factor of each face takes the maximum value of the discontinuity detection factors of the left and right adjacent grid cells, that is, \(\alpha\) f = max(\(\alpha\) i , \(\alpha\) j ), where \(\alpha\) f is the discontinuity detection factor of the \(f\)-th face, and \(\alpha\) i , \(\alpha\) j are the discontinuity detection factors of the \(i\)-th and \(j\)-th grid cells adjacent to the left and right of the \(f\)-th face.

[0069] When \(\alpha\) f < \(\alpha\) thrWhen there is a discontinuity in the flow field, the dissipative term of the Roe scheme remains unchanged; when α f > α thr it is considered that the local flow field is smooth. At this time, the dissipative term of the Roe scheme is multiplied by a coefficient β, usually taking a value of 0 < β < 1. Thus, two parameters α thr and β are introduced in the adaptive dissipation of the Roe scheme. Generally speaking, these two parameters need to be given different values according to different cases. In particular, the setting of the discontinuity detection threshold directly determines whether the detector judges the local flow field as discontinuous or smooth. It can be summarized from a large number of cases that α thr = 0.5 is suitable for most cases. And the parameter β determines the dissipation magnitude of the flow field in the smooth region, and according to experience, it can usually be set to 0.1. To sum up, after adopting the discontinuity detection function of the present invention, the dissipative term coefficient at the smooth part of the flow field is less than 1, while the dissipative term coefficient at the flow field discontinuity is equal to 1, thus realizing the adaptive dissipation of the Roe scheme.

[0070] Next, briefly introduce the Roe scheme:

[0071]

[0072] In the above formula, F represents the value of the numerical flux, and the subscript i represents the i-th grid cell, represents the face of the i-th grid cell; the second term is the dissipative term of the Roe scheme, and respectively represent the k-th eigenvalue and eigenvector of the Roe average matrix, N is the number of eigenvalues, ΔU = U i+1 - U i represents the change of the conserved quantity U, and U i represents the conserved quantity on the i-th grid cell; as mentioned above, this scheme hopes that the dissipative term can adaptively change with the local flow field information. Therefore, only a coefficient β needs to be added before the dissipative term to control the size of the dissipative term:

[0073]

[0074] And the size of the coefficient β depends on the local flow field information. Near the shock wave, β should be 1.0 to increase the dissipation and improve the stability of numerical calculation, while in the smooth region of the flow field, the value of β can be less than 1.0 to increase the resolution of the flow field solution.

[0075] In summary, after adopting the discontinuous detection function of the present invention, the dissipation term coefficient at the smooth part of the flow field is less than 1, while the dissipation term coefficient at the discontinuous part of the flow field is equal to 1, thus realizing the adaptive dissipation of the numerical format. Further, an adaptive dissipation numerical format is used to solve the discrete control equations until the final iteration converges and the lift and drag coefficients, surface pressure coefficient distribution, and other physical quantity parameters of the flow field around the aircraft are output.

[0076] The present invention enables the Roe format with adaptive dissipation to dynamically adjust the numerical dissipation characteristics in different flow regions. It can increase dissipation in complex flow regions to maintain numerical stability and reduce dissipation in smooth regions to improve simulation accuracy, providing efficient and reliable aerodynamic analysis for the numerical simulation of complex flows of aircraft at high Reynolds numbers. The implementation of this method is relatively simple. Compared with the original Roe format, using this method can significantly reduce the dependence of the program on the user's experience, improve the robustness of the program, and while maintaining numerical stability, reduce numerical dissipation to enhance the ability to capture small-scale flow structures.

[0077] Example:

[0078] The interaction between supersonic cross-flow and high-speed jet in a high-speed aircraft will generate complex shock structures and vortex distributions, including shock-shock interference, shock-boundary layer interference, and complex vortex structures such as horseshoe vortices, hovering vortices, and counter-rotating vortex pairs (CVP), which can be used to test the effectiveness of the adaptive dissipation numerical format. The calculation conditions of the example are as follows: free-stream Reynolds number Re = 2.1*10 6 , Mach number Ma = 3, temperature T = 103.2K, pressure P = 0.1949×10 5 Pa; lateral jet: Mach number Ma = 1.0, temperature T = 223K, total pressure ratio Pr = 200. For comparative verification, the dissipation term coefficient β of the original method is 1; in the improved adaptive dissipation method, β = 0.1 in the smooth region and β = 1 in the discontinuous region. In addition, the discontinuous detection threshold α thr = 0.5.

[0079] Figure 1 The algorithm flow chart of the method of the present application is given. Figure 2 It is the cloud map of the pressure coefficient distribution on the symmetry plane of the current example obtained by the adaptive method. It can be clearly seen from the figure that there are complex shock structures in the flow field, and the effectiveness of the adaptive method can be tested through the drastic changes in the local flow field. Figure 3 It gives the cloud map of the discontinuous detection factor distribution on the symmetry plane. It can be clearly seen from the figure that the values of the discontinuous detection factor are relatively small at the discontinuous part of the flow field, while the values of the discontinuous detection factor are relatively large at the smooth part of the flow field. Therefore, the discontinuous detection function proposed by the method of the present application can effectively detect the discontinuity in the flow field.

[0080] Figures 4 - 5 They are the vortex structures at the same moment calculated by the original Roe format and the adaptive dissipation Roe format. In this example, the Rortex method is used to identify the vortex structures in the flow field, and the pressure coefficient is used for coloring. It can be clearly seen from the figure that at the same moment, the vortex structures calculated by the adaptive dissipation Roe format are very obvious. The horseshoe vortices in the main separation region, the large-scale vortex structures in the jet wake region, and the obvious counter-rotating vortex pair (CVP) and other typical flow field structures can be clearly seen. However, fewer vortex structures are calculated by using the original Roe format. Except for the weak horseshoe vortex structure that can be seen, almost no counter-rotating vortices can be seen. This shows that the original Roe format brings greater dissipation in the smooth region of the flow field, so that the typical vortex structures in the flow field are dissipated. The adaptive dissipation Roe format of the method of the present application can effectively reduce the dissipation in the smooth region and capture the small-scale structures in the flow field, especially the vortex structures in the jet wake region.

[0081] The above embodiments are only used to illustrate the technical solutions of the present application, rather than to limit them; although the present application has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that: they can still modify the technical solutions recorded in the foregoing embodiments, or perform equivalent replacements on some of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present application, and should all be included in the protection scope of the present application.

Claims

1. A method for adaptively controlling the dissipation of a computational fluid dynamics flux format, characterized in that: include: A fluid mechanics solution for an aircraft to be subjected to high Reynolds number flow simulation is obtained by solving the flow rate using a discrete control equation group at each iteration step to obtain the value of a physical quantity on a grid unit in the flow field of the next iteration step, and based on the value of the physical quantity, a left and right value of each face contained in each grid unit is determined, and based on the left and right values ​​of each face, a left and right temperature value and a left and right pressure value of each face are determined, and the positive and negative of the left and right temperature value and the left and right pressure value of each face are judged to ensure that they are positive values, and on this basis, a left and right Mach number of each face is determined; The left and right pressure gradient of each surface is calculated according to the pressure values ​​of the left and right grid cells of each surface; wherein the left and right pressure gradient of the surface is equal to the pressure gradient at the centers of the left and right grid cells adjacent to the surface; Based on the left-right pressure gradient, the left-right pressure value, and the left-right Mach number of each surface, the compressibility factor on each surface is calculated, and the discontinuity detection factor of each grid cell is determined based on the compressibility factor; The discontinuity detection threshold is introduced. When the flow rate of the discrete control equations is calculated using the Roe format, the discontinuity detection factor of the grid unit is converted into the discontinuity detection factor of each surface on the grid unit. The discontinuity detection factor of each surface and the discontinuity detection threshold are used to determine whether there is an interruption in the flow field.

2. The method for adaptively controlling the dissipation of a computational fluid dynamics flux format according to claim 1, characterized in that: After the discrete control equations are used to iterate and solve the flow rate in the n-1th step, the value of the physical quantity on the grid unit in the flow field in the nth iteration step can be obtained: The parameters ρ, u, v, w, and E represent density, velocity in the x, y, and z directions, and energy, respectively. The superscript n indicates the nth iteration step, and the subscript i indicates the i-th grid unit. Based on the value of the physical quantity, calculate the left and right values ​​of each face contained in each grid unit The superscripts L and R of the parameters represent that the parameter is the left value calculated by the grid cells on the left side of the surface and the right value calculated by the grid cells on the right side of the surface respectively; the subscripts i and j represent the i-th and j-th grid cells adjacent to the current calculated surface on the left and right sides, and the double subscripts i and j represent that the parameter is the physical quantity on the surface between grid cells i and j; Among them, the temperature of And left and right pressure values The formulas used in the calculation are: In the above formula, Indicates that the parameter in the brackets is the physical quantity in the left value of the face. Indicates that the parameters in the brackets are the physical quantities in the right value of the surface; γ is the specific heat ratio of the gas, U represents the local velocity, 3. The method for adaptively controlling the dissipation of computational fluid dynamics flux format according to claim 1, characterized in that: The positive and negative values ​​of the left and right temperature values ​​and the left and right pressure values ​​of each surface are judged to ensure that they are positive values, including: Check whether the left and right pressure values ​​of the surface are negative. If so, the pressure takes the preset value ε, otherwise it remains unchanged. For the left and right temperature values ​​of the surface, if they are negative, the absolute value is taken, otherwise it remains unchanged.

4. The method for adaptively controlling the dissipation of a computational fluid dynamics flux format according to claim 2, characterized in that: The calculation process of the left and right Mach numbers of each surface is as follows: First, we need to calculate the left and right sound speed values ​​of the surface based on the temperature of the left and right sides of the surface Where R is the gas constant, and the left and right Mach numbers of the surface are calculated Represent the local velocities of the left and right sides of the surface respectively; then according to the normal vector n of the surface i,j Calculate the left and right normal Mach numbers of the surface and the left and right tangential Mach numbers 5. The method for adaptively controlling the dissipation of computational fluid dynamics flux format according to claim 3, characterized in that: The calculating the left and right pressure gradients of each surface according to the pressure values ​​of the left and right grid cells of each surface includes: For the second-order finite volume method, the calculation method of the pressure gradient at the grid center is consistent with the method of calculating the gradient of the conserved quantity, based on the pressure value on the i-th grid cell at the n-th iteration step The pressure gradient at the center of the i-th grid cell is calculated using the Green-Gauss method or the least squares method.

6. The method for adaptively controlling the dissipation of a computational fluid dynamics flux format according to claim 4, characterized in that: Calculate the compression factor α on the surface at the nth iteration step ij : In the above formula, α ij Indicates the compression factor on the currently calculated surface, which is the surface between the i-th and j-th grid cells adjacent to each other on the left and right sides; p i ,p j Represent the pressure values ​​on the i-th and j-th grid cells respectively, Represents the pressure gradient at the center of the i-th and j-th grid cells, Ω i ,Ω j Represent the volumes of the i-th and j-th grid cells respectively, ||·|| represents the modulo operation, and ∈ is the preset pressure threshold.

7. The method for adaptively controlling the dissipation of a computational fluid dynamics flux format according to claim 6, characterized in that: Calculate the discontinuity detection factor α for each grid cell i : In the above formula, α i represents the discontinuity detection factor of the i-th grid unit, N f Represents the number of faces in the ith grid cell.

8. The method for adaptively controlling the dissipation of computational fluid dynamics flux format according to claim 1, characterized in that: The discontinuity detection factor on each face adopts the maximum value of the discontinuity detection factors of the two adjacent left and right grid cells; When α f <α thr When α f >α thr When , the local flow field is considered to be smooth. At this time, the Roe format dissipation term is multiplied by a coefficient β, which is usually taken as 0<β<1; where α f is the discontinuity detection factor of the fth surface, α thr is the discontinuity detection threshold.

9. A terminal device comprising a processor, a memory and a computer program stored in the memory; characterized in that: When the processor executes the computer program, the method for adaptively controlling the dissipation of the computational fluid dynamics flux format according to any one of claims 1 to 8 is implemented.

10. A computer-readable storage medium, wherein a computer program is stored in the medium; characterized in that: When the computer program is executed by a processor, the method for adaptively controlling the dissipation of a computational fluid dynamics flux format according to any one of claims 1 to 8 is implemented.

Images