A method for accelerating convergence of numerical simulation of aerodynamics of an aircraft

Abstract

The application discloses an accelerating convergence method for aerodynamic numerical simulation of an aircraft. The method only needs to add a function for judging and correcting the variation of physical quantities in the current step after each iteration, and the original program is slightly changed. The method realizes the target of dynamically adjusting the local CFL number according to the local flow field variation, reduces the experience dependence of the program on users, and improves the robustness and calculation efficiency of the numerical simulation process.

Description

Technical Field

[0001] This invention belongs to the field of aircraft technology, specifically relating to an accelerated convergence method for aerodynamic numerical simulation of aircraft. Background Technology

[0002] Computational Fluid Dynamics (CFD) is a technique that uses numerical methods and algorithms to analyze and solve fluid flow problems. It is widely used in the aerospace field and plays an irreplaceable role in aircraft aerodynamic design. During flight, aircraft flow often exhibits complex, multi-dimensional, wide-velocity flows, including turbulence, separation, and reattachment. These flows are characterized by inter-dimensional interference and the coexistence of high- and low-speed flows, significantly increasing the difficulty of numerical simulation. Furthermore, wind tunnel testing equipment is costly to operate and has limited ability to simulate real flow fields, especially for multiphysics coupling problems and flow separation issues. Therefore, using numerical simulation to predict and analyze aerodynamic characteristics in aircraft design (such as lift and drag prediction, torque prediction, aerodynamic shape optimization, and dynamic characteristic analysis) can effectively reduce experimental costs and time, and improve the accuracy and reliability of aerodynamic design. To accurately and quickly simulate these flow problems, high-precision, high-fidelity numerical simulation methods and large amounts of high-quality meshes can be used, but these significantly increase the computational load and cost. Therefore, strategies are needed to accelerate the numerical solution process while ensuring its stability. In aerodynamic numerical simulation algorithms, the Courant-Friedrichs-Lewy number (CFL) is a crucial concept. It describes the relationship between the time step and the spatial step, ensuring the stability and accuracy of the numerical calculation and, to a certain extent, determining the efficiency of the simulation. A large CFL means that information travels a greater distance across the grid within each time step, potentially leading to instability in the numerical solution. Conversely, a small CFL, while improving the stability of the solution, increases computation time and cost. The physical meaning of the CFL is that it represents the ratio of the distance fluid information travels across the grid within a time step to the grid cell size. In explicit time-progression methods, the CFL must be less than or equal to 1 to guarantee numerical stability, as explicit methods require that the information travel no more than one grid cell within each time step. In implicit time-progression methods, the CFL number is generally not required to meet a certain condition. However, practical experience shows that the CFL number cannot be infinitely large. This mainly depends on the local flow field and mesh. A larger local CFL number results in a longer local time step, which can accelerate convergence and save computation time. Therefore, in practical applications, an appropriate CFL number needs to be selected based on the specific problem and the characteristics of the numerical method. This presents a dilemma in choosing the CFL number. Generally, it is desirable to accelerate the convergence of numerical iterations by providing a large CFL number, but a large CFL number will increase the local time step, leading to divergence in the flow field iteration calculation.Because of this characteristic of the CFL number, users generally need to manually adjust the CFL number according to the convergence of the iterations during numerical simulation in order to balance the efficiency and stability of the numerical simulation, which often depends heavily on the user's experience. Summary of the Invention

[0003] To overcome the shortcomings of existing technologies, this invention provides an accelerated convergence method for aerodynamic numerical simulation of aircraft. It only requires adding a function after each iteration to judge and correct the changes in the physical quantities of the current step, with minimal changes to the original program. This method achieves the goal of dynamically adjusting the local CFL number according to local flow field changes, reduces the program's dependence on the user's experience, and improves the robustness and computational efficiency of the numerical simulation process.

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

[0005] Step 1: After the current iteration step is completed, examine the change in physical quantities on each grid.

[0006] Step 1-1: After iteratively solving the discrete control equations using numerical simulation in the nth step, obtain the changes in physical quantities on grid cell i in the flow field in the nth step. ,in These represent density, velocity in the x, y, and z directions, and total energy, respectively. Indicates the amount of change in the following variable;

[0007] Step 1-2: Based on the change in physical quantity at grid cell i, determine whether there is a non-numerical quantity NaN in all grid cells. If the result shows that at least one grid cell i has a change in physical quantity containing NaN, then the current iteration is considered to have failed, and the CFL number of all grid cells is reduced. Meanwhile, the physical quantity of the current iteration step is no longer updated. The current iteration step is solved again using a new CFL number. If the change of the physical quantity of the nth iteration step is still NaN, the CFL number is reduced until the maximum number of times the CFL number can be continuously reduced, k, or the program fails to solve and exits.

[0008] Step 2: Check whether the density, pressure, and internal energy of the grid cells in step n+1 are negative;

[0009] Step 2-1: After the check in Step 1 passes, there are no NaN cells in any of the mesh elements in the flow field. At this point, update the physical quantities of all mesh elements in step n+1. ; The physical quantity at step n+1 based on the grid cell The pressure on each grid cell is obtained by solving the problem. and internal energy Check the density on each grid cell. ,pressure and internal energy If negative values ​​exist in the density, pressure, or internal energy of grid cell i, then the CFL number of the current grid cell is reduced, and the current grid cell is set to a dormant state. For grid cells without negative values, their CFL number remains unchanged, and the current grid cell is set to active state. ;

[0010] Step 2-2: If the density, pressure, or internal energy of at least one grid cell in the flow field has a negative value, the current flow field will not be updated; that is, the physical quantity from step n will be used. Repeat the solution for the (n+1)th iteration. If, after decreasing the CFL number for the (n+1)th iteration, the grid cells in the flow field still have negative density, pressure, or internal energy, then multiply the changes in these grid cell physical quantities by a relaxation factor less than 1 to update the physical quantities for the (n+1)th iteration. ;

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

[0012] The physical quantity at step n based on the grid cell The density on each grid cell is obtained by solving. ,pressure and the (n+1)th step physical quantity based on the grid cell The density on each grid cell is obtained by solving. ,pressure Determine whether the relative changes in density and pressure in the (n+1)th iteration step exceed the set maximum range, i.e. or If at least one of the two conditions is not met, the change in the physical quantity 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 conditions are met, the relaxation factor before the change in the physical quantity is 1.0, and the current mesh cell is set to active state. ;

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

[0014] Set the status flag of all grid cells to active state. If the active state indicator of grid cell i reaches the preset value m, i.e. Increase the CFL number of the current grid cell by Simultaneously, the relaxation factor of the change in physical quantities of all mesh elements is set to 1, and the physical quantities are updated to... .

[0015] A computer program that causes a computer to execute the above-described accelerated convergence method.

[0016] An electronic device includes: 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 to enable the electronic device to perform the above-described accelerated convergence method.

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

[0018] A chip includes a processor for retrieving and running a computer program from a memory, causing a device on which the chip is mounted to perform the aforementioned accelerated convergence method.

[0019] A computer program product includes a computer storage medium storing a computer program, the computer program including instructions executable by at least one processor, which, when executed by the at least one processor, implement the above-described accelerated convergence method.

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

[0021] This invention establishes a dynamic CFL number adjustment method based on local flow field information, thereby achieving accelerated convergence in aircraft aerodynamic numerical simulations. Combined with the flow field solver HUNS3D, a program for dynamically adjusting the CFL number was developed. Users need to set the minimum allowed CFL number for the relevant calculation example in the input file of the HUNS3D solver program. Maximum number of CFLs Relaxation factor CFL number reduction ratio CFL number magnification ratio The active state preset value m, the maximum number of times the density and CFL number continuously decrease k, and the maximum relative change in pressure. This method is relatively simple to implement. Compared with the original method of fixed CFL number, this method can significantly reduce the program's dependence on the user's experience, improve the program's robustness, and greatly accelerate the program's convergence, thus shortening the simulation time of aircraft aerodynamic numerical simulation. Attached Figure Description

[0022] Figure 1 This is a flowchart of the method of the present invention;

[0023] Figure 2This is a cloud map showing the distribution of pressure coefficient on the upper surface of the wing calculated according to an embodiment of the present invention;

[0024] Figure 3 This is a cloud map of the pressure coefficient distribution on the upper surface of the wing calculated using the original constant CFL number method;

[0025] Figure 4 This is a cloud map showing the distribution of the pressure coefficient on the lower surface of the wing calculated according to an embodiment of the present invention;

[0026] Figure 5 This is a cloud map of the pressure coefficient distribution on the lower surface of the wing calculated using the original constant CFL number method;

[0027] Figure 6 These are the residual convergence curves of the original fixed CFL number method and the method of this invention;

[0028] Figure 7 These are the lift coefficient convergence curves of the original constant CFL number method and the method of this invention. Detailed Implementation

[0029] The present invention will be further described below with reference to the accompanying drawings and embodiments.

[0030] like Figure 1 As shown, this embodiment provides a method for dynamically adjusting the local CFL number based on local flow field changes, reducing the program's dependence on user experience and improving the robustness and computational efficiency of the CFD program. The method includes the following:

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

[0032] In practical applications, after performing the nth iterative solution of the discrete control equations using numerical simulation, the changes in physical quantities on grid cell i in the flow field at the nth step can be obtained. ,in These represent density, velocity 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 solutions, 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 example. Maximum number of CFLs Relaxation factor CFL number reduction ratio CFL number magnification ratio The active state preset value m, the maximum number of consecutive decreases in CFL number k, and the maximum relative change in density and pressure. Based on the applicant's experience using HUNS3D, for most aircraft aerodynamic simulation problems, the above parameters can be set sequentially as follows: , , , , m=5, k=5 Furthermore, when the iteration step n=1, the CFL number of each grid cell in the flow field is the same, which is the initial CFL number. Generally speaking, in order to ensure the stability of numerical calculation in the initial calculation, the initial CFL number should be a small value, such as 1.0 (the user can decide according to the working conditions of the example and the grid quality).

[0034] Next, we need to determine whether the current iteration has failed based on the changes in the physical quantities in each grid cell obtained from the nth iteration. The criterion is to detect the changes in all the physical quantities to be solved in all grid cells. If NaN or infinity exists, and the change in physical quantity on at least one grid cell in the determination result contains NaN, then the current iteration is determined to have failed, and the result is multiplied by the CFL number reduction ratio. To reduce the number of CFLs in all grid cells It's important to explain why we need to reduce not only the CFL number of the grid cells containing NaN, but also the CFL number of other grid cells. The main reason is that the HUNS3D solver uses the LUSGS implicit time-progression solution method. When a grid cell in the flow field contains NaN, the problem often isn't just with that particular grid cell. Instead, it's because disturbances from other areas are amplified infinitely after being solved iteratively by LUSGS, ultimately causing NaN values ​​in the physical quantities of the problematic grid cell. The cause of the problem is that some grid cells in the flow field have excessively high CFL numbers, but the physical quantities of these grid cells themselves don't contain NaN values ​​and therefore won't be detected. So, simply reducing the CFL number of the current grid cell cannot solve the problem. A more robust approach is to reduce the CFL number of all grid cells to ensure stable iterative solutions for the flow field. Meanwhile, the physical quantity of the current step is no longer updated. The current physical step is solved again using the new reduced CFL number. If the change of the physical quantity of the nth iteration step is still NaN, the CFL number is reduced until the maximum number of consecutive reductions allowed by the program is reached, k, or the program fails to solve and exits.

[0035] Step 2: Check whether there are negative values ​​for density, pressure, and internal energy on the grid cell in step n+1.

[0036] Once the check in step 1 passes, there is no NaN in any of the mesh elements in the flow field. At this point, the physical values ​​of all mesh elements in step n+1 are updated. In practical applications, the physical quantity to be solved is determined based on each grid cell. The formulas used to calculate the pressure and internal energy on the current mesh cell are as follows: , Then check the density on each grid cell. ,pressure and internal energy If a negative value exists, the CFL number of the grid cell remains unchanged, and the current grid cell is set to active status. If the density, pressure, or internal energy in grid cell i is negative, the CFL number of the current grid cell is reduced. The formula for reducing the CFL number is as follows: At the same time, set all grid cells to sleep mode. Furthermore, since the current flow field is not being updated, the physical quantities of the mesh elements need to be changed from the values ​​updated in step n+1. Back Then, the flow field unsolved variables are obtained using the new CFL number in step n. Then, solve the problem again for the (n+1)th iteration step based on the previous solution.

[0037] In addition to density, this invention also incorporates pressure and internal energy as criteria. This is primarily because actual aerodynamic numerical simulations of aircraft reveal that pressure and internal energy are the first to diverge. If non-physical solutions can be predicted in advance, disturbances during the iteration process can be suppressed by reducing the CFL number. Furthermore, the HUNS3D solver uses the second-order finite volume method to solve the flow control equations. For numerical simulations of transonic and supersonic flows in aircraft, shock waves may exist in the flow field. As a strongly discontinuous flow region, shock waves are the most common numerical divergence region in numerical simulations. Because using second-order precision solutions requires numerically calculating the gradient of the flow field near the shock wave, this often introduces significant errors. Even with limiters in the program, errors from gradient calculation can still 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 computational stability. According to the relationship before and after the shock wave, momentum and energy are conserved in these regions. If momentum or energy is used as the main criterion, the divergence characteristics in these regions cannot be effectively captured. Conversely, density and pressure change drastically 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 regions in the flow field.

[0038] When decreasing the CFL number to continue solving in the nth iteration step, it is necessary to ensure that the CFL number of each grid cell is greater than the minimum allowed CFL number set by the program, i.e. If, in the next iteration, some grid cells in the flow field still have negative density, pressure, or internal energy, then the CFL number of the current grid cell continues to decrease according to the CFL number reduction criterion, ensuring it is greater than the minimum CFL number set by the program, until the maximum number of consecutive decreases in the CFL number allowed by the program, k, is reached. If, at this point, some grid cells in the flow field still have negative density, pressure, or internal energy, then the changes in the physical quantities of these grid cells are directly multiplied by a relaxation factor less than 1 to update the physical quantities of the (n+1)th step. And perform the iterative solution in the (n+1)th step.

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

[0040] Once steps 1 and 2 pass the checks, it means that the updated physical quantities in the flow field do not contain NaN (NaN) values, nor do they contain grid cells with negative pressure, density, or internal energy. Generally, the relaxation factor in iterative solutions reduces the amount of change advanced in each iteration, thus making the iteration process more stable. In numerical simulations of aircraft aerodynamics, if the CFL number is too large, even if the change in the current iteration step meets the requirements of steps 1 and 2, it may still lead to an excessively large change, introducing significant errors that diverge in subsequent iterations. Therefore, in specific applications, it is necessary to further limit the changes in physical quantities to enhance the stability of the numerical simulation. The nth step physical quantity based on the grid cells... The density on each grid cell is obtained by solving. ,pressure and the (n+1)th step physical quantity based on the grid cell The density on each grid cell is obtained by solving. ,pressure Then use the formula. and Obtain the changes in density and pressure at step n, where the pressure is calculated using the following formula: Then, it is determined whether the relative changes in density and pressure exceed the program's set values. or If at least one of the two conditions is not met, then the change in the physical quantity of the current mesh element is multiplied by a relaxation factor z less than 1. ,Right now

[0041]

[0042] As mentioned earlier, to prevent iterative instability caused by excessively large changes in each iteration, a relaxation factor is introduced during the iteration process for these unsatisfactory mesh cells. Generally, the introduction of a relaxation factor reduces the convergence speed, so the relaxation factor should not be set too low; it is generally recommended to... Simultaneously, for these grid cells with large variations, the current grid cell is set to a dormant state. It should be noted that, due to the relatively large changes in the physical quantities of these mesh elements, the CFL numbers of these mesh elements will be frozen, meaning they will neither increase nor decrease in the next iteration. If all conditions are met, the relaxation factor before the change in physical quantities is set to 1.0, and the current mesh element is set to an active state. ;

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

[0044] Once steps 1, 2, and 3 have all passed the checks, it means that the updated physical quantities in the flow field do not contain NaN, nor do they contain any grid cells with negative pressure, density, or internal energy. Furthermore, the relative changes in pressure and density of all grid cells do not exceed the set maximum values. At this point, the status flags of all grid cells are set to active. It should be noted that the active state of the grid cells introduced in this invention is based on numerical simulation experience. In practice, it was found that frequently alternating increases and decreases in the CFL number of grid cells can lead to sawtooth-shaped fluctuations in the CFL number, causing computational instability. Therefore, when the CFL number decreases, the grid cell is set to a dormant state for the next m steps. Even if the grid cell simultaneously meets the criteria of steps 1, 2, and 3, the CFL number of the current grid cell will not increase until the active state indicator of grid cell i reaches a preset value m. Then the CFL number of the current grid cell can be increased to... For CFL number growth factor Based on actual numerical simulation experience, it is recommended to set the value to [value to be inserted here]. If the value is too small, convergence will be slow; if it is too large, computational stability issues may arise. Furthermore, the increased number of CFLs must still meet the program's setting of not exceeding the maximum number of CFLs, i.e.

[0045] Furthermore, it should be noted that when a mesh element simultaneously satisfies the criteria in steps 1, 2, and 3, regardless of whether the mesh element's state is active or dormant, the relaxation factor of the changes in the physical quantities of all mesh elements is set to 1, and the physical quantities of the mesh elements are updated. .

[0046] Example:

