CN113792432A - Flow field calculation method based on improved FVM-LBFS method - Google Patents
Flow field calculation method based on improved FVM-LBFS method Download PDFInfo
- Publication number
- CN113792432A CN113792432A CN202111078887.2A CN202111078887A CN113792432A CN 113792432 A CN113792432 A CN 113792432A CN 202111078887 A CN202111078887 A CN 202111078887A CN 113792432 A CN113792432 A CN 113792432A
- Authority
- CN
- China
- Prior art keywords
- equation
- flux
- lbfs
- interface
- fvm
- 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 121
- 238000004364 calculation method Methods 0.000 title claims abstract description 43
- 230000004907 flux Effects 0.000 claims abstract description 60
- 230000008569 process Effects 0.000 claims abstract description 20
- 230000035939 shock Effects 0.000 claims abstract description 13
- 238000012937 correction Methods 0.000 claims abstract description 7
- 238000004458 analytical method Methods 0.000 claims abstract description 6
- 239000006185 dispersion Substances 0.000 claims abstract description 5
- 238000005315 distribution function Methods 0.000 claims description 39
- 239000002245 particle Substances 0.000 claims description 28
- 230000006870 function Effects 0.000 claims description 24
- 238000009795 derivation Methods 0.000 claims description 11
- 238000009826 distribution Methods 0.000 claims description 11
- 230000014509 gene expression Effects 0.000 claims description 10
- 239000000853 adhesive Substances 0.000 claims description 7
- 230000001070 adhesive effect Effects 0.000 claims description 7
- 239000012530 fluid Substances 0.000 claims description 7
- 150000001875 compounds Chemical class 0.000 claims description 5
- 230000005012 migration Effects 0.000 claims description 5
- 238000013508 migration Methods 0.000 claims description 5
- 239000000126 substance Substances 0.000 claims description 4
- 238000004519 manufacturing process Methods 0.000 claims description 2
- 239000000203 mixture Substances 0.000 claims description 2
- 238000005381 potential energy Methods 0.000 claims description 2
- 230000007704 transition Effects 0.000 claims description 2
- 238000006467 substitution reaction Methods 0.000 claims 1
- 238000004088 simulation Methods 0.000 abstract description 18
- 238000012795 verification Methods 0.000 abstract description 2
- 238000011160 research Methods 0.000 description 6
- 230000008901 benefit Effects 0.000 description 4
- 238000011161 development Methods 0.000 description 2
- 238000010438 heat treatment Methods 0.000 description 2
- 238000000329 molecular dynamics simulation Methods 0.000 description 2
- 229910001220 stainless steel Inorganic materials 0.000 description 2
- 239000010935 stainless steel Substances 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- CLOMYZFHNHFSIQ-UHFFFAOYSA-N clonixin Chemical compound CC1=C(Cl)C=CC=C1NC1=NC=CC=C1C(O)=O CLOMYZFHNHFSIQ-UHFFFAOYSA-N 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000010534 mechanism of action Effects 0.000 description 1
- 230000010355 oscillation Effects 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000012546 transfer Methods 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
-
- 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
- G06F17/12—Simultaneous equations, e.g. systems of linear equations
-
- 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/08—Thermal analysis or thermal optimisation
-
- 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)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Theoretical Computer Science (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Computational Mathematics (AREA)
- Data Mining & Analysis (AREA)
- General Engineering & Computer Science (AREA)
- Geometry (AREA)
- Evolutionary Computation (AREA)
- Computer Hardware Design (AREA)
- Operations Research (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention discloses a flow field calculation method based on an improved FVM-LBFS method, and belongs to the field of aircraft pneumatic calculation. The method combines a traditional finite volume method with an LBM method, adopts the finite volume method to carry out space dispersion on a macroscopic N-S control equation, utilizes the inviscid flux of a local solution reconstruction unit interface of a one-dimensional compressible lattice Boltzmann model, and has simple and efficient calculation process. The improved FVM-LBFS method realizes the accurate control of the viscosity of the inviscid flux value in the existing LBFS method by adopting a newly proposed improved switch control function for controlling the pressure and the temperature simultaneously, and introduces a grid unit slenderness ratio correction coefficient to expand the application range of the non-uniform grid. Through a series of verification and analysis of numerical simulation examples, the improved FVM-LBFS method can simultaneously stably capture the flow of complex strong shock waves and accurately predict the aerodynamic thermal parameters in the numerical simulation of hypersonic flow.
Description
Technical Field
The invention discloses a new flow field calculation method based on an improved FVM-LBFS method, and belongs to the field of aircraft pneumatic calculation.
Background
In the later period of the 20 th century and the 80 s, along with the rapid development of computer technology, the running speed of a computer is faster and faster, and a numerical simulation method for fluid mechanics is developed from a mesoscopic level and is gradually started. The Lattice Boltzman Method (LBM), one of the typical representatives, has received considerable attention from researchers since 1988, and has been rapidly developed in a period of approximately 20 years. Different from an N-S equation based on continuity assumption, the lattice Boltzman method is a mesoscale fluid computational mechanics method between a fluid micro-molecular dynamics calculation method and a macroscopic continuity equation calculation method, so that compared with the traditional CFD method, the lattice Boltzman method is not limited by continuous medium assumption and can realize cross-scale flow simulation.
The lattice Boltzman equation is obtained by dispersing the Boltzman equation in a velocity space, a physical space and time, and describes the process of a particle moving to an adjacent point in the next time step at a certain velocity. One advantage of the lattice Boltzmann method is that: the method solves the difficulty of solving nonlinear partial differential equation sets (such as Euler equation sets, N-S equation sets and the like) by the traditional method, and is completed once during model establishment, so that only lattice Boltzmann equation simple linear equations or equation sets need to be processed in numerical simulation, and compared with a macroscopic method, the number and the form of the equations are greatly simplified. The solving process of the standard lattice Boltzmann method only comprises two parts, namely a simple linear migration process and a simple collision process, only algebraic operation is needed, and a pressure Poisson equation does not need to be solved, so that the program is very simple to implement. Given the many unique advantages of the lattice Boltzmann method, it has evolved into a new computational tool that simulates complex flow and physical phenomena. Therefore, in many research fields, such as compressible/incompressible flow, hypersonic flow, multiphase flow, porous media flow, chemical reaction flow, lean gas flow, etc., the numerical method based on the lattice Boltzman model is well applied, and the physical mechanism of action of many complex flow problems is revealed.
The implementation of the Lattice Boltzmann method is mainly based on two parts of a Lattice Boltzmann Model (Lattice Boltzmann Model) and a corresponding numerical method, a discrete speed Model and a corresponding balanced distribution function form in the Lattice Boltzmann Model are key, and a large amount of research work is done by researchers at home and abroad at present. The traditional lattice Boltzmann model is mostly built for incompressible flows, such as the DdQm series basic model proposed by Qian et al, 1992 (d denotes the spatial dimension and m denotes the number of discrete velocities). In order to expand the application of the LBM method in compressible flow, a compressible lattice Boltzmann model is also continuously proposed, such as a multi-speed model, a specific heat ratio adjustable model, a coupled double-distribution model and the like. Based on these models, the LBM method achieves numerical simulations of Ma >0.3 flows, but many lattice Boltzmann models are currently only applicable to low mach number compressible flows (flows with significant changes in aerodynamic density due to thermodynamic processes). This is because these models are mostly built based on the taylor series expansion of Maxwellian distribution functions, and the truncation error is directly related to mach number, resulting in unsuitability for high mach number flows. In order to overcome the problem, Qu et al uses a circular function to replace a maxwell distribution function, establishes a series of one-dimensional and two-dimensional models, such as a D1Q5L2 model and a D2Q13L2 model, and successfully applies the models to numerical simulation of a shock tube, a double-mach-number reflection and other strong shock wave supersonic flow fields. Deltar also derived some two-dimensional compressible lattice Boltzmann models, such as the non-split seven velocity model, the 4+3 split velocity model. However, these models contain many free parameters and the lattice velocity needs to be artificially set, which leads to poor versatility in application to complex high mach number flow numerical simulations. In view of the above, Shu and its research team provide a non-free parameter lattice Boltzmann model derivation platform based on a moment conservation form through theoretical derivation, and a series of non-free parameter lattice Boltzmann models for solving compressible non-viscous flow are obtained through the platform, such as a D1Q3 model, a D1Q4 model, a D2Q8 model, a D2Q12 model and the like. By utilizing the models, macroscopic equations can be locally and accurately recovered, and a one-dimensional compressible inviscid flow Riemann flux solver is established at the local unit interface.
Due to the requirement of grid uniformity, the traditional LBM method has very high requirement on the uniformity of a computational grid, and is not suitable for the computation of complex shapes. Moreover, for high Reynolds number flows, the computational grid is very dense due to the limitation of migration time step, resulting in a significant reduction in computational efficiency. Therefore, when numerical simulation of compressible flow is performed using the lattice Boltzmann method, most calculations are performed using a method based on the Discrete Velocity Boltzmann Equation (DVBE) to update the particle distribution function in time. Many scholars such as Qu, q.li, Ubertini and Stiebler have done a lot of research work on the DVBE method. However, it should be noted that in the DVBE method, the number of particle distribution functions is much greater than the number of conservative variables, and the form of the distribution functions is very complex, which leads to a complicated direct solving process of the DVBE method, and the calculation efficiency is much lower than that of the conventional CFD method. To develop a more efficient compressible flow solver with LBM properties, Shu and his research team proposed a Lattice Boltzmann Flux Solver (LBFS) based on a finite volume method. It is well known that the Finite Volume Method (FVM) has its unique advantages in CFD numerical simulation: the FVM starts from a conservation equation in an integral form, the physical concept of the derivation process is clear, and the conservation characteristic of the physical quantity is kept in each control body; and the method has good flexibility for computing grids (structural/non-structural grids), and is suitable for computing complex geometric shapes. Therefore, in the LBFS method, the complex process of directly solving the DVBE equation is avoided, the macroscopic control equation is discretized by utilizing the advantages of the traditional FVM, and the non-viscous flux of the interface of the local solution reconstruction unit of the one-dimensional compressible lattice Boltzmann model is adopted. Numerical calculation results show that the LBFS method has good positivity and high efficiency in the strong shock wave compressible flow simulation. In the conventional LBFS method, in order to capture strong discontinuous physical flow such as strong shock waves, numerical dissipation is introduced by using an unbalanced part of a unit interface distribution function to ensure the stability of calculation. However, the amount of numerical dissipation directly introduced in this way is too large to seriously disturb the true solution of smooth regions such as boundary layers, and is not acceptable particularly for accurately predicting the heat flow density of the hypersonic flow field. Therefore, how to ensure that the LBFS method can stably capture the strong shock wave flow and accurately predict the aerodynamic heat in the hypersonic flow at the same time needs to be further improved, which is also the main content of the research of the present invention.
Aiming at the defect that the existing numerical method based on the lattice Boltzmann model simulates high Mach flow, the invention provides an improved FVM-LBFS method for numerical simulation of hypersonic flow and provides a detailed derivation process. The method adopts a finite volume method to carry out space dispersion on a macroscopic N-S control equation, utilizes a compressible LB model to construct a local one-dimensional Riemann flux solver, and has a simple and efficient calculation process. The improved FVM-LBFS method realizes the accurate control of the viscosity of the inviscid flux value in the existing LBFS method by adopting a newly proposed improved switch control function for controlling the pressure and the temperature simultaneously, and introduces a grid unit slenderness ratio correction coefficient to expand the application range of the non-uniform grid. And finally, selecting a series of high-speed numerical value examples to verify and analyze the reliability and accuracy of the provided calculation method.
Disclosure of Invention
The invention mainly adopts a numerical calculation means, and develops a new flow field calculation method based on an improved FVM-LBFS method around the hypersonic pneumatic heating problem.
Aiming at the problem that the high Reynolds number flow is difficult to simulate by the existing numerical method based on the lattice Boltzmann model, an improved FVM-LBFS method is provided for numerical simulation of hypersonic flow. The method combines a traditional finite volume method with an LBM method, adopts the finite volume method to carry out space dispersion on a macroscopic N-S control equation, utilizes the inviscid flux of a local solution reconstruction unit interface of a one-dimensional compressible lattice Boltzmann model, and has simple and efficient calculation process. The improved FVM-LBFS method realizes the accurate control of the viscosity of the inviscid flux value in the existing LBFS method by adopting a newly proposed improved switch control function for controlling the pressure and the temperature simultaneously, and introduces a grid unit slenderness ratio correction coefficient to expand the application range of the non-uniform grid. Through a series of verification and analysis of numerical simulation examples, the improved FVM-LBFS method can simultaneously and stably capture the flow of complex strong shock waves and accurately predict the aerodynamic thermal parameters in the numerical simulation of the hypersonic flow, and the application range of the numerical method based on the lattice Boltzmann model in the field of the numerical calculation of the hypersonic flow is well expanded.
The technical scheme of the invention is as follows:
a flow field calculation method based on an improved FVM-LBFS method comprises the following steps:
the average Reynolds N-S (RANS) equation is used as the control equation, and the integral form (passive term) is as follows,
in the formula: w is a conservative value; fcNo adhesive flux; fvIs the viscous flux; w, Fc、FvIs defined as:
wherein ρ, p and T represent densities,Pressure, temperature, and satisfy the ideal gas state equation p ═ ρ RT, R is the ideal gas constant. u, v, w are velocities of the fluid in the body in three directions respectively, nx、ny、nzRespectively, the outer normal unit vector component of the control body unit surface dS, where V is nxu+nyv+nzw is the surface normal velocity, τijIs a friction stress tensor component, k is a heat conduction coefficient, E is p/rho (gamma-1) is the energy of the fluid in the control body, H is the total enthalpy, and gamma is the specific heat ratio; h is E + p/rho.
From the above, it can be seen from equation (1) that the key step in solving it is to solve for the inviscid flux FcAnd a viscous flux FvAnd (4) calculating. The invention adopts the LBFS method based on the non-free parameter D1Q4 model to calculate the non-adhesive flux FcConsidering the physical property of the viscous flux, the viscous flux F is calculated by using a central difference methodv. The solution of the LBFS method to the inviscid flux F will be described and derived in detail belowcThe calculation process of (2).
The LBFS method is established based on a lattice Boltzmann model, and in order to accurately simulate hypersonic flow, reasonable selection of a compressible LB model is very important for the LBFS method. The invention selects a D1Q4 model with non-free parameters to be used in the LBFS method of the invention, and the derivation process of the D1Q4 model is briefly explained as follows:
to recover the one-dimensional Euler equation by using the Maxwell' S distribution function to obtain the relationship of the moments required to recover the Euler equation or the N-S equation, the conservation of momentum relationship can be simplified to the form of a summation,
in the formula (I), the compound is shown in the specification,
as a function of the equilibrium distribution of the particles in the i-th direction, ξiThe particle velocities in the respective directions, p and u being the density and velocity, respectively, of the one-dimensional macroscopic flow, where
Indicating a particular velocity of the moving particles. When the derivation of the lattice Boltzmann model is performed, the particle velocity distribution shape of the LB model needs to be determined first. It can be seen that the compressible D1Q4 model consists of two particle velocities (D)1，d2) And four equilibrium distribution functions
Composition, there are thus a total of 6 unknown parameters to be claimed. Since the equation (1) only contains 4 independent equations, in order to solve the 6 unknowns, the following two independent equations need to be introduced
Equations (3) and (4) are two higher order conservation of momentum relationships, and these two additional higher order moment relationships can further constrain the compressible D1Q4 model, so that it can reflect the physical flow properties more realistically, and especially has better positivity for the computation of the strong shock wave problem in hypersonic flow.
Thus, simultaneous equations (2) - (4), the following form of the system of equations can be obtained,
solving equation set (5) can obtain 4 equilibrium distribution functions and 2 discrete particle velocity expressions in the compressible D1Q4 model as
FVM-LBFS establishes a connection between an N-S equation and a non-free parameter D1Q4 model, and according to a physical conservation relation (2) for restoring an Euler equation, a conservation variable related to a normal direction speed
And no adhesive flux FcCan be calculated by the following formula,
in the formula, xiiThe discrete particle phase velocity in the ith direction is shown as xi in the one-dimensional non-free parameter D1Q4 model1＝d1，ξ2＝-d1，ξ3＝d2，ξ4＝-d2。VnIs the interface normal velocity, e is the unit mass internal energy;
then represents a moment relation vector, whose expression is
In addition, epIs the potential energy of particles, fi(r, t) represents the particle distribution function of the ith direction at the spatial position r at the time t, and the whole conservation variable W and the non-viscous flux F can be respectively obtained by combining the above equationscConservation variable with respect to normal direction velocity
And no adhesive flux
The relationship between them is as follows
Wherein the content of the first and second substances,
are respectively the corresponding items 1, 2, 3, u in the formula (8)τIs tangential direction velocity, uτx、uτy、uτzAre each uτComponent in three directions,
Are respectively the corresponding items 1, 2 and 3 in (9);
in the improved FVM-LBFS method, the core part is derivation and solution of the inviscid flux, the invention starts from Chapman-Enskog development analysis, a balanced part and an unbalanced part of a particle distribution function on a unit interface are simultaneously reserved in the LBFS method, the LBFS method is only used for calculating the inviscid flux of the unit interface, the contribution of the unbalanced part of the distribution function on the unit interface to the flux can be regarded as artificial numerical stickiness, the control of the numerical stickiness needs to be accurately controlled according to the physical characteristics of a flow field, and the derivation process of inviscid flux calculation is introduced as follows:
in the present LBFS, the lattice Boltzmann model is used to build a local Riemann flux solver at the cell interface to compute the inviscid flux. A non-viscous flux with respect to normal velocity at the cell interface can be obtained
The expression is as follows:
in the formula, the symbol "+" represents a numerical solution, the subscript "j + 1/2" represents the cell interface position, the upper subscript I represents the contribution of the equilibrium distribution function at the cell interface to the no-stick flux, and the upper subscript II represents the contribution of the equilibrium distribution function at the points around the cell interface to the no-stick flux. Further, by incorporating the contribution of tangential directional velocity to the inviscid flux into equation (12), the complete expression of inviscid flux at the cell interface can be written as follows
δ t is the particle migration time step, and the present invention treats the contribution of the non-equilibrium portion of the distribution function as the numerical viscosity, τ0Can be used as a weight coefficient of numerical viscosity;
the density, normal momentum and normal velocity at the cell interface contribute to energy production as follows
From the compatibility condition, the calculation of the conservation variable is independent of the non-equilibrium portion of the particle distribution function. Thus, the unbalanced portion of the distribution function in equation (14) can be written as,
it is possible to obtain,
as can be seen from equation (16), the conservation variable at the cell interface can be represented by the equilibrium distribution function gi(-ξiδ t, t- δ t). Further, we can get the expression of the conservation variable on the cell interface,
in the formula (I), the compound is shown in the specification,
respectively representing the balance distribution functions of the left and right of the unit interface;
since the one-dimensional non-free parametric LB model is applied in equation (17), only the conservation variable related to the normal velocity on the cell interface can be obtained. In order to calculate the contribution to tangential velocity on the cell interface, the present invention employs a feasible approximation method, as shown below,
in the formula (I), the compound is shown in the specification,
the tangential velocities on the cell interface left side, the cell interface right side, and the cell interface, respectively. Wherein the content of the first and second substances,
and
the difference between the total speed and the normal speed of the left side and the right side of the reconstructed unit interface can be calculated respectively. Through simultaneous equations, the conservation variable rho on the unit interface can be directly obtained through calculation*，
And p*. Further, the calculated conservation variable rho*，
And p*Substituting equations (6) - (7) can calculate the equilibrium distribution function g on the cell interfacei(0, t). Function g is distributed with equilibriumi(-ξiδ t, t- δ t) and gi(0, t) is determined and the corresponding inviscid flux at the cell interface can be calculated explicitly.
In order to accurately control the magnitude of the numerical viscosity generated by the unbalanced portion of the cell interface distribution function in the whole calculation domain, a global switch function is required to be constructed to calculate the numerical viscosity weight coefficient tau0The magnitude of the value of (c).
Since the accuracy of the heat flow calculation in hypersonic flow depends on the accurate capture of the shock and temperature boundary layers, which are very sensitive to the numerical viscosity in the calculation. According to the theory of molecular dynamics, dimensionless relaxation times τ0Is both pressure and temperature dependent, and therefore the numerical viscosity weight coefficient τ in the domain is calculated for more accurate control0The invention simultaneously introduces the temperature value and the pressure value of the control body units on the left and the right sides of the unit interface to the weight coefficient tau0The value of (c) is controlled. Based on the above theoretical analysis, a new improved switch control function is provided to calculate tau0The value of (b) is defined as follows,
wherein, tanh (x) is hyperbolic tangent function in hyperbolic function, and the function value of tanh (x) can be in [0, 1%]And smooth transition between them. p is a radical ofLAnd pRRespectively the pressure values at the left and right sides of the unit interface, and C is a relaxation factor constant with a value range of [1,100 ]]，pL，pRAnd TL，TRRespectively are pressure value and temperature value at left and right sides of the unit interface, delta is correction coefficient of length-to-thickness ratio of grid unit, omegaL、ΩRThe volumes of the control body units on the left side and the right side of the unit interface are respectively,
then the average volume of the left and right control volume units, LL,RIs the distance between the centers of the left and right control units. N is a radical offLAnd NfRRespectively representing the surface numbers of the grid cells on the left side and the grid cells on the right side of the cell interface,
respectively, the weighting coefficients on the ith plane.
The specific solving process is as follows:
(1) calculating the derivative of the conservation variable of each grid unit, and obtaining the conservation variables on two sides of the unit interface through linear reconstruction;
(2) substituting the conservation variables on both sides of the cell interface into equations (6) - (7) to calculate the equilibrium distribution function
And
(3) by equation (17)(18) The conservation variable at the cell interface is calculated, and further, the original variable ρ is calculated using equation (17)*，
And p*；
(4) Calculating a numerical viscosity weight coefficient τ using equation set (19)0Then the total non-viscous flux at the cell interface is calculated by equation (13)
(5) Calculating F of viscous flux at cell interfacev；
(6) Solving an N-S equation by a time advancing method, wherein a new conservation variable value of the center of the grid unit at the next time step can be calculated in the step;
(7) and (4) repeating the steps (1) - (7) until a convergence solution meeting the condition is obtained.
The invention has the following beneficial effects:
1. the invention discloses a compressible lattice Boltzmann model for hypersonic velocity aerodynamic heat calculation
A non-free parametric one-dimensional D1Q4 model is used. The distribution function of the model is simple in structure and high in calculation efficiency, and the particle speed is determined through a high-order moment of momentum relation without manually giving free parameters, so that the model is more reasonable.
2. The invention discloses a novel grid effect-based numerical value viscosity control function
The invention simultaneously introduces the temperature value and the pressure value of the control body units on the left side and the right side of the unit interface to the weight coefficient tau0The value of (c) is controlled. Based on the above theoretical analysis, a new improved switch control function is provided to calculate tau0The value of (a) is,
3. the invention discloses a new flow field calculation method based on an improved FVM-LBFS method
In order to expand the application range of a lattice Boltzmann model-based numerical method in a hypersonic flow area, the invention provides an improved FVM-LBFS method for numerical simulation of hypersonic flow. The method adopts a finite volume method to carry out space dispersion on a macroscopic N-S control equation, adopts an LBFS method based on a non-free parameter D1Q4 model to calculate the inviscid flux, and adopts a second-order central difference method to calculate the viscosity flux. The newly proposed improved switch control function for controlling pressure and temperature simultaneously realizes the accurate control of the viscosity of the non-viscous flux numerical value in the existing LBFS method, and introduces a grid unit slenderness ratio correction coefficient to expand the application range of the non-uniform grid.
Drawings
FIG. 1 is a computational grid.
Fig. 2(a) is a pressure cloud of the flow field.
Fig. 2(b) is a cloud of flow field temperatures.
FIG. 2(c) is a cloud graph of the flow field numerical viscosity control function.
Fig. 3(a) cylinder surface pressure distribution.
FIG. 3(b) cylindrical surface heat flow distribution.
FIG. 4 is a schematic diagram of interface migration.
Detailed Description
The invention is described below with reference to application examples.
1. Technical parameters
Wieting. a. R (Allen R W1987 NASA TM-100484) performed an aerodynamic heating test of the leading edge of stainless steel round tubes in an 8-foot high enthalpy wind tunnel at NASA lanley research center, usa in 1987 (Wieting a R, Holden M S1987 AIAA 22nd Thermophysics Conference honolu Hawaii June 8-10,1987), which has been used many times to verify the accuracy of the heat transfer coupling calculations for flow field structures. The invention also selects infinite-length stainless steel round tube models with the same conditions, and the size of the round tube is the outer diameter RoutThe incoming flow condition calculation parameters are given in table 1, 0.0381 m.
Table 1 incoming flow condition parameters
2. Model mesh
For the present example, choose320 x 320 non-uniform structure grid, and the height of the first layer grid on the near wall surface is determined according to the Reynolds number of the local grid, the Reynolds number of the incoming flow of the embodiment is 1.835 x 105In order to calculate the Reynolds number Re of the first layer of grid in the wall surface normal directioncell1.835, meets the requirement of accurate capture of thermal boundary layer in the literature on the reynolds number of the grid. FIG. 1 shows Ma∞8.03 hypersonic cylindrical streaming computational grid schematic.
3. Numerical simulation results
The pressure, temperature and switching control function cloud images of the flow field of the hypersonic cylindrical flow, which are calculated by adopting the improved LBFS method of the invention and have the Ma of 8.03, are respectively shown in the figure. From fig. 2(a) and 2(b) we can see that the improved LBFS method can stably and clearly capture the strong shock wave region in front of the cylinder without emitting any significant numerical oscillations that are present. In addition, as can be seen from fig. 2(c), the improved switching control function of the present invention can play a good role in controlling the numerical viscosity in hypersonic flow, and the numerical viscosity weight coefficient τ in smooth flow field0The value is close to zero, and tau around the strong shock wave region0The value is close to 1, therefore, τ0The distribution of values is consistent with the expected theory, and the design target of the control function is achieved.
The calculated pressure (upper) and heat flow (lower) distributions (normalized by stagnation values) along the surface of the cylinder are shown separately. In addition, experimental data distribution, calculation results in AUSM + format and calculation results in pneumatic BGK format of Xu et al are respectively shown in the figure, and it can be observed that the pressure and heat flow distribution predicted by the improved LBFS method of the present invention is well matched with the reference values. In addition, table 2 shows the pressure and heat flux at the stagnation point calculated using the LBFS method compared to the results calculated in references xu. As can be seen from the table, the results of the modified LBFS calculation fit well with the reference values.
The data comparison shows that for numerical simulation of hypersonic flow problems, the improved LBFS method can simultaneously and stably capture strong shock waves and accurately calculate the flow information of the thin thermal boundary layer.
TABLE 1 comparison of stagnation pressure to Heat flow value
Claims (7)
1. The flow field calculation method based on the improved FVM-LBFS method is characterized in that the flow field calculation method combines a finite volume method with an LBM method, a macroscopic N-S control equation is subjected to spatial dispersion by adopting the finite volume method, the inviscid flux of a unit interface is reconstructed by utilizing a local solution of a one-dimensional compressible lattice Boltzmann model, the improved FVM-LBFS method realizes accurate control of inviscid flux numerical value viscosity by adopting an improved switch control function for controlling pressure and temperature simultaneously, and a grid unit slenderness ratio correction coefficient is introduced to expand the application range of the non-uniform grid.
2. The improved FVM-LBFS method of claim 1, wherein an average Reynolds N-S (RANS) equation is used as a control equation, and an integral form of the equation is as follows:
in the formula: w is a conservative value; fcNo adhesive flux; fvIs the viscous flux; w, Fc、FvIs defined as:
where ρ, p and T represent density, pressure, temperature, respectively, and are fullThe equation of state p of ideal gas is rho RT, and R is an ideal gas constant; u, v, w are velocities of the fluid in the body in three directions respectively, nx、ny、nzRespectively, the outer normal unit vector component of the control body unit surface dS, where V is nxu+nyv+nzw is the surface normal velocity, τijIs a friction stress tensor component, k is a heat conduction coefficient, E is p/rho (gamma-1) is the energy of the fluid in the control body, and gamma is a specific heat ratio; h is total enthalpy, and H is E + p/rho;
solving equation (1) consists in solving the inviscid flux FcAnd a viscous flux FvCalculating (1); computing inviscid flux F by adopting LBFS method based on non-free parameter D1Q4 modelcThe central difference method is adopted to calculate the viscous flux Fv。
3. The improved FVM-LBFS method based flow field calculation method of claim 2, wherein the derivation process of D1Q4 model is:
to recover the one-dimensional Euler equation by using the Maxwell' S distribution function to obtain the relationship of the moments required to recover the Euler equation or the N-S equation, the conservation of momentum relationship can be simplified to the form of a summation,
in the formula, fieqAs a function of the equilibrium distribution of the particles in the i-th direction, ξiThe particle velocities in the respective directions, p and u being the density and velocity, respectively, of the one-dimensional macroscopic flow, where
Represents a specific velocity of the moving particle; when the lattice Boltzmann model is deduced, firstly, the particle velocity distribution form of an LB model needs to be determined; it can be seen that the compressible D1Q4 model consists of two particle velocities (D)1，d2) And four equilibrium distribution functions
Composition, so there are 6 unknown parameters to be solved in total; since equation (1) contains only 4 independent equations, to solve these 6 unknowns, two more independent equations need to be introduced as follows:
equations (3) and (4) are two higher-order momentum conservation relations, and the two additional higher-order moment relation equations can further constrain the compressible D1Q4 model, so that the physical flow property can be reflected more truly, and particularly, the calculation of the strong shock wave problem in hypersonic flow has better positivity;
thus, simultaneous equations (2) - (4), the following form of the system of equations can be obtained,
solving equation set (5) can obtain 4 equilibrium distribution functions and 2 discrete particle velocity expressions in the compressible D1Q4 model as
4. The improved FVM-LBFS method of calculating a flow field according to claim 3, wherein the FVM-LBFS establishes a connection between an N-S equation and a non-free parameter D1Q4 model, and a conservation variable with respect to a normal direction velocity is obtained according to a physical conservation relation (2) restoring the Euler equation
And no adhesive flux
Can be calculated by the following formula,
in the formula, xiiThe discrete particle phase velocity in the ith direction is shown as xi in the one-dimensional non-free parameter D1Q4 model1＝d1，ξ2＝-d1，ξ3＝d2，ξ4＝-d2；VnIs the normal velocity of the interface, e is the internal energy per unit mass,
then represents a moment relation vector, whose expression is
In addition, epIs the potential energy of particles, fi(r, t) represents the particle distribution function of the ith direction at the spatial position r at the time t, and the whole conservation variable W and the non-viscous flux F can be respectively obtained by combining the above equationscConservation variable with respect to normal direction velocity
And no adhesive flux
The relationship between them is as follows:
wherein the content of the first and second substances,
are respectively the corresponding items 1, 2, 3, u in the formula (8)τIs tangential direction velocity, uτx、uτy、uτzAre each uτComponent in three directions,
Items 1, 2 and 3 in (9), respectively.
5. The flow field calculation method based on the improved FVM-LBFS method of claim 4, wherein in the derivation and solution process of the non-viscous flux, based on Chapman-Enskog expansion analysis, the balance part and the non-balance part of the particle distribution function on the unit interface are simultaneously reserved in the LBFS method, the LBFS method is only used for calculating the non-viscous flux of the unit interface, the contribution of the non-balance part of the distribution function on the unit interface to the flux is regarded as artificial numerical viscosity, the control of the numerical viscosity needs to be accurately controlled according to the physical characteristics of the flow field, and the derivation process of the non-viscous flux calculation is as follows:
in LBFS, the lattice Boltzmann model is used to build a local Riemann flux solver at the cell interface to compute the inviscid flux; a non-viscous flux with respect to normal velocity at the cell interface can be obtained
The expression is as follows:
in the formula, the symbol "+" represents a numerical solution, and the subscript "j + 1/2" represents a cell interface position; the upper corner mark I represents the contribution of the balanced distribution function at the unit interface to the inviscid flux, and the upper corner mark II represents the contribution of the balanced distribution function at the points around the unit interface to the inviscid flux;
the contribution of tangential directional velocity to the inviscid flux is incorporated into equation (12), and the complete expression of inviscid flux at the cell interface is written as follows:
δ t is the particle migration time step, and the contribution of the non-equilibrium portion of the distribution function is considered as the numerical viscosity, τ0Can be used as a weight coefficient of numerical viscosity;
the density, normal momentum and normal velocity at the cell interface contribute to energy production as follows
According to the compatibility condition, the calculation of the conservation variable does not depend on the non-equilibrium part of the particle distribution function; thus, the unbalanced portion of the distribution function in equation (14) can be written as,
it is possible to obtain,
from equation (16), the conservation variable at the cell interface can be represented by the equilibrium distribution function gi(-ξiDelta t, t-delta t) is calculated; thereby obtaining a conservation variable expression on the unit interface,
in the formula (I), the compound is shown in the specification,
respectively representing the balance distribution functions of the left and right of the unit interface;
since the one-dimensional non-free parameter LB model is applied in equation (17), only the conservation variable related to the normal velocity on the cell interface can be obtained; to calculate the contribution to tangential velocity on the cell interface, an approximation method is used, as shown below,
in the formula (I), the compound is shown in the specification,
the tangential velocities on the left side of the unit interface, the right side of the unit interface and the unit interface are respectively; wherein the content of the first and second substances,
and
respectively calculating the difference value between the total speed and the normal speed of the left side and the right side of the reconstructed unit interface; directly calculating to obtain a conservation variable rho on a unit interface through a simultaneous equation*，
And p*；
The calculated conservation variable rho*，
And p*Substituting equations (6) - (7) to calculate the equilibrium distribution function g on the cell interfacei(0, t); function g is distributed with equilibriumi(-ξiδ t, t- δ t) and gi(0, t) is determined and the corresponding inviscid flux at the cell interface can be calculated explicitly.
6. The improved FVM-LBFS method of claim 5, wherein the temperature and pressure values of the control body units introduced to the left and right sides of the unit interface at the same time are set to the weighting coefficient tau0The value of (c) is controlled; a new and improved switching control function is proposed to calculate tau0The value of (b) is defined as follows,
wherein, tanh (x) is hyperbolic tangent function in hyperbolic function, and the function value of tanh (x) can be in [0, 1%]Smooth transition between the two; c is a relaxation factor constant with a value range of [1,100 ]]，pL，pRAnd TL，TRRespectively are pressure value and temperature value at left and right sides of the unit interface, delta is correction coefficient of length-to-thickness ratio of grid unit, omegaL、ΩRThe volumes of the control body units on the left side and the right side of the unit interface are respectively,
then the average volume of the left and right control volume units, LL,RThe distance is the center of the left control unit and the right control unit; n is a radical offLAnd NfRRespectively representing the surface numbers of the grid cells on the left side and the grid cells on the right side of the cell interface,
respectively, the weighting coefficients on the ith plane.
7. The improved FVM-LBFS method-based flow field calculation method of claim 6, wherein the specific solving process is as follows:
(1) calculating the derivative of the conservation variable of each grid unit, and obtaining the conservation variables on two sides of the unit interface through linear reconstruction;
(2) substitution of conservation variables on both sides of cell interfaceEquations (6) - (7), the equilibrium distribution function is calculated
And
(3) the conservation variable at the cell interface is calculated by equations (17) - (18), and further, the original variable ρ is calculated by equation (17)*，
And p*；
(4) Calculating a numerical viscosity weight coefficient τ using equation set (19)0Then the total non-viscous flux at the cell interface is calculated by equation (13)
(5) Calculating F of viscous flux at cell interfacev；
(6) Solving an N-S equation by a time advancing method, wherein a new conservation variable value of the center of the grid unit at the next time step can be calculated in the step;
(7) and (4) repeating the steps (1) - (7) until a convergence solution meeting the condition is obtained.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111078887.2A CN113792432A (en) | 2021-09-15 | 2021-09-15 | Flow field calculation method based on improved FVM-LBFS method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111078887.2A CN113792432A (en) | 2021-09-15 | 2021-09-15 | Flow field calculation method based on improved FVM-LBFS method |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN113792432A true CN113792432A (en) | 2021-12-14 |
Family
ID=79183429
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111078887.2A Pending CN113792432A (en) | 2021-09-15 | 2021-09-15 | Flow field calculation method based on improved FVM-LBFS method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113792432A (en) |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114218833A (en) * | 2021-12-16 | 2022-03-22 | 西北工业大学太仓长三角研究院 | Method and system for predicting performance of internal flow field of secondary light gas gun |
| CN115238397A (en) * | 2022-09-15 | 2022-10-25 | 中国人民解放军国防科技大学 | Method and device for calculating thermal environment of hypersonic aircraft and computer equipment |
| CN116070554A (en) * | 2023-04-06 | 2023-05-05 | 中国人民解放军国防科技大学 | Hypersonic aircraft aerodynamic heat load calculation method, device and equipment |
| CN116882322A (en) * | 2023-09-06 | 2023-10-13 | 中国空气动力研究与发展中心计算空气动力研究所 | Calculation method and device for non-sticky flux, terminal equipment and storage medium |
| CN117393062A (en) * | 2023-12-13 | 2024-01-12 | 上海交通大学四川研究院 | Simulation method for rigid chemical reaction flow rollback self-adaptive semi-hidden semi-explicit coupling time |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20080177511A1 (en) * | 2007-01-04 | 2008-07-24 | Honda Motor Co., Ltd. | Method and system for simulating flow of fluid around a body |
| US20100023276A1 (en) * | 2008-07-22 | 2010-01-28 | General Electric Company | System and method for assessing fluid dynamics |
| CN104573232A (en) * | 2015-01-06 | 2015-04-29 | 浙江理工大学 | Method for determining offset of distribution blade inlet based on energy gradient theory |
| CN107729691A (en) * | 2017-11-13 | 2018-02-23 | 西北工业大学 | A kind of gas flow characteristic method for numerical simulation of thin continuum one |
| US20190258764A1 (en) * | 2018-02-20 | 2019-08-22 | Dassault Systemes Simulia Corp. | Lattice Boltzmann Based Solver for High Speed Flows |
| US20200394277A1 (en) * | 2019-06-11 | 2020-12-17 | Dassault Systemes Simulia Corp. | Computer simulation of physical fluids on irregular spatial grids stabilized for explicit numerical diffusion problems |
-
2021
- 2021-09-15 CN CN202111078887.2A patent/CN113792432A/en active Pending
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20080177511A1 (en) * | 2007-01-04 | 2008-07-24 | Honda Motor Co., Ltd. | Method and system for simulating flow of fluid around a body |
| US20100023276A1 (en) * | 2008-07-22 | 2010-01-28 | General Electric Company | System and method for assessing fluid dynamics |
| CN104573232A (en) * | 2015-01-06 | 2015-04-29 | 浙江理工大学 | Method for determining offset of distribution blade inlet based on energy gradient theory |
| CN107729691A (en) * | 2017-11-13 | 2018-02-23 | 西北工业大学 | A kind of gas flow characteristic method for numerical simulation of thin continuum one |
| US20190258764A1 (en) * | 2018-02-20 | 2019-08-22 | Dassault Systemes Simulia Corp. | Lattice Boltzmann Based Solver for High Speed Flows |
| US20200394277A1 (en) * | 2019-06-11 | 2020-12-17 | Dassault Systemes Simulia Corp. | Computer simulation of physical fluids on irregular spatial grids stabilized for explicit numerical diffusion problems |
Non-Patent Citations (1)
| Title |
|---|
| JIAWEI LI 等: ""A hybrid lattice Boltzmann flux solver for integrated hypersonic fluid-thermal-structural analysis", CHINESE JOURNAL OF AERONAUTICS, vol. 32, no. 9, pages 2295 - 2312 * |
Cited By (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114218833A (en) * | 2021-12-16 | 2022-03-22 | 西北工业大学太仓长三角研究院 | Method and system for predicting performance of internal flow field of secondary light gas gun |
| CN114218833B (en) * | 2021-12-16 | 2023-11-10 | 西北工业大学太仓长三角研究院 | Method and system for predicting performance of flow field in secondary light gas gun |
| CN115238397A (en) * | 2022-09-15 | 2022-10-25 | 中国人民解放军国防科技大学 | Method and device for calculating thermal environment of hypersonic aircraft and computer equipment |
| CN116070554A (en) * | 2023-04-06 | 2023-05-05 | 中国人民解放军国防科技大学 | Hypersonic aircraft aerodynamic heat load calculation method, device and equipment |
| CN116882322A (en) * | 2023-09-06 | 2023-10-13 | 中国空气动力研究与发展中心计算空气动力研究所 | Calculation method and device for non-sticky flux, terminal equipment and storage medium |
| CN116882322B (en) * | 2023-09-06 | 2024-02-13 | 中国空气动力研究与发展中心计算空气动力研究所 | Calculation method and device for non-sticky flux, terminal equipment and storage medium |
| 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 |
|---|---|---|
| CN113792432A (en) | Flow field calculation method based on improved FVM-LBFS method | |
| Afshari et al. | On numerical methods; optimization of CFD solution to evaluate fluid flow around a sample object at low Re numbers | |
| Cao et al. | Implicit high-order gas kinetic scheme for turbulence simulation | |
| Madaliev et al. | Comparison of finite-difference schemes for the first order wave equation problem | |
| CN112069689B (en) | Simulation method and system for fuel atomization characteristic of aircraft engine | |
| CN104461677A (en) | Virtual thermal test method based on CFD and FEM | |
| Radenac | Validation of a 3D ice accretion tool on swept wings of the SUNSET2 program | |
| CN105808954A (en) | Periodic unsteady flow field prediction method suitable for CFD numerical simulation | |
| Benim et al. | Modelling turbulent flow past a circular cylinder by RANS, URANS, LES and DES | |
| CN115587551B (en) | Multi-scale prediction method for ablation behavior of heat-proof structure of hypersonic aircraft | |
| CN115544675B (en) | Multi-scale prediction method for surface catalytic properties of heat-proof material of hypersonic aircraft | |
| CN114168796A (en) | Method for establishing high-altitude aerodynamic database of aircraft | |
| Wei | The development and application of CFD technology in mechanical engineering | |
| Chen et al. | A rotated lattice Boltzmann flux solver with improved stability for the simulation of compressible flows with intense shock waves at high Mach number | |
| Niu | Advection upwinding splitting method to solve a compressible two‐fluid model | |
| Yang et al. | Gas kinetic flux solver based high-order finite-volume method for simulation of two-dimensional compressible flows | |
| Tongqing et al. | CFD/CSD-based flutter prediction method for experimental models in a transonic wind tunnel with porous wall | |
| Lee et al. | Convergence characteristics of upwind method for modified artificial compressibility method | |
| CN114611423A (en) | Three-dimensional multiphase compressible fluid-solid coupling rapid calculation method | |
| Athavale et al. | Application of an unstructured grid solution methodology to turbomachinery flows | |
| Ma et al. | High-order flux reconstruction thermal lattice Boltzmann flux solver for simulation of incompressible thermal flows | |
| CN114021499B (en) | Aircraft heat protection structure heat conduction calculation method based on FVM-TLBFS method | |
| Gutović et al. | CFD Analysis of Pressure Losses in Flat-Oval Duct Fittings. | |
| Taddei et al. | CFD-based analysis of multistage throughflow surfaces with incidence | |
| Gleize et al. | RANS simulations on TMR test cases and M6 wing with the Onera elsA flow solver |
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 |