CN113792432B - 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
- CN113792432B CN113792432B CN202111078887.2A CN202111078887A CN113792432B CN 113792432 B CN113792432 B CN 113792432B CN 202111078887 A CN202111078887 A CN 202111078887A CN 113792432 B CN113792432 B CN 113792432B
- Authority
- CN
- China
- Prior art keywords
- flux
- equation
- interface
- velocity
- lbfs
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 126
- 238000004364 calculation method Methods 0.000 title claims abstract description 47
- 230000004907 flux Effects 0.000 claims abstract description 64
- 230000008569 process Effects 0.000 claims abstract description 24
- 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 5
- 238000005315 distribution function Methods 0.000 claims description 43
- 239000002245 particle Substances 0.000 claims description 28
- 230000014509 gene expression Effects 0.000 claims description 10
- 238000009795 derivation Methods 0.000 claims description 7
- 239000012530 fluid Substances 0.000 claims description 7
- 238000011161 development Methods 0.000 claims description 5
- 238000013508 migration Methods 0.000 claims description 5
- 230000005012 migration Effects 0.000 claims description 5
- 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
- 238000004088 simulation Methods 0.000 abstract description 16
- 239000006185 dispersion Substances 0.000 abstract description 4
- 238000012795 verification Methods 0.000 abstract description 2
- 238000011160 research Methods 0.000 description 6
- 230000008901 benefit Effects 0.000 description 4
- 238000010586 diagram Methods 0.000 description 2
- 238000010438 heat treatment Methods 0.000 description 2
- 229910001220 stainless steel Inorganic materials 0.000 description 2
- 239000010935 stainless steel Substances 0.000 description 2
- 238000012546 transfer Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000006243 chemical reaction Methods 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
- 238000013461 design Methods 0.000 description 1
- 230000009977 dual effect Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000000329 molecular dynamics simulation Methods 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
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 aerodynamic calculation of aircrafts. The method combines the traditional finite volume method with the LBM method, adopts the finite volume method to carry out space dispersion on a macroscopic N-S control equation, and utilizes the non-sticky flux of a local deconstructed unit interface of a one-dimensional compressible lattice Boltzmann model, so that the calculation process is simple and efficient. The improved FVM-LBFS method realizes the accurate control of the non-viscous flux numerical viscosity in the existing LBFS method by adopting the newly proposed improved switch control function for simultaneously controlling the pressure and the temperature, and introduces the grid cell slenderness ratio correction coefficient to expand the application range of the non-uniform grid. Through a series of numerical simulation example verification and analysis, the improved FVM-LBFS method can simultaneously stably capture complex strong shock wave flow and accurately predict aerodynamic heat parameters in hypersonic flow numerical simulation.
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 aerodynamic calculation of aircrafts.
Background
In the later 80 s of the 20 th century, along with the rapid development of computer technology, the running speed of a computer is faster and faster, and the development of a numerical simulation method of fluid mechanics starts to rise gradually from a mesoscopic level. The lattice Boltzman method (Lattice Boltzmann Method, LBM), one of which is typically represented, has received extensive attention from researchers since 1988, and has evolved rapidly over a period of nearly 20 years. Unlike the N-S equation based on continuity hypothesis, the lattice Boltzman method is a mesoscale fluid computational mechanics method between the fluid micromolecular dynamics computation method and the macroscopic continuity equation computation method, so compared with the traditional CFD method, the lattice Boltzman method is not limited by the continuous medium hypothesis, and can realize the trans-scale flow simulation.
The lattice Boltzman equation is a discrete representation of the Boltzman equation in velocity space, physical space and time, which describes the process of moving a particle 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 a nonlinear partial differential equation set (such as an Euler equation set, an N-S equation set and the like) by a traditional method, and is completed once when a model is established, so that only a lattice Boltzmann equation simple linear equation or equation set needs 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 a simple linear migration process and a collision process, only algebraic operation is needed, and a pressure Poisson equation is not needed to be solved, so that the program is very simple to realize. Given the many unique advantages of the lattice Boltzmann method, it has evolved as 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 medium flow, chemical reaction flow, thin gas flow, etc., numerical methods based on lattice Boltzman model are well applied and the physical mechanism of many complex flow problems is also 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, and a discrete speed model and a corresponding balanced distribution function form in the lattice Boltzmann model are key, so that a great deal of research work is performed by researchers at home and abroad at present. The traditional lattice Boltzmann model is mostly built for incompressible flows, such as the DdQm series of basic models proposed by Qian et al in 1992 (d represents the spatial dimension and m represents 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 enables numerical modeling 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, with truncation errors directly related to mach numbers, resulting in inapplicability to high mach number flows. To overcome this problem, qu et al have used a round function instead of maxwell's distribution function to build a series of one-dimensional and two-dimensional models, such as D1Q5L2 model, D2Q13L2 model, and successfully applied these models to numerical simulations of strong shock supersonic flow fields such as shock tubes, dual mach number reflections, and the like. Dellar also derive two-dimensional compressible lattice Boltzmann models, such as a seven-speed model without splitting, a 4+3 split speed model. However, these models contain many free parameters, and the lattice velocity needs to be artificially given, which results in poor versatility of application in complex high mach number flow numerical simulations. In view of this, shu and its research team derive theoretically, a non-free parameter lattice Boltzmann model derivation platform based on moment conservation is provided, 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. The macroscopic equation can be recovered locally and accurately by using the models, and a one-dimensional compressible non-viscous flow Riemann flux solver is established at a local unit interface.
Due to the requirement of grid uniformity, the traditional LBM method has very high requirement on the uniformity of the calculation grid, and is not suitable for calculating the complex appearance. Moreover, for high Reynolds number flows, the computational grid can be very dense due to the limitations of migration time steps, resulting in a substantial reduction in computational efficiency. Thus, when using the lattice Boltzmann method for numerical modeling of compressible flows, most use a discrete velocity Boltzmann equation (Discrete Velocity Boltzmann Equation, DVBE) based method calculation to update the particle distribution function in time. Many scholars, qu, q.li, ubertini, stiebler, etc., have made a great deal of research work on the DVBE method. However, it should be noted that in the DVBE method, the number of particle distribution functions is far greater than the number of conservation variables, and the distribution function form is very complex, so that the direct solving process of the DVBE method is complicated, and the calculation efficiency is far worse than that of the traditional 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 of an integral form, the physical concept of the deduction process is clear, and the conservation characteristic of the physical quantity is maintained in each control body; and has good flexibility for calculation grids (structural/non-structural grids), and is suitable for calculation of complex geometric shapes. Therefore, in the LBFS method, the complex process of directly solving DVBE equations is avoided, the advantage of the traditional FVM is utilized to discrete a macroscopic control equation, and the non-sticky flux of a unit interface is reconstructed by adopting the local solution of a one-dimensional compressible lattice Boltzmann model. Numerical calculation results show that the LBFS method has good positive value and high efficiency in strong shock wave compressible flow simulation. In the existing LBFS method, in order to capture strong discontinuous physical flow such as strong shock waves, the unbalanced part of a unit interface distribution function is utilized to introduce numerical dissipation so as to ensure the stability of calculation. However, the value dissipation directly introduced in this way is too large to seriously disturb the true solution of smooth areas such as boundary layers, and is unacceptable in particular for accurately predicting the heat flux density of hypersonic flow fields. Therefore, how to ensure that the LBFS method can simultaneously stably capture the strong shock wave flow and accurately predict the aerodynamic heat in hypersonic flow is also required to be further perfected, which is also the main content of the research of the 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 deduction 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 simple and efficient calculation process. The improved FVM-LBFS method realizes the accurate control of the non-viscous flux numerical viscosity in the existing LBFS method by adopting the newly proposed improved switch control function for simultaneously controlling the pressure and the temperature, and introduces the grid cell slenderness ratio correction coefficient to expand the application range of the non-uniform grid. Finally, a series of high-speed numerical examples are selected to verify and analyze the reliability and accuracy of the proposed 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 existing numerical method based on the lattice Boltzmann model is difficult to simulate high Reynolds number flow, an improved FVM-LBFS method is provided for numerical simulation of hypersonic flow. The method combines the traditional finite volume method with the LBM method, adopts the finite volume method to carry out space dispersion on a macroscopic N-S control equation, and utilizes the non-sticky flux of a local deconstructed unit interface of a one-dimensional compressible lattice Boltzmann model, so that the calculation process is simple and efficient. The improved FVM-LBFS method realizes the accurate control of the non-viscous flux numerical viscosity in the existing LBFS method by adopting the newly proposed improved switch control function for simultaneously controlling the pressure and the temperature, and introduces the grid cell slenderness ratio correction coefficient to expand the application range of the non-uniform grid. Through a series of numerical simulation example verification and analysis, the improved FVM-LBFS method can simultaneously and stably capture complex strong shock wave flow and accurately predict aerodynamic heat parameters in hypersonic flow numerical simulation, and well expands the application range of the numerical method based on the lattice Boltzmann model in the hypersonic flow numerical calculation field.
The technical scheme of the invention is as follows:
The flow field calculation method based on the improved FVM-LBFS method comprises the following steps:
the average reynolds N-S (RANS) equation is used as the control equation, whose integral form (passive term) is as follows,
Wherein: w is a constant; f c is non-stick flux; f v is viscous flux; w, F c、Fv is defined as:
Where ρ, p and T represent density, pressure, temperature, respectively, and satisfying the ideal gas state equation p=ρrt, R is the ideal gas constant. u, V, w are the velocities of the fluid in the control body in three directions, n x、ny、nz is the external normal unit vector component of the control body unit surface dS, v=n xu+nyv+nz w is the surface normal velocity, τ ij is the friction stress tensor component, k is the heat transfer coefficient, e=p/ρ (γ -1) is the energy of the fluid in the control body, H is the total enthalpy, γ is the specific heat ratio; h=e+p/p.
From the above, it can be seen from equation (1) that the key step in solving it is the calculation of the non-tacky flux F c and the tacky flux F v. The invention calculates the non-viscous flux F c by using an LBFS method based on a non-free parameter D1Q4 model, and calculates the viscous flux F v by using a center difference method in consideration of the physical properties of the viscous flux. The calculation of the LBFS method to solve for the non-stick flux F c will be described and deduced in detail below.
The LBFS method is built 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 uses a D1Q4 model with non-free parameters for the LBFS method of the invention, and the following briefly describes the derivation process of the D1Q4 model:
In order to obtain the moment relation of each order needed for restoring the Euler equation or the N-S equation by using the Maxwell distribution function, the momentum conservation relation can be simplified into a summation form,
In the method, in the process of the invention,For the particle equilibrium distribution function in the ith direction, ζ i is the particle velocity in the corresponding direction, ρ and u are the density and velocity of one-dimensional macroscopic flow, respectively, where/>Indicating the particular velocity of the moving particles. In the derivation of the lattice Boltzmann model, it is first necessary to determine the particle velocity distribution pattern of the LB model. It can be seen that the compressible D1Q4 model consists of two particle velocities (D 1,d2) and four equilibrium distribution functions/> The composition, therefore, there are a total of 6 unknown parameters to be solved. Since equation (1) contains only 4 independent equations, to solve these 6 unknowns, the following two independent equations need to be introduced again
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 properties can be reflected more truly, and the calculation of the strong shock problem in hypersonic flow has better positive value.
Thus, simultaneous equations (2) - (4) can be obtained in the form of a set of equations,
Solving equation set (5) can obtain 4 equilibrium distribution functions and 2 discrete particle velocity expressions in the compressible D1Q4 model as
FVM-LBFS constructs a link between the N-S equation and the non-free parameter D1Q4 model, based on the physical conservation relation (2) of recovering Euler' S equation, conservation variable for normal direction velocityAnd the tack-free flux F c can be calculated by the following formula,
In the formula, xi i represents the discrete particle phase velocity in the i-th direction, ξ1＝d1,ξ2＝-d1,ξ3＝d2,ξ4＝-d2.Vn is the interface normal velocity in a one-dimensional non-free parameter D1Q4 model, and e is the unit mass internal energy; then the moment relation vector is represented, the expression is
In addition, e p is the potential energy of particles, F i (r, t) represents the particle distribution function of the ith direction at the space position r at the time t, and the whole conservation variable W, the non-sticky flux F c and the conservation variable about the normal direction speed can be respectively obtained by combining the equationsAnd non-sticking flux/>The connection between them is as follows
Wherein,The corresponding 1 st, 2 nd and 3 rd items in the formula (8) are respectively shown, u τ is tangential direction speed, u τx、uτy、uτz is the component of u τ in three directions respectively,/>The corresponding items 1, 2 and 3 in the item (9) respectively;
In the improved FVM-LBFS method, the core part is derivation and solution of non-sticky flux, the balanced part and the unbalanced part of a particle distribution function on a unit interface are reserved in the LBFS method at the same time, the LBFS method is only used for calculating the non-sticky 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 viscosity, the control of the numerical viscosity needs to be accurately controlled according to the physical characteristics of a flow field, and the derivation process of the non-sticky flux calculation is described below:
In the LBFS of the present invention, a lattice Boltzmann model is used to build a local Riemann flux solver at the cell interface to calculate the non-stick flux. Can obtain non-sticking flux at the unit interface with respect to normal velocity The expression is as follows:
In the formula, the symbol "×" represents a numerical solution, the subscript "j+1/2" represents a unit interface position, the upper corner mark I represents a contribution of a balanced distribution function at the unit interface to the non-stick flux, and the upper corner mark II represents a contribution of a balanced distribution function at points around the unit interface to the non-stick flux. Further, incorporating the contribution of tangential direction velocity to non-stick flux into equation (12), the complete expression of non-stick flux at the cell interface can be written as follows
Δt is a particle migration time step, and the contribution of the unbalanced part of the distribution function is regarded as numerical viscosity, and τ 0 can be used as a weight coefficient of the numerical viscosity;
The contributions to energy production from density, normal momentum and normal velocity at the cell interface are as follows
From the compatibility conditions, the calculation of the conservation variable is independent of the unbalanced 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 a solution that,
From equation (16), it can be seen that the conservation variable at the cell interface can be calculated by the equilibrium distribution function g i(-ξi δt, t—δt. Further, we can derive the conservation variable expression on the cell interface,
In the method, in the process of the invention,Respectively representing the balanced 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 at the cell interface can be obtained. In order to calculate the contribution of tangential velocity at the cell interface, the present invention employs a viable approximation method, as follows,
In the method, in the process of the invention,Tangential velocities on the left side of the cell interface, the right side of the cell interface, and the cell interface, respectively. Wherein/>And/>The difference between the total velocity and the normal velocity at the left and right sides of the reconstructed unit interface can be calculated. Through simultaneous equations, we can directly calculate the conservation variable ρ *,/>, on the cell interfaceAnd p *. Further, the calculated conservation variable ρ *,/>And p * are substituted into equations (6) - (7) to calculate the equilibrium distribution function g i (0, t) at the cell interface. With the equilibrium distribution functions g i(-ξi δt, t- δt and g i (0, t) determined, the corresponding non-stick flux at the cell interface can be explicitly calculated.
In order to accurately control the magnitude of the numerical sticky generated by the unbalanced portion of the unit interface distribution function within the overall computational domain, a global switching function is constructed to calculate the magnitude of the numerical sticky weight coefficient τ 0.
Since the accuracy of the calculation of heat flow in hypersonic flow depends on the exact capture of shock and temperature boundary layers, which in turn are very sensitive to the numerical viscosity in the calculation. According to the molecular dynamics theory, the dimensionless relaxation time tau 0 is related to pressure and temperature, so that in order to more accurately control the value of the numerical viscosity weight coefficient tau 0 in a calculation domain, the invention controls the value of the weight coefficient tau 0 by simultaneously introducing the temperature value and the pressure value of a control body unit at the left side and the right side of a unit interface. Based on the above theoretical analysis, a new and improved switch control function is proposed to calculate the value of τ 0, defined as follows,
Where tan h (x) is the hyperbolic tangent function in the hyperbolic function, the tan h (x) function value may be smoothly excessive between [0,1 ]. p L and p R are respectively the pressure values of the left and right sides of the cell interface, C is a relaxation factor constant, the range of the values is [1,100], p L,pR and T L,TR are respectively the pressure values and the temperature values of the left and right sides of the cell interface, delta is a grid cell slenderness ratio correction coefficient, omega L、ΩR is the volume of a control body unit at the left and right sides of the cell interface,The average volume of the left and right control units is L L,R, the distance between the centers of the left and right control units. N fL and N fR represent the number of faces of the left and right grid cells, respectively, of the cell interface,/>The weight coefficients on the i-th face respectively.
The specific solving process is as follows:
(1) Calculating the derivative of the conservation variable of each grid unit, and obtaining the conservation variable at two sides of the unit interface through linear reconstruction;
(2) Substituting the conservation variable at two sides of the unit interface into equations (6) - (7) to calculate the equilibrium distribution function And
(3) The conservation variable on 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 τ 0 by using equation (19), and then calculating the total non-viscous flux at the cell interface by using equation (13)
(5) F v calculating the viscous flux at the cell interface;
(6) Solving an N-S equation by a time advancing method, wherein the new conservation variable value of the center of the grid unit of the next time step can be calculated;
(7) 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 aerodynamic heat calculation
A non-free parameter one-dimensional D1Q4 model is used. The distribution function of the model has simple structure and high calculation efficiency, and the particle speed is determined by a higher-order moment-of-momentum relation without manually giving free parameters, so that the model is more reasonable.
2. The invention discloses a novel numerical value viscosity control function based on grid effect
The invention controls the value of the weight coefficient tau 0 by simultaneously introducing the temperature value and the pressure value of the control body units at the left side and the right side of the unit interface. Based on the above theoretical analysis, a new and improved switch control function is proposed to calculate the value of τ 0,
3. The invention discloses a new flow field calculation method based on an improved FVM-LBFS method
The invention provides an improved FVM-LBFS method for numerical simulation of hypersonic flow in order to expand the application range of a numerical method based on a lattice Boltzmann model in hypersonic flow fields. 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 non-viscous flux, and adopts a second-order center difference method to calculate viscous flux. The precise control of the non-viscous flux numerical value viscosity in the existing LBFS method is realized through the newly-proposed improved switch control function for simultaneously controlling the pressure and the temperature, and the grid cell slenderness ratio correction coefficient is introduced to expand the application range of the non-uniform grid.
Drawings
FIG. 1 is a computational grid.
Fig. 2 (a) is a flow field pressure cloud.
Fig. 2 (b) is a flow field temperature cloud.
Fig. 2 (c) is a flow field numerical viscosity control function cloud.
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 present invention will be described below with reference to examples of application.
1. Technical parameters
The 8 foot high enthalpy wind tunnel in the NASA orchid research center in the united states in 1987 was completed by wietting. A. R (Allen R W1987 NASA TM-100484) for the stainless steel round tube leading edge aerodynamic heating test (WIETING A R, holden M S1987AIAA 22nd Thermophysics Conference Honolulu Hawaii June 8-10,1987) which has been used many times to verify the accuracy of the flow field structure heat transfer coupling calculations. The invention also selects an infinitely long stainless steel circular tube model with identical conditions, the circular tube size is that of an external diameter R out =0.0381 m, and the calculation parameters of the incoming flow conditions are given in a table 1.
TABLE 1 incoming flow Condition parameters
2. Model mesh
For the calculation example, 320×320 non-uniform structural grids are selected, the height of the first grid layer near the wall surface is determined according to the local grid Reynolds number, the incoming flow Reynolds number of the calculation example is 1.835 ×10 5, the Reynolds number Re cell = 1.835 of the first grid layer near the wall surface is selected conveniently for calculation, and the requirement of accurate capture of the Reynolds number of the grid by a thermal boundary layer in the literature is met. Fig. 1 is a schematic diagram of a hypersonic cylinder bypass flow calculation grid with Ma ∞ =8.03.
3. Numerical simulation results
The flow field pressure, temperature and switch control function cloud charts of Ma=8.03 hypersonic cylindrical flow calculated by the improved LBFS method are respectively shown in the figures. From fig. 2 (a) and 2 (b), we can see that the improved LBFS method can stably and clearly capture the strong shock region in front of the cylinder without any significant numerical oscillation. In addition, as can be seen from fig. 2 (c), the improved switch control function of the present invention can play a good role in controlling numerical viscosity in hypersonic flow, the numerical viscosity weight coefficient τ 0 is close to zero in a smooth river basin, and the value τ 0 is close to 1 around a strong shock wave region, so that the distribution of the values τ 0 is consistent with the expected theory, and the design goal of the control function is achieved.
The calculated pressure (up) and heat flow (down) profiles along the cylindrical surface (normalized by the dwell point values) of the present invention are shown. In addition, experimental data distribution, calculation results adopting an AUSM+ format and calculation results adopting an aerodynamic BGK format by Xu et al are also respectively provided in the graph, and the pressure and heat flow distribution predicted by the improved LBFS method of the invention can be observed to be better matched with the reference values. In addition, table 2 shows a comparison of the pressure and heat flux at the stagnation point calculated using the LBFS method with the results calculated in references xu.et al and Yang.et al. It can be seen from the table that the results of the improved LBFS calculation agree 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 flow information of a thin thermal boundary layer.
Table 1 comparison of stagnation pressure and Heat flow values
Claims (1)
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 spatially discretized by adopting the finite volume method, the non-sticky flux of a local solution reconstruction unit interface of a one-dimensional compressible grid Boltzmann model is utilized, the improved FVM-LBFS method realizes accurate control of the numerical viscosity of the non-sticky flux by adopting an improved switch control function for simultaneously controlling pressure and temperature, and a grid unit slenderness ratio correction coefficient is introduced to expand the application range of a non-uniform grid;
the average Reynolds N-S (RANS) equation is used as the control equation, and the integral form is as follows:
Wherein: w is a constant; f c is non-stick flux; f v is viscous flux; w, F c、Fv is defined as:
Wherein ρ, p and T represent density, pressure, temperature, respectively, and satisfy an ideal gas state equation p=ρrt, R being an ideal gas constant; u, V, w are the velocities of the fluid in the control body in three directions, n x、ny、nz is the external normal unit vector component of the control body unit surface dS, v=n xu+nyv+nz w is the surface normal velocity, τ ij is the friction stress tensor component, k is the heat conduction coefficient, e=p/ρ (γ -1) is the energy of the fluid in the control body, γ is the specific heat ratio; h is total enthalpy, h=e+p/p;
Solving equation (1) consists in the calculation of a non-tacky flux F c and a tacky flux F v; calculating a non-viscous flux F c by adopting an LBFS method based on a non-free parameter D1Q4 model, and calculating a viscous flux F v by adopting a center difference method;
the derivation process of the D1Q4 model is as follows:
In order to obtain the moment relation of each order needed for restoring the Euler equation or the N-S equation by using the Maxwell distribution function, the momentum conservation relation can be simplified into a summation form,
Wherein f i eq is the particle balance distribution function in the ith direction, ζ i is the particle velocity in the corresponding direction, ρ and u are the density and velocity of one-dimensional macroscopic flow, respectively, in the formulaRepresenting a particular velocity of the moving particle; when the derivation of the lattice Boltzmann model is carried out, firstly, the particle velocity distribution form of the 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, therefore there are a total of 6 unknown parameters to be solved; since equation (1) contains only 4 independent equations, to solve these 6 unknowns, the following two independent equations need to be introduced:
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 properties can be reflected more truly, and the calculation of the strong shock problem in hypersonic flow has better positive value;
Thus, simultaneous equations (2) - (4) can be obtained in the form of a set of equations,
Solving equation set (5) can obtain 4 equilibrium distribution functions and 2 discrete particle velocity expressions in the compressible D1Q4 model as
FVM-LBFS constructs a link between the N-S equation and the non-free parameter D1Q4 model, based on the physical conservation relation (2) of recovering Euler' S equation, conservation variable for normal direction velocityAnd non-sticking flux/>It can be calculated by the following formula,
Wherein xi i represents the discrete particle phase velocity in the i-th direction, ξ1＝d1,ξ2＝-d1,ξ3＝d2,ξ4＝-d2;Vn is the interface normal velocity in the one-dimensional non-free parameter D1Q4 model, e is the unit mass internal energy,Then the moment relation vector is represented, the expression is
In addition, e p is the potential energy of particles, F i (r, t) represents the particle distribution function of the ith direction at the space position r at the time t, and the whole conservation variable W, the non-sticky flux F c and the conservation variable about the normal direction speed can be respectively obtained by combining the equationsAnd non-sticking flux/>The links between are as follows:
Wherein, The corresponding 1,2,3 items in the formula (8), u τ is the tangential direction velocity, u τx、uτy、uτz is the components of u τ in three directions, respectively,/>The corresponding items 1, 2 and 3 in the item (9) respectively;
In the development and solving process of the non-sticking flux, starting from the development analysis of Chapman-Enskog, a balanced part and an unbalanced part of a particle distribution function on a unit interface are reserved in an LBFS method, the LBFS method is only used for calculating the non-sticking flux of the unit interface, the contribution of the unbalanced 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 a flow field, and the following development process of the non-sticking flux calculation is carried out:
In LBFS, a lattice Boltzmann model is used to build a local Riemann flux solver at the cell interface to calculate the non-stick flux; can obtain non-sticking flux at the unit interface with respect to normal velocity The expression is as follows:
Wherein, the symbol ". Times." 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 equilibrium distribution function at the unit interface to the non-sticky flux, and the upper corner mark II represents the contribution of the equilibrium distribution function at the points around the unit interface to the non-sticky flux;
The contribution of tangential direction velocity to non-stick flux is incorporated into equation (12), and the complete expression of non-stick flux at the cell interface is written as follows:
δt is a particle migration time step, the contribution of the unbalanced part of the distribution function is regarded as numerical viscosity, and τ 0 can be used as a weight coefficient of the numerical viscosity;
The contributions to energy production from density, normal momentum and normal velocity at the cell interface are as follows
From the compatibility conditions, the calculation of the conservation variable is independent of the unbalanced 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 a solution that,
From equation (16), the conservation variable at the cell interface can be calculated by the equilibrium distribution function g i(-ξi δt, t—δt; thereby obtaining the conservation variable expression on the unit interface,
In the method, in the process of the invention,Respectively representing the balanced 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; in order to calculate the contribution of tangential velocity at the cell interface, an approximation method is used, as follows,
In the method, in the process of the invention,Tangential velocity at the left side of the cell interface, the right side of the cell interface, and the cell interface, respectively; wherein/>And/>Respectively calculating the difference value between the total speed and the normal speed at the left side and the right side of the reconstructed unit interface; directly calculating to obtain conservation variable rho on the unit interface through simultaneous equations And p *;
The calculated conservation variable ρ *, And p * are substituted into equations (6) - (7) to calculate a balanced distribution function g i (0, t) on the cell interface; with the equilibrium distribution functions g i(-ξi δt, t- δt, and g i (0, t) determined, the corresponding non-stick flux at the cell interface can be explicitly calculated;
The temperature value and the pressure value of the control body units simultaneously introduced into the left side and the right side of the unit interface are controlled to the value of the weight coefficient tau 0; a new and improved switch control function is proposed to calculate the value of tau 0, defined as follows,
Wherein, tan h (x) is a hyperbolic tangent function in the hyperbolic function, and the tan h (x) function value can be smoothly transited between [0,1 ]; c is a relaxation factor constant, the value range is [1,100], p L,pR and T L,TR are respectively the pressure value and the temperature value of the left side and the right side of the cell interface, delta is a grid cell slenderness ratio correction coefficient, omega L、ΩR is the volume of a control body cell of the left side and the right side of the cell interface,Then the average volume of the left and right control units, L L,R is the distance between the centers of the left and right control units; n fL and N fR represent the number of faces of the left and right grid cells, respectively, of the cell interface,/>Respectively the weight coefficients on the ith surface;
the specific solving process is as follows:
(1) Calculating the derivative of the conservation variable of each grid unit, and obtaining the conservation variable at two sides of the unit interface through linear reconstruction;
(2) Substituting the conservation variable at two sides of the unit interface into equations (6) - (7) to calculate the equilibrium distribution function And
(3) The conservation variable on 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 τ 0 by using equation (19), and then calculating the total non-viscous flux at the cell interface by using equation (13)
(5) F v calculating the viscous flux at the cell interface;
(6) Solving an N-S equation by a time advancing method, wherein the new conservation variable value of the center of the grid unit of the next time step can be calculated;
(7) 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 CN113792432B (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 CN113792432B (en) | 2021-09-15 | 2021-09-15 | Flow field calculation method based on improved FVM-LBFS method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113792432A CN113792432A (en) | 2021-12-14 |
| CN113792432B true CN113792432B (en) | 2024-06-18 |
Family
ID=79183429
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111078887.2A Active CN113792432B (en) | 2021-09-15 | 2021-09-15 | Flow field calculation method based on improved FVM-LBFS method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113792432B (en) |
Families Citing this family (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114218833B (en) * | 2021-12-16 | 2023-11-10 | 西北工业大学太仓长三角研究院 | Method and system for predicting performance of flow field in secondary light gas gun |
| CN115178723B (en) * | 2022-07-07 | 2024-08-23 | 苏州数算软云科技有限公司 | Calculation method of die air pressure in high-pressure casting process |
| CN115238397B (en) * | 2022-09-15 | 2022-12-02 | 中国人民解放军国防科技大学 | Method and device for calculating thermal environment of hypersonic aircraft and computer equipment |
| CN116070554B (en) * | 2023-04-06 | 2023-06-09 | 中国人民解放军国防科技大学 | Hypersonic aircraft aerodynamic heat load calculation method, device and equipment |
| CN116882322B (en) * | 2023-09-06 | 2024-02-13 | 中国空气动力研究与发展中心计算空气动力研究所 | Calculation method and device for non-sticky flux, terminal equipment and storage medium |
| CN117393062B (en) * | 2023-12-13 | 2024-02-23 | 上海交通大学四川研究院 | Simulation method for rigid chemical reaction flow rollback self-adaptive semi-hidden semi-explicit coupling time |
Family Cites Families (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7921002B2 (en) * | 2007-01-04 | 2011-04-05 | Honda Motor Co., Ltd. | Method and system for simulating flow of fluid around a body |
| US8577626B2 (en) * | 2008-07-22 | 2013-11-05 | General Electric Company | System and method for assessing fluid dynamics |
| CN104573232B (en) * | 2015-01-06 | 2018-02-16 | 浙江理工大学 | Method is determined based on the theoretical splitterr vanes inlet offset degree of energy gradient |
| CN107729691B (en) * | 2017-11-13 | 2020-05-22 | 西北工业大学 | Rarefied continuous unified gas flow characteristic numerical simulation method |
| US12118279B2 (en) * | 2018-02-20 | 2024-10-15 | Dassault Systemes Americas Corp. | Lattice Boltzmann based solver for high speed flows |
| US11645433B2 (en) * | 2019-06-11 | 2023-05-09 | 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/CN113792432B/en active Active
Non-Patent Citations (1)
| Title |
|---|
| "A hybrid lattice Boltzmann flux solver for integrated hypersonic fluid-thermal-structural analysis;jiawei li 等;Chinese Journal of Aeronautics;第32卷(第9期);2295-2312 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113792432A (en) | 2021-12-14 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113792432B (en) | Flow field calculation method based on improved FVM-LBFS method | |
| CN104461677B (en) | A kind of virtual thermal test method based on CFD and FEM technologies | |
| Madaliev et al. | Comparison of finite-difference schemes for the first order wave equation problem | |
| Cao et al. | Implicit high-order gas kinetic scheme for turbulence simulation | |
| Su et al. | Low-speed flow simulation by the gas-kinetic scheme | |
| Yang et al. | A three-dimensional explicit sphere function-based gas-kinetic flux solver for simulation of inviscid compressible flows | |
| CN105808954A (en) | Periodic unsteady flow field prediction method suitable for CFD numerical simulation | |
| CN107423511A (en) | Meet to immerse border implicit iterative solving method without sliding boundary condition and the condition of continuity | |
| Namkoong et al. | Numerical analysis of two-dimensional motion of a freely falling circular cylinder in an infinite fluid | |
| Li et al. | An efficient implementation of aeroelastic tailoring based on efficient computational fluid dynamics-based reduced order model | |
| Rodi et al. | ERCOFTAC workshop on data bases and testing of calculation methods for turbulent flows | |
| CN114611423A (en) | Three-dimensional multiphase compressible fluid-solid coupling rapid calculation method | |
| Sharif et al. | Numerical study of turbulent natural convection in a side-heated square cavityat various angles of inclination | |
| CN114021499B (en) | Aircraft heat protection structure heat conduction calculation method based on FVM-TLBFS method | |
| Yang et al. | Gas kinetic flux solver based high-order finite-volume method for simulation of two-dimensional compressible flows | |
| Patel et al. | A MAC scheme in boundary-fitted curvilinear coordinates | |
| Doolan et al. | Modeling mass entrainment in a quasi-one-dimensional shock tube code | |
| Jia-hui et al. | Dynamic simulation of nozzle structure based on thermal-fluid-solid coupling analysis | |
| Pan et al. | High-order gas-kinetic scheme in curvilinear coordinates for the Euler and Navier-Stokes solutions | |
| Taddei et al. | CFD-based analysis of multistage throughflow surfaces with incidence | |
| Verma et al. | GeoSteam. Net: steam transport simulator for geothermal pipeline network | |
| Regan et al. | Development of a Linear Stirling System Model with Varying Heat Input | |
| Scha¨ fer | Coupled fluid-solid problems: Survey on numerical approaches and applications | |
| Cheng et al. | Review and Assessment of Turbulence Transition Models | |
| Haoui | Finite volume analysis of a supersonic non-equilibrium flow around an axisymmetric blunt body |
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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |