CN115329689A - High-efficiency calculation method for complex turbulent flow based on pseudo-unsteady time propulsion - Google Patents
High-efficiency calculation method for complex turbulent flow based on pseudo-unsteady time propulsion Download PDFInfo
- Publication number
- CN115329689A CN115329689A CN202210846089.8A CN202210846089A CN115329689A CN 115329689 A CN115329689 A CN 115329689A CN 202210846089 A CN202210846089 A CN 202210846089A CN 115329689 A CN115329689 A CN 115329689A
- Authority
- CN
- China
- Prior art keywords
- time
- equation
- grid
- unsteady
- calculation
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/10—Geometric CAD
- G06F30/15—Vehicle, aircraft or watercraft design
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/22—Design optimisation, verification or simulation using Petri net models
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/08—Fluids
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Geometry (AREA)
- Mathematical Optimization (AREA)
- General Engineering & Computer Science (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- Evolutionary Computation (AREA)
- Computer Hardware Design (AREA)
- Data Mining & Analysis (AREA)
- Algebra (AREA)
- Automation & Control Theory (AREA)
- Operations Research (AREA)
- Aviation & Aerospace Engineering (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Computing Systems (AREA)
- Fluid Mechanics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Feedback Control In General (AREA)
Abstract
The invention discloses a high-efficiency calculation method for complex turbulent flow based on pseudo-unsteady time propulsion, which comprises the following steps: step 1, establishing a control equation of unsteady turbulence simulation; step 2, determining a solution area, dividing grid units in the area, and assigning initial values of a flow field to all the grid units; step 3, utilizing the current t n Time grid sheetComputing the non-viscous flux and the viscous flux of all grid units by using the metavariable, and summing the non-viscous flux and the viscous flux as the current t n A variable time derivative value of a flow field of a grid unit at a moment; step 4, utilizing the current t n Solving the time derivative value of the flow field variable at the moment n+1 The value of the flow field variable at the moment; step 5, judging the next step t n+1 Whether the time is greater than or equal to the calculation termination t end At the moment, if the condition is met, turning to the step 6; otherwise, enabling n = n +1, and returning to the step 2; and 6, outputting the final result. The method can effectively ensure the diagonal dominance property of the matrix and ensure the stability of numerical calculation, thereby improving the calculation efficiency and reducing the calculation period and the calculation cost.
Description
[ technical field ] A method for producing a semiconductor device
The invention belongs to the field of computational mathematics and computational fluid mechanics, and particularly relates to a high-efficiency computation method for complex turbulent flow based on pseudo-unsteady time propulsion.
[ background of the invention ]
As a complex flow phenomenon in nature and engineering, turbulent flow often exists in aircraft bypass flow in various complex flow forms, such as bottom flow, jet flow interference, separation flow, shock wave/boundary layer interference and the like, and has important influence on flow characteristics of bottom resistance, flow field vortex system structure, air inlet passage flow quality, engine fuel mixing and combustion and the like of the aircraft. In the process of subsonic to hypersonic flight, particularly in a large attack angle state, the side edge of an aircraft body, local body protrusions, a control surface and the like can induce large-scale flow separation, and a complex separation vortex structure is formed. Separation vortices have strong unsteady characteristics and tend to have a significant impact on aircraft performance. Therefore, accurate realization of fine simulations of unsteady complex turbulent flows is critical to the aerodynamic design, control systems, and propulsion system design of aircraft.
However, it is very difficult to achieve high accuracy and low cost of complex turbulent flow prediction. The current Unsteady turbulence high-precision Simulation method comprises an Unsteady Reynolds average Navier-Stokes (Unsteady Reynolds-averaged Navier-Stokes, URANS) method, a Large vortex Simulation (LES) method, a hybrid RANS/LES method, a Direct Numerical Simulation (DNS) method and the like, wherein the calculation amount of the LES and the DNS method is too Large, and the current method does not have the actual application capability of engineering complex appearance. On the other hand, practical experience shows that the calculation efficiency and convergence efficiency of the extremely complex turbulence are problems to be solved. Researches show that in order to accurately calculate unsteady complex turbulent flow with high Reynolds (Reynolds, re) number, a grid with a large slenderness ratio needs to be adopted near a wall surface in the conventional numerical calculation, and the maximum time step allowed by calculation is limited by the extremely low wall surface normal grid size, so that the time efficiency requirement of engineering design is difficult to meet.
In recent years, much work has been done domestically and abroad in improving the calculation efficiency of numerical methods. Current time integration schemes can be classified into explicit and implicit according to a boosting method, and can be classified into spatio-temporal decoupling dispersion (also called semi-dispersion) and spatio-temporal coupling dispersion (also called full dispersion) according to a spatio-temporal dispersion method. The unsteady turbulence high-precision simulation mostly adopts space-time decoupling dispersion, and has the advantages that only the space dispersion can be considered when the format is constructed, and a mature ordinary differential equation solution, such as an explicit multistage Runge-Kutta (RK) method, is directly applied in time. The traditional explicit format has a small time step, and besides developing an explicit large time step format, researchers have conducted beneficial exploration on implicit formats suitable for unsteady computation in recent years. However, in general, unsteady calculation has many difficulties in theoretical research, numerical calculation research and the like, and becomes one of the most challenging leading-edge research hotspots of fluid mechanics. Compared with the steady flow time calculation method with long development history and maturity, the unsteady flow time calculation method is still immature in development and has many problems to be further researched, wherein the difficulties and research focuses mainly on the problems of time discrete format with second-order precision, unsteady mass calculation and verification and confirmation of unsteady flow calculation results.
The method mainly aims at the problems of large unsteady calculation amount and low efficiency of the traditional URANS and mixed RANS/LES methods, the calculation grids are layered according to a pseudo-unsteady time advancing idea, different sub-iteration time step lengths are adopted by different layers of grids, and the calculation efficiency of turbulent flow simulation is improved by a strategy of reducing time dispersion precision and time step length. Compared with the existing numerical simulation method, the method used by the invention is constructed on the basis of a double-time-step discrete method according to the ideas of implicit time advance and layered calculation, inherits the advantages of good stability of implicit time format and no limitation of time step, and can effectively improve the unsteady calculation efficiency when facing to the grid with large wall surface slenderness ratio used in turbulence simulation. In addition, the diagonal dominance property of the matrix is effectively ensured due to the reduction of the time step, so that the stability of numerical simulation is further ensured. Therefore, on the premise of ensuring that the calculation precision is not obviously influenced, the method effectively improves the calculation efficiency of the complex unsteady turbulence and reduces the calculation cost.
[ summary of the invention ]
The technical problem to be solved by the invention is as follows: the method overcomes the defects of the prior art and provides a time advancing method suitable for unsteady simulation of the complex turbulent flow. Aiming at URANS and mixed RANS/LES methods, a large-length-to-thin-ratio grid is used near the wall surface of an aircraft, and based on a pseudo-time propulsion idea, the grid layering and time step reduction are utilized, so that the diagonal dominance property of a matrix is effectively ensured, the numerical calculation stability is ensured, the calculation efficiency is improved, and the calculation period and the calculation cost are reduced.
The technical scheme adopted by the invention for solving the technical problems is as follows: a high-efficiency calculation method for complex turbulent flow based on pseudo-unsteady time propulsion comprises the following steps:
in the formula (1)
The partial derivative of the variable is calculated, Q is a conservative variable, F, G and H are non-viscous fluxes in x, y and z directions, F μ G μ And H μ The viscous flux in the three directions x, y and z is shown, the lower subscript mu represents viscosity, t is time term, x, y and z are space terms, e is energy, tau is viscous stress, and q is heat flow. In the formula (2), X represents a variable of the turbulence model equation, P is a generation (production) source term of the turbulence model equation, epsilon is a destruction (destruction) source term of the turbulence model equation, and D is a diffusion (diffusion) term of the turbulence model equation.
And 2, determining a solving area according to the turbulence flow to be simulated, dividing grid units in the area, and assigning initial values of a flow field, including density rho, speed u, v, w, pressure p and initial values of a turbulence model equation variable X, to all the grid units.
Step 3, utilizing the current t n The variation of the grid unit at the moment, calculating the non-adhesive flux and the adhesive flux of all the grid units, and then summing the non-adhesive flux and the adhesive flux to obtain the current t n A variable time derivative value of a flow field of a grid unit at a moment;
step 4, utilizing the current t n Solving the time derivative value of the flow field variable at the moment n+1 The value of the flow field variable at the moment is changed, turbulent flow viscosity information in a main flow equation is updated, and variable information of the main flow equation in a turbulent flow model equation is synchronously updated;
step 5, judging the next step t n+1 Whether the time is greater than or equal to the calculation termination t end At the moment, if the condition is met, turning to the step 6; if not, enabling n = n +1, and returning to the step 2;
and 6, outputting the final result.
Said step 4 utilizes the current t n Solving the time derivative value of the flow field variable at the moment n+1 The variable value of the flow field at the moment adopts an unsteady implicit time advance calculation method, and the expression is
In the formula, R represents the summation result of the inviscid flux and the viscous flux obtained by calculation in the step 3, the superscript m represents the iteration step of the non-constant internal iteration process, and the superscripts n and n-1 represent the current t n Time and previous step t n-1 The calculation result of the flow field variable at the moment, delta t represents t n+1 Time and t n Time step size of time instant, i.e. Δ t = t n+1 -t n The specific size is calculated by a given number of CFLs.
Furthermore, the unsteady implicit time advance calculation method adopts pseudo unsteady time advance solution, the iteration step number of the unsteady internal iteration process is taken as 1, and the time step delta t is set as the time step depending on the grid layering result
Wherein i represents the mesh level and k represents the number of iteration steps in the ith mesh. Solving equation writing of pseudo-unsteady time-marching method
In the formula
Representing the flux residual of the mainstream equation in the process of any iteration of the ith grid,
the Jacobian matrix representing the mainstream equation can be decomposed into a nonnegative matrix A of eigenvalues according to the eigenvalues of the matrix + And eigenvalue non-positive matrix A - Form of summation
Meanwhile, on the premise that the time step is small enough, two steps are usedThe increments of the solutions for 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
Thus, at any grid level i, the solved equations are reduced to pseudo-unsteady implicit discrete equations.
Furthermore, the pseudo-unsteady time advancing method adopts grid hierarchical computation and uses time step length on the grid of the level i
Grid layering should ensure that grids with similar slenderness ratios are at the same level. To determine the mesh hierarchy, one may artificially select an appropriate direction from three directions of the structural mesh (in the case of three-dimensional computation), compute the mesh slenderness ratio along that direction, and classify meshes with similar slenderness ratios into the same level along that direction. In the invention, the slenderness ratio of the first layer of grid close to the wall surface of the aircraft and the outermost layer of grid far away from the wall surface is firstly calculated, the slenderness ratio of the first layer of grid close to the wall surface and the outermost layer of grid far away from the wall surface are different by a multiple of n, so that the growth ratio of each layer can be roughly determined as
Wherein the number of grid layers is I. Considering the complexity of the actual grid structure, the growth ratio can be adjusted appropriately, and the situation that the number of grids in a certain level is too large or too small is avoided. Then, all meshes on level I (I =1, \8230;, I) are given a feature scale L i A typical characteristic scale may be
L i ＝{L 1 ＝0.01,..L i ＝0.1...,L I ＝1} (7)
t n+1 Time and t n The time step delta t of the moment is calculated by giving the CFL number, and for the unsteady calculation, the time steps of all grid units need to be unified into the minimum time step in the flow field grid unit, so that the time step delta t of the moment is calculated by giving the CFL number, and the time step delta t of the moment is calculated by giving the CFL number to the unsteady calculation
Where j traverses all grid cells, V j Represents the volume of a grid cell and,
represents the maximum eigenvalue of the grid cell,
representing the grid cell area vector, c diff V and v diff Is a characteristic coefficient related to viscosity. Using time steps on a grid of level i
And corresponding number of inner iteration steps k may pass through a given feature scale L i And the minimum time step Δ t
Furthermore, the pseudo-unsteady time propulsion method solves an implicit discrete equation (6) on each level of grid, and according to the windward principle, backward and forward differences are respectively carried out on terms of the equation (6) to obtain
To ensure convergence performance, constructing an approximate Jacobian matrix ensures diagonal dominance properties. Specifically, it is possible to construct
Thus A is + -A - Becomes a scalar quantity, and avoids matrix fitting operation. Overlapping according to the hierarchical order of grids in the solving processFirst, the 1 st grid is iterated for 1/L 1 Then iterating the second stage for 1/L 2 And until the I level of the last iteration, the iteration number is 1. And after the iteration of all the grids is finished, updating flow field variables in a main flow equation and a turbulence model equation.
Furthermore, the pseudo-unsteady time propulsion method solves the implicit discrete equation of the turbulence model, which is different from the mainstream equation, and the turbulence model equation does not adopt hierarchical grid solution, namely, all grids are directly solved
The subscript X in equation (12) represents the turbulence model equation variables. It can be seen that the update of the turbulence model equation variable is consistent with the update of the level I of the main flow equation, so that the solution of the turbulence model equation is prevented from occupying more computing resources, and the updated turbulence model variable is used for updating the turbulence viscosity in the main flow equation.
Further, the pseudo-unsteady time-marching method uses single-step LU diagonal scanning in implicit iterations, single-step LU diagonal scanning in inner iterations of each level mesh for mainstream equations, and single-step LU diagonal scanning in iterations of all meshes for turbulence model equations. Without loss of generality, taking the iteration of the mainstream equation on the ith level grid as an example, a single-step forward scanning is used to obtain an internal iteration intermediate solution delta Q *
Matrix A to be calculated in the forward scan + The updated solution Δ Q obtained during the scanning is obtained from equation (11) * Can be directly stored in the right term vector
In addition, the memory overhead is further saved.
End of forward scanThen, backward scanning is started to obtain Δ Q n
Matrix A to be calculated in backward scanning - The updated solution delta Q of the iteration is obtained after the backward scanning is finished according to the formula (11) n And the calculation is used for the next iterative calculation.
Compared with the prior art, the invention has the advantages that:
1. the high-efficiency calculation method for the complex turbulent flow based on the pseudo-unsteady time propulsion can use a layered iteration method according to the calculated mesh slenderness ratio, and effectively avoids the reduction of the whole unsteady calculation efficiency caused by the over-small time step of the mesh region with the large slenderness ratio.
2. The pseudo-unsteady implicit time advancing method adopted by the invention inherits the characteristics of good implicit discrete stability and no limitation of time step length.
3. Aiming at a mainstream equation and a turbulence model equation, the method adopts the implicit time dispersion in the same form and has a unified calculation method, compared with the traditional unsteady implicit propulsion method, the method avoids the influence that the sub-iteration can not be effectively converged, and is more suitable for RANS equation-based calculation methods such as URANS, mixed RANS/LES and the like.
4. The invention uses the method of reducing the time discrete precision and the smaller time step length, so that the discrete matrix has better diagonal dominance characteristic, and the calculation stability can be further enhanced.
5. The single-step LU diagonal scanning used in the implicit iteration can avoid storage block diagonal matrixes and solution vectors at the previous step, saves memory overhead and is easy to realize programming.
[ description of the drawings ]
FIG. 1 is a flowchart of the overall process of the present invention.
Fig. 2a and b are partial flow charts of the complex turbulent flow high efficiency calculation method based on pseudo-unsteady time advance proposed by the method of the invention.
FIG. 3 illustrates a grid pattern near the airfoil surface in an embodiment of the present invention.
FIG. 4 is a cloud of transient vortex viscosity distributions of simulation results of an embodiment of the present invention.
FIG. 5 shows the comparison of the predicted lift resistance characteristics of the method of the present invention and the conventional two-time-step method.
[ detailed description ] A
The invention is described in further detail below with reference to the use flow chart of the method:
as shown in FIG. 1, the invention relates to a high-efficiency calculation method for complex turbulent flow based on pseudo-unsteady time propulsion, which comprises the following steps:
In the formula (1)
The partial derivative of the solved variable is shown, Q is a conservative variable, F, G and H are inviscid flux in the x, y and z directions, F μ G μ And H μ The viscous flux in the three directions x, y and z is shown, the lower subscript mu represents viscosity, t is time term, x, y and z are space terms, e is energy, tau is viscous stress and q is heat flow. In the formula (2), X represents a variable of a turbulence model equation, P is a generation (production) source term of the turbulence model equation, epsilon is a destruction (destruction) source term of the turbulence model equation, and D is a diffusion (diffusion) term of the turbulence model equation.
Step 2, determining a solving area according to the turbulent flow to be simulated, dividing grid units in the area, and assigning initial values of a flow field to all the grid units, wherein the initial values comprise density rho, speed u, speed v, pressure p and initial value of a turbulent model equation variable X;
step 3, utilizing the current t n The variation of the grid unit at the moment, calculating the non-adhesive flux and the adhesive flux of all the grid units, and then summing the non-adhesive flux and the adhesive flux to obtain the current t n A variable time derivative value of a flow field of a grid unit at a moment;
as shown in fig. 2a and b, the following steps are embodied in the method for calculating the high efficiency of the complex turbulent flow based on the pseudo-unsteady time advance:
step 4.1, use the current t n Time derivative value solving t of flow field variable of moment n+1 The variable value of the flow field at the moment adopts a calculation method based on unsteady implicit time advance, and the expression is
In the formula, R represents the summation result of the inviscid flux and the viscous flux obtained by calculation in the 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 calculation result of the flow field variable at the moment, delta t represents t n+1 Time t and n time step size of time instant, i.e. Δ t = t n+1 -t n The specific size is calculated by giving the number of CFLs;
and 4.2, carrying out grid layering, manually selecting a proper direction from three directions (under the three-dimensional calculation condition) of the structural grid, calculating the slenderness ratio of the grid along the direction, and classifying the grids with similar slenderness ratios into the same level along the direction. Noting the number of mesh levels I, all the meshes at level I (I =1, \8230;, I) are given a characteristic dimension L i A typical characteristic scale may be
L i ＝{L 1 ＝0.01,..L i ＝0.1...,L I ＝1} (4)
Step 4.3, calculate t n+1 Time t and n the time step at time instant Δ t. The calculation is carried out by giving the CFL number, and for the unsteady calculation, the time step length of all grid units needs to be unified into the minimum time step length in the flow field grid unit, so the time step length of all grid units needs to be unified into the minimum time step length in the flow field grid unit
Where j traverses all grid cells, V j Represents the volume of a grid cell and,
represents the maximum eigenvalue of the grid cell,
representing the grid cell area vector, c diff V and v diff Is a characteristic coefficient related to viscosity.
Step 4.4, calculate the time step Δ t on the grid of level i i k And corresponding inner iteration step number k, wherein i represents the grid level, and k represents the number of iteration steps in the ith layer of grid. By giving a characteristic dimension L i And the minimum time step Δ t
Step 4.5, adopting pseudo-unsteady time to advance and solve equation (3), taking the iteration step number of the iteration process in the unsteady time as 1, and setting the time step delta t as the time step depending on the grid layering result
Solving equation writing of pseudo-unsteady time-marching method
In the formula
Representing the flux residual of the mainstream equation during any iteration of the ith layer of the grid,
the Jacobian matrix representing the mainstream equation can be decomposed into a nonnegative matrix A of eigenvalues according to the eigenvalues of the matrix + And eigenvalue non-positive matrix A - Form of summation
Meanwhile, it is assumed that the increment of the solution at two adjacent time steps is approximately equal on the premise that the time step is sufficiently small, i.e., Δ Q = Q n+1 -Q n ≈Q n -Q n-1 Therefore, equation (8) can be further written as
Thus, at any grid level i, the solved equations are reduced to pseudo-unsteady implicit discrete equations.
Step 4.6, solving an implicit discrete equation (9) on each level grid, and respectively carrying out backward and forward difference on each item of the equation (9) according to the windward principle to obtain
Constructing an approximate Jacobian matrix then ensures the diagonal dominance property. Specifically, it is possible to construct
At this time A + -A - Becomes a scalar quantity, and avoids matrix inversion operation.
Step 4.7, starting pseudo-unsteady implicit iteration on the grid of the level 1, wherein the iteration times k =1/L 1 The step of time is
Single step LU diagonal scanning is used in the implicit iteration and single step LU diagonal scanning is used in the inner iteration of the level 1 grid. The matrix A to be calculated in the forward scanning is obtained from the formula (11) + Starting the forward scan, an inner iterative intermediate solution Δ Q is obtained from equation (12) *
Updated solution Δ Q obtained during forward scan * Can be directly stored in the right term vector
In addition, the memory overhead is further saved.
After the forward scanning is finished, the matrix A required to be calculated in the backward scanning is obtained by the formula (11) - Starting the backward scan, solving equation (13)
Obtaining the updated solution delta Q of the current internal iteration n And performing iteration calculation on the initial value in the next step of the 1 st level grid, and continuing the iteration until the iteration number reaches k =1/L 1 。
And then repeating the step 4.7, finishing the pseudo-unsteady implicit iteration on all the hierarchical grids, and updating flow field variables in the mainstream equation and the turbulence model equation after finishing the iteration of all the hierarchical grids.
And 4.8, solving a turbulence model equation (14) on all grids, wherein the iteration number is 1, and the time step is delta t.
In the implicit iteration, single-step LU diagonal scanning is used, and the matrix required to be calculated in the forward scanning of the turbulence model equation is obtained by the formula (8)
Starting the forward scan, an inner iterative intermediate solution Δ X is obtained from equation (15) *
Updated solution Δ X obtained during forward scan * Can be directly stored in the right end item R X In addition, the memory overhead is further saved.
After the forward scanning is finished, the matrix required to be calculated in the backward scanning is obtained by the formula (11)
Starting the backward scan, Δ X is obtained from equation (16) n
After the backward scanning is finished, obtaining the equation updating solution delta X of the turbulence model n The turbulence model variables and the turbulence viscosity in the mainstream equations are updated.
Step 5, judging the next step t n+1 Whether the time is greater than or equal to the calculation termination t end At the moment, if yes, go to step 6; if not, enabling n = n +1, and returning to the step 3;
and 6, outputting a final result.
Example 1: airfoil separation turbulent flow aerodynamic performance prediction
The evaluation method of the embodiment develops unsteady turbulence numerical simulation prediction performance. The embodiments describe the airfoil having a large area of flow separation at the upper surface under high angle of attack flow conditions, resulting in periodic vortex shedding. The calculation area is [ -20,20] × [ -10,10], and the calculation is carried out by adopting a structured/unstructured mixed grid, and the grid near the airfoil surface is shown in figure 3. The incoming flow conditions are mach number Ma =0.2, static temperature T =300K, angle of attack α =10 °, initial conditions are given as the value from the incoming flow. The calculation was performed according to the above procedure.
Fig. 4 is a cloud chart of the transient vortex viscosity distribution of the simulation result of the present embodiment. As can be seen in FIG. 4, the present invention clearly captures the flow separation and vortex shedding phenomena, demonstrating that the present invention is able to accurately depict unsteady turbulent flow phenomena. Fig. 5 is a comparison of the results of the prediction of lift resistance characteristics by the method of the present invention and the conventional dual time-step method, and it can be seen that the prediction results of the method of the present invention and the conventional unsteady flow calculation method are basically consistent, but the calculation efficiency of the present invention is significantly improved, and in the embodiment, the calculation efficiency of the present invention is improved by about 80% compared with the calculation efficiency of the conventional dual time-step method, which indicates that the present invention can realize high-efficiency simulation of unsteady turbulent flow.
Portions of the invention not disclosed in detail are well within the skill of the art.
Although illustrative embodiments of the present invention have been described above to facilitate the understanding of the present invention by those skilled in the art, it should be understood that the present invention is not limited to the scope of the embodiments, and various changes may be made apparent to those skilled in the art as long as they are within the spirit and scope of the present invention as defined and defined by the appended claims, and all matters of the invention which utilize the inventive concepts are protected.
Claims (8)
1. A high-efficiency calculation method for complex turbulent flow based on pseudo-unsteady time propulsion is characterized by comprising the following steps: the method comprises the following steps:
step 1, establishing a control equation of unsteady turbulence simulation, wherein the control equation adopted by the method is a compressible flow Reynolds average Navier-Stokes (RANS) equation, the equation is written into a main flow equation and turbulence model equation decoupling form, the main flow equation is written into a formula (1), and the turbulence model equation is written into a formula (2):
in the formula (1)
The partial derivative of the solved variable is shown, Q is a conservative variable, F, G and H are inviscid flux in the x, y and z directions, F μ G μ And H μ The flux of viscosity in the x, y and z directions is shown, the lower corner mark mu represents viscosity, t is a time term, x, y and z are space terms, e is energy, tau is viscosity stress, and q is heat flow; in the formula (2), X represents a variable of a turbulence model equation, P is a generation source term of the turbulence model equation, epsilon is a destruction source term of the turbulence model equation, and D is a diffusion term of the turbulence model equation;
step 2, determining a solving area according to the turbulent flow to be simulated, dividing grid units in the area, and assigning initial values of a flow field to all the grid units, wherein the initial values comprise density rho, speed u, speed v, pressure p and initial value of a turbulent model equation variable X;
step 3, utilizing the current t n The variation of the grid unit at the moment, calculating the non-adhesive flux and the adhesive flux of all the grid units, and then summing the non-adhesive flux and the adhesive flux to obtain the current t n A variable time derivative value of a flow field of a grid unit at a moment;
step 4, utilizing the current t n Time derivative value calculation of a current field variable at a timeSolution of the next step t n+1 The value of the flow field variable at the moment is changed, turbulent flow viscosity information in a main flow equation is updated, and variable information of the main flow equation in a turbulent flow model equation is synchronously updated;
step 5, judging the next step t n+1 Whether the time is greater than or equal to the calculation termination t end At the moment, if yes, go to step 6; if not, enabling n = n +1, and returning to the step 2;
and 6, outputting the final result.
2. The pseudo-unsteady time-marching-based high-efficiency computation method of complex turbulent flow according to claim 1, characterized in that:
said step 4 utilizing the current t n Solving the time derivative value of the flow field variable at the moment n+1 The variable value of the flow field at the moment adopts a calculation method based on unsteady implicit time advance, and the expression is
Wherein R represents the summation result of the inviscid flux and the viscous flux obtained by the calculation in the 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 last step t n-1 The calculation result of the flow field variable at the moment, delta t represents t n+1 Time t and n time step size of time instant, i.e. Δ t = t n+1 -t n The specific size is calculated by a given number of CFLs.
3. The pseudo-unsteady time-marching-based high-efficiency computation method of complex turbulent flow according to claim 2, characterized in that: the unsteady implicit time advance calculation method adopts pseudo unsteady time advance solution, the iteration step number of the unsteady internal iteration process is taken as 1, and the time step delta t is set as the time step depending on the grid layering result
Wherein i represents a grid level, and k represents the number of iteration steps in the ith grid; writing of solving equations for pseudo-unsteady time-marching methods
In the formula
Representing the flux residual of the mainstream equation during any iteration of the ith layer of the grid,
the Jacobian matrix representing the mainstream equation can be decomposed into a nonnegative matrix A of eigenvalues according to the eigenvalues of the matrix + And eigenvalue non-positive matrix A - Form of summation
At the same time, it is assumed that the increments of the solutions at two adjacent time steps are approximately equal, i.e. Δ Q = Q, on the premise that the time steps are sufficiently small n+1 -Q n ≈Q n -Q n-1 Therefore, equation (5) can be further written as
Thus, at any grid level i, the solved equations are reduced to pseudo-unsteady implicit discrete equations.
4. The pseudo-unsteady time-based complex turbulent flow high efficiency computation method of claim 3, characterized in that: the pseudo-unsteady time advancing method adopts grid hierarchical calculation and uses time step length on a grid of a hierarchy i
The grid layering ensures that grids with similar slenderness ratios are in the same level; to determine the mesh hierarchy, an appropriate direction is selected from the directions of the structural meshes, mesh slenderness ratios are calculated in the appropriate direction, and meshes having similar slenderness ratios are classified into the same level in the appropriate direction.
5. The pseudo-unsteady time-marching-based complex turbulent flow high efficiency computation method of claim 4, characterized in that: the mesh slenderness ratio is calculated, firstly, the slenderness ratio of a first layer of mesh close to the wall surface of the aircraft and an outermost layer of mesh far away from the wall surface are calculated, the difference between the slenderness ratio and the outermost layer of mesh is n, and therefore the growth ratio of each layer can be determined as
Wherein the number of grid layers is I; then, all meshes on level I (I =1, \8230;, I) are given a feature scale L i A typical characteristic scale may be
L i ＝{L 1 ＝0.01,..L i ＝0.1...,L I ＝1} (7)
t n+1 Time t and n the time step delta t of the moment is calculated by giving the CFL number, and for the unsteady calculation, the time steps of all grid units need to be unified into the minimum time step in the flow field grid unit, so that the time step delta t of the moment is calculated by giving the CFL number, and the time step delta t of the moment is calculated by giving the CFL number to the unsteady calculation
Where j traverses all grid cells, V j Representing the volume of a grid cell,
represents the maximum eigenvalue of the grid cell,
representing the grid cell area vector, c diff V and v diff Is a characteristic coefficient related to viscosity; using time steps on a grid of level i
And the corresponding number of inner iteration steps k through a given feature scale L i And the minimum time step Δ t
6. The pseudo-unsteady time-marching-based complex turbulent flow high efficiency computation method of claim 3, characterized in that: the pseudo-unsteady time propulsion method solves an implicit discrete equation (6) on each level of grid, and according to the windward principle, backward and forward differences are respectively carried out on each item of the equation (6) to obtain
To ensure convergence performance, an approximate Jacobian matrix is constructed to ensure diagonal dominance: specifically, the method comprises the following steps:
thus A is + -A - The scalar quantity is formed, and matrix fitting operation is avoided; iteration is carried out according to the hierarchical order of grids in the solving process, the 1 st level grid is iterated firstly, and the iteration times are 1/L 1 Then iterating the second stage for 1/L 2 Until the I level of the last iteration, the iteration number is 1; and after the iteration of all the grids is finished, updating flow field variables in a main flow equation and a turbulence model equation.
7. The pseudo-unsteady time-based complex turbulent flow high efficiency computation method of claim 3, characterized in that: the method for solving the implicit discrete equation of the turbulence model by the pseudo-unsteady time propulsion method is different from the mainstream equation, and the turbulence model equation does not adopt hierarchical grid solution, namely, all grids are directly solved
The subscript X in equation (12) represents the turbulence model equation variables; the update of the turbulence model equation variable is consistent with the update of the I level of the main flow equation, so that the solution of the turbulence model equation is prevented from occupying more computing resources, and the updated turbulence model variable is used for updating the turbulence viscosity in the main flow equation.
8. The pseudo-unsteady time-marching-based complex turbulent flow high efficiency computation method of claim 3, characterized in that: the pseudo-unsteady time advancing method uses single-step LU diagonal scanning in implicit iteration, uses single-step LU diagonal scanning in inner iteration of each level grid aiming at a mainstream equation, and uses single-step LU diagonal scanning in iteration of all grids aiming at a turbulence model equation; taking the iteration of the mainstream equation on the ith level grid as an example, the single-step forward scanning is used to obtain the intermediate solution delta Q of the internal iteration *
Matrix A to be calculated in the forward scan + The updated solution Δ Q obtained during the scanning is obtained from equation (11) * Directly stored in the right term vector
Performing the following steps;
end of forward scanThen, backward scanning is started to obtain Δ Q n
Matrix A to be calculated in backward scan - The updated solution delta Q of the iteration is obtained after the backward scanning is finished according to the formula (11) n And the calculation is used for the next iterative calculation.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210846089.8A CN115329689A (en) | 2022-07-05 | 2022-07-05 | High-efficiency calculation method for complex turbulent flow based on pseudo-unsteady time propulsion |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210846089.8A CN115329689A (en) | 2022-07-05 | 2022-07-05 | High-efficiency calculation method for complex turbulent flow based on pseudo-unsteady time propulsion |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN115329689A true CN115329689A (en) | 2022-11-11 |
Family
ID=83918080
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210846089.8A Pending CN115329689A (en) | 2022-07-05 | 2022-07-05 | High-efficiency calculation method for complex turbulent flow based on pseudo-unsteady time propulsion |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115329689A (en) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116227388A (en) * | 2023-04-21 | 2023-06-06 | 中国空气动力研究与发展中心计算空气动力研究所 | Dynamic adjustment method, system, equipment and medium for CFL number of high-ultra-flow simulation |
| CN116244853A (en) * | 2023-02-22 | 2023-06-09 | 北京航空航天大学 | Harmonic balance method accelerating convergence method for predicting unsteady flow of turbine |
| CN117393062A (en) * | 2023-12-13 | 2024-01-12 | 上海交通大学四川研究院 | Simulation method for rigid chemical reaction flow rollback self-adaptive semi-hidden semi-explicit coupling time |
-
2022
- 2022-07-05 CN CN202210846089.8A patent/CN115329689A/en active Pending
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116244853A (en) * | 2023-02-22 | 2023-06-09 | 北京航空航天大学 | Harmonic balance method accelerating convergence method for predicting unsteady flow of turbine |
| CN116244853B (en) * | 2023-02-22 | 2023-10-20 | 北京航空航天大学 | Harmonic balance method accelerating convergence method for predicting unsteady flow of turbine |
| CN116227388A (en) * | 2023-04-21 | 2023-06-06 | 中国空气动力研究与发展中心计算空气动力研究所 | Dynamic adjustment method, system, equipment and medium for CFL number of high-ultra-flow simulation |
| CN116227388B (en) * | 2023-04-21 | 2023-07-07 | 中国空气动力研究与发展中心计算空气动力研究所 | Dynamic adjustment method, system, equipment and medium for CFL number of high-ultra-flow simulation |
| CN117393062A (en) * | 2023-12-13 | 2024-01-12 | 上海交通大学四川研究院 | Simulation method for rigid chemical reaction flow rollback self-adaptive semi-hidden semi-explicit coupling time |
| CN117393062B (en) * | 2023-12-13 | 2024-02-23 | 上海交通大学四川研究院 | Simulation method for rigid chemical reaction flow rollback self-adaptive semi-hidden semi-explicit coupling time |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN115329689A (en) | High-efficiency calculation method for complex turbulent flow based on pseudo-unsteady time propulsion | |
| Ito et al. | Unstructured Mesh Generation for Viscous Flow Computations. | |
| Pirzadeh | An adaptive unstructured grid method by grid subdivision, local remeshing, and grid movement | |
| CN112016167A (en) | Aircraft aerodynamic shape design method and system based on simulation and optimization coupling | |
| CN110309571B (en) | Wing body fusion underwater glider external shape optimization method based on radial basis function model | |
| CN111079228A (en) | Pneumatic shape optimization method based on flow field prediction | |
| CN105718634A (en) | Airfoil robust optimization design method based on non-probability interval analysis model | |
| Leifsson et al. | Variable-fidelity aerodynamic shape optimization | |
| CN114065662B (en) | Method suitable for rapidly predicting airfoil flow field with variable grid topology | |
| CN110276131B (en) | Wing body fusion underwater glider appearance optimization method based on polynomial response surface model | |
| CN107766620A (en) | A kind of Aerodynamic Heating structural optimization method based on reduced-order model | |
| CN110543677A (en) | vortex characteristic driven rotational turbulence PANS model | |
| CN116628854A (en) | Wing section aerodynamic characteristic prediction method, system, electronic equipment and storage medium | |
| Murayama et al. | Japan Aerospace Exploration Agency Studies for the Fifth AIAA Drag Prediction Workshop | |
| Hongtao et al. | Numerical simualtion research on the transonic aeroelasticity of a highaspect-ratio wing | |
| Barrett et al. | Airfoil shape design and optimization using multifidelity analysis and embedded inverse design | |
| Iuliano et al. | Design of a supersonic natural laminar flow wing-body | |
| Ekelschot et al. | Enabling metric-based mesh adaptation for advanced compressible flow simulations using US3D | |
| Chi et al. | Hybrid inverse/optimization design approach for transonic natural-laminar-flow airfoils | |
| Ghoreyshi et al. | Grid Quality and Resolution Effects on the Aerodynamic Modeling of Parachute Canopies | |
| CN117077297B (en) | Natural laminar flow nacelle pneumatic robust optimization design method based on double-layer proxy model | |
| Chen et al. | Efficient multidisciplinary aerodynamic optimization design based on discrete adjoint method | |
| Namura et al. | Multipoint design of vortex generators on a swept infinite-wing under cruise and critical condition | |
| CN117077298B (en) | Aircraft robust optimization design method based on gradient enhancement random Co-Kriging model | |
| CN116702331A (en) | Airship comprehensive optimization method based on neural network |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |