CN114021499B - Aircraft heat protection structure heat conduction calculation method based on FVM-TLBFS method - Google Patents
Aircraft heat protection structure heat conduction calculation method based on FVM-TLBFS method Download PDFInfo
- Publication number
- CN114021499B CN114021499B CN202111303732.4A CN202111303732A CN114021499B CN 114021499 B CN114021499 B CN 114021499B CN 202111303732 A CN202111303732 A CN 202111303732A CN 114021499 B CN114021499 B CN 114021499B
- Authority
- CN
- China
- Prior art keywords
- equation
- distribution function
- heat transfer
- calculation
- model
- 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 101
- 238000004364 calculation method Methods 0.000 title claims abstract description 86
- 238000012546 transfer Methods 0.000 claims abstract description 52
- 230000004907 flux Effects 0.000 claims abstract description 31
- 238000005315 distribution function Methods 0.000 claims description 59
- 238000009826 distribution Methods 0.000 claims description 27
- 239000002245 particle Substances 0.000 claims description 22
- 239000007787 solid Substances 0.000 claims description 17
- 230000008569 process Effects 0.000 claims description 14
- 238000004088 simulation Methods 0.000 claims description 10
- 238000004458 analytical method Methods 0.000 claims description 9
- 238000009795 derivation Methods 0.000 claims description 9
- 238000013508 migration Methods 0.000 claims description 7
- 230000005012 migration Effects 0.000 claims description 7
- 230000014509 gene expression Effects 0.000 claims description 5
- 238000009792 diffusion process Methods 0.000 claims description 4
- 239000000463 material Substances 0.000 claims description 4
- 239000006185 dispersion Substances 0.000 claims description 3
- 238000013213 extrapolation Methods 0.000 claims description 3
- 238000011084 recovery Methods 0.000 claims description 2
- CLOMYZFHNHFSIQ-UHFFFAOYSA-N clonixin Chemical compound CC1=C(Cl)C=CC=C1NC1=NC=CC=C1C(O)=O CLOMYZFHNHFSIQ-UHFFFAOYSA-N 0.000 claims 1
- 230000008646 thermal stress Effects 0.000 abstract description 2
- 238000012795 verification Methods 0.000 abstract 1
- 239000012530 fluid Substances 0.000 description 8
- 238000011160 research Methods 0.000 description 7
- 238000013459 approach Methods 0.000 description 4
- 230000035939 shock Effects 0.000 description 4
- 230000009286 beneficial effect Effects 0.000 description 3
- 230000008901 benefit Effects 0.000 description 3
- 238000013461 design Methods 0.000 description 3
- 238000011161 development Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000010438 heat treatment Methods 0.000 description 2
- 230000007246 mechanism Effects 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 230000000007 visual effect Effects 0.000 description 2
- 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
- 238000010586 diagram Methods 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 230000009977 dual effect Effects 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000012552 review Methods 0.000 description 1
- 230000002277 temperature effect Effects 0.000 description 1
Abstract
The invention discloses an aircraft heat protection structure heat conduction calculation method based on a FVM-TLBFS method, which is applied to two-dimensional and three-dimensional structure heat conduction calculation for the first time. Aiming at the heat transfer and thermal stress/strain characteristic problems of the thermal protection structure, the existing high-efficiency LBE model is utilized to popularize the heat transfer and strain characteristic problems into structural heat transfer calculation, a two-dimensional and three-dimensional TLBFS flux solver is successfully constructed and used for solving the numerical flux of a structural heat transfer equation, and corresponding numerical calculation example verification is carried out. The method for solving the structural heat conduction is simple and efficient, and can be suitable for structural heat conduction calculation of complex geometric shapes. Meanwhile, the traditional first-order heat transfer precision based on the finite volume method can be improved to second-order precision, and method support is provided for the establishment of a flow field and structure heat transfer integrated calculation method.
Description
Technical Field
The invention relates to a thermal conduction calculation method for a two-dimensional and three-dimensional structure of an aircraft thermal protection structure based on a FVM-TLBFS method, which is suitable for calculating the thermal conduction of the two-dimensional and three-dimensional structure.
Background
The hypersonic aircraft can generate serious pneumatic heating in the long-time flight process in the near space, so that the temperature rise in the structure is obviously increased, the structural safety is seriously endangered, and the design of the heat-proof structure of the aircraft faces great challenges. In the thermal protection problem, the heat transfer of the structure is a basic physical phenomenon, and the thermal safety problem of the thermal structure needs to be considered with great importance: limitation of the surface temperature or heating amount of the structure; the internal temperature of the structure does not exceed the working temperature of the material and the load operation temperature range; the structure can still keep enough strength after the temperature is increased, and the structural damage is avoided; the structural displacement deformation is kept in an allowable range, and the deformation does not influence the normal operation of the movable parts (the full-motion tail wing, the movable hatch cover and the like). The hypersonic flow-solid-thermal multi-field coupling problem analysis must consider the influence of temperature effect, and the effective prediction of structural heat transfer and thermal stress/strain characteristics is an important premise and basis for heat-proof structural design and structural safety assessment. Therefore, the numerical simulation research of the heat transfer characteristics of the heat protection structure is developed, and the analysis and summary of the heat transfer mechanism and law of the heat protection structure are very important for the heat protection design of the aircraft.
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. Compared with macroscopic and microscopic computational fluid mechanics methods, the lattice Boltzman method has the characteristics of clear physical meaning, simple fluid interaction description, simple boundary condition processing, easy program implementation, good parallelism and high model robustness, and is currently considered as an effective means for describing fluid movement and processing engineering problems. 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, for example the DdQm series of basic models (d represents the spatial dimension, m represents the number of discrete velocities) proposed by Qian et al (Qian Y H, D.D' Humi res, P.Lallemannd. Europhysics Letters,199217 6 479484) in 1992. 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 (Qu K, shu C, chew Y t. Physical Review E200775036706) 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, etc. Dellar (Dellar, paul J.progress in Computational Fluid Dynamics, an International Journal 2008884) also derived 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 (Yang L M, shu C, wu j. Computers & Fluids 201379190199) have given a non-free parameter lattice Boltzmann model derivation platform based on the form of moment conservation by theoretical derivation, by which a series of non-free parameter lattice Boltzmann models, such as D1Q3 model, D1Q4 model, D2Q8 model, D2Q12 model, etc., are obtained that solve for compressible, non-viscous flow. 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. Numerous scholars 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 (Yang L M, shu C, wu j. Computers & MATHEMATICS WITH Applications 20167120692081) 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.
Based on the analysis of the research results, the LBM method is found to be widely explored and applied only in hydrodynamic calculation, and the heat transfer calculation in a solid structure is not involved. Meanwhile, through practical analysis, it is found that high-quality structural grids are difficult to produce due to the complexity of a real thermal protection structure and the irregularity of the shape, and great difficulties and challenges are brought to numerical simulation. The finite volume method has good applicability to structural grids and non-structural grids, has excellent geometric flexibility, is suitable for simulating real complex shapes, and meanwhile, the data structure is beneficial to grid self-adaption, so that the heat transfer calculation method based on the finite volume method highlights unique advantages. Therefore, the patent aims to provide the two-dimensional and three-dimensional structure heat conduction calculation method for the aircraft heat protection structure based on the FVM-TLBFS method based on the current heat protection calculation heat transfer problem. The patent promotes the LBE method to structural heat conduction calculation, provides a FVM-TLBFS method, and compared with the TLBM method, the TLBFS method has higher efficiency, can be applied to non-uniform grids, and can be applied to structural heat protection calculation with irregular shapes, so that a hypersonic heat protection structural heat transfer calculation technology is formed. The derivation of the calculated structural heat transfer by the FVM-TLBFS method will be described in detail.
Disclosure of Invention
The invention mainly adopts a numerical calculation means to explore and establish a two-dimensional and three-dimensional structure heat conduction calculation method of the aircraft heat protection structure based on the FVM-TLBFS method.
The technical scheme of the invention is as follows:
The heat conduction calculation method of the two-dimensional and three-dimensional structure of the aircraft heat protection structure based on the FVM-TLBFS method aims at the heat transfer characteristic problem of the heat protection structure, and the LBE method is popularized to the structure heat conduction calculation. The TLBFS method has higher efficiency and can be applied to non-uniform grids. The thermal protection structure is spatially discretized by adopting a FVM method, and the numerical flux of a unit interface is reconstructed by utilizing the local solution of a thermal Lattice Boltzmann Equation (LBE) model, and the D3Q6 model is applied to a TLBFS method for solving the structural thermal conduction equation. The limited volume method based on the grid-core format is adopted for space dispersion, and the display Runge-Kutta method is adopted for time propulsion.
The heat conduction calculation method of the two-dimensional structure of the aircraft heat protection structure specifically comprises the following steps:
1) Structural heat conduction equation derivation based on Boltzmann equation
For isotropic materials, the structural heat transfer control equation can be written in the form of an integral as follows:
Wherein ρ s is the solid structure density, c s is the constant pressure specific heat capacity, T is the solid temperature, k s is the solid heat transfer coefficient, T is the temperature gradient, n is the unit normal vector, S is the control body surface area, Ω is the control body volume, and Q is the heat flow source term.
For incompressible flow heat transfer problems, the heat transfer equation can also be written as follows
Where u is the structural heat transfer equation, κ is the thermal diffusivity, and v is the gradient operator.
When the heat flux term is zero (q=0), equation (1) is equivalent to equation (2) for the structural heat transfer equation (u=0). The relationship between the thermal conductivity k s in equation (1) and the thermal diffusivity k in equation (2) is as follows,
To reduce macroscopic heat transfer equation (2), we achieve this by constructing the Boltzmann equation as follows:
Where h is a temperature distribution function, h eq is an equilibrium state where the temperature distribution function h approaches by particle collisions within the collision time scale τ κ, ζ= (ζ 1,ξ2) is a particle velocity in the particle velocity space, and v h is a temperature distribution function gradient.
2) Thermal lattice Boltzmann model D2Q4 model
The D2Q4 model will be described in detail below:
with the discrete lattice LBE (lattice Boltzmann equation) model, boltzmann equation (4) can be rewritten as follows:
in equation (4), the particle velocity ζ may be written as e α. In two-dimensional heat transfer equation calculations, the D2Q4-LBE model may be applied to discrete grid velocity spaces. In the D2Q4 model, the equilibrium distribution function and discrete particle velocity can be calculated by the following equation:
eα＝(cos[π(α-1)/2],sin[π(α-1)/2])
α＝1,2,3,4 (7)
Through Chapman-Enskog expansion analysis, the relationship between the distribution function of the energy process and the temperature and normal flux can be obtained as follows:
In the method, in the process of the invention, Is a unit interface normal conservation variable,/>For the unit interface normal flux, u 1 is the one-dimensional macroscopic velocity, e α,1 is the first component of particle velocity in the local coordinate system. In addition, the thermal diffusivity, κ, and the time-to-collision scale, τ κ, in equation (9) can be written as follows:
τκ＝2κ (10)
3) Calculation of energy square Cheng Tongliang
For the heat conduction equation (energy equation), the temperature distribution function of the cell interface can also be expressed as
Where h α (0, t) is the distribution function at the x=0 position,For a balanced distribution function at the x=0 position,For an unbalanced distribution function at the x=0 position,/>Is dimensionless collision time.
Substituting equation (11) into equation (9) can obtain the normal energy flux of the cell interface as follows,
We can see that the normal energy fluxAlso consists of two parts: one is F s,I, representing the interface equilibrium distribution function/>The contribution produced; second, F s,II represents the equilibrium distribution function/>, around the cell interfaceIs a contribution of (a). F s,I is directly related to macroscopic flow velocity. When u=0, F s,I is also zero, so for solid structure heat transfer (u=0), only the second partial flux F s,II needs to be calculated.
When F s,II is calculated, the equilibrium distribution function of the unit interface near points should be obtainedFor any variable/>The equilibrium distribution function family at the point of approach of the cell interface can be calculated by the following equation,
In the method, in the process of the invention,And/>Respectively represent the left side and the right side of the unit interface/>Value sum/>Gradient values of (a). /(I)And/>Then represent left and right cell/>, respectivelyValue sum/>Average of gradients of (i.e./>)And/>
Simultaneous equations (12) and (13), for two-dimensional heat transfer equation calculation, can be found using the D2Q4 model
Further, substituting formula (14) into equation (12) while letting u=0, a complete two-dimensional solid heat transfer energy equation flux calculation expression can be obtained as
The last undetermined parameter in the above equation is thus the migration time step δt, which is calculated based on the principle of avoiding numerical extrapolation by the following equation,
Where Δl and Δr are the shortest side lengths of the left and right units of the interface, respectively.
4) Calculation flow
Having described the entire derivation and calculation of the FVM-TLBFS method in detail above, the specific steps implemented by the solver are given below 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 variables at the two sides of the unit interface into equations (6), (7) and (13) to calculate a balanced distribution function;
(3) Calculating a conservation variable on the unit interface through an equation (8);
(4) Calculating the energy flux by equations (9), (14), (15);
(6) Calculating a migration time step δt by equation (16);
(7) Solving a control equation (1) 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;
(8) Repeating the steps (1) - (7) until a temperature distribution convergence solution meeting the condition is obtained.
The three-dimensional structure heat conduction calculation method of the aircraft heat protection structure specifically comprises the following steps:
1) Structural heat conduction equation derivation based on Boltzmann equation
For isotropic materials, the structural heat transfer control equation can be written in the form of an integral as follows:
Wherein ρ s is the solid structure density, c s is the constant pressure specific heat capacity, T is the solid temperature, k s is the solid heat transfer coefficient, T is the temperature gradient, n is the unit normal vector, S is the control body surface area, Ω is the control body volume, and Q is the heat flow source term.
For incompressible flow heat transfer problems, the heat transfer equation can also be written as follows
Where u is the structural heat transfer equation, κ is the thermal diffusivity, and v is the gradient operator.
When the heat flux term is zero (q=0), equation (1) is equivalent to equation (2) for the structural heat transfer equation (u=0). The relationship between the thermal conductivity k s in equation (1) and the thermal diffusivity k in equation (2) is as follows,
To reduce macroscopic heat transfer equation (2), this is achieved by constructing the Boltzmann equation as follows:
where h is a temperature distribution function, h eq is an equilibrium state where the temperature distribution function h approaches by particle collisions within the collision time scale τ κ, ζ= (ζ 1,ξ2) is the particle velocity in the particle velocity space, and v h is a temperature distribution function gradient.
2) Thermal lattice Boltzmann model D3Q6 model
The D3Q6 model will be described in detail below:
with the discrete lattice LBE (lattice Boltzmann equation) model, boltzmann equation (4) can be rewritten as follows:
In equation (4), the particle velocity ζ may be written as e α. In three-dimensional heat transfer equation calculations, the D3Q6-LBE model may be applied to discrete grid velocity spaces. In the D3Q6 model, the equilibrium distribution function and discrete particle velocity are calculated by the following equation:
eα＝(±1,0,0),(0,±1,0),(0,0,±1),α＝1,2,3,4,5,6
the distribution function, time derivative and spatial derivative were developed by Chapman-Enskog expansion analysis into the following form
In the method, in the process of the invention,For the equilibrium distribution function ε is a small parameter proportional to Knudsen number, h α is the distribution function,/>Develop a distribution function for epsilon-order,/>Develop a distribution function for ε 2 level,/>Is the time derivative,/>For the t 1-order time derivative where t1=εt,/>The time derivative is the t2 stage, where t2=εt 2,▽1 and. Substituting equation (7) into equation (5) yields the following three expressions:
Wherein O (epsilon 0) is a epsilon 0 level truncation error, O (epsilon 1) is a epsilon 1 level truncation error, and O (epsilon 2) is a epsilon 2 level truncation error.
Summing equation (9) and equation (10) with respect to index α by the D3Q6 model can result in
According to the Boussinesq approximation, pi (1) can be expressed as
Here, ma is mach number. Simultaneous equations (11) and (12) can be obtained
By comparison, the equation thermal diffusivity, κ, was found to relate to the time-to-collision scale, τ κ, as follows:
τκ＝3κ (14)
In addition, the error term of the recovery energy equation is O (τ κMa2▽2 T). For the simulation of incompressible flow, this error term is ignored.
From this, the distribution function of the energy equation is obtained as follows, with respect to temperature, normal flux:
In the method, in the process of the invention, Is a unit interface normal conservation variable,/>For the unit interface normal flux, u 1 is the one-dimensional macroscopic velocity, x 1 is the one-dimensional direction, and e α,1 is the first component of the particle velocity in the local coordinate system.
3) Calculation of energy square Cheng Tongliang
For the heat conduction equation (energy equation), the temperature distribution function of the cell interface can also be expressed as
Where h α (0, t) is the distribution function at the x=0 position,For a balanced distribution function at the x=0 position,For an unbalanced distribution function at the x=0 position,/>Is dimensionless collision time.
Substituting equation (17) into equation (16) can obtain the normal energy flux of the cell interface as follows,
It follows that the normal energy fluxAlso consists of two parts: one is F s,I, representing the interface equilibrium distribution function/>The contribution produced; second, F s,II represents the equilibrium distribution function/>, around the cell interfaceIs a contribution of (a). Since F s,I is directly related to macroscopic flow velocity. When u=0, F s,I is also zero, so for solid structure heat transfer (u=0), only the second partial flux F s,II needs to be calculated.
When F s,II is calculated, the equilibrium distribution function of the unit interface near points should be obtainedFor any variable/>The equilibrium distribution function family at the point of approach of the cell interface can be calculated by the following equation,
In the method, in the process of the invention,And/>Respectively represent the left side and the right side of the unit interface/>Value sum/>Gradient values of (a). /(I)And/>Then represent left and right cell/>, respectivelyValue sum/>Average of gradients of (i.e./>)And/>
Simultaneous equations (18) and (19), for three-dimensional heat transfer equation calculations, can be found using the D3Q6 model
Further, substituting formula (20) into equation (19) while letting u=0, a complete three-dimensional solid heat transfer energy equation flux calculation expression can be obtained as
The last undetermined parameter in the above equation is thus the migration time step δt, which is calculated based on the principle of avoiding numerical extrapolation by the following equation,
Where Δl and Δr are the shortest side lengths of the left and right units of the interface, respectively.
4) Calculation flow
Having described the entire derivation and calculation of the FVM-TLBFS method in detail above, the specific steps implemented by the solver are given below 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 conservation variables at two sides of the unit interface into equations (6) - (19), and calculating a balanced distribution function;
(3) Calculating conservation variables on the cell interfaces through equations (15) - (17);
(4) Calculating the energy flux by equations (18), (20), (21);
(6) Calculating a migration time step δt by equation (22);
(7) Solving a control equation (1) 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;
(8) Repeating the steps (1) - (7) until a temperature distribution convergence solution meeting the condition is obtained.
The invention has the following beneficial effects: the FVM-TLBFS method is simple and efficient in solving structure heat conduction, can be suitable for complex and real complex geometric shapes, is beneficial to grid self-adaption, and can improve the first-order heat transfer calculation accuracy to second-order based on the traditional finite volume method. Theoretical basis and method support are provided for the establishment of a follow-up flow field and structural heat transfer integrated calculation method.
Drawings
FIG. 1 (a) shows a plate temperature distribution analysis solution.
Fig. 1 (b) is a grid of flat plate structures.
Fig. 1 (c) shows the calculation result of the grid of the flat panel structure.
FIG. 1 (d) is a flat plate unstructured grid.
Fig. 1 (e) shows the calculation result of the flat plate unstructured grid.
Fig. 2 (a) is a comparison of y=0.4m cross-sectional temperature distribution results.
Fig. 2 (b) is a comparison of y=0.8m cross-sectional temperature distribution results.
Fig. 3 is a model of the D2Q4 model particle velocity spatial distribution.
Fig. 4 (a) is a cloud chart of the temperature distribution of the outer surface of the round tube obtained by adopting a least square method at the time t=0.1 s.
Fig. 4 (b) is a cloud chart of the temperature distribution of the outer surface of the round tube obtained by the gaussian-grime method at time t=0.1 s.
Fig. 4 (c) is a cloud chart of the temperature distribution of the inner surface of the circular tube obtained by using a least square method at the time t=0.1 s.
Fig. 4 (d) is a cloud chart of the temperature distribution of the inner surface of the round tube obtained by the gaussian-grime method at time t=0.1 s.
Fig. 5 (a) shows a temperature distribution curve of the cross section of the round tube x=0.25 m at time t=0.1 s.
Fig. 5 (b) shows a temperature distribution curve of the cross section of the round tube x=0.30m at time t=0.1s.
Fig. 5 (c) shows a temperature distribution curve of the cross section of the round tube x=0.35 m at time t=0.1 s.
Fig. 5 (d) shows a temperature distribution curve of the cross section of the round tube x=0.45 m at time t=0.1s.
Fig. 6 is a model of the D3Q6 model particle velocity spatial distribution.
Detailed Description
The present invention will be described below with reference to examples of application.
Example 1: two-dimensional rectangular domain steady state heat conduction:
The two-dimensional rectangular flat plate heat conduction calculation example is a basic calculation example for verifying the accuracy of a heat conduction calculation method. Assuming that the length and width of the two-dimensional plate are a and b, respectively, the temperature boundary conditions are as follows,
Under the condition of no heat source, the temperature distribution of the example has analytical solution,
T(x,y)＝T1(x,y)+T2(x,y) (2)
Wherein,
For ease of calculation, the geometric dimensions of the two-dimensional rectangular flat plate are set to a=1m and b=1m. The plate thermal characteristic parameters are: t 0＝1K,ρs＝1kg/m3,cs＝1J/(kg·K),ks = 1W/(m·k). For this example, calculations were performed using unstructured grids and uniformly structured grids of N x×Ny = 20 x 20 units, respectively, and the gradient calculation was performed using the gaussian-green method.
Fig. 1 (a) to 1 (e) show the temperature distribution calculated with the structured grid and the unstructured grid, respectively, compared with the exact solution. As can be seen from the comparison of the figures, the calculation temperature cloud image of the structural grid and the non-structural grid is almost the same as the analysis solution result. For a more visual observation of the plate temperature distribution changes, fig. 2 (a) and 2 (b) show temperature distribution curves along the cross section (y=0.4 m and y=0.8 m) in the rectangular plate region, respectively. As can be clearly seen from the comparison result, the calculation result is better matched with the analytic solution no matter the structural grid or the non-structural grid is adopted. Further, in order to evaluate the accuracy of the calculation result, the following two error calculation formulas are adopted for evaluation,
Wherein, T num and T exact are the calculated value and the analytical solution respectively,The average Error value and the maximum Error value are respectively defined as Error. Table 1 lists two kinds of grid calculation error results, and it can be seen that the two kinds of grid calculation results have very small errors, but the accuracy of the non-structural grid calculation results is slightly lower than that of the uniform structural grid, because the quadratic diffusion term of the quadrangle orthogonal structural grid is zero, the quadratic diffusion term of the triangular non-structural grid is not zero, the calculation of the quadratic diffusion term causes errors, and therefore, the calculation accuracy is lower than that of the orthogonalization structural grid.
Table 1 two grid computing errors
In addition, in order to study the influence of the grid size on the calculation accuracy, 4 sets of uniform structural grids with different cell sizes were selected for calculation test (N x×Ny =20×20, 40×40, 60×60, 80×80). Two points are randomly selected A, B in the computational domain to compare the error values between the exact solution and the numerical solutions obtained with different size grids. The point A coordinates are 0.9872,0.8123 and are positioned in boundary units of the computing grid; the coordinates of point B are (0.8875,0.1875) located in the interior cells of the computational grid. Table 2 shows the calculated values for A, B points at different cell size grids versus the exact solution. It is clear from the table that progressive refinement of the grid can effectively reduce the error caused by the grid, and in addition, when the grid size is N x×Ny =20×20, the maximum error of the calculation result is also within the effective range, which indicates that the accuracy of the calculation method is reliable.
Table 2 comparison of different size grid calculations
Example 2: three-dimensional circular tube unsteady state heat conduction
To verify the accuracy of the unsteady structure thermal conduction calculation method based on the FVM-TLBFS method in time advance, a three-dimensional round tube unsteady state thermal conduction example in literature (Shang Liying. Solution and application of unsteady state thermal conduction problem on unstructured grid [ D ]. University of south-jing theory, 2006) was selected for testing. The three-dimensional circular tube has the following dimensions: outer diameter R out =1.0m, inner diameter R in =0.5m, length l=0.5m. The calculation conditions of the round tube thermal boundary are as follows: the upper end and the lower end of the round tube are at fixed temperature, the temperature of the upper end is 10 ℃, and the temperature of the lower end is 0 ℃; the outer wall of the circular tube continuously releases heat outwards, the heat flux density is-10W/m 2, the outer wall of the circular tube continuously absorbs heat outwards, and the heat flux density is 50W/m 2. The initial temperature of the round tube is T 0 =0 ℃, and the thermal characteristic parameter is ρ s＝1kg/m3,cs＝1J/(kg·K),ks =1W/(m.K). The calculation grid of the calculation example is a uniform structural grid of N x×Ny×Nz = 40 x 15 x 20, and the actual physical time step value of unsteady calculation is deltat = 0.01s.
Fig. 4 (a) to 4 (d) show calculation cloud diagrams of internal and external temperature distribution of a three-dimensional circular tube at the time Δt=0.1 s. In addition, temperature distribution results obtained by calculation of two gradient calculation methods of a Gaussian-Green method and a least square method are also respectively listed in the graph, and the temperature distribution calculated by the FVM-TLBFS method based on the two gradient calculation methods is found to be almost the same by visual comparison of temperature distribution cloud images. In order to more accurately verify the calculation results, fig. 5 (a) to 5 (d) show comparison results of temperature distribution curves of different sections along the X-axis direction, and it can be found that the temperature distributions obtained by the two gradient calculation methods of the gaussian-green method and the least square method are substantially identical and all agree well with the reference values.
In summary, the developed FVM-TLBFS calculation method was also used to calculate the unsteady solid structure heat conduction with reliable time-advancing accuracy.
Claims (1)
1. The aircraft heat protection structure heat conduction calculation method based on the FVM-TLBFS method is characterized in that the LBE method is popularized to structure heat conduction calculation aiming at the heat transfer characteristic problem of the heat protection structure; TLBFS method is applied to non-uniform grids; performing space dispersion on the thermal protection structure by adopting a FVM method, reconstructing the numerical flux of a unit interface by utilizing the local solution of the thermal lattice Boltzmann equation model, and applying a D3Q6 model to a TLBFS method for solving a three-dimensional structure thermal conduction equation; performing space dispersion by adopting a limited volume method based on a grid format, and performing time propulsion by adopting a display Runge-Kutta method;
the three-dimensional structure heat conduction calculation method of the aircraft heat protection structure specifically comprises the following steps:
1) Structural heat conduction equation derivation based on Boltzmann equation
For isotropic materials, the structural heat transfer control equation can be written in the form of an integral as follows:
wherein ρ s is the solid structure density, c s is the constant pressure specific heat capacity, T is the solid temperature, k s is the solid heat transfer coefficient, The temperature gradient is represented by n, the unit normal vector is represented by S, the surface area of the control body is represented by omega, the volume of the control body is represented by Q, and the heat flow source item is represented by Q;
for the incompressible flow heat transfer problem, the heat transfer equation is written as follows
Where u is the structural heat transfer equation, κ is the thermal diffusivity,Is a gradient operator;
When the heat flux term is zero, q=0, for structural heat transfer equation u=0, equation (1) is equivalent to equation (2); the relationship between the thermal conductivity k s in equation (1) and the thermal diffusivity k in equation (2) is as follows,
To reduce macroscopic heat transfer equation (2), this is achieved by constructing the Boltzmann equation as follows:
Where h is the temperature distribution function, h eq is the equilibrium state of the temperature distribution function h within the collision time scale τ κ, which is approximated by particle collisions, ζ= (ζ 1,ξ2) is the particle velocity in the particle velocity space, Gradient as a temperature distribution function;
2) Thermal lattice Boltzmann model D3Q6 model
The D3Q6 model will be described in detail below:
since a discrete lattice LBE model is employed, boltzmann equation (4) is rewritten as follows:
In equation (4), the particle velocity ζ becomes e α; in two-dimensional heat transfer equation computation, the D2Q4-LBE model can be applied to discrete lattice velocity spaces; in the D3Q6 model, the equilibrium distribution function and discrete particle velocity are calculated by the following equation:
eα＝(±1,0,0),(0,±1,0),(0,0,±1),α＝1,2,3,4,5,6
The distribution function, time derivative and space derivative are developed into the following form by the analysis of binary diffusion number Chapman-Enskog expansion
In the method, in the process of the invention,For the equilibrium distribution function ε is a small parameter proportional to Knudsen number, h α is the distribution function,/>Develop a distribution function for epsilon-order,/>Develop a distribution function for ε 2 level,/>Is the time derivative,/>For the t 1-order time derivative where t1=εt,/>For the t 2-order time derivative where t2=εt 2,/>And/>Are spatial derivatives; substituting equation (7) into equation (5) yields the following three expressions:
Wherein O (epsilon 0) is a epsilon 0 level truncation error, O (epsilon 1) is a epsilon 1 level truncation error, and O (epsilon 2) is a epsilon 2 level truncation error;
Summing equation (9) and equation (10) with respect to index α by the D3Q6 model, results in
According to the Boussinesq approximation, pi (1) is denoted as
Here, ma is mach number; simultaneous equations (11) and (12) to obtain
By comparison, the equation thermal diffusivity, κ, was found to relate to the time-to-collision scale, τ κ, as follows:
τκ＝3κ (14)
in addition, the error term of the recovery energy equation is For simulations of incompressible flows, this error term is ignored;
from this, the distribution function of the energy equation is obtained as follows, with respect to temperature, normal flux:
In the method, in the process of the invention, Is a unit interface normal conservation variable,/>For unit interface normal flux, u 1 is the macroscopic velocity in one dimension, x 1 is the one dimension, e α,1 is the first component of particle velocity in the local coordinate system;
3) Calculation of energy square Cheng Tongliang
For the heat conduction equation, the temperature distribution function of the cell interface is expressed as
Where h α (0, t) is the distribution function at the x=0 position,For a balanced distribution function at the x=0 position,/>For an unbalanced distribution function at the x=0 position,/>Is dimensionless collision time;
Substituting equation (17) into equation (16) to obtain the normal energy flux of the cell interface as follows,
It follows that the normal energy fluxAlso consists of two parts: one is F s,I, representing the interface equilibrium distribution functionThe contribution produced; second, F s,II represents the equilibrium distribution function/>, around the cell interfaceContribution of (2); since F s,I is directly related to macroscopic flow velocity; when u=0, F s,I is also zero, so for solid structure heat transfer u=0, only the second partial flux F s,II needs to be calculated;
When F s,II is calculated, the equilibrium distribution function of the unit interface near points should be obtained For any variable/>The equilibrium distribution function family at the point near the cell interface is calculated by the following equation,
In the method, in the process of the invention,And/>Respectively represent the left side and the right side of the unit interface/>Value sum/>Gradient values of (2); /(I)And/>Then represent left and right cell/>, respectivelyValue sum/>Average of gradients of (i.e./>)And/>
Simultaneous equations (18) and (19), for three-dimensional heat transfer equation calculation, the D3Q6 model is used to calculate
Substituting formula (20) into equation (19) and letting u=0 to obtain the complete three-dimensional solid heat transfer energy equation flux calculation expression as
The last undetermined parameter in the above equation is thus the migration time step δt, which is calculated based on the principle of avoiding numerical extrapolation by the following equation,
Wherein Deltal and Deltar are the shortest side lengths of the left and right units of the interface respectively;
4) Calculation flow
Having described the entire derivation and calculation of the FVM-TLBFS method in detail above, the specific steps implemented by the solver are given below 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 conservation variables at two sides of the unit interface into equations (6) - (19), and calculating a balanced distribution function;
(3) Calculating conservation variables on the cell interfaces through equations (15) - (17);
(4) Calculating the energy flux by equations (18), (20), (21);
(6) Calculating a migration time step δt by equation (22);
(7) Solving a control equation (1) 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;
(8) Repeating the steps (1) - (7) until a temperature distribution convergence solution meeting the condition is obtained.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111303732.4A CN114021499B (en) | 2021-11-05 | Aircraft heat protection structure heat conduction calculation method based on FVM-TLBFS method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111303732.4A CN114021499B (en) | 2021-11-05 | Aircraft heat protection structure heat conduction calculation method based on FVM-TLBFS method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114021499A CN114021499A (en) | 2022-02-08 |
| CN114021499B true CN114021499B (en) | 2024-06-11 |
Family
ID=
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103823976A (en) * | 2014-02-25 | 2014-05-28 | 北京农业信息技术研究中心 | Photo-thermal environment calculating method of solar greenhouse |
| CN107808065A (en) * | 2017-11-23 | 2018-03-16 | 南京航空航天大学 | The solid hot quick calculation method of 3 D complex profile high-speed aircraft stream |
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103823976A (en) * | 2014-02-25 | 2014-05-28 | 北京农业信息技术研究中心 | Photo-thermal environment calculating method of solar greenhouse |
| CN107808065A (en) * | 2017-11-23 | 2018-03-16 | 南京航空航天大学 | The solid hot quick calculation method of 3 D complex profile high-speed aircraft stream |
Non-Patent Citations (2)
| Title |
|---|
| 结构隔绝内-外空间的流-热-固耦合仿真应用;郭亮;王一丁;黄立;谢云恺;童明波;空气动力学学报;20150415;第33卷(第2期);全文 * |
| 高超声速全机外形气动加热与结构传热快速计算方法;李佳伟;王江峰;程克明;伍贻兆;;空气动力学学报;20191215(第06期);全文 * |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Thomas et al. | Geometric conservation law and its application to flow computations on moving grids | |
| Mellor et al. | A survey of the mean turbulent field closure models. | |
| Jameson et al. | Finite volume solution of the two-dimensional Euler equations on a regular triangular mesh | |
| CN108763683B (en) | New WENO format construction method under trigonometric function framework | |
| CN113792432A (en) | Flow field calculation method based on improved FVM-LBFS method | |
| Zhao et al. | Aeroelastic response of rocket nozzles to asymmetric thrust loading | |
| Ye et al. | Numerical investigation on the aerothermoelastic deformation of the hypersonic wing | |
| Yang et al. | Aileron buzz simulation using an implicit multiblock aeroelastic solver | |
| CN115587551B (en) | Multi-scale prediction method for ablation behavior of heat-proof structure of hypersonic aircraft | |
| CN104281730A (en) | Great-rotating-deformation plate shell structure dynamic response finite element analysis method | |
| Huang et al. | Numerical studies of static aeroelastic effects on grid fin aerodynamic performances | |
| An et al. | Reduced-order extrapolation spectral-finite difference scheme based on POD method and error estimation for three-dimensional parabolic equation | |
| Haider et al. | Mathematical analysis of flow passing through a rectangular nozzle | |
| Liu et al. | Modeling and experimental study on free vibration of plates with curved edges | |
| Kamali et al. | Development and validation of a High-Fidelity Aero-Thermo-Elastic analysis capability | |
| CN114021499B (en) | Aircraft heat protection structure heat conduction calculation method based on FVM-TLBFS method | |
| Wen et al. | Numerical and theoretical study on the varying speed impact of wedge bodies on a water surface | |
| Tongqing et al. | CFD/CSD-based flutter prediction method for experimental models in a transonic wind tunnel with porous wall | |
| Yang et al. | Gas kinetic flux solver based high-order finite-volume method for simulation of two-dimensional compressible flows | |
| CN114021499A (en) | Method for calculating heat conduction of aircraft thermal protection structure based on FVM-TLBFS method | |
| Li et al. | An engineering method of aerothermodynamic environments prediction for complex reentry configurations | |
| Patel et al. | A MAC scheme in boundary-fitted curvilinear coordinates | |
| Šekutkovski et al. | A partitioned solution approach for the fluid–structure interaction of thin-walled structures and high-Reynolds number flows using RANS and hybrid RANS–LES turbulence models | |
| Byun et al. | Wing-body aeroelasticity using finite-difference fluid/finite element structural equations on parallel computers | |
| Athavale et al. | Application of an unstructured grid solution methodology to turbomachinery flows |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant |