Accelerated convergence method for aerodynamic numerical simulation of aircraft

Abstract

The invention discloses an accelerated convergence method for aerodynamic numerical simulation of an aircraft, which only needs to add a function for judging and correcting the variable quantity of the physical quantity of the current step after each iteration, and the original program is slightly changed. According to the method, the target of dynamically adjusting the local CFL number according to the local flow field change is achieved, the dependence of a program on the experience of a user is reduced, and the robustness and the calculation efficiency of the numerical simulation process are improved.

Description

Technical Field

[0001] The invention belongs to the technical field of aircraft, and in particular relates to an accelerated convergence method for aircraft aerodynamic numerical simulation. Background Art

[0002] Computational fluid dynamics (CFD) is a technology that uses numerical methods and algorithms to analyze and solve fluid flow problems. It is widely used in the field of aerospace and plays an irreplaceable role in the aerodynamic design of aircraft. During the flight of an aircraft, the flow often presents complex flows in a multi-dimensional and wide speed range, including turbulence, separation, reattachment and other phenomena. This flow has the characteristics of mutual interference between dimensions and the coexistence of high and low speed flows, which significantly increases the difficulty of numerical simulation. In addition, due to the high operating cost of wind tunnel experimental equipment and limited ability to simulate real flow fields, especially for multi-physics field coupling problems and flow separation problems, the use of numerical simulation methods to predict and analyze aerodynamic characteristics problems in the aircraft design process (such as lift, drag, torque prediction, aerodynamic shape optimization, dynamic characteristics analysis of aircraft, etc.) can not only effectively reduce experimental costs and time, but also improve the accuracy and reliability of aerodynamic design. In order to accurately and quickly simulate these flow problems, high-precision and high-fidelity numerical simulation methods and a large number of high-quality grids can be used, but this will significantly increase the amount of calculation and cost of numerical simulation. To this end, it is necessary to adopt strategies in the numerical simulation process to accelerate the numerical solution process, but at the same time, it is necessary to ensure the stability of the numerical solution. In the aerodynamic numerical simulation algorithm, the CFL number (Courant-Friedrichs-Lewy number, also known as Courant number) is a very important concept. It is used to describe the relationship between the time step and the space step to ensure the stability and accuracy of the numerical calculation, and to a certain extent determines the efficiency of the numerical simulation. When the CFL number is large, it means that the information propagates a long distance on the grid in each time step, which may lead to the instability of the numerical solution; when the CFL number is small, although it can improve the stability of the numerical solution, it will increase the calculation time and cost. The physical meaning of the CFL number is that it represents the ratio of the distance that the fluid information propagates on the grid to the size of the grid unit in one time step. In the explicit time marching method, the CFL number must be less than or equal to 1 to ensure the stability of the numerical calculation, because the explicit method requires that the distance of information propagation in each time step cannot exceed one grid unit. In the implicit time marching method, the CFL number is generally not required to meet a certain condition. However, from practical experience, the CFL number cannot be infinitely large. This mainly depends on the local flow field and grid. If the local CFL number is larger, the local time step is longer, which can accelerate convergence and save calculation time. Therefore, in practical applications, it is necessary to select a suitable CFL value based on the characteristics of the specific problem and numerical method. This brings trouble to the selection of the CFL number. Generally speaking, it is hoped that a large CFL number can be given to accelerate the convergence of numerical iterations, but a larger CFL number will increase the local time step, resulting in divergence of flow field iterative calculations.Due to this characteristic of the CFL number, in the numerical simulation process, the user generally needs to manually adjust the CFL number according to the iterative convergence situation to balance the efficiency and stability of the numerical simulation, which is often very dependent on the user's experience. Summary of the invention

[0003] In order to overcome the shortcomings of the prior art, the present invention provides an accelerated convergence method for numerical simulation of aircraft aerodynamics. It only needs to add a function to judge and correct the change of the physical quantity of the current step after each iteration, and the original program needs only a small change. The method achieves the goal of dynamically adjusting the local CFL number according to the local flow field changes, reduces the program's dependence on user experience, and improves the robustness and computational efficiency of the numerical simulation process.

[0004] The technical solution adopted by the present invention to solve the technical problem is as follows:

[0005] Step 1: After the current iteration step is completed, check the change of the physical quantity on each grid;

[0006] Step 1-1: After the discrete control equations are solved iteratively for the nth step using the numerical simulation method, the change in the physical quantity on the grid unit i in the flow field at the nth step is obtained. ,in They represent density, speed in the x, y, and z directions, and total energy respectively; Indicates the change of the following variable;

[0007] Step 1-2: Based on the change in the physical quantity on grid cell i, determine whether there is a non-numeric value NaN in all grid cells. If the change in the physical quantity on at least one grid cell i in the determination result contains NaN, the current iteration is determined to have failed, and the CFL number of all grid cells is reduced. , and the physical quantity of the current iteration step is no longer updated, and the new CFL number is used to advance the solution of the current iteration step again. If the change in the physical quantity of the nth iteration step calculated next time is still NaN, the CFL number is continuously reduced until the maximum number of times k that the program allows the CFL number to be reduced continuously is reached or the program fails to solve and exits;

[0008] Step 2: Check whether the density, pressure, and internal energy on the grid cell of step n+1 have negative values;

[0009] Step 2-1: When the test in step 1 passes, there is no NaN in all the grid cells in the flow field. At this time, update the physical quantities of all the grid cells in step n+1 ; Based on the physical quantity of the n+1th step on the grid unit Solve to get the pressure on each grid element and internal energy , check the density on each grid cell ,pressure and internal energy Whether there is a negative value. If the density, pressure or internal energy in grid cell i has a negative value, the CFL number of the current grid cell is reduced and the current grid cell is set to sleep state. For grid cells without negative values, the CFL number remains unchanged and the current grid cell is set to active. ;

[0010] Step 2-2: If the density, pressure or internal energy of at least one grid cell in the flow field is negative, the current flow field is not updated, that is, the physical quantity of step n is used. Re-solve the n+1th iteration step; when the CFL number is reduced for the kth time in the n+1th iteration, if there are still negative density, pressure or internal energy in the grid cells in the flow field, then directly multiply the physical quantity changes of these grid cells by a relaxation factor less than 1 to update the physical quantity of the n+1th step ;

[0011] Step 3: Check whether the change of the physical quantity on the grid unit in step n is too large;

[0012] Physical quantities based on the nth step on the grid cells Solve to get the density on each grid cell ,pressure And the n+1th step physical quantity based on the grid unit Solve to get the density on each grid cell ,pressure , to determine whether the relative changes in density and pressure at the n+1th iteration step exceed the set maximum range, that is, or If at least one of the two conditions is not satisfied, the physical quantity change of the current grid cell is multiplied by a relaxation factor less than 1, and the current grid cell is set to a dormant state. If both are satisfied, the relaxation factor before the change of the physical quantity is 1.0, and the current grid cell is set to active ;

[0013] Step 4: Check whether the number of CFLs on the grid cells can be increased;

[0014] Set the status flags of all grid cells to active ; If the active state flag of grid unit i reaches the preset value m, that is , increase the CFL number of the current grid cell to , and at the same time set the relaxation factor of the change of the physical quantity of all grid cells to 1 and update the physical quantity to .

[0015] A computer program, wherein the computer program enables a computer to execute the above-mentioned accelerated convergence method.

[0016] An electronic device comprises: a processor and a memory; the memory is used to store a computer program, and the processor is used to execute the computer program stored in the memory, so that the electronic device executes the above-mentioned accelerated convergence method.

[0017] A computer-readable storage medium stores a computer program, which implements the above-mentioned accelerated convergence method when executed by a processor.

[0018] A chip includes: a processor, which is used to call and run a computer program from a memory, so that a device equipped with the chip executes the above-mentioned accelerated convergence method.

[0019] A computer program product comprises a computer storage medium storing a computer program, wherein the computer program comprises instructions executable by at least one processor, and the accelerated convergence method is implemented when the instructions are executed by the at least one processor.

[0020] The beneficial effects of the present invention are as follows:

[0021] The present invention establishes a dynamic CFL number adjustment method based on local flow field information, thereby achieving accelerated convergence of aircraft aerodynamic numerical simulation. Combined with the flow field solver HUNS3D, a related program for dynamically adjusting the CFL number is written. The user needs to set the minimum CFL number allowed by the relevant example in the input file of the HUSN3D solver. , Maximum number of CFL , relaxation factor 、CFL number reduction ratio , CFL number magnification ratio , active state preset value m, density, maximum number of consecutive decreases of CFL number k, and relative maximum change of pressure Etc. The implementation of this method is relatively simple. Compared with the original method of setting the CFL number, this method can greatly reduce the program's dependence on the user's experience, improve the program's robustness, greatly accelerate the program's convergence, and shorten the simulation time of the aircraft aerodynamic numerical simulation. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] Figure 1 is a flow chart of the method of the present invention;

[0023] Figure 2is a cloud diagram of the pressure coefficient distribution on the wing upper surface calculated by an embodiment of the present invention;

[0024] Figure 3 It is the distribution cloud diagram of the pressure coefficient on the wing surface calculated by the original fixed CFL number method;

[0025] Figure 4 is a cloud diagram of the pressure coefficient distribution on the lower surface of the wing calculated in an embodiment of the present invention;

[0026] Figure 5 It is the pressure coefficient distribution cloud diagram of the wing lower surface calculated by the original fixed CFL number method;

[0027] Figure 6 It is the residual convergence curve of the original fixed CFL number method and the method of the present invention;

[0028] Figure 7 It is the lift coefficient convergence curve of the original fixed CFL number method and the method of the present invention. DETAILED DESCRIPTION

[0029] The present invention is further described below in conjunction with the accompanying drawings and embodiments.

[0030] like Figure 1 As shown, this embodiment provides a method for dynamically adjusting the local CFL number according to the local flow field change, which reduces the program's dependence on the user's experience and improves the robustness and computational efficiency of the CFD program, including the following contents:

[0031] Step 1: After the current iteration step is completed, check the change of the physical quantity on each grid.

[0032] In practical applications, when the discrete control equations are iteratively solved by numerical simulation methods, the change in the physical quantity on the grid unit i in the flow field at the nth step can be obtained. ,in They represent density, speed in the x, y, and z directions, and total energy respectively.

[0033] 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 initial CFL number, the following parameters related to the dynamic adjustment of the CFL number also need to be set in the input file of the flow field solver: the minimum CFL number allowed in the current case , Maximum number of CFL , relaxation factor 、CFL number reduction ratio , CFL number magnification ratio , active state preset value m, maximum number of consecutive reductions of CFL number k, and relative maximum changes in density and pressure Etc. According to the applicant's experience in using HUNS3D, for most aircraft aerodynamic simulation problems, the above parameters can be set in sequence as follows: , , , , , m=5, k=5, In addition, when the iteration step n=1, the CFL number of each grid unit in the flow field is the same as the initial CFL number. Generally speaking, in order to ensure the stability of numerical calculations during the initial calculation, the initial CFL number should be a smaller value, such as 1.0 (the user can decide based on the working conditions of the example and the quality of the grid).

[0034] Next, we need to judge whether the current iteration has failed based on the change in the physical quantity on each grid unit obtained in the nth iteration. The judgment criterion is to detect the change in all physical quantities to be solved in all grid units. Whether there is NaN or infinity. If the change of the physical quantity on at least one grid cell in the judgment result contains NaN, the current iteration is judged to have failed, and the reduction ratio is multiplied by the CFL number. To reduce the CFL number of all grid cells . It should be explained here why we should not only reduce the CFL number of the grid cells with NaN, but also reduce the CFL number of other grid cells. The main reason is that the HUNS3D solver adopts the LUSGS implicit time advancement solution method. When NaN appears in a grid cell in the flow field, the grid cell with the problem is often not just the current grid cell, but because the disturbance generated in other places is infinitely magnified after the LUSGS iterative solution, which eventually leads to the appearance of NaN in the physical quantity of the current grid cell with the problem. The reason for the problem is that the CFL number of some grid cells in the flow field is too large, but the values ​​of the physical quantities of these grid cells themselves will not appear NaN and will not be detected. Therefore, only reducing the CFL number of the current grid cell cannot solve the problem. A more robust solution is to reduce the CFL number of all grid cells so that the flow field can be solved stably by iterative solution. At the same time, the physical quantity of the current step is no longer updated, and the new reduced CFL number is used to re-advance the solution of the current physical step. If the change in the physical quantity of the nth iteration step calculated next time is still NaN, the CFL number is continuously reduced until the maximum number of times k set by the program to allow continuous reduction is reached or the program fails to exit;

[0035] Step 2: Check whether the density, pressure, and internal energy on the grid cells in step n+1 have negative values.

[0036] When the test in step 1 passes, there is no NaN in all the grid cells in the flow field. At this time, the physical quantity values ​​of all the grid cells in step n+1 are updated to obtain In practical applications, according to the physical quantity to be solved on each grid unit To solve the pressure and internal energy on the current grid unit, the calculation formulas used are , . Then check the density on each grid cell ,pressure and internal energy Whether there is a negative value. If there is no negative value grid cell, its CFL number remains unchanged and the current grid cell is set to active state. ; If the density, pressure or internal energy in grid cell i has a negative value, the CFL number of the current grid cell is reduced. The calculation formula for reducing the CFL number is: , and set all grid cells to sleep mode In addition, the current flow field is not updated, so the physical quantity of the grid unit needs to be changed from the updated value of the n+1th step Back to Then use the new CFL number to solve the flow field variables in the nth step On this basis, the solution of the n+1th iteration step is repeated again.

[0037] In addition to selecting density, the present invention also adds pressure and internal energy as part of the discrimination criteria. The main reason is that according to the actual numerical simulation of aircraft aerodynamics, it is found that pressure and internal energy are the first to diverge. If it is possible to predict in advance whether they will have non-physical solutions, the disturbance in the iteration process can be suppressed in advance by reducing the CFL number. In addition, the HUNS3D solver uses the second-order finite volume method to solve the flow control equation. For the numerical simulation of aircraft transonic flow and supersonic flow, shock waves may exist in the flow field. As a strong discontinuous flow area, shock waves are the most common numerical divergence areas in the numerical simulation process. Because when using the second-order accuracy solution, it is necessary to use a numerical method to solve the gradient of the flow field near the shock wave, which often brings a large error. Even if a limiter is used in the program, the error caused by the gradient solution will cause the flow field near the shock wave to diverge. Therefore, it is generally necessary to reduce the CFL number near the shock wave to improve the stability of the calculation. According to the relationship between the shock wave and the surrounding area, momentum and energy are conserved in these areas. If momentum or energy is used as the main criterion, the divergence characteristics of these areas cannot be effectively captured. On the contrary, density and pressure change dramatically before and after the shock wave, and non-physical solutions are prone to occur during the iteration process. Therefore, using these variables can effectively capture the numerical divergence area in the flow field.

[0038] When reducing the CFL number to continue solving the nth iteration step, it is necessary to ensure that the CFL number of each grid cell is greater than the minimum CFL number allowed by the program, that is, If some grid cells in the flow field still have negative density, pressure or internal energy in the next iteration, the CFL number of the current grid cell will continue to be reduced according to the CFL number reduction criterion, and it will be guaranteed to be greater than the minimum CFL number set by the program until the maximum number of continuous reductions k allowed by the program is reached. At this time, if some grid cells in the flow field still have negative density, pressure or internal energy, then the change in the physical quantity of these grid cells is directly multiplied by a relaxation factor less than 1 to update the physical quantity of the n+1th step. , and perform iterative solution in the n+1th step.

[0039] Step 3: Check whether the change of the physical quantity on the grid unit in step n is too large;

[0040] When the checks of step 1 and step 2 are passed, it means that there is no NaN in the updated physical quantities in the flow field, and there is no grid unit with negative pressure, density or internal energy. Generally speaking, the role of the relaxation factor in the iterative solution process is to reduce the amount of change promoted by each iteration so as to make the iterative process more stable. In the numerical simulation of aircraft aerodynamics, if the CFL number is too large, even if the amount of change in the current iteration step meets the requirements of step 1 and step 2, it may cause the current amount of change to be too large, introduce large errors, and diverge in subsequent iterations. Therefore, in specific application scenarios, it is necessary to further limit the amount of change in physical quantities to enhance the stability of numerical simulation. Based on the physical quantity of the nth step on the grid unit Solve to get the density on each grid cell ,pressure And the n+1th step physical quantity based on the grid unit Solve to get the density on each grid cell ,pressure Then use the formula and Obtain the changes in density and pressure at step n, where the pressure is calculated as , and then judge whether the relative changes of density and pressure exceed the set value of the program, that is, or If at least one of the two is not satisfied, the change in the physical quantity of the current grid unit is multiplied by a relaxation factor z less than 1. ,Right now

[0041]

[0042] As mentioned above, in order to prevent the iteration from being unstable due to the large amount of change each time, a relaxation factor is introduced during the iteration process for those mesh units that do not meet the requirements. Generally speaking, the introduction of a relaxation factor will reduce the convergence speed, so the relaxation factor should not be set too low. It is generally recommended to At the same time, for those grid cells with large changes, the current grid cells are set to dormant state. It should be noted that, due to the large relative changes in the physical quantities of these grid cells, the CFL numbers of these grid cells will be frozen, that is, the CFL numbers of these grid cells will not be increased or decreased in the next iteration step. If all conditions are met, the relaxation factor before the change in the physical quantity is 1.0, and the current grid cell is set to active. ;

[0043] Step 4: Check whether the number of CFLs on the grid cells can be increased;

[0044] When the tests of step 1, step 2 and step 3 are all passed, it means that there is no NaN in the updated physical quantity in the flow field, and there is no grid cell with negative pressure, density or internal energy. At the same time, the relative changes of pressure and density of all grid cells do not exceed the set maximum value. At this time, the status flags of all grid cells are set to active . It should be noted that the active state of the grid unit introduced in the present invention is obtained based on numerical simulation experience. In actual operation, it is found that if the CFL number of the grid unit is frequently increased and decreased alternately, it will cause a sawtooth fluctuation of the CFL number, which will cause calculation instability. Therefore, when the CFL number decreases, the state of the grid unit is set to a dormant state in the next m steps. At this time, even if the grid unit meets the criteria of step 1, step 2 and step 3 at the same time, the CFL number of the current grid unit will not increase. Only when the active state flag of the grid unit i reaches the preset value m, that is, , then the CFL number of the current grid cell can be increased to For the CFL number growth factor According to the actual numerical simulation experience, the recommended value is If the value is too small, the convergence will be slow, and if it is too large, there will be computational stability problems. And the increased CFL number still needs to meet the program setting of not exceeding the maximum CFL number, that is,

[0045] In addition, it should be noted that when a grid cell satisfies the criteria of steps 1, 2, and 3 at the same time, no matter whether the state flag of the grid cell is active or dormant, the relaxation factor of the change amount of the physical quantity of all grid cells is set to 1 and the physical quantity of the grid cell is updated to .

[0046] Example:

[0047] The M6 ​​wing is a CFD calculation model and can also be used to test the performance of the CFD solver. In most calculation results, the incoming flow Mach number is 0.84 and the Reynolds number is 11.7*10 6 , with an initial condition of an angle of attack of 3.06°. A typical feature of the M6 ​​wing under this initial condition is that there is a λ-shaped shock wave on the upper wing surface, with a phenomenon of two shock waves "merging". For such a transonic calculation example, it can be used to better test the effectiveness of the method of the present invention. For comparative verification, the original method uses a fixed CFL number, and the full-field CFL number is 10.0; the parameters of the dynamic CFL number method are set as follows: , , , , , m=5, k=5, , in addition, the initial value of the CFL number is set to 10.0.

[0048] Figure 2 , Figure 3 The pressure coefficient distribution cloud diagram of the upper wing surface obtained by the two methods is given. It can be clearly seen from the figure that the surface flow fields calculated by the two methods are almost the same, and the lambda-shaped shock wave on the upper wing surface is accurately captured. Figure 4 , Figure 5 The lower surface pressure coefficient distribution cloud diagrams calculated by the two methods are almost consistent with each other, indicating that the method of the present invention and the original CFL number method have consistent calculation accuracy. Figure 6 The residual convergence curves of the two methods are given. It can be seen that the method of the present invention can accelerate the convergence of the residual faster than the original method. At the same time, the residual value finally obtained is also significantly smaller than the original method. The residual of the original method is still -6 after 120,000 steps of calculation, while the residual of the method of the present invention is reduced to -12, showing better convergence. In addition Figure 7 The convergence curves of the lift coefficients obtained by the two calculation methods are given. It can be clearly seen that the lift coefficient of the method of the present invention has basically converged at about 10,000 steps, while the original method can barely converge after 120,000 steps of calculation, which shows that the method of the present invention has achieved the effect of accelerated convergence and can significantly improve the calculation efficiency.

