Disclosure of Invention
In order to solve the problems in the prior art, meet the numerical simulation requirement of high precision and high efficiency and simultaneously meet the physical practical situation, the invention synthesizes the three curved surface processing formats based on the practical numerical simulation calculation requirement, provides a novel curved surface boundary optimization processing method suitable for LBM numerical simulation, and can be used for developing a calculation platform meeting the numerical simulation of complex flow fields of high precision and high efficiency. Based on the method, corresponding airfoil bypass numerical simulation and engine inner cavity inner flow numerical simulation methods are further provided, the actual application fields of high-performance airfoil design, new generation engine design and the like are carried out, the precision and efficiency of airfoil bypass numerical simulation and engine inner cavity inner flow numerical simulation are improved, and the development and innovation of the subject are promoted.
The technical scheme of the invention is as follows:
a processing method for curved surface boundary in lattice Boltzmann method numerical simulation solves a point f edge of a point f in an object surface boundary through the following processDistribution function of direction:
Step1, point f is followedThe distribution function of the direction is divided into equilibrium distribution functionsAnd an unbalanced distribution function:
Wherein the equilibrium distribution functionBy a distribution function of equilibrium statesObtaining a calculation formula;
step 2, solving the unbalanced distribution function of the point f by adopting the following steps:
Step 2.1 integration of information of points ff and w using linear interpolation
Wherein the method comprises the steps ofIs the embedding depth;
step 2.2, synthesizing the information of points fff and ff by means of linear extrapolation
Step 2.3, obtaining the average linear internal and external interpolation result:
wherein the information of w is calculated by the following procedure:
The distribution function of the point w is divided into an unbalanced state and an equilibrium state
And is determined by object plane slip-free condition limitationObtainingWherein the equilibrium distribution function of the w pointBy a distribution function of equilibrium statesObtaining a calculation formula;
simultaneously using linear internal and external interpolation to solve for the distribution function of the point w
Thereby obtainingWherein the point b is alongThe distribution function of the direction is followed by the f point after collision migrationThe distribution function of the direction is migrated from the point f alongThe distribution function of the direction is followed by the ff point after collision migrationThe distribution function of the direction is migrated from the ff point alongThe distribution function of the direction is followed by the fff point after collision migrationThe distribution function of the direction is migrated;
Thereby obtaining the following steps:
And then obtain
。
Further, in step 1, an equilibrium distribution function is calculatedWhen the density and speed of the point f are interpolated by the points around the point f:
。
Further, in step 3, the velocity and density of the point w are interpolated by using the information of the points ff and f, and then the equilibrium distribution function of the point w is obtained by calculating the equilibrium distribution function。
Advantageous effects
The processing method for the curved surface boundary in the numerical simulation of the lattice Boltzmann method can meet the curved surface boundary processing problem encountered in the numerical simulation of an LBM complex flow field and a complex appearance, is more accurate than the three existing curved surface boundary processing formats, and has higher convergence speed due to more accurate calculation in each step, thereby realizing high precision and high efficiency of the numerical simulation and obtaining better application effects in the practical application of airfoil design, engine design and the like.
Additional aspects and advantages of the invention will be set forth in part in the description which follows, and in part will be obvious from the description, or may be learned by practice of the invention.
Detailed Description
Aiming at the problem of processing the boundary of the curved surface object in complex shape numerical simulation based on Lattice Boltzmann Method (LBM) algorithm, the invention provides a novel curved surface boundary processing method, which not only can process the external flow fields such as airfoil bypass and the like in actual engineering calculation, but also is suitable for numerical simulation of the internal flow field of the cavity with the curved surface physical boundary.
The new curved surface boundary processing method is described below by taking cylindrical bypass numerical simulation, NACA0012 airfoil turbulent flow field numerical simulation, SD7003 airfoil turbulent flow field numerical simulation and U cavity internal flow field numerical simulation as examples respectively.
Example 1, cylindrical bypass numerical simulation:
when the cylindrical bypass numerical simulation is carried out, after CFD conventional operations such as a cylindrical three-dimensional model, grid division and the like are established, a distribution function is solved for points in the boundary of a cylindrical object plane through the following processes:
as shown in figure two, along with point fDistribution function of directionSolving for example, the points fff, ff and f are fluid points, and after each evolution, the information of the points fff and ff is updated. And point f is alongDistribution function of directionUnknown because it is migrated from virtual point b, which is a virtual point, there is no distribution function. We can therefore construct the f-point edge only from the information of points fff, ff and wDistribution function of direction。
First, the point f is followedThe distribution function of the direction is divided into equilibrium distribution functionsAnd an unbalanced distribution functionTo be processed.
(1)
Wherein the equilibrium distribution functionBy a distribution function of equilibrium statesThe calculation formula is obtained, so that we only need to know the density and the speed of the point f, and the density and the speed of the point f can be interpolated by the points around the point f in the method
(Black triangle marked dot)
Thus we can obtain the equilibrium distribution function of point fNext we find the non-equilibrium distribution function of the solution point f。
First we use the information of the linear interpolation integration points ff and w
(2)
Wherein the method comprises the steps ofTo embed depth, then the information of the linear extrapolation combining points fff and ff is utilized
(3)
Then average the result of the linear interpolation and the linear interpolation
(4)
Except for those in equation (4) after each evolutionThis is otherwise known, so we need to solve next。
First, we know that the distribution function of the point w can be divided into two parts of non-equilibrium state and equilibrium state
(5)
Or (b)
(6)
At the same time we know that there is no slip condition limit due to object plane
(7)
Thus, the formula (6) becomes
(8)
Since the equilibrium state distribution function of the w point can pass through the equilibrium state distribution functionThe calculation formula obtains that the speed and the density of the point w are interpolated by only using the information of the points ff and f, if the object plane is considered, the speed of the object plane is 0, and if the object plane is not stationary, the speed of the point w can be determined only according to specific conditions. Then we again use the linear interpolation inside and outside to solve the distribution function of the point w
(9)
(10)
Comprehensive formulas (9) and (10)
(11)
Wherein the point b is alongThe distribution function of the direction is followed by the f point after collision migrationThe distribution function of the direction is migrated from the same way that the f point is alongThe distribution function of the direction is followed by the ff point after collision migrationThe distribution function of the direction is migrated from the ff point alongThe distribution function of the direction is followed by the fff point after collision migrationThe distribution function of the direction is migrated. Bringing formula (11) into formula (8)
(12)
Substituting equation (12) into equation (4) we obtain
(13)
To this point f edgeDistribution function of directionAnd after the solution is finished, the information of the points fff, ff, f and w is comprehensively used to construct the unbalanced distribution function of the points, and meanwhile, the comprehensive information of the red circle punctuation is used to construct the balanced distribution function.
Compared with the Chen format, the invention distinguishes the unbalanced state and balanced state distribution functions when constructing the distribution functions. And the new format adopts linear interpolation value and extrapolation value to comprehensively consider the information and influence of the peripheral points. Compared with the Mei format and the Guo format, more reasonable surrounding points are selected as interpolation information sources of destination point speeds and densities in constructing the equilibrium state distribution, and speeds and densities of virtual points (which are not in agreement with physical reality) are not forcedly constructed. B point edge in the formatDistribution function of directionBut one concept for data exchange (assignment). Meanwhile, when the unbalanced distribution function is calculated, only the information of the point f is used as an interpolation information source by the Mei format and the Guo format, and the internal and external linear interpolation is selected to be more reasonable in the format, and meanwhile, the comprehensive information of peripheral multiple points is considered.
Table 1. Comparison of the results of the calculation of the values of the cylindrical bypass flow using the invention with other data (re=40)
Table 1 shows a comparison of flow data for steady flow at re=40 for cylindrical wrap flow using the novel surface boundary processing format proposed by the present invention with other results. It can be seen from the table that the calculation accuracy of the novel curved surface boundary processing format is superior to that of the other three popular curved surface boundary processing formats. Because the new format is different from other formats and is divided into two parts of pre-processing and post-construction, the time required for each iteration step is longer than that of other formats, but because the new format itself has high precision and is closer to a theoretical solution after each iteration step, convergence can be accelerated, so that the time consumed by the new format is basically the same as that of other three formats as a whole.
Fig. 3 shows a flow chart of steady flow at cylindrical bypass re=40, and two stable attached vortices can be clearly observed.
TABLE 2 comparison of Hopf flow bifurcation with other results using New Format cylindrical bypass flow
Table 2 shows predictions of cylindrical bypass Hopf flow bifurcation points for different documents and three popular formats in comparison with new format calculations. It can be seen that the calculation results of the new format are very consistent with those of other documents, indicating that the calculation accuracy of the new format is higher.
Table 3. Comparison of the calculation results of the cylindrical wrap flow using the new format with other literature (re=3900).
Fig. 4 and table 3 show the calculation results of the unsteady periodic flow state using the cylindrical wrap flow of the new format. Figure 4 (a) depicts a plot of lift coefficient versus time,,,,,A complete cycle is shown, with the instantaneous flow diagram for cylindrical wrap re=3900 shown in turn for each time node diagram 4 (b), (c), (d), (e), (f).
Example 2, a numerical simulation of the turbulent flow field of a naca0012 airfoil, comprising two states re=3000, angle of attack of 0 ° and re=100000, angle of attack of 0 °;
Example 3, turbulent flow field numerical simulation of an sd7003 airfoil, state re=60000, angle of attack 14 °;
in the embodiment 4,U, the numerical simulation of the flow field in the cavity is carried out, and the state is re=1000;
In examples 2,3 and 4, the distribution function was solved by the same method as in example 1 for the points in the boundary of the object plane of the NACA0012 airfoil, SD7003 airfoil and U-shaped cavity, and the obtained numerical simulation results are shown in FIG. 5 to FIG. 7.
Fig. 5 shows a numerical simulation of the turbulent flow field of an NACA0012 airfoil, re=3000, upper surface pressure coefficient at 0 ° angle of attack, flow chart and pressure cloud. Fig. 6 depicts a numerical simulation of the turbulent flow field of an NACA0012 airfoil, re=100000, coefficient of drag rise at 0 ° angle of attack, flow chart and pressure cloud. Fig. 7 shows a numerical simulation of the turbulent flow field of the SD7003 airfoil, re=60000, the flow chart at an angle of attack of 14 ° and the flow chart of the flow field in the U-cavity is also given (re=1000). It can be seen that the curved surface boundary format provided by the invention is suitable for numerical simulation of the wing profile bypass flow field and the U-shaped cavity inner flow field.
Although embodiments of the present invention have been shown and described above, it will be understood that the above embodiments are illustrative and not to be construed as limiting the invention, and that variations, modifications, alternatives, and variations may be made in the above embodiments by those skilled in the art without departing from the spirit and principles of the invention.