High efficiency calculation method for complex turbulent flow based on pseudo-non-steady time marching

Abstract

The application discloses a high-efficiency calculation method for complex turbulent flow based on pseudo-non-steady time advancement, which comprises the following steps: step 1, establishing control equations for non-steady turbulent flow simulation; step 2, determining a solution region, dividing grid units in the region, and assigning initial values of flow fields to all the grid units; step 3, using grid unit variables at a current t n time to calculate non-viscous flux and viscous flux of all the grid units, and summing the fluxes to obtain a time derivative value of the grid unit flow field variables at the current t n time; step 4, using the time derivative value of the flow field variables at the current t n time to solve flow field variable values at a next t n+1 time; step 5, judging whether the next t n+1 time is greater than or equal to a calculation termination t end time; if yes, the method proceeds to step 6; otherwise, n is set to n+1, and the method returns to step 2; and step 6, outputting final results. The method can effectively ensure the diagonal dominance of a matrix, ensure the stability of numerical calculation, improve the calculation efficiency, and reduce the calculation period and the calculation cost.

Description

[Technical Field]

[0001] This invention belongs to the fields of computational mathematics and computational fluid dynamics, and specifically relates to a high-efficiency calculation method for complex turbulent flows based on pseudo-unsteady time propagation. [Background Technology]

[0002] As a complex flow phenomenon in nature and engineering, turbulent flow often exists in various complex forms around aircraft, such as bottom flow, jet interference, separated flow, and shock wave / boundary layer interference. It significantly impacts the aircraft's bottom drag, vortex structure, inlet flow quality, and engine fuel mixing and combustion characteristics. During subsonic to hypersonic flight, especially at high angles of attack, the aircraft's side edges, local fuselage protrusions, and control surfaces induce large-scale flow separation, forming complex separated vortex structures. These separated vortices exhibit strong unsteady characteristics and often significantly affect aircraft performance. Therefore, accurate and detailed simulation of unsteady complex turbulent flows is crucial for the aerodynamic design, control system design, and propulsion system design of aircraft.

[0003] However, achieving high-precision, low-cost prediction of complex turbulent flows is extremely difficult. Current high-precision simulation methods for unsteady turbulence include Unsteady Reynolds-averaged Navier-Stokes (URANS), Large Eddy Simulation (LES), hybrid RANS / LES, and Direct Numerical Simulation (DNS). Among these, LES and DNS methods are computationally too intensive and currently lack practical application capabilities for complex engineering shapes. Furthermore, practical experience shows that computational efficiency and convergence efficiency for unsteady complex turbulence are also pressing issues. Research indicates that to accurately calculate high Reynolds (Re) number unsteady complex turbulent flows, current numerical calculations require large aspect ratio meshes near the wall. The extremely low mesh size in the wall normal direction limits the maximum allowable time step, making it difficult to meet the time efficiency requirements of engineering designs.

[0004] In recent years, much work has been done both domestically and internationally to improve the computational efficiency of numerical methods. Current time integration schemes can be categorized into explicit and implicit schemes based on their implementation method, and into spatiotemporally decoupled discretization (also known as semi-discretization) and spatiotemporally coupled discretization (also known as full discretization) based on their spatiotemporal discretization method. Most of the aforementioned high-precision simulations of unsteady turbulence employ spatiotemporally decoupled discretization, which has the advantage of considering only spatial discretization during scheme construction while directly applying mature ordinary differential equation solution methods in time, such as the explicit multi-stage Runge-Kutta (RK) method. Traditional explicit schemes have relatively small time steps. Besides developing explicit schemes with large time steps, researchers have also made valuable explorations in recent years regarding implicit schemes applicable to unsteady computation. However, overall, unsteady computation still faces many difficulties in both theoretical and numerical research, making it one of the most challenging frontier research topics in fluid mechanics. Compared to the more mature and established methods for calculating the time of steady flow, the methods for calculating the time of unsteady flow are still in their infancy and many issues require further research. The main difficulties and research focus are the second-order precision time discretization scheme, the massive computation of unsteady flow, and the verification and confirmation of the calculation results of unsteady flow.

[0005] This invention primarily addresses the problems of high computational cost and low efficiency in traditional URANS and hybrid RANS / LES methods for unsteady flow. Based on the pseudo-unsteady time-stepping concept, it layers the computational grid, with different sub-iteration time steps for different grid layers. By reducing the time discretization precision and decreasing the time step size, the computational efficiency of turbulent flow simulation is improved. Compared to existing numerical simulation methods, this invention is constructed based on the implicit time-stepping and layered computation approach, building upon the dual-time-step discretization method. It inherits the advantages of implicit time schemes, such as good stability and independence from time step limitations, effectively improving unsteady computational efficiency when dealing with grids with large aspect ratios used in turbulence simulations. Furthermore, by reducing the time step size, the diagonal dominance of the matrix is ​​effectively guaranteed, further ensuring the stability of the numerical simulation. Therefore, this invention effectively improves the computational efficiency of complex unsteady turbulence while reducing computational cost, without significantly affecting computational accuracy. [Summary of the Invention]

[0006] The technical problem this invention aims to solve is to overcome the shortcomings of existing technologies and provide a time-progression method suitable for unsteady simulation of complex turbulent flows. For URANS and hybrid RANS / LES methods, this invention considers using a high aspect ratio mesh near the aircraft wall. Based on the pseudo-time-progression concept, it effectively ensures the diagonal dominance of the matrix by utilizing mesh layering and reducing the time step, thus guaranteeing numerical computation stability, thereby improving computational efficiency and reducing computational cycle and cost.

[0007] The technical solution adopted by this invention to solve the above-mentioned technical problems is as follows: a high-efficiency calculation method for complex turbulent flow based on pseudo-unsteady time propagation, comprising the following steps:

[0008] Step 1: Establish the governing equations for unsteady turbulence simulation. The governing equations adopted in this invention are the Reynolds-averaged Navier-Stokes (RANS) equations for compressible flow. The equations are written in a decoupled form of the mainstream equations and the turbulence model equations. The mainstream equations are written as equation (1), and the turbulence model equations are written as equation (2).

[0009]

[0010]

[0011] In formula (1) This represents the partial derivatives of the variables, where Q is a conserved variable, and F, G, and H are the inviscid fluxes in the x, y, and z directions, respectively. μ G μ and H μ Let μ be the viscous flux in the x, y, and z directions, where μ represents viscosity, t is the time term, x, y, and z are the spatial terms, e is the energy, τ is the viscous stress, and q is the heat flux. In equation (2), X represents the variables of the turbulence model equation, P is the production source term of the turbulence model equation, ε is the destruction source term of the turbulence model equation, and D is the diffusion term of the turbulence model equation.

[0012] Step 2: Determine the solution domain based on the turbulent flow to be simulated, divide the domain into grid cells, and assign initial flow field values ​​to all grid cells, including density ρ, velocity u, v, w, pressure p, and initial values ​​of the variable X in the turbulence model equation.

[0013] Step 3, using the current t n Given the grid cell variables at time t, calculate the inviscid and viscous fluxes for all grid cells, then sum the inviscid and viscous fluxes to obtain the current flux t. n The time derivative of the flow field variables in the grid cells at time t;

[0014] Step 4, using the current t n The next step is to solve the time derivative of the flow field variables at time t. n+1 The flow field variable values ​​at time t, and update the turbulent viscosity information in the mainstream equation, and synchronously update the variable information of the mainstream equation in the turbulence model equation;

[0015] Step 5, determine the next step t n+1 Is the time greater than or equal to the time at which the calculation terminates (t)? endIf the condition is met, proceed to step 6; otherwise, let n = n + 1 and return to step 2.

[0016] Step 6: Output the final result.

[0017] Step 4 utilizes the current t n The next step is to solve the time derivative of the flow field variables at time t. n+1 The flow field variable values ​​at time t are calculated using an unsteady implicit time-propagation method, expressed as follows:

[0018]

[0019] In the formula, R represents the sum of the inviscid flux and viscous flux calculated in step 3, the superscript m represents the iteration step of the unsteady internal iteration process, and the superscripts n and n-1 represent the current t. n Time and previous step t n-1 The result of the flow field variable calculation at time t, where Δt represents t n+1 Time and t n The size of the time step at time t, i.e., Δt = t n+1 -t n Its specific size is calculated based on a given number of CFLs.

[0020] Furthermore, the unsteady implicit time-progression calculation method employs pseudo-unsteady time-progression, setting the number of iterations in the unsteady internal iteration process to 1, and setting the time step Δt to depend on the time step of the mesh layering results. Where i represents the grid level, and k represents the number of iterations within the i-th grid level. The solution equations for the pseudo-unsteady time-progression method are written as follows:

[0021]

[0022] In the formula This represents the flux residual of the mainstream equation during any iteration of the i-th grid layer. The Jacobian matrix representing the main equation can be decomposed into a non-negative eigenvalue matrix A based on its eigenvalues. + and the non-positive eigenvalue matrix A - Summation form

[0023]

[0024] Meanwhile, assuming that the time step is sufficiently small, the increments of the solutions at two adjacent time steps are approximately equal, i.e., ΔQ = Q. n+1 -Q n ≈Q n -Q n-1 Therefore, equation (5) can be further written as

[0025]

[0026] Therefore, at any grid level i, the equations to be solved are reduced to pseudo-unsteady implicit discrete equations.

[0027] Furthermore, the pseudo-unsteady time-progression method employs grid-based hierarchical computation, using a time step size on the grid at level i. Mesh layering should ensure that meshes with similar slenderness ratios are grouped into the same layer. To determine mesh layering, a suitable direction can be manually selected from the three directions of the structural mesh (in the case of 3D calculation), the slenderness ratio of the mesh is calculated along that direction, and meshes with similar slenderness ratios are grouped into the same layer along that direction. In this invention, the slenderness ratios of the first layer of meshes near the aircraft wall and the outermost layer of meshes far from the wall are first calculated. The difference between the two is a multiple of n, therefore the growth rate of each layer can be roughly determined as follows: The number of grid layers is I. Considering the complexity of actual grid structures, the growth ratio can be adjusted appropriately to avoid the situation where the number of grid layers in a certain layer is too large or too small. Then, all grid layers in layer i (i = 1, ..., I) are given a feature scale L. i Typical feature scales can be

[0028] L i ={L1=0.01,..L i =0.1...,L I =1} (7)

[0029] t n+1 Time and t n The time step Δt is calculated using a given CFL number. For unsteady calculations, the time step of all grid cells needs to be unified to the smallest time step in the flow field grid cells.

[0030]

[0031] In the formula, j traverses all mesh cells, and V j Represents the volume of the mesh cell. Represents the largest eigenvalue of the grid cell. c represents the area vector of the grid cells. diff and ν diff These are characteristic coefficients related to viscosity. The time step is used on the grid at level i. The corresponding number of inner iterations k can be obtained by giving the feature scale L. i And the minimum time step Δt are obtained

[0032]

[0033] Furthermore, the pseudo-unsteady time-progression method solves the implicit discrete equation (6) at each level of the grid. According to the upwind principle, backward and forward differences can be performed on each term of equation (6) to obtain...

[0034]

[0035] To ensure convergence performance, an approximate Jacobian matrix is ​​constructed to guarantee the diagonal dominance property. Specifically, it can be constructed as follows:

[0036]

[0037] Therefore A + -A - By converting the data to scalars, matrix fitting operations are avoided. The solution process iterates according to the grid hierarchy, first iterating through the first level (1 / L1), then the second level (1 / L2), and so on, until finally iterating through the first level (1). After all grid levels have been iterated, the flow field variables in the main flow equations and the turbulence model equations are updated.