[0047] The M6 ​​wing is a computational benchmark for CFD and can also be used to test the performance of CFD solvers. Most calculations use an incoming Mach number of 0.84 and a Reynolds number of 11.7 × 10⁻⁶. 6 The initial assumption is an angle of attack of 3.06°. A typical characteristic of the M6 ​​wing under this initial assumption is a λ-shaped shock wave on the upper wing surface, exhibiting a phenomenon of two shock waves "merging". This transonic example can be used to effectively test the effectiveness of the method of this invention. For comparative verification, the original method uses a fixed CFL number, with a full-field CFL number of 10.0; the parameter settings for the dynamic CFL number method are 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 contour maps of the upper surface of the wing obtained by two methods are presented. It can be clearly seen from the figures that the surface flow fields calculated by the two methods are almost identical, and the λ-shaped shock wave of the upper wing surface is accurately captured. Figure 4 , Figure 5 The diagrams show the distribution contours of the lower surface pressure coefficients calculated by the two methods. The results obtained by the two methods are almost identical, indicating that the method of the present invention has the same calculation accuracy as the original fixed CFL number method. Figure 6 The residual convergence curves of the two methods are presented. It can be seen that the method of the present invention can accelerate the convergence of residuals faster than the original method, and the final residual value obtained is also significantly smaller than that of the original method. After 120,000 calculation steps, the residual of the original method is still at the power of -6, while the residual of the method of the present invention is reduced to the power of -12, showing better convergence. In addition, Figure 7 The convergence curves of the lift coefficient 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 after about 10,000 steps, while the original method can barely converge after 120,000 steps. This shows that the method of the present invention has achieved the effect of accelerating convergence and can significantly improve the calculation efficiency.

Claims

1. An accelerated convergence method for numerical simulation of aircraft aerodynamics, characterized in that, Includes the following steps: Step 1: After the current iteration step is completed, examine the change in physical quantities on each grid. Step 1-1: After iteratively solving the discrete control equations using numerical simulation in the nth step, obtain the changes in physical quantities on grid cell i in the flow field in the nth step. ,in These represent density, velocity in the x, y, and z directions, and total energy, respectively. Indicates the amount of change in the following variable; Step 1-2: Based on the change in physical quantity at grid cell i, determine whether there is a non-numerical quantity NaN in all grid cells. If the result shows that at least one grid cell i has a change in physical quantity containing NaN, then the current iteration is considered to have failed, and the CFL number of all grid cells is reduced. Meanwhile, the physical quantity of the current iteration step is no longer updated. The current iteration step is solved again using a new CFL number. If the change of the physical quantity of the nth iteration step is still NaN, the CFL number is reduced until the maximum number of times the CFL number can be continuously reduced, k, or the program fails to solve and exits. Step 2: Check whether the density, pressure, and internal energy of the grid cells in step n+1 are negative; Step 2-1: After the check in Step 1 passes, there are no NaN cells in any of the mesh elements in the flow field. At this point, update the physical quantities of all mesh elements in step n+1. ; The physical quantity at step n+1 based on the grid cell The pressure on each grid cell is obtained by solving the problem. and internal energy Check the density on each grid cell. ,pressure and internal energy If negative values ​​exist in the density, pressure, or internal energy of grid cell i, then the CFL number of the current grid cell is reduced, and the current grid cell is set to a dormant state. For grid cells without negative values, their CFL number remains unchanged, and the current grid cell is set to active state. ; Step 2-2: If the density, pressure, or internal energy of at least one grid cell in the flow field has a negative value, the current flow field will not be updated; that is, the physical quantity from step n will be used. Repeat the solution for the (n+1)th iteration. If, after decreasing the CFL number for the (n+1)th iteration, the grid cells in the flow field still have negative density, pressure, or internal energy, then multiply the changes in these grid cell physical quantities by a relaxation factor less than 1 to update the physical quantities for the (n+1)th iteration. ; Step 3: Check whether the change in physical quantities on the grid cell in step n is too large; The physical quantity at step n based on the grid cell The density on each grid cell is obtained by solving. ,pressure and the (n+1)th step physical quantity based on the grid cell The density on each grid cell is obtained by solving. ,pressure Determine whether the relative changes in density and pressure in the (n+1)th iteration step exceed the set maximum range, i.e. or If at least one of the two conditions is not met, the change in the physical quantity 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 conditions are met, the relaxation factor before the change in the physical quantity is 1.0, and the current mesh cell is set to active state. ; Step 4: Check if the number of CFLs on the mesh cells has been increased; Set the status flag of all grid cells to active state. If the active state indicator of grid cell i reaches the preset value m, i.e. Increase the CFL number of the current grid cell by Simultaneously, the relaxation factor of the change in physical quantities of all mesh elements is set to 1, and the physical quantities are updated to... .

2. 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 to cause the electronic device to perform the method as described in claim 1.

3. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the method as described in claim 1.

4. A chip, characterized in that, include: A processor for retrieving and running a computer program from memory, causing a device on which the chip is mounted to perform the method as described in claim 1.

5. A computer program product, characterized in that, The computer program product includes a computer storage medium storing a computer program, the computer program including instructions executable by at least one processor, which, when executed by the at least one processor, implement the method as described in claim 1.

Images