CN112632874A - Optimization method and system for numerical simulation of helicopter flow field - Google Patents
Optimization method and system for numerical simulation of helicopter flow field Download PDFInfo
- Publication number
- CN112632874A CN112632874A CN202011632168.6A CN202011632168A CN112632874A CN 112632874 A CN112632874 A CN 112632874A CN 202011632168 A CN202011632168 A CN 202011632168A CN 112632874 A CN112632874 A CN 112632874A
- Authority
- CN
- China
- Prior art keywords
- matrix
- flow field
- grid
- airflow
- data
- 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
- 238000000034 method Methods 0.000 title claims abstract description 57
- 238000004088 simulation Methods 0.000 title claims abstract description 26
- 238000005457 optimization Methods 0.000 title claims abstract description 19
- 239000011159 matrix material Substances 0.000 claims abstract description 126
- 230000008569 process Effects 0.000 claims abstract description 20
- 238000007781 pre-processing Methods 0.000 claims abstract description 9
- 238000000354 decomposition reaction Methods 0.000 claims description 23
- 238000004364 calculation method Methods 0.000 claims description 21
- 238000004422 calculation algorithm Methods 0.000 claims description 19
- 230000006835 compression Effects 0.000 claims description 10
- 238000007906 compression Methods 0.000 claims description 10
- 239000013598 vector Substances 0.000 claims description 9
- 239000012530 fluid Substances 0.000 claims description 5
- 238000012937 correction Methods 0.000 claims description 4
- 230000004907 flux Effects 0.000 claims description 4
- 238000006467 substitution reaction Methods 0.000 claims description 4
- 230000017105 transposition Effects 0.000 claims description 4
- 238000012545 processing Methods 0.000 claims description 3
- 238000004590 computer program Methods 0.000 claims description 2
- 230000002194 synthesizing effect Effects 0.000 claims description 2
- 238000010276 construction Methods 0.000 claims 1
- 238000013461 design Methods 0.000 abstract description 4
- 238000004458 analytical method Methods 0.000 abstract description 3
- 230000008859 change Effects 0.000 description 3
- 239000007789 gas Substances 0.000 description 3
- 238000004134 energy conservation Methods 0.000 description 2
- 230000008521 reorganization Effects 0.000 description 2
- FGUUSXIOTUKUDN-IBGZPJMESA-N C1(=CC=CC=C1)N1C2=C(NC([C@H](C1)NC=1OC(=NN=1)C1=CC=CC=C1)=O)C=CC=C2 Chemical compound C1(=CC=CC=C1)N1C2=C(NC([C@H](C1)NC=1OC(=NN=1)C1=CC=CC=C1)=O)C=CC=C2 FGUUSXIOTUKUDN-IBGZPJMESA-N 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000008030 elimination Effects 0.000 description 1
- 238000003379 elimination reaction Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
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
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- 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
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Computer Hardware Design (AREA)
- General Engineering & Computer Science (AREA)
- Geometry (AREA)
- Evolutionary Computation (AREA)
- Fluid Mechanics (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Physics (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Computing Systems (AREA)
- Algebra (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention discloses an optimization method and an optimization system for numerical simulation of a helicopter flow field, which can accelerate numerical simulation analysis of the airflow state of a helicopter in the flight process and provide data reference for aerodynamic design of the helicopter. The whole system comprises a data acquisition module, a data preprocessing module, a flow field solving module and a grid boundary data updating module. The invention mainly optimizes the performance of the flow field solving module and can improve the solving speed. The flow field solving module mainly calculates the amount of the algebraic equation. The randomization treatment of the SRBT and the application of the low-rank matrix are combined, so that the operation times are reduced to the maximum extent. The method greatly improves the solving speed of the flow field on the premise of ensuring certain accuracy, and accelerates the updating speed of the numerical simulation of the flow field.
Description
Technical Field
The invention belongs to the technical field of aerodynamic force of helicopters, and relates to an optimization method and system based on numerical simulation of a helicopter flow field.
Background
The aerodynamic force of the whole helicopter is an important means for checking the flight capability of the helicopter. There are three main approaches to the research of aerodynamic problems: theoretical analysis (based on flow simplification assumptions giving an analytical solution to the problem under study), experimental simulation (ground based experiments as a means of study), and numerical calculations. The numerical simulation method is characterized in that a basic law (mass conservation law, Newton's second law, energy conservation law, gas state equation and the like) and related constitutive equations and state equations form a partial differential equation set for describing the change rate of the air state (density, speed and pressure) along with time, then integral and differential terms in the equations are approximately expressed into discrete algebraic forms, so that control equations in the integral forms are converted into algebraic equation sets, and finally the algebraic equation sets are solved through a computer, so that the numerical solution of a flow field on time and space is obtained. The invention improves the performance of the whole numerical simulation mainly by optimizing the solution algebraic equation. The algebraic equation coefficient matrix may be represented by a sparse matrix, which we compute using the multi-wave front Method. In the multiple wavefront method, the calculation of the wavefront matrix is a dense matrix calculation problem, and since the calculation task is divided into the calculation of a plurality of wavefront matrices, the calculation process involving the wavefront matrices becomes a core step. The direct solver of a large linear equation set is widely used due to its good numerical robustness, reliability and high availability. However, it also has the disadvantages of high computational complexity and high memory usage. This also limits the scope of application of the direct solution method to large scale problems (matrices with hundreds of millions of unknowns). The technology aims to reduce the time consumption of algebraic equation solution on the premise of ensuring robustness, calculation performance and high availability, and can improve the performance of numerical simulation of the whole helicopter flow field.
Disclosure of Invention
An object of the present invention is to solve the above problems and provide an optimization system for numerical simulation of a helicopter flow field, which can accelerate numerical simulation analysis of airflow state during helicopter flight and provide data reference for aerodynamic design of a helicopter. The whole process is divided into the following four modules: the device comprises a data acquisition module, a data preprocessing module, a flow field solving module and a grid boundary data updating module. The invention mainly optimizes the performance of the flow field solving module and can improve the solving speed. The flow field solving module mainly calculates the amount of the algebraic equation. By combining the randomization process of the srbt (symmetric Random transmit transform) and the application of the low rank matrix, the operation times are reduced to the maximum extent. Under the data layout of tile (tile) transformation, the calculation task granularity is reduced, the parallelism degree is improved, a low-rank matrix is introduced, the existing ACA (adaptive Cross application) algorithm is improved, the defect that the ACA algorithm has too many calculation iteration times when some tile blocks with poor low-rank characteristics are compressed is overcome, the operation number is greatly reduced, the compression operation can be completed under the condition of a low-rank threshold value required by the application, and then the OpenMP task constraint is used for parallelization, so that the calculation speed is greatly improved on a platform with a shared memory. And finally, the accuracy of the final result is improved by using an iterative refinement algorithm, the solving speed of the flow field is greatly improved on the premise of ensuring certain accuracy, and the updating speed of the numerical simulation of the flow field is accelerated.
The system comprises a data acquisition module, a data preprocessing module, a flow field solving module and a grid boundary data updating module.
The data acquisition module is used for acquiring geometrical parameters of a helicopter body and physical parameter information of airflow in the flight process, wherein the physical parameter information of the airflow comprises incoming flow speed, attack angle and the like in the flight process of the helicopter;
the data preprocessing module is used for constructing an airflow mesh topological structure, wherein the airflow mesh topological structure is subjected to mesh division according to an overlapping mesh method, and each mesh contains physical parameter information of airflow at the current position;
the flow field solving module is used for carrying out numerical solution on the flow field in each grid of the airflow grid topological structure to obtain flow field data;
the grid boundary data updating module is used for updating the grid nodes of the airflow grid topological structure according to the flow field data obtained by the flow field solving module, so that the purpose of numerical simulation of the flow field is achieved.
Another object of the present invention is to provide a method for optimizing numerical simulation of a helicopter flow field.
The method comprises the following steps of (1) obtaining geometrical parameters of a helicopter body and physical airflow parameter information in a flight process, wherein the physical airflow parameter information comprises incoming flow speed, attack angle and the like in the flight process of a helicopter;
step (2), constructing an airflow mesh topological structure according to the data acquired in the step (1), wherein the airflow mesh topological structure is subjected to mesh division according to an overlapping mesh method, and each mesh contains physical parameter information of the airflow at the current position;
step (3), each grid in the airflow grid topological structure carries out independent flow field calculation to obtain grid boundary data of flow field spatial distribution so as to provide a numerical value required by flow field updating;
and (4) updating the topological structure of the airflow grid according to the grid boundary data of the flow field spatial distribution in the step (3).
And (5) repeating the steps (2), (3) and (4) until the flow field is converged.
It is a further object of the present invention to provide a computer-readable storage medium having stored thereon a computer program which, when executed in a computer, causes the computer to perform the above-mentioned method.
It is a further object of the present invention to provide a computing device comprising a memory having stored therein executable code and a processor that, when executing the executable code, implements the method described above.
The technical scheme provided by the invention has the following beneficial effects:
on the premise of ensuring the accuracy requirement, the method can greatly improve the solving speed of the flow field, improve the prediction efficiency of the flow field and provide data reference for the aerodynamic design of the helicopter.
Drawings
FIG. 1 is a block diagram of the system of the present invention;
FIG. 2 is a flow chart of matrix decomposition by the LDLT tile algorithm;
FIGS. 3(a), (b) are a column-first data layout and a tile block data layout, respectively;
FIG. 4 is a data dependency of LDLT decomposition on the realization flow and parallelism of tasks in a solver;
FIGS. 5(a) - (b) are full rank matrix data format and low rank matrix data format before and after compression, respectively;
FIG. 6 is a flow chart of an improved ACA algorithm;
FIG. 7 is a flow chart of the mesh boundary data update module.
Detailed Description
Embodiments of the invention are described in further detail below with reference to the accompanying drawings:
the flow field spatial distribution obtained by the optimization method based on the numerical simulation of the helicopter flow field provides data support for the optimization design of the aerodynamic appearance of the aircraft;
as shown in fig. 1, the optimization system based on the numerical simulation of the helicopter flow field in the present invention includes a data acquisition module, a data preprocessing module, a flow field solving module, and a grid boundary data updating module.
The data acquisition module is used for acquiring geometrical parameters of a helicopter body and physical parameter information of airflow in the flight process, wherein the physical parameter information of the airflow comprises incoming flow speed, attack angle and the like in the flight process of the helicopter;
the data preprocessing module is used for constructing an airflow mesh topological structure, wherein the airflow mesh topological structure is subjected to mesh division according to an overlapping mesh method, and each mesh contains physical parameter information of airflow at the current position;
the flow field solving module is used for carrying out numerical solution on the flow field in each grid of the airflow grid topological structure to obtain flow field data;
the grid boundary data updating module is used for updating the grid nodes of the airflow grid topological structure according to the flow field data obtained by the flow field solving module, so that the purpose of numerical simulation of the flow field is achieved.
The optimization method based on the numerical simulation of the helicopter flow field specifically comprises the following steps:
the method comprises the following steps of (1) obtaining geometrical parameters of a helicopter body and physical airflow parameter information in a flight process, wherein the physical airflow parameter information comprises incoming flow speed, attack angle and the like in the flight process of a helicopter;
step (2), constructing an airflow grid topological structure according to the data acquired in the step (1); the airflow grid topological structure is subjected to grid division according to an overlapping grid method, and each grid contains physical parameter information of airflow at the current position;
the airflow grid topological structure is established by an overlapping grid method, and specifically comprises the following steps:
(1) dividing the flow field into a plurality of regions having overlapping portions;
(2) respectively generating a grid structure on each area;
(3) each area is independently solved;
(4) and carrying out flow field information exchange at the overlapped part through inter-grid interpolation.
Step (3), each grid in the airflow grid topological structure carries out independent flow field calculation to obtain grid boundary data of flow field spatial distribution so as to provide a numerical value required by flow field updating;
3-1 air flow follows mainly four laws: kinematic aspects, kinetic aspects, thermodynamic aspects and physical and chemical properties of gases, including the laws of mass conservation, newton's second law, law of energy conservation, entropy equations, equations of gas states, and the like. Based on the above basic law and related constitutive equations and state equations, a partial differential equation set describing the time rate of change of the air state (density, velocity, pressure) can be formed, i.e., a Navier-Stocks equation, which is referred to as an N-S equation for short. Establishing a three-dimensional depressible Reynolds average N-S equation on each grid in the airflow grid topological structure:
wherein W is a conservative variable, F (W) and FvConvective flux and viscous flux, respectively. The partial differential equation set on each grid and the final result of the N-S equation can be approximated to a linear algebraic equation set, and the linear equation set is as follows:
Ax＝b (1)
the matrix dimension of the coefficient matrix A is n, nb is the block size, and x is the equation unknown quantity;
the coefficient matrix A is composed of data blocks with the size of nb × nb, and total nt × nt data blocks;
extracting a coefficient matrix A of an algebraic equation;
3-2 randomizing the coefficient matrix A by SRBT method, and processing the matrix ArTransforming into a tile (tile) shape;
Ara random matrix after coefficient matrix a randomization:
Ar＝UAUT (4)
wherein T represents transposition, U represents random matrix, and A after randomizationrThe matrix size of (a) is still equal to the original matrix a;
the random matrix U is a recursive butterfly matrix:
U＝Ud×...×U1 (5)
wherein d represents deepDegree, usually less than or equal to 2, Uq(1. ltoreq. q. ltoreq. d) represents a block diagonal matrix generating the random matrix U:
wherein B ispRepresenting a butterfly matrix of size n/2q-1(1≤p≤2q-1) Remember n1＝n/2q-1Then the size is n1×n1The butterfly matrix of (a) is defined as:
wherein R ispAnd SpRespectively a random diagonal matrix and a nonsingular matrix with the size of n1/2×n1/2；RpThe value of the middle element is represented by eρ/10Is generated, ρ is
Random values within a range, thus the butterfly matrix BpThe element range of (B) is as follows:
e-1/20≤rp≤e1/20 (8)
tile-shaped means that a matrix is divided into grid units, and one tile block is used as a unit for sequential storage, so that the task granularity is refined, and the use efficiency of the cache is improved;
fig. 3(a) is a data layout of prior matrix column priority storage, which is only stored according to an actual structure of data, and there is no relevance of extracting associated data in an actual use process, and there is no good cache reuse, so that when data is accessed, the cache hit rate is low, and thus the access performance is reduced. Fig. 3(b) is a tile-shaped data layout, and since a single tile block is used as a unit in the matrix operation process, the tile-shaped data layout can improve the hit rate of the cache and make full use of the fast access characteristic of the cache. Will matrix ArFrom the layout format shown in FIG. 3(a) to that shown in FIG. 3(b)The layout format shown.
3-3 decomposing tile (tile) like matrix A by multi-wavefront methodrExtracting a wave front matrix;
the multi-wave front method reorganizes the operation related to sparse decomposition, the reorganization is based on the structure of a elimination tree, the reorganization converts the task of decomposing a large sparse matrix into a series of subtasks, and each subtask only comprises the local decomposition of a smaller dense wave front matrix Q;
3-4, performing LDLT decomposition on the wavefront matrix Q by adopting an improved LDLT tile algorithm, and parallelizing by utilizing OpenMP task constraint to obtain a lower triangular matrix L and a diagonal matrix D;
fig. 2 describes the whole decomposition process using the modified LDLT tile algorithm, and fig. 4 takes nt as 3 as an example, i.e. assuming that the matrix is divided into 3 × 3 blocks, the task relevance and execution order in the outermost one cycle process are in the whole decomposition process. In the figure, the factor step is the following step b), the solve step is the following step c), the compress step is the following step d), and the update step is the following step f). The arrows in the figure indicate the direct sequential dependency of the steps and also illustrate how the decomposition process can be performed in parallel.
The LDLT decomposition decomposes the matrix Q into:
Q＝LDLT (10)
where each element in the L, D matrix represents a tile of nb x nb. T refers to transposing the matrix.
The improved LDLT tile algorithm is specifically as follows:
a) initializing k to 1;
b) wave front matrix Q diagonal block Q by LDLT algorithmkkDecomposing to obtain a lower triangular matrix LkkAnd diagonal matrix Dkk：
Wherein k is more than or equal to 1 and less than or equal to nt, and T represents transposition;
c) obtaining the lower triangular matrix L according to the abovekkAnd diagonal matrix DkkCalculating a matrix L of all tile blocks in the kth column of the wavefront matrix Qik；
Wkk＝LkkDkk (15)
Wherein k is not less than i not more than nt.
d) Adopting an improved ACA algorithm to compress the tile blocks;
the format change of data before and after compression is described in fig. 5, the matrix after low-rank compression has the characteristics of low storage space occupancy rate, low calculation operand and the like, and for the step e) with the largest calculation amount, the compression processing is performed before the step e), so that the calculation performance of the step e) can be greatly improved.
The flow of the improved ACA algorithm implementation is depicted in fig. 6. Each cycle generates a group of vectors (a row vector and a column vector which are respectively merged into the matrix X and the matrix Y) according to the ACA algorithm traditional mode to approximately express the value of the whole matrix, and the cycle is ended and the compression step is ended until the approximation precision reaches a low-rank threshold or the cycle number reaches a set cycle threshold. The accuracy of the approximation can be calculated by the two-norm of the difference of the low rank matrix and the original matrix. When a matrix with poor low-rank property is encountered in the conventional ACA algorithm, the situations of excessive cycle times, overlong operation time and excessive generated vector groups can occur, and the subsequent calculation optimization is not obvious. We therefore set a threshold for the number of cycles to limit this, and when the number of cycles reaches the threshold, we will stop the compression and output the result, and then the operation involving the matrix block will take the original matrix value instead of the compressed matrix.
Therefore, the compression algorithm is as follows, the lower triangular matrix L of all tile blocks of the k column of step c) isikMatrix L is divided according to ACA algorithm traditional modeikDecomposing the matrix into X, Y low-rank matrix; if the precision reaches a low-rank threshold value, or the number of vector groups in the low-rank matrix reaches a cycle threshold value, the cycle is ended, otherwise, a new vector group is generated and added to X, Y, and a new X, Y low-rank matrix is generated.
Wherein k +1 is more than or equal to i and less than or equal to nt.
e) Update Qij
Judging whether the rank of the low-rank matrix X, Y after the compression in the step d) is more than or equal to a cycle threshold value, if so, using an original full-rank matrix LikUpdate Qij：
If not, substituting equation (17) into equation (18), updating Q with the low rank matrixij：
Wherein, i is more than or equal to k +1 and less than or equal to nt, j is more than or equal to k +1 and less than or equal to i
f) And judging whether k under the current iteration is nt, if so, finishing the decomposition of all tile blocks, otherwise, updating k to k +1, and returning to the step b).
3-5, obtaining the solution of the final equation set through forward substitution and reverse substitution;
obtaining the decomposition result of the wave front matrix Q from the previous step
Q＝LDLT (10)
And synthesizing the partial decomposition results of the series of wavefront matrixes Q to obtain the decomposition result of the original coefficient matrix A.
The decomposition result of the original coefficient matrix A is obtained by the subsequent calculation of the multi-wave front method
Wherein L ismergeAnd DmergeAre respectively the matrix ALDLTThe lower triangular matrix and the diagonal matrix of the decomposition result.
According to the formula (4), it can be seen that
Ax＝UTbU (21)
With reference to formulas (20) to (21), it can be seen that:
the intermediate variable y is found, then the solution x is found by substituting (22) into equation (1):
x＝Uy (23)
3-6, one iteration optimization is used to improve the accuracy of the solution.
3-6-1 calculating the residual r according to the solution x given by equation (23):
r＝b-Ax (24)
3-6-2 calculating the residual r from the solution x obtained from equation (24):
Ad＝r (25)
where d is a correction value.
3-6-3 update correction value
x′＝x+d (26)
Where m represents the number of iterations.
Because the previously obtained solution has a certain accuracy and the time cost for executing the iterative refinement is considered to be higher, the iterative refinement operation is executed only once, namely, the steps 3-6-1 to 3-6-3 are executed only once, so that x' is the solved solution.
And obtaining grid boundary data of flow field spatial distribution according to the solution x' after iterative optimization updating.
Step (4), updating the topological structure of the airflow grid according to the grid boundary data of the flow field spatial distribution in the step (3); as fig. 7 specifically shows:
4-1 dynamic mesh generation (mesh node update);
4-2, updating geometric data (calculating the normal speed of the control surface according to a geometric conservation law);
4-3 establishing an overlapping grid relation;
4-4 outer iteration sub-trellis interpolation.
And (5) repeating the operations of grid establishment, flow field calculation in the grid and grid updating in the non-constant-field fluid calculation until the flow field is converged.
Therefore, the aim of numerical simulation of the helicopter flow field can be achieved through the data obtained after the flow field is solved and the updating of the grid boundary data.
Claims (7)
1. An optimization method based on computational fluid dynamics numerical simulation is characterized by comprising the following steps:
the method comprises the following steps of (1) obtaining geometrical parameters of a helicopter body and physical parameter information of airflow in a flight process;
step (2), constructing an airflow grid topological structure according to the data acquired in the step (1); wherein each grid contains airflow physical parameter information for a current location;
step (3), each grid in the airflow grid topological structure carries out independent flow field calculation to obtain grid boundary data of flow field spatial distribution;
3-1, establishing a three-dimensional depressible Reynolds average N-S equation on each grid in the airflow grid topological structure:
wherein W is a conservative variable, F (W) and FvConvective flux and viscous flux, respectively;
the partial differential equation set on each grid and the final result of the N-S equation can be approximated to a linear algebraic equation set, and the linear equation set is as follows:
Ax＝b (1)
the matrix dimension of the coefficient matrix A is n, nb is the block size, and x is the equation unknown quantity;
the coefficient matrix A is composed of data blocks with the size of nb × nb, and total nt × nt data blocks;
extracting a coefficient matrix A of an algebraic equation;
3-2 randomizing the coefficient matrix A by SRBT method, and processing the matrix ArTransforming into a tile (tile) shape;
3-3 decomposing tile (tile) like matrix A by multi-wavefront methodrExtracting a wavefront matrix Q;
3-4, performing LDLT decomposition on the wavefront matrix Q by adopting an improved LDLT tile algorithm, and parallelizing by utilizing OpenMP task constraint to obtain a lower triangular matrix L and a diagonal matrix D;
the LDLT decomposition decomposes the matrix Q into:
Q＝LDLT (10)
wherein L, D each element in the matrix represents a tile block of nb by nb; t refers to transposing the matrix;
the improved LDLT tile algorithm is specifically as follows:
a) initializing k to 1;
b) wave front matrix Q diagonal block Q by LDLT algorithmkkDecomposing to obtain a lower triangular matrix LkkAnd diagonal matrix Dkk：
Wherein k is more than or equal to 1 and less than or equal to nt, and T represents transposition;
c) obtaining the lower triangular matrix L according to the abovekkAnd diagonal matrix DkkCalculating a matrix L of all tile blocks in the kth column of the wavefront matrix Qik；
Wkk＝LkkDkk (15)
Wherein k is not less than i and not more than nt;
d) adopting an improved ACA algorithm to compress the tile blocks;
lower triangle for all tile blocks in the k column of step c)Matrix LikDecomposing the matrix into X and Y low-rank matrixes; if the precision reaches a low-rank threshold value, or the number of vector groups in the low-rank matrix reaches a cycle threshold value, ending the cycle, otherwise, returning to generate a new vector group, adding the new vector group to X and Y, and generating a new X and Y low-rank matrix;
wherein k +1 is more than or equal to i and less than or equal to nt;
e) update Qij
Judging whether the rank of the low-rank matrix X, Y after the compression in the step d) is more than or equal to a cycle threshold value, if so, using an original full-rank matrix LikUpdate Qij：
If not, substituting equation (17) into equation (18), updating Q with the low rank matrixij：
Wherein k +1 is more than or equal to i and less than or equal to nt, and k +1 is more than or equal to j and less than or equal to i;
f) judging whether k under the current iteration is nt, if so, finishing the decomposition of all tile blocks, otherwise, updating k to k +1, and returning to the step b);
3-5, obtaining the solution of the final equation set (1) through forward substitution and reverse substitution;
3-6 use one iteration to optimize and promote the accuracy of the solution:
3-6-1 calculating the residual r according to the solution x given by equation (23):
r＝b-Ax (24)
3-6-2 calculating the residual r from the solution x obtained from equation (24):
Ad＝r (25)
wherein d is a correction value;
3-6-3 update correction value
x′＝x+d (26)
Wherein m represents the number of iterations;
obtaining grid boundary data of flow field spatial distribution according to the solution x' after iterative optimization updating;
step (4), updating the topological structure of the airflow grid according to the grid boundary data of the flow field spatial distribution in the step (3);
and (5) repeating the steps (2), (3) and (4) until the flow field is converged.
2. The optimization method based on computational fluid dynamics numerical simulation of claim 1, wherein the physical parameter information of the air flow in step (1) comprises the incoming flow speed, the attack angle and the like during the construction of the helicopter flight.
3. The optimization method based on computational fluid dynamics numerical simulation according to claim 1 or 2, wherein the coefficient matrix A randomization preprocessing of the step (3-2) is specifically:
Ara random matrix after coefficient matrix a randomization:
Ar＝UAUT (11)
wherein T represents transposition, U represents random matrix, and A after randomizationrThe matrix size of (a) is still equal to the original matrix A;
the random matrix U is a recursive butterfly matrix:
U＝Ud×...×U1 (12)
wherein d represents depth, Uq(1. ltoreq. q. ltoreq. d) represents a block diagonal matrix generating the random matrix U:
wherein B ispRepresenting a butterfly matrix of size n/2q-1(1≤p≤2q-1) Remember n1＝n/2q-1(ii) a Then the size is n1×n1The butterfly matrix of (a) is defined as:
wherein R ispAnd SpRespectively representing a random diagonal matrix and a nonsingular matrix with the size of n1/2×n1/2；RpThe value of the middle element is represented by eρ/10Is generated, ρ is
Random values within a range, thus the butterfly matrix BpThe element range of (B) is as follows:
e-1/20≤rp≤e1/20 (15)。
4. the optimization method based on computational fluid dynamics numerical simulation according to claim 3, wherein the steps (3-5) are specifically:
obtaining the decomposition result of the wave front matrix Q from the step (3-4)
Q＝LDLT (16)
Synthesizing the partial decomposition results of a series of wave front matrixes Q to obtain the decomposition result of the original coefficient matrix A;
the LDL of the original coefficient matrix A is obtained by subsequent calculation of a multi-wave-front methodTThe result of decomposition is
Wherein L ismergeAnd DmergeRespectively, matrix A LDLTA lower triangular matrix and a diagonal matrix of the decomposition result;
according to the formula (11):
Ax＝UTbU (18)
with reference to formulas (17) to (18), it can be seen that:
the intermediate variable y is found, then equation (8) is substituted into (19) to find the solution x:
x＝Uy (20)。
5. an optimization system for numerical simulation of a helicopter flow field is characterized by comprising a data acquisition module, a data preprocessing module, a flow field solving module and a grid boundary data updating module.
The data acquisition module is used for acquiring the geometrical parameters of the helicopter body and the physical parameter information of airflow in the flight process;
the data preprocessing module is used for constructing an airflow mesh topological structure; the airflow grid topological structure is subjected to grid division according to an overlapping grid method, and each grid contains airflow physical parameter information of the current position;
the flow field solving module is used for carrying out numerical solution on the flow field in each grid of the airflow grid topological structure to obtain flow field data;
the grid boundary data updating module is used for updating the grid nodes of the airflow grid topological structure according to the flow field data obtained by the flow field solving module, so that the purpose of numerical simulation of the flow field is achieved.
6. A computer-readable storage medium, on which a computer program is stored which, when executed in a computer, causes the computer to carry out the method of any one of claims 1-4.
7. A computing device comprising a memory having executable code stored therein and a processor that, when executing the executable code, implements the method of any of claims 1-4.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011632168.6A CN112632874A (en) | 2020-12-31 | 2020-12-31 | Optimization method and system for numerical simulation of helicopter flow field |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011632168.6A CN112632874A (en) | 2020-12-31 | 2020-12-31 | Optimization method and system for numerical simulation of helicopter flow field |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN112632874A true CN112632874A (en) | 2021-04-09 |
Family
ID=75289878
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202011632168.6A Pending CN112632874A (en) | 2020-12-31 | 2020-12-31 | Optimization method and system for numerical simulation of helicopter flow field |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112632874A (en) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113392472A (en) * | 2021-08-17 | 2021-09-14 | 北京航空航天大学 | OpenMP parallel disturbance domain updating method for aircraft aerodynamic characteristic simulation |
| CN114925627A (en) * | 2022-05-12 | 2022-08-19 | 南京航空航天大学 | Helicopter flow field numerical simulation system and method based on graphic processor |
| CN115700572A (en) * | 2021-07-23 | 2023-02-07 | 合肥本源量子计算科技有限责任公司 | Quantum solving method and device for iteration of grid equation and physical equation |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108920811A (en) * | 2018-06-28 | 2018-11-30 | 南京航空航天大学 | A kind of analogy method and system for helicopter flight simulation |
| CN109088776A (en) * | 2018-09-12 | 2018-12-25 | 西安交通大学 | A kind of parallel C FD calculating optimization method Chong Die with communication based on supercomputer |
| CN110096838A (en) * | 2019-05-16 | 2019-08-06 | 杭州电子科技大学 | A kind of helicopter flow field numerical value Parallel Implicit method for solving based on N-S equation |
-
2020
- 2020-12-31 CN CN202011632168.6A patent/CN112632874A/en active Pending
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108920811A (en) * | 2018-06-28 | 2018-11-30 | 南京航空航天大学 | A kind of analogy method and system for helicopter flight simulation |
| CN109088776A (en) * | 2018-09-12 | 2018-12-25 | 西安交通大学 | A kind of parallel C FD calculating optimization method Chong Die with communication based on supercomputer |
| CN110096838A (en) * | 2019-05-16 | 2019-08-06 | 杭州电子科技大学 | A kind of helicopter flow field numerical value Parallel Implicit method for solving based on N-S equation |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115700572A (en) * | 2021-07-23 | 2023-02-07 | 合肥本源量子计算科技有限责任公司 | Quantum solving method and device for iteration of grid equation and physical equation |
| CN113392472A (en) * | 2021-08-17 | 2021-09-14 | 北京航空航天大学 | OpenMP parallel disturbance domain updating method for aircraft aerodynamic characteristic simulation |
| CN113392472B (en) * | 2021-08-17 | 2021-11-09 | 北京航空航天大学 | OpenMP parallel disturbance domain updating method for aircraft aerodynamic characteristic simulation |
| CN114925627A (en) * | 2022-05-12 | 2022-08-19 | 南京航空航天大学 | Helicopter flow field numerical simulation system and method based on graphic processor |
| CN114925627B (en) * | 2022-05-12 | 2024-03-15 | 南京航空航天大学 | Helicopter flow field numerical simulation system and method based on graphic processor |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112632874A (en) | Optimization method and system for numerical simulation of helicopter flow field | |
| CN110110413B (en) | Structural topology optimization method based on material field reduction progression expansion | |
| JP6784780B2 (en) | How to build a probabilistic model for large-scale renewable energy data | |
| CN105808829A (en) | CPU+GPU heterogeneous parallel computing based natural frequency characteristic analysis method for turbomachinery blade | |
| CN106570204A (en) | Method for analyzing static strength characteristics of turbomachinery blade based on CPU+GPU heterogeneous parallel computing | |
| CN103345580B (en) | Based on the parallel CFD method of lattice Boltzmann method | |
| CN107403466A (en) | Ultra-large unstrctured grid generation method based on overall situation encryption | |
| CN112597610B (en) | Optimization method, device and equipment for lightweight design of mechanical arm structure | |
| CN116436012B (en) | FPGA-based power flow calculation system and method | |
| CN101980182A (en) | Matrix operation-based parallel computing method | |
| JP5316433B2 (en) | Optimization processing program, method and apparatus | |
| CN111489447A (en) | Right-angle grid adaptive modeling method suitable for lattice Boltzmann method | |
| CN111079326B (en) | Two-dimensional anisotropic grid cell measurement tensor field smoothing method | |
| CN114676522B (en) | Pneumatic shape optimization design method, system and equipment integrating GAN and migration learning | |
| CN114925627A (en) | Helicopter flow field numerical simulation system and method based on graphic processor | |
| CN109840348B (en) | Triple acceleration topology optimization method | |
| Watanabe et al. | Optimizing memory layout of hyperplane ordering for vector supercomputer SX-Aurora TSUBASA | |
| Kulikov | The numerical modeling of the collapse of molecular cloud on adaptive nested mesh | |
| CN113239591A (en) | DCU cluster-oriented large-scale finite element grid parallel partitioning method and device | |
| CN117421703A (en) | Depth sign regression accelerator and depth sign regression method | |
| CN109948253A (en) | The GPU accelerated method of thin plate mesh free Galerkin Constructional Modal Analysis | |
| CN116303219A (en) | Grid file acquisition method and device and electronic equipment | |
| CN110021339B (en) | Cluster parallel computing acceleration method based on protein folding calculation protein structure | |
| CN113609585B (en) | Multi-level model construction method and device for ship power system and electronic equipment | |
| Gao et al. | A multi-level parallel tie-dye algorithm for auto-CFD |
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 |