[0038] Furthermore, the pseudo-unsteady time-progression method for solving the implicit discrete equations of the turbulence model differs from the mainstream equations. The turbulence model equations are solved directly over all grids, without using a hierarchical grid.

[0039]

[0040] In equation (12), the subscript X represents the variables of the turbulence model equation. It can be seen that the update of the variables of the turbulence model equation is consistent with the update of the first level of the mainstream equation. This avoids the solution of the turbulence model equation from consuming more computational resources. The updated turbulence model variables are used to update the turbulence viscosity in the mainstream equation.

[0041] Furthermore, the pseudo-unsteady time-progression method uses a single-step LU diagonal scan in the implicit iteration. For the mainstream equation, a single-step LU diagonal scan is used in the inner iteration of each grid level; for the turbulence model equation, a single-step LU diagonal scan is used in the iteration of all grids. Without loss of generality, taking the iteration of the mainstream equation on the i-th grid level as an example, a single-step forward scan is used to obtain the intermediate solution ΔQ in the inner iteration. *

[0042]

[0043] Matrix A that needs to be calculated in the forward scan + From formula (11), the updated solution ΔQ obtained during the scanning process is... * It can be directly stored in the right-hand term vector This further reduces memory consumption.

[0044] After the forward scan is completed, the backward scan begins to obtain ΔQ. n

[0045]

[0046] Matrix A that needs to be calculated in the backward scan - From formula (11), after the backward scan is completed, the updated solution ΔQ for this iteration is obtained. n This is used for the next iterative calculation.

[0047] The advantages of this invention compared to the prior art are:

[0048] 1. The high-efficiency calculation method for complex turbulent flow based on pseudo-unsteady time propagation used in this invention can use a hierarchical iteration method according to the aspect ratio of the calculation grid, which can effectively avoid the overall unsteady calculation efficiency reduction caused by the excessively small time step in the grid region with a large aspect ratio.

[0049] 2. The pseudo-unsteady implicit time-progression method adopted in this invention inherits the characteristics of good implicit discrete stability and no limitation on time step size.

[0050] 3. This invention employs the same implicit time discretization for both mainstream equations and turbulence model equations, resulting in a unified computational method. Compared to traditional unsteady implicit propagation methods, it avoids the ineffective convergence of sub-iterations and is more suitable for computational methods based on RANS equations, such as URANS and hybrid RANS / LES.

[0051] 4. This invention uses a method of reducing time discretization precision and a smaller time step size, so that the discrete matrix has better diagonal dominance characteristics, which can further enhance computational stability.

[0052] 5. The single-step LU diagonal scan used in the implicit iteration of this invention can avoid storing the block diagonal matrix and the solution vector of the previous step, saving memory overhead and being easier to implement in programming. [Attached Image Description]

[0053] Figure 1 This is a flowchart of the overall process of the method of the present invention.

[0054] Figure 2a Figures 1 and 2b are partial flowcharts of the high-efficiency calculation method for complex turbulent flow based on pseudo-unsteady time propagation proposed in this invention.

[0055] Figure 3 The image shown is a mesh diagram near the airfoil surface in an embodiment of the present invention.

[0056] Figure 4 The image shown is a cloud map of the instantaneous eddy viscosity distribution of the simulation results in an embodiment of the present invention.

[0057] Figure 5 The figure shows a comparison of the results of predicting lift and drag characteristics using the method of this invention and the traditional dual-time-step method.

Detailed Implementation Methods

[0058] The invention will be further described in detail below with reference to the flowchart of the method:

[0059] like Figure 1 As shown, this invention relates to a high-efficiency calculation method for complex turbulent flows based on pseudo-unsteady time propagation, and the specific implementation steps are as follows:

[0060] Step 1: Based on the unsteady turbulent flow to be simulated, select a turbulence model and establish the RANS equations for unsteady turbulence simulation. Write the equations in a decoupled form of the mainstream equation and the turbulence model equation. The mainstream equation is written as equation (1), and the turbulence model equation is written as equation (2).

[0061]

[0062]

[0063] In formula (1) This represents the partial derivatives of the variables, where Q is a conserved variable, and F, G, and H are the inviscid fluxes in the x, y, and z directions, respectively. μ G μ and H μ Let μ be the viscous flux in the x, y, and z directions, where μ represents viscosity, t is the time term, x, y, and z are the spatial terms, e is the energy, τ is the viscous stress, and q is the heat flux. In equation (2), X represents the variables of the turbulence model equation, P is the production source term of the turbulence model equation, ε is the destruction source term of the turbulence model equation, and D is the diffusion term of the turbulence model equation.

[0064] Step 2: Determine the solution domain based on the turbulent flow to be simulated, divide the domain into grid cells, and assign initial flow field values ​​to all grid cells, including density ρ, velocity u, v, w, pressure p, and initial values ​​of the variable X in the turbulence model equation.

[0065] Step 3, using the current t n Given the grid cell variables at time t, calculate the inviscid and viscous fluxes for all grid cells, then sum the inviscid and viscous fluxes to obtain the current flux t. n The time derivative of the flow field variables in the grid cells at time t;

[0066] like Figure 2a As shown in Figures 1 and 2, the following are the specific implementation steps of a high-efficiency calculation method for complex turbulent flows based on pseudo-unsteady time propagation:

[0067] Step 4.1, using the current t n Numerical solution of the time derivative of the flow field variables at time t n+1 The flow field variable values ​​at time t are calculated using an unsteady implicit time-propagation method, expressed as follows:

[0068]

[0069] In the formula, R represents the sum of the inviscid flux and viscous flux calculated in step 3, the superscript m represents the iteration step of the unsteady internal iteration process, and the superscripts n and n-1 represent the current t. n Time and previous step t n-1 The result of the flow field variable calculation at time t, where Δt represents t n+1 Time and t n The size of the time step at time t, i.e., Δt = t n+1 -t n Its specific size is calculated based on a given number of CFLs;

[0070] Step 4.2 involves mesh layering. A suitable direction can be manually selected from the three directions of the structural mesh (in the case of 3D calculation). The mesh aspect ratio is calculated along this direction, and meshes with similar aspect ratios are grouped into the same level along this direction. Let I be the number of mesh layers, and let L be the characteristic scale of all meshes at level i (i = 1, ..., I). i Typical feature scales can be

[0071] L i ={L1=0.01,..L i =0.1...,L I =1} (4)

[0072] Step 4.3, calculate t n+1 Time and t n The time step Δt is given at time t. Calculations are performed using a given CFL number. For unsteady calculations, the time step of all grid cells needs to be unified to the smallest time step in the flow field grid cells.

[0073]

[0074] In the formula, j traverses all mesh cells, and V j Represents the volume of the mesh cell. Represents the largest eigenvalue of the grid cell. c represents the area vector of the grid cells. diff and ν diff These are characteristic coefficients related to viscosity.

[0075] Step 4.4, calculate the time step Δt on the grid at level i. i kAnd the corresponding number of inner iterations k, where i represents the grid level and k represents the number of iterations within the i-th grid level. Given a feature scale L... i And the minimum time step Δt are obtained

[0076]

[0077] Step 4.5: Solve equation (3) using pseudo-unsteady time-progression. Set the number of iterations in the unsteady internal iteration process to 1, and set the time step Δt to the time step that depends on the grid layering results. Writing the solution equations for the pseudo-unsteady time-progression method

[0078]

[0079] In the formula This represents the flux residual of the mainstream equation during any iteration of the i-th grid layer. The Jacobian matrix representing the main equation can be decomposed into a non-negative eigenvalue matrix A based on its eigenvalues. + and the non-positive eigenvalue matrix A - Summation form

[0080]

[0081] Meanwhile, assuming that the time step is sufficiently small, the increments of the solutions at two adjacent time steps are approximately equal, i.e., ΔQ = Q. n+1 -Q n ≈Q n -Q n-1 Therefore, equation (8) can be further written as

[0082]

[0083] Therefore, at any grid level i, the equations to be solved are reduced to pseudo-unsteady implicit discrete equations.

[0084] Step 4.6: Solve the implicit discrete equation (9) on each level of the grid. According to the upwind principle, the terms of equation (9) can be obtained by performing backward and forward differences respectively.

[0085]

[0086] Then, an approximate Jacobian matrix is ​​constructed to ensure the diagonal dominance property. Specifically, it can be constructed as follows:

[0087]

[0088] At this time, A + -A - By making it a scalar, the matrix inversion operation is avoided.

[0089] Step 4.7: Begin pseudo-unsteady implicit iteration at level 1, with iteration number k = 1 / L1 and time step size of [missing value]. In the implicit iteration, a single-step LU diagonal scan is used, and in the inner iteration of the first-level mesh, a single-step LU diagonal scan is also used. The matrix A to be calculated in the forward scan is obtained from equation (11). + The forward scan begins, and the intermediate solution ΔQ of the inner iteration is obtained from equation (12). *

[0090]

[0091] The updated solution ΔQ obtained during the forward scan process * It can be directly stored in the right-hand term vector This further reduces memory consumption.

[0092] After the forward scan is completed, the matrix A to be calculated in the backward scan is obtained by formula (11). - Then begin the backward scan to solve equation (13).

[0093]

[0094] The updated solution ΔQ for this inner iteration is obtained. n This is used to calculate the initial value for the next iteration of the first-level grid, and the iteration continues until the number of iterations reaches k = 1 / L1.

[0095] Then repeat step 4.7 to complete the pseudo-unsteady implicit iteration on all grid levels. After all grid levels have been iterated, update the flow field variables in the mainstream equation and the turbulence model equation.

[0096] Step 4.8: Solve the turbulence model equation (14) on all grids, with an iteration count of 1 and a time step of Δt.

[0097]

[0098] In implicit iteration, a single-step LU diagonal scan is used, and the matrix to be calculated in the forward scan of the turbulence model equation is obtained from formula (8). The forward scan begins, and the intermediate solution ΔX of the inner iteration is obtained from equation (15). *

[0099]

[0100] The updated solution ΔX obtained during the forward scan process * It can be directly stored in the right-hand item R. X This further reduces memory consumption.

[0101] After the forward scan is completed, the matrix to be calculated in the backward scan is obtained by formula (11). Start the backward scan, and obtain ΔX from equation (16). n

[0102]

[0103] After the backward scan is completed, the updated solution ΔX of the turbulence model equations is obtained. n Update the turbulence model variables and the turbulent viscosity in the mainstream equations.

[0104] Step 5, determine the next step t n+1 Is the time greater than or equal to the time at which the calculation terminates (t)? end If the condition is met, proceed to step 6; otherwise, let n = n + 1 and return to step 3.

[0105] Step 6: Output the final result.

[0106] Example 1: Prediction of Aerodynamic Performance of Airfoil Separated Turbulent Flow

[0107] This embodiment evaluates the performance of the method in predicting unsteady turbulent flow using numerical simulation. The embodiment describes an airfoil under high angle-of-attack flow conditions where large-scale separation flow occurs on its upper surface, generating periodic vortex shedding. The computational domain is [-20, 20] × [-10, 10], using a structured / unstructured hybrid mesh. The mesh near the airfoil surface is shown in the figure. Figure 3 The incoming flow conditions are Mach number Ma = 0.2, static temperature T = 300K, angle of attack α = 10°, and initial conditions are given as free incoming flow values. Calculations are performed according to the steps described above.

[0108] Figure 4 This is a cloud map showing the instantaneous eddy viscosity distribution of the simulation results in this embodiment. Figure 4 As can be seen, the present invention clearly captures the flow separation and eddy shedding phenomena, indicating that the present invention can accurately characterize unsteady turbulent flow phenomena. Figure 5 The results of predicting lift and drag characteristics using the method of this invention and the traditional two-time-step method are compared. It can be seen that the prediction results of this invention and the traditional unsteady calculation method are basically consistent. However, the calculation efficiency of this invention is significantly improved. In the embodiment, the calculation efficiency of this invention is improved by about 80% compared with the traditional two-time-step method, indicating that this invention can achieve high-efficiency simulation of unsteady turbulent flow.

[0109] The parts of this invention not disclosed in detail are well-known technologies in the field.

[0110] Although the illustrative specific embodiments of the present invention have been described above to enable those skilled in the art to understand the invention, it should be understood that the invention is not limited to the scope of the specific embodiments. For those skilled in the art, various changes are obvious as long as they are within the spirit and scope of the invention as defined and determined by the appended claims, and all inventions utilizing the concept of the present invention are protected.

Claims

1. A high-efficiency calculation method for complex turbulent flows based on pseudo-unsteady time propagation, characterized in that: The method includes the following steps: Step 1: Establish the governing equations for unsteady turbulence simulation. The governing equations adopted are the Reynolds-averaged Navier-Stokes (RANS) equations for compressible flow. The equations are written in a decoupled form of the mainstream equations and the turbulence model equations. The mainstream equations are written as equation (1), and the turbulence model equations are written as equation (2). （1） （2） In formula (1) This indicates finding the partial derivative of the variable. As a conserved variable, , and for x , y , z Inviscid flux in three directions, and for x , y , z Viscous flux in three directions, subscript Represents viscosity. t For time terms, x , y , z For spatial terms, e For energy, For viscous stress, q For heat flow; in equation (2) Representing the variables in the turbulence model equations, It is the source term for generating the turbulence model equations. It is the source term that is destroyed in the turbulence model equation. D It is the diffusion term in the turbulence model equation; Step 2: Based on the turbulent flow to be simulated, determine the solution domain, divide the domain into grid cells, and assign initial flow field values, including density, to all grid cells. ,speed ,pressure and turbulence model equation variables X Initial value; Step 3, utilize the current t n Given the grid cell variables at each time step, calculate the inviscid and viscous fluxes for all grid cells, then sum the inviscid and viscous fluxes to obtain the current flux. t n The time derivative of the flow field variables in the grid cells at time t; Step 4, utilize the current Next step: Solving the time derivative of the flow field variables at time t. The flow field variable values ​​at time t, and update the turbulent viscosity information in the mainstream equation, and synchronously update the variable information of the mainstream equation in the turbulence model equation; Step 5, determine the next step. Is the time greater than or equal to the time at which the calculation terminates? If the condition is met, proceed to step 6; otherwise, let... Return to step 2; Step 6: Output the final result; Step 4 utilizes the current Next step: Solving the time derivative of the flow field variables at time t. The flow field variable values ​​at time t are calculated using an unsteady implicit time-propagation method, expressed as: （3） In the formula R This represents the sum of the inviscid and viscous fluxes calculated in step 3, indicated by the superscript. m The superscript represents the iteration step in an unsteady internal iterative process. n and n -1 indicates the current Time and the previous step Calculation results of flow field variables at time points. express Time and The size of the time step at time t, i.e. Its specific size is determined by a given... CFL Calculate the numbers; The unsteady implicit time-progression calculation method employs pseudo-unsteady time-progression, taking the number of iterations in the unsteady internal iteration process as 1, and the time step size as... Set to depend on the time step of the grid layering results ,in i Represents grid hierarchy, k Representing the i Number of iteration steps within a layer grid; writing the solution equations for the pseudo-unsteady time-progression method: （4） In the formula The dominant equation is in the th... i Flux residuals during arbitrary iterations of a layered mesh The Jacobian matrix representing the main equation is decomposed into non-negative eigenvalue matrices based on the matrix's eigenvalues. and eigenvalue non-positive matrices The form of summation: （5） Meanwhile, assuming that the time step is sufficiently small, the increments of the solutions at two adjacent time steps are approximately equal, i.e. Therefore, equation (5) can be further written as: （6） Therefore, at any grid level i The equations to be solved are reduced to pseudo-unsteady implicit discrete equations.

2. The high-efficiency calculation method for complex turbulent flow based on pseudo-unsteady time propagation as described in claim 1, characterized in that: The pseudo-unsteady time-progression method employs grid-based hierarchical computation at each level. i Using time steps on the grid Mesh layering should ensure that meshes with similar aspect ratios are in the same layer. To determine the mesh layering, a suitable direction is selected from the directions of the structural mesh, the mesh slenderness ratio is calculated along this suitable direction, and meshes with similar slenderness ratios are grouped into the same level along this suitable direction.

3. The high-efficiency calculation method for complex turbulent flow based on pseudo-unsteady time propagation as described in claim 2, characterized in that: The calculation of the mesh slenderness ratio first involves calculating the slenderness ratio of the first layer of mesh closest to the aircraft wall and the outermost layer of mesh furthest from the wall; the difference between the two is a multiple of [missing information]. Therefore, the growth rate for each level is determined as follows: The number of grid layers is I ; Then, at the level i All grids on a given feature scale L i ,in i =1,…, I Typical feature scales are: （7） Time and Time step of a moment By giving CFL For unsteady calculations, the time step of all grid cells needs to be unified to the smallest time step in the flow field grid cells. Therefore: （8） In the formula j Traverse all grid cells, Represents the volume of the mesh cell. Represents the largest eigenvalue of the grid cell. Represents the area vector of the grid cells. and These are characteristic coefficients related to viscosity; At the level i Using time steps on the grid and the corresponding number of inner iterations k By giving feature scale L i and minimum time step get: （9）。 4. The high-efficiency calculation method for complex turbulent flow based on pseudo-unsteady time propagation according to claim 1, characterized in that: The pseudo-unsteady time-progression method solves the implicit discrete equation (6) at each level of the grid. According to the upwind principle, the terms of equation (6) are obtained by backward and forward difference respectively: （10） To ensure convergence performance, an approximate Jacobian matrix is ​​constructed to guarantee the diagonal dominance property: specifically, it is: (11) therefore To avoid matrix inversion operations, the solution is converted to scalars; iterations are performed according to the grid hierarchy during the solution process, starting with the first-level grid, and the number of iterations is specified. Then iterate through the second level, the number of iterations... Until the final iteration I Level, iteration count is 1; After all levels of mesh iteration are completed, the flow field variables in the mainstream equation and the turbulence model equation are updated.

5. The high-efficiency calculation method for complex turbulent flow based on pseudo-unsteady time propagation according to claim 1, characterized in that: The pseudo-unsteady time-progression method solves the implicit discrete equations of the turbulence model. Unlike the mainstream equations, the turbulence model equations are solved directly over all grids, without using a hierarchical grid. （12） In equation (12) the subscript X Representing the variables of the turbulence model equations; the update of the variables of the turbulence model equations and the first mainstream equation. I The updates at each level are kept consistent, thus avoiding the need for solving the turbulence model equations to consume more computational resources. The updated turbulence model variables are used to update the turbulence viscosity in the mainstream equations.

6. The high-efficiency calculation method for complex turbulent flow based on pseudo-unsteady time propagation according to claim 4, characterized in that: The pseudo-unsteady time-progression method uses single-step LU diagonal scanning in implicit iteration. For the mainstream equation, single-step LU diagonal scanning is used in the inner iteration of each grid level; for the turbulence model equation, single-step LU diagonal scanning is used in the iteration of all grids. The mainstream equation is in the... i Iteration on the hierarchical grid, using a single-step forward scan, yields intermediate solutions for the inner iterations. : （13） The matrix that needs to be calculated in the forward scan The updated solution obtained during the scanning process is derived from formula (11). Stored directly in the right-hand term vector middle; After the forward scan is completed, the backward scan begins. : （14） The matrix that needs to be calculated in the backward scan From formula (11), after the backward scan is completed, the updated solution for this iteration is obtained. This is used for the next iterative calculation.

Images