Claims

1. A method for accelerating convergence of numerical simulation of aircraft aerodynamics, characterized in that: The steps include: Step 1: After the current iteration step is completed, check the change of the physical quantity on each grid; Step 1-1: After the discrete control equations are solved iteratively for the nth step using the numerical simulation method, the change in the physical quantity on the grid unit i in the flow field at the nth step is obtained. ,in They represent density, speed in the x, y, and z directions, and total energy respectively; Indicates the change of the following variable; Step 1-2: Based on the change in the physical quantity on grid cell i, determine whether there is a non-numeric value NaN in all grid cells. If the change in the physical quantity on at least one grid cell i in the determination result contains NaN, the current iteration is determined to have failed, and the CFL number of all grid cells is reduced. , and the physical quantity of the current iteration step is no longer updated, and the new CFL number is used to advance the solution of the current iteration step again. If the change in the physical quantity of the nth iteration step calculated next time is still NaN, the CFL number is continuously reduced until the maximum number of times k that the program allows the CFL number to be reduced continuously is reached or the program fails to solve and exits; Step 2: Check whether the density, pressure, and internal energy on the grid cell of step n+1 have negative values; Step 2-1: When the test in step 1 passes, there is no NaN in all the grid cells in the flow field. At this time, update the physical quantities of all the grid cells in step n+1 ; Based on the physical quantity of the n+1th step on the grid unit Solve to get the pressure on each grid element and internal energy , check the density on each grid cell ,pressure and internal energy Whether there is a negative value. If the density, pressure or internal energy in grid cell i has a negative value, the CFL number of the current grid cell is reduced and the current grid cell is set to sleep state. For grid cells without negative values, the CFL number remains unchanged and the current grid cell is set to active. ; Step 2-2: If the density, pressure or internal energy of at least one grid cell in the flow field is negative, the current flow field is not updated, that is, the physical quantity of step n is used. Re-solve the n+1th iteration step; when the CFL number is reduced for the kth time in the n+1th iteration, if there are still negative density, pressure or internal energy in the grid cells in the flow field, then directly multiply the physical quantity changes of these grid cells by a relaxation factor less than 1 to update the physical quantity of the n+1th step ; Step 3: Check whether the change of the physical quantity on the grid unit in step n is too large; Physical quantities based on the nth step on the grid cells Solve to get the density on each grid cell ,pressure And the n+1th step physical quantity based on the grid unit Solve to get the density on each grid cell ,pressure , to determine whether the relative changes in density and pressure at the n+1th iteration step exceed the set maximum range, that is, or If at least one of the two conditions is not satisfied, the physical quantity change of the current grid cell is multiplied by a relaxation factor less than 1, and the current grid cell is set to a dormant state. If both are satisfied, the relaxation factor before the change of the physical quantity is 1.0, and the current grid cell is set to active ; Step 4: Check whether the number of CFLs on the grid cells can be increased; Set the status flags of all grid cells to active ; If the active state flag of grid unit i reaches the preset value m, that is , increase the CFL number of the current grid cell to , and at the same time set the relaxation factor of the change of the physical quantity of all grid cells to 1 and update the physical quantity to .

2. A computer program, characterized in that The computer program enables a computer to execute the method as claimed in claim 1.

3. An electronic device, characterized in that: include: Processor and memory; The memory is used to store a computer program, and the processor is used to execute the computer program stored in the memory, so that the electronic device executes the method as claimed in claim 1.

4. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the method as claimed in claim 1 is implemented.

5. A chip, characterized in that: include: A processor, used to call and run a computer program from a memory, so that a device equipped with the chip executes the method as claimed in claim 1.

6. A computer program product, characterized in that The computer program product comprises a computer storage medium storing a computer program, wherein the computer program comprises instructions executable by at least one processor, and when the instructions are executed by the at least one processor, the method according to claim 1 is implemented.

Images