CN104461677B - A kind of virtual thermal test method based on CFD and FEM technologies - Google Patents
A kind of virtual thermal test method based on CFD and FEM technologies Download PDFInfo
- Publication number
- CN104461677B CN104461677B CN201410601982.XA CN201410601982A CN104461677B CN 104461677 B CN104461677 B CN 104461677B CN 201410601982 A CN201410601982 A CN 201410601982A CN 104461677 B CN104461677 B CN 104461677B
- Authority
- CN
- China
- Prior art keywords
- mrow
- msub
- test
- mfrac
- virtual
- 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
- 238000005516 engineering process Methods 0.000 title claims abstract description 22
- 238000010998 test method Methods 0.000 title claims abstract description 11
- 238000012360 testing method Methods 0.000 claims abstract description 108
- 239000012530 fluid Substances 0.000 claims abstract description 39
- 238000000034 method Methods 0.000 claims abstract description 26
- 230000004907 flux Effects 0.000 claims description 35
- 230000008878 coupling Effects 0.000 claims description 30
- 238000010168 coupling process Methods 0.000 claims description 30
- 238000005859 coupling reaction Methods 0.000 claims description 30
- 239000007787 solid Substances 0.000 claims description 25
- 238000012546 transfer Methods 0.000 claims description 23
- 238000004364 calculation method Methods 0.000 claims description 22
- 239000007921 spray Substances 0.000 claims description 12
- 238000007664 blowing Methods 0.000 claims description 11
- 238000006243 chemical reaction Methods 0.000 claims description 9
- 238000004458 analytical method Methods 0.000 claims description 7
- 239000000463 material Substances 0.000 claims description 4
- 239000013598 vector Substances 0.000 claims description 4
- 239000000853 adhesive Substances 0.000 claims description 3
- 230000001070 adhesive effect Effects 0.000 claims description 3
- 238000006073 displacement reaction Methods 0.000 claims description 3
- 238000010494 dissociation reaction Methods 0.000 claims description 3
- 230000005593 dissociations Effects 0.000 claims description 3
- 238000004088 simulation Methods 0.000 abstract description 4
- 238000012795 verification Methods 0.000 abstract description 2
- 238000007711 solidification Methods 0.000 abstract 1
- 230000008023 solidification Effects 0.000 abstract 1
- 239000007789 gas Substances 0.000 description 19
- 238000013461 design Methods 0.000 description 5
- 230000008569 process Effects 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 230000005540 biological transmission Effects 0.000 description 3
- 238000012937 correction Methods 0.000 description 3
- 238000011161 development Methods 0.000 description 3
- 239000011159 matrix material Substances 0.000 description 3
- 230000003595 spectral effect Effects 0.000 description 3
- 230000003197 catalytic effect Effects 0.000 description 2
- 239000006185 dispersion Substances 0.000 description 2
- 238000010891 electric arc Methods 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 238000010438 heat treatment Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 230000035939 shock Effects 0.000 description 2
- 239000000126 substance Substances 0.000 description 2
- 238000012935 Averaging Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 239000000470 constituent Substances 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000010924 continuous production Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 230000005284 excitation Effects 0.000 description 1
- 238000013213 extrapolation Methods 0.000 description 1
- 238000009413 insulation Methods 0.000 description 1
- 230000008520 organization Effects 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 230000001131 transforming effect Effects 0.000 description 1
- 230000001052 transient effect Effects 0.000 description 1
- 238000011144 upstream manufacturing Methods 0.000 description 1
Landscapes
- Aerodynamic Tests, Hydrodynamic Tests, Wind Tunnels, And Water Tanks (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
A kind of virtual thermal test method based on CFD and FEM technologies, arc tunnel dummy model is set up using Fluid Mechanics Computation method, realizes trystate thermal environment accurate simulation and amendment；Testpieces model is set up using finite element method, is set up including typical model parameter, testpieces model boundary is accurately set up, and determines that test model being capable of accurate simulation actual tests part；Two kinds of model combine analogs are realized using fluid structurecoupling technology, test model amendment and solidification are completed using Modifying model and model verification technique.All kinds of heat structure virtual tests are completed using this method, can significantly shorten the test period, testing expenses are reduced.
Description
Technical Field
The invention relates to a virtual thermal test method based on CFD and FEM technologies, and belongs to the related fields of hypersonic aerodynamic thermal environment, thermal protection design and thermal structure wind tunnel test technologies.
Background
The hypersonic aircraft thermal environment has the characteristics of high enthalpy value, long time, large total heat load and the like, the structural scheme of the whole aircraft thermal protection system is complex, the heat-proof and heat-insulating performance is met, and the integrity of the system is ensured, so that the heat-proof and heat-insulating test assessment of the thermal protection system is significant. Because the thermal assessment test of the thermal protection structure has the characteristics of high test difficulty, long test period and high test expenditure, in order to ensure the project development progress and further improve the design and evaluation level of the hypersonic thermal protection system, a thermal protection test technology and a numerical simulation technology are necessary to be combined to develop a virtual test technology of the thermal protection structure.
In the design of the thermal protection system, aerodynamic thermal environment numerical value and test research are firstly carried out, the scheme design of the thermal protection system is carried out on the basis of definite environmental conditions, and the typical test piece arc wind tunnel thermal assessment test is carried out on the basis of the scheme design. At present, no relevant data of a thermal structure virtual test technology aiming at a high-ultrasonic quick-heating protection system is found at home and abroad. However, the implementation of the thermal assessment test faces great difficulty, and particularly comprises the following three reasons. Firstly, the test difficulty is large: the method is limited to wind tunnel capability, and the assessment indexes cannot be realized for large-size structures, and the test has great risk; secondly, the test period is long: the test period of more than one year is usually needed when the organization completes one heat examination test; thirdly, the test cost is high: due to the characteristics of the thermal protection component of the hypersonic aircraft, the thermal assessment test of the hypersonic aircraft usually reaches the limit of the domestic test level, the consumption is huge, and therefore the test cost is very high.
Although pneumatic heating and structural thermal response have become an important direction of research. However, currently in many research efforts, the two are separated. In fact, these two problems are not isolated, but a unified continuous process. Energy exchange occurs between the structure and the flow field at the object plane through heat flow, and the heat flow calculation accuracy can directly influence the heat transfer calculation result, which is particularly important for the coupling calculation of the heat transfer between the flow field and the structure. In the traditional heat transfer calculation, a flow field is usually calculated firstly, the surface heat flow distribution is given through the flow field calculation, then the heat flow subjected to hot wall correction is loaded on the surface of a solid, and then the temperature distribution in the solid is calculated. However, when the temperature distribution in the structure is calculated by adopting the method, the temperature distribution of the structure is often overestimated to different degrees, so that the calculation method of the flow field/structure coupling heat transfer is generated.
Disclosure of Invention
The technical problem to be solved by the invention is as follows: the defects of the prior art are overcome, and the thermal protection structure virtual test method based on the CFD & FEM technology is provided. The problems of long test and assessment period, high cost, high test technology difficulty and low maturity of the conventional thermal protection structure are solved.
The technical solution of the invention is as follows:
a virtual thermal test method based on CFD and FEM technologies comprises the following steps:
(1) calculating the thermal environment of the hypersonic aircraft, specifically: the method comprises the following steps of carrying out regional division on the outer surface of the hypersonic aircraft, dividing the hypersonic aircraft into an aircraft head spherical region and an aircraft fuselage plane region, and calculating the thermal environment of the aircraft spherical region through a formula:
in the formula: q is heat flux density, ω is 0.52, ρswThe density of the air on the wall surface of the aircraft; mu.sswThe viscosity coefficient of the wall gas of the aircraft is Pr equal to 0.71 and Le equal to 1.0-2.0, rhosGas density at stagnation point, musViscosity coefficient of gas as stagnation point, hDIs the enthalpy of dissociation, hsIs the stagnation enthalpy, hswIs wall enthalpy;for the velocity gradient at the stagnation point, obtained from the modified newton's equation, the calculation formula is:Rois the radius of curvature, psIs the stagnation pressure, p∞Is the far field pressure;
for an aircraft fuselage plan area, the thermal environment calculation is performed by the following formula: when the aircraft surface flow field is in a laminar state,when the aircraft surface flow field is in a turbulent state,in the formula:reynolds number estimated under reference conditions, ρ is far-field density, UeFor reference speed, x is the reference length, gcμ 32.174 reference viscosity coefficient, τm＝1，hrTo recover enthalpy;
(2) selecting a block-shaped area on a surface thermal protection layer of the hypersonic aircraft as a test piece, and carrying out geometric modeling on the test piece in a three-dimensional modeling tool;
(3) determining a wind tunnel used for carrying out a virtual thermal test, wherein the wind tunnel needs to meet the following requirements: a. the diameter of the wind tunnel spray pipe is larger than 1.2 times of the envelope size of the test piece, and b, the maximum value of the heat flow of the outlet of the wind tunnel spray pipe is larger than 1.2 times of the heat flow obtained by calculating the thermal environment of the aircraft in the step (1);
when the selected position of the test piece is the aircraft head, the maximum value of the outlet heat flow of the wind tunnel spray pipe is 1.2 times larger than the heat flow of the aircraft spherical area in the step (1); when the selected position of the test piece is the plane area of the aircraft fuselage, the maximum value of the outlet heat flow of the wind tunnel spray pipe is more than 1.2 times of the heat flow of the thermal environment of the plane area of the aircraft fuselage in the step (1);
(4) establishing a virtual arc wind tunnel CFD model in a computational fluid dynamics tool, and solving a three-dimensional constant compressible Reynolds average NS equation of complete gas of thermochemical reaction through the virtual arc wind tunnel CFD model under a Cartesian coordinate system, wherein the equation is as follows:
whereinIs a conservation variable, t is time，Is the non-adhesive flux in the three directions of x, y and z,obtaining a wall surface heat flux density q for viscous flux vectors in the x, y and z directions;
(5) establishing a virtual test piece FEM heat transfer model in a finite element analysis tool according to a geometric model and a material of the test piece, and calculating through the following heat transfer equation to obtain the wall temperature T of the test piecesw；
The heat transfer equation is:
where ρ is the structure effective density, c is the effective specific heat, T is the temperature, kx,ky,kzEffective heat transfer coefficients in the x, y and z directions;
(6) by the formulaThe wall surface temperature TswPerforming interpolation as the boundary condition for solving the wall surface heat flux density q in the step (4), wherein Ti,j,1Temperature, T, interpolated for fluid boundary (i, j,1) cellsi+1,j,1,Ti1,j,1,Ti,j+1,1,Ti,j-1,1The temperatures obtained by projecting the structural units on the fluid units respectively;
by the formulaInterpolating the wall heat flux density q to obtain the wall temperature T in the step (5)swIn which q isn,snRespectively the heat flux density over the fluid cell n and the cell area,sirespectively representing the heat flux density and the unit area on the structural unit i;
adopting the structural heat transfer characteristic time as the time step of the coupling calculation, and repeatedly and circularly iterating to solve the virtual heat test fluid-solid coupling model obtained in the step (4) and the step (5);
(7) performing a virtual thermal test, namely a wind tunnel blowing test, according to the virtual thermal test fluid-solid coupling model obtained in the step (6), so as to obtain the time-space distribution of the temperature, the time-space distribution of the displacement and the time-space distribution of the stress of the test piece;
(8) comparing the result of the real test piece subjected to the blowing test in the real wind tunnel with the result obtained by the virtual blowing test, and if the comparison result is within a preset range, determining that the virtual thermal test fluid-solid coupling model is accurate; and if the comparison result is out of the preset range, the virtual thermal test fluid-solid coupling model is deemed to be inaccurate, the virtual thermal test fluid-solid coupling model is corrected, and the step (4) is returned.
The computational fluid dynamics tool is CFD + + or FASTRAN.
The finite element analysis tool is ANSYS or ABAQUS.
The method for correcting the fluid-solid coupling model of the virtual thermal test comprises the following specific steps:
(4.1) correcting the CFD model and adjusting the fluid grid scale;
(4.2) correcting the FEM heat transfer model of the test piece, and adjusting the scale of the structural grid;
(4.3) using the formulaThe heat flux density q at the fluid and structure interface is modified to ensure conservation of energy at the fluid and structure interface.
Compared with the prior art, the invention has the beneficial effects that:
(1) the invention provides a novel virtual arc wind tunnel blowing test method, which combines a CFD (computational fluid dynamics) technology and an FEM (finite element modeling) method, realizes the coupling heat transfer of a fluid model and a structural model through a fluid-solid coupling technology, establishes a virtual thermal test fluid-solid coupling model to simulate the state of a real arc wind tunnel blowing test, reduces the risk of failure of the real arc wind tunnel test, shortens the model development period and greatly reduces the project development expenditure;
(2) the invention can numerically simulate the thermochemical reaction complete gas, and can adopt 7-component and 11-component models according to the situation to simulate the blowing test situation of the electric arc wind tunnel more truly;
(3) the invention establishes the interpolation method of wall surface temperature and wall surface heat flow transmission between the fluid unit and the structural unit, can carry out fluid-solid coupling calculation under the condition that the fluid grid is not matched with the structural grid, greatly improves the efficiency of structural heat transmission, and pushes the fluid-solid coupling heat transmission calculation to engineering application from theory;
drawings
FIG. 1 is a flow chart of the present invention;
FIG. 2 is a schematic view of a fluid grid and a structural grid
Fig. 3 is a graph comparing a test result of a typical thermal protection assembly with a virtual test result, wherein fig. 3(a) is a virtual blowing test result, and fig. 3(b) is a real arc wind tunnel test result.
Detailed Description
The following describes embodiments of the present invention in further detail with reference to the accompanying drawings.
The virtual test method for the thermal protection structure can effectively solve the problems of large difficulty, long period and high cost of an actual test. In the establishing process of the method provided by the invention, a CFD technology is used for establishing a virtual wind tunnel model, a FEM method is used for establishing a virtual test piece model, a fluid-solid coupling technology is used for realizing the coupling of the two models, and the models are corrected by using the evaluation result of the thermal protection test to form various thermal virtual test databases.
As shown in fig. 1, the present invention provides a virtual thermal test method based on CFD and FEM technologies, which includes the following steps:
(1) calculating the thermal environment of the hypersonic aircraft, specifically: the method comprises the following steps of carrying out regional division on the outer surface of the hypersonic aircraft, dividing the hypersonic aircraft into an aircraft head spherical region and an aircraft fuselage plane region, and calculating the thermal environment of the aircraft spherical region through a formula:
in the formula: q is heat flux density, ω is 0.52, ρswThe density of the air on the wall surface of the aircraft; mu.sswThe viscosity coefficient of the wall gas of the aircraft is Pr equal to 0.71 and Le equal to 1.0-2.0, rhosGas density at stagnation point, musViscosity coefficient of gas as stagnation point, hDIs the enthalpy of dissociation, hsIs the stagnation enthalpy, hswIs wall enthalpy;for the velocity gradient at the stagnation point, obtained from the modified newton's equation, the calculation formula is:Rois the radius of curvature, psIs the stagnation pressure, p∞Is the far field pressure;
for an aircraft fuselage plan area, the thermal environment calculation is performed by the following formula: when the aircraft surface flow field is in a laminar state,when the aircraft surface flow field is in a turbulent state,in the formula:is the Reynolds number, p, estimated under the reference condition*For far field density, UeFor reference speed, x is the reference length, gcμ 32.174 reference viscosity coefficient, τm＝1，hrTo recover enthalpy;
(2) selecting a block-shaped area on a surface thermal protection layer of the hypersonic aircraft as a test piece, and carrying out geometric modeling on the test piece in a three-dimensional modeling tool;
(3) determining the test working condition, the test condition and the wind tunnel used for carrying out the virtual heat test according to the wall surface heat flux density obtained in the step (1), wherein the wind tunnel needs to meet the following requirements: a. the diameter of the wind tunnel spray pipe is larger than 1.2 times of the envelope size of the test piece, and b, the maximum value of the heat flow of the outlet of the wind tunnel spray pipe is larger than 1.2 times of the heat flow obtained by calculating the thermal environment of the aircraft in the step (1);
when the selected position of the test piece is the aircraft head, the maximum value of the outlet heat flow of the wind tunnel spray pipe is 1.2 times larger than the heat flow of the aircraft spherical area in the step (1); when the selected position of the test piece is the plane area of the aircraft fuselage, the maximum value of the outlet heat flow of the wind tunnel spray pipe is more than 1.2 times of the heat flow of the thermal environment of the plane area of the aircraft fuselage in the step (1);
(4) dividing fluid grids in a computational fluid dynamics tool, establishing a virtual arc wind tunnel CFD model, and solving a three-dimensional highly-constantly compressible Reynolds average NS equation of complete gas of thermochemical reaction through the virtual arc wind tunnel CFD model under a Cartesian coordinate system, wherein the equation is as follows:
wherein,
the viscous stress terms are respectively:
for a complete gas, there are equations of state p ═ ρ RT, h ═ CpT
Total energy per unit gas of
For isotropic fluids, thermal conductivity
In the formulaIs a conservative variable, t is time,is the non-adhesive flux in the three directions of x, y and z,is the viscous flux vector in the three directions x, y and z, u, v and w are the velocities in the three directions x, y and z, CpThe specific heat capacity at constant pressure is obtained, the specific heat ratio gamma is 1.4, and R is 287;
by the followingThe method solves the partial differential equation set to obtain the wall surface heat flux densityn is the wall normal;
a gas model adopted
In the electric arc wind tunnel test process, because the gas temperature is extremely high, chemical reaction occurs, so a thermochemical reaction complete gas model is adopted in the equation. The multi-component mixed gas is composed of various thermal complete gases, and the total internal energy of unit mass of the gas is composed of translational energy, rotational energy, vibration energy and electronic excitation energy in various energy modes. Assuming a two-temperature model, the energy mode of each component i is as follows:
internal energy value and component internal energy of mixed gas:
enthalpy value h and component enthalpy h of mixed gasi：
Wherein, YiIs the percentage of each component i, cp,iThe constant pressure specific heat capacity of the component i,is the chemical enthalpy (enthalpy of formation).
The study mainly focuses on the reaction between air components under high temperature conditions, and for pure air considering only N, O two elements, 7-component and 11-component models and a finite-rate chemical reaction model are respectively adopted.
b value solving method
Spatial dispersion
The spatial discrete format has a great influence on the calculation accuracy and stability of the flow field, and the research adopts the FDS format of Roe which is verified by a great deal of engineering practice to carry out spatial discrete. The construction method according to the Roe format comprises the following steps:
definition ofThen
Represents a pair QL,QRThe conservation variable obtained by calculation by using the Roe averaging method is as follows:
defining:then there is
According to the Jacobian matrixCan split the artificial dissipation term into the sum of three groups of vectors, namely:
wherein:
Δ()＝()R-()L
the FDS format of Roe is an approximate Riemann's solution that has a natural high resolution to the discontinuity problem and is therefore well suited to solving shear-like motion of viscous fluids within the boundary layer. However, because the Roe format is essentially a kind of approximate fitting for transforming the nonlinear problem into the linear Riemannian problem, it is difficult to completely and accurately describe all the characteristics of the nonlinear problem, and sometimes even a wrong physical description is generated, for example, when the hypersonic expansion occurs, the Roe format calculates the expansion shock wave. The reason why such a non-physical phenomenon occurs is that it is difficult for the Roe format to correctly determine the propagation direction of a wave in the case where the eigenvalue tends to zero due to lack of limitation of the entropy condition. For this reason, entropy correction needs to be artificially introduced to dissipate the non-physical expansion shock wave into an expansion sector so that the expansion sector meets the entropy condition. Here the anisotropic Muller-type entropy correction format is used:
wherein sigmanIs the normal spectral radius, στ1,στ2In the form of the tangential spectral radius,the empirical constants are typically taken to be 0.1-0.4.
Time dispersion
For the finite volume semi-discrete format, implicit format processing is adopted for the inviscid flux, and explicit processing is adopted for the sticky flux, which can be written as follows:
where Ω is the unit volume, defining a Jacobian matrix without viscous flux:
then there is
Substituting into the left end of the equation, and finishing to obtain
Wherein
Splitting the non-viscous flux Jacobian matrix according to the windward principle according to the positive and negative eigenvalues,
(AΔQ)i+1,J,K＝(A+ΔQ)I,J,K+(A-ΔQ)I+1,J,K
(AΔQ)i,J,K＝(A+ΔQ)I-1,J,K+(A-ΔQ)I,J,K
(BΔQ)I,j+1,K＝(B+ΔQ)I,J,K+(B-ΔQ)I,J+1,K
(BΔQ)I,j,K＝(B+ΔQ)I,J-1,K+(B-ΔQ)I,J,K
(CΔQ)I,J,k+1＝(C+ΔQ)I,J,K+(C-ΔQ)I,J,K+1
(CΔQ)I,J,k＝(C+ΔQ)I,J,K-1+(C-ΔQ)I,J,K
definition of
Wherein
The maximum of the absolute values of the eigenvalues of the matrices a, B, C, respectively, i.e. the spectral radius, where χ is 1.01, such equation can be further formulated as:
and is also provided with
A+-A-＝rAI B+-B-＝rBI C+-C-＝rCI
Note the book
Then there are: (N + U + L) Δ Qn＝RHS
Decomposed from the approximate LU: (N + L) N-1(N+U)ΔQn＝RHS
The original equation can then be solved in the following two-step format:
c initial and boundary conditions
Initial conditions
The initial flow field is a convergent flow field when the temperature of the wall surface is equal to 300K;
boundary of entrance
Giving a free incoming flow condition according to actual conditions at an upstream boundary;
outlet boundary
Adopting a flow variable extrapolation method at a downstream boundary;
boundary of object plane
The object surface temperature is the temperature of the wall surface of the junction of the fluid and the structure, and the normal pressure gradient is 0. When taking the complete catalytic wall condition, approximately assuming that no atomic component is contained in the chemical constituent, i.e. the mass fraction of the given incoming air flow; when non-catalytic wall conditions are adopted, the object plane material is assumed to have no influence on chemical reaction, namely the normal gradient of the component mass fraction is 0.
(5) Establishing a virtual test piece FEM heat transfer model in a finite element analysis tool according to a geometric model and a material of the test piece, and calculating through the following heat transfer equation to obtain the wall temperature T of the test piecesw；
The heat transfer equation is:
where ρ is the structure effective density, c is the effective specific heat, T is the temperature, kx,ky,kzEffective heat transfer coefficients in the x, y and z directions;
initial and boundary conditions:
the initial conditions were: t (x, y, z) & gtnon & gtt＝0＝300K；
Boundary conditions: calculating the heat flux density of the fluid and structure interface in the step (4) and interpolating the heat flux density on the structure model
Upper heat flux density q*，
Taking heat insulation wall surfaces from other wall surfaces of the structural model;
(6) by the formulaThe wall surface temperature TswPerforming interpolation as the boundary condition for solving the wall surface heat flux density q in the step (4), wherein Ti,j,1Temperature, T, interpolated for fluid boundary (i, j,1) cellsi+1,j,1,Ti1,j,1,Ti,j+1,1,Ti,j-1,1The temperatures obtained by projecting the structural units on the fluid units are respectively shown in fig. 2, which is a schematic diagram of the corresponding relationship between the structural units and the fluid units, wherein the dense background grid is a fluid grid;
by the formulaInterpolating the wall heat flux density q to obtain the wall temperature T in the step (5)swIn which q isn,snRespectively the heat flux density over the fluid cell n and the cell area,siare respectively a structural unitHeat flux density and cell area over i;
the invention adopts the structure heat transfer characteristic time delta t as the time step of the coupling calculation, and the flow field is considered to be transient stable within each time step of the structure calculation. The specific method comprises the following steps:
a. obtaining a steady-state flow field by taking the structure boundary temperature at the time 0 as a boundary condition;
b. under the condition of the heat flow distribution of the flow field and the object plane at the moment 0, solving a structural heat conduction control equation to obtain the temperature distribution of the structure within the time period of 0-delta t;
c. taking the structure boundary temperature at the delta t moment as a boundary condition, and obtaining a stable flow field at the delta t moment;
d. and under the condition of the heat flow distribution of the flow field and the object plane at the time of delta t, solving a heat conduction control equation to obtain the temperature distribution of the structure within the time period of delta t-2 delta t, and then repeating the above process for calculation until the simulation blowing test is finished.
The method is a calculation process of the virtual thermal test fluid-solid coupling model;
(7) performing a virtual thermal test, namely simulating a wind tunnel blowing test, according to the virtual thermal test fluid-solid coupling model obtained in the step (6), so as to obtain the time-space distribution of the temperature, the time-space distribution of the displacement and the time-space distribution of the stress of the test piece;
(8) comparing the result of the real test piece in the real wind tunnel with the result obtained by the virtual blowing test, and if the comparison result is within a preset range (generally, the range is specified to be within 10%), considering that the virtual thermal test fluid-solid coupling model is accurate; and if the comparison result is out of the preset range, the virtual thermal test fluid-solid coupling model is deemed to be inaccurate, the virtual thermal test fluid-solid coupling model is corrected, and the step (4) is returned. Wherein, fig. 3 is the comparison condition between the virtual blowing verification example result and the real blowing test, fig. 3(a) is the virtual blowing test result, fig. 3(b) is the real arc wind tunnel test result, the abscissa axis is the test time, and the ordinate axis is the temperature of the central point of the test piece (model), the comparison proves that the invention can simulate the arc wind tunnel blowing test more accurately.
The computational fluid dynamics tool is CFD + + or FASTRAN.
The finite element analysis tool is ANSYS or ABAQUS.
The method for correcting the fluid-solid coupling model of the virtual thermal test comprises the following specific steps:
(a) correcting the CFD model, and adjusting the fluid grid dimension;
(b) correcting the FEM heat transfer model of the test piece, and adjusting the scale of the structural grid;
(c) by the formulaThe heat flux density q at the fluid and structure interface is modified to ensure conservation of energy at the fluid and structure interface.
Claims (4)
1. A virtual thermal test method based on CFD and FEM technologies is characterized by comprising the following steps:
(1) calculating the thermal environment of the hypersonic aircraft, specifically: the method comprises the following steps of carrying out regional division on the outer surface of the hypersonic aircraft, dividing the hypersonic aircraft into an aircraft head spherical region and an aircraft fuselage plane region, and calculating the thermal environment of the aircraft spherical region through a formula:
<mrow> <mi>q</mi> <mo>=</mo> <mn>0.763</mn> <msup> <mi>Pr</mi> <mrow> <mo>-</mo> <mn>0.6</mn> </mrow> </msup> <msup> <mrow> <mo>(</mo> <msub> <mi>&rho;</mi> <mrow> <mi>s</mi> <mi>w</mi> </mrow> </msub> <msub> <mi>&mu;</mi> <mrow> <mi>s</mi> <mi>w</mi> </mrow> </msub> <mo>)</mo> </mrow> <mn>0.1</mn> </msup> <msup> <mrow> <mo>(</mo> <msub> <mi>&rho;</mi> <mi>s</mi> </msub> <msub> <mi>&mu;</mi> <mi>s</mi> </msub> <mo>)</mo> </mrow> <mn>0.4</mn> </msup> <msqrt> <msub> <mrow> <mo>(</mo> <mfrac> <mrow> <msub> <mi>du</mi> <mi>e</mi> </msub> </mrow> <mrow> <mi>d</mi> <mi>x</mi> </mrow> </mfrac> <mo>)</mo> </mrow> <mi>s</mi> </msub> </msqrt> <mo>&lsqb;</mo> <mn>1</mn> <mo>+</mo> <mrow> <mo>(</mo> <msup> <mi>Le</mi> <mi>&omega;</mi> </msup> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mfrac> <msub> <mi>h</mi> <mi>D</mi> </msub> <msub> <mi>h</mi> <mi>s</mi> </msub> </mfrac> <mo>&rsqb;</mo> <mrow> <mo>(</mo> <msub> <mi>h</mi> <mi>s</mi> </msub> <mo>-</mo> <msub> <mi>h</mi> <mrow> <mi>s</mi> <mi>w</mi> </mrow> </msub> <mo>)</mo> </mrow> </mrow>
in the formula: q is heat flux density, ω is 0.52, ρswThe density of the air on the wall surface of the aircraft; mu.sswThe viscosity coefficient of the aircraft wall gas, Pr ═ 0.71, Le ∈ [1.0, 2.0%]，ρsGas density at stagnation point, musViscosity coefficient of gas as stagnation point, hDIs the enthalpy of dissociation, hsIs the stagnation enthalpy, hswIs wall enthalpy;for the velocity gradient at the stagnation point, obtained from the modified newton's equation, the calculation formula is:
Rois the radius of curvature, psIs the stagnation pressure, p∞Is the far field pressure;
for an aircraft fuselage plan area, the thermal environment calculation is performed by the following formula:
when the aircraft surface flow field is in a laminar state,
when the aircraft surface flow field is in a turbulent state,
in the formula:reynolds number estimated under reference conditions, ρ is far-field density, UeFor reference speed, x is the reference length, gcμ 32.174 reference viscosity coefficient, τm＝1，hrTo recover enthalpy;
(2) selecting a block-shaped area on a surface thermal protection layer of the hypersonic aircraft as a test piece, and carrying out geometric modeling on the test piece in a three-dimensional modeling tool;
(3) determining a wind tunnel used for carrying out a virtual thermal test, wherein the wind tunnel needs to meet the following requirements: a. the diameter of the wind tunnel spray pipe is larger than 1.2 times of the envelope size of the test piece, and b, the maximum value of the heat flow of the outlet of the wind tunnel spray pipe is larger than 1.2 times of the heat flow obtained by calculating the thermal environment of the aircraft in the step (1);
when the selected position of the test piece is the aircraft head, the maximum value of the outlet heat flow of the wind tunnel spray pipe is 1.2 times larger than the heat flow of the aircraft spherical area in the step (1); when the selected position of the test piece is the plane area of the aircraft fuselage, the maximum value of the outlet heat flow of the wind tunnel spray pipe is more than 1.2 times of the heat flow of the thermal environment of the plane area of the aircraft fuselage in the step (1);
(4) establishing a virtual arc wind tunnel CFD model in a computational fluid dynamics tool, and solving a three-dimensional constant compressible Reynolds average NS equation of complete gas of thermochemical reaction through the virtual arc wind tunnel CFD model under a Cartesian coordinate system, wherein the equation is as follows:
<mrow> <mfrac> <mrow> <mo>&part;</mo> <mover> <mi>Q</mi> <mo>&RightArrow;</mo> </mover> </mrow> <mrow> <mo>&part;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mfrac> <mrow> <mo>&part;</mo> <mover> <mi>F</mi> <mo>&RightArrow;</mo> </mover> </mrow> <mrow> <mo>&part;</mo> <mi>x</mi> </mrow> </mfrac> <mo>+</mo> <mfrac> <mrow> <mo>&part;</mo> <mover> <mi>G</mi> <mo>&RightArrow;</mo> </mover> </mrow> <mrow> <mo>&part;</mo> <mi>y</mi> </mrow> </mfrac> <mo>+</mo> <mfrac> <mrow> <mo>&part;</mo> <mover> <mi>H</mi> <mo>&RightArrow;</mo> </mover> </mrow> <mrow> <mo>&part;</mo> <mi>z</mi> </mrow> </mfrac> <mo>=</mo> <mfrac> <mrow> <mo>&part;</mo> <msub> <mover> <mi>F</mi> <mo>&RightArrow;</mo> </mover> <mi>v</mi> </msub> </mrow> <mrow> <mo>&part;</mo> <mi>x</mi> </mrow> </mfrac> <mo>+</mo> <mfrac> <mrow> <mo>&part;</mo> <msub> <mover> <mi>G</mi> <mo>&RightArrow;</mo> </mover> <mi>v</mi> </msub> </mrow> <mrow> <mo>&part;</mo> <mi>y</mi> </mrow> </mfrac> <mo>+</mo> <mfrac> <mrow> <mo>&part;</mo> <msub> <mover> <mi>H</mi> <mo>&RightArrow;</mo> </mover> <mi>v</mi> </msub> </mrow> <mrow> <mo>&part;</mo> <mi>z</mi> </mrow> </mfrac> </mrow>
whereinIs a conservative variable, t is time,is the non-adhesive flux in the three directions of x, y and z,obtaining a wall surface heat flux density q for viscous flux vectors in the x, y and z directions;
(5) establishing a virtual test piece FEM heat transfer model in a finite element analysis tool according to a geometric model and a material of the test piece, and calculating through the following heat transfer equation to obtain the wall temperature T of the test piecesw；
The heat transfer equation is:
where ρ is the structure effective density, c is the effective specific heat, T is the temperature, kx,ky,kzEffective heat transfer coefficients in the x, y and z directions;
(6) by the formulaThe wall surface temperature TswPerforming interpolation as the boundary condition for solving the wall surface heat flux density q in the step (4), wherein Ti,j,1Temperature, T, interpolated for fluid boundary (i, j,1) cellsi+1,j,1,Ti-1,j,1,Ti,j+1,1,Ti,j-1,1Are respectively provided withThe temperature obtained by projecting the structural unit on the fluid unit;
by the formulaInterpolating the wall heat flux density q to obtain the wall temperature T in the step (5)swIn which q isn,snRespectively the heat flux density over the fluid cell n and the cell area,sirespectively representing the heat flux density and the unit area on the structural unit i;
adopting the structural heat transfer characteristic time as the time step of the coupling calculation, and repeatedly and circularly iterating to solve the virtual heat test fluid-solid coupling model obtained in the step (4) and the step (5);
(7) performing a virtual thermal test, namely a wind tunnel blowing test, according to the virtual thermal test fluid-solid coupling model obtained in the step (6), so as to obtain the time-space distribution of the temperature, the time-space distribution of the displacement and the time-space distribution of the stress of the test piece;
(8) comparing the result of the real test piece subjected to the blowing test in the real wind tunnel with the result obtained by the virtual blowing test, and if the comparison result is within a preset range, determining that the virtual thermal test fluid-solid coupling model is accurate; and if the comparison result is out of the preset range, the virtual thermal test fluid-solid coupling model is deemed to be inaccurate, the virtual thermal test fluid-solid coupling model is corrected, and the step (4) is returned.
2. The virtual thermal test method based on CFD and FEM technologies as claimed in claim 1, wherein: the computational fluid dynamics tool is CFD + + or FASTRAN.
3. The virtual thermal test method based on CFD and FEM technologies as claimed in claim 1, wherein: the finite element analysis tool is ANSYS or ABAQUS.
4. The virtual thermal test method based on CFD and FEM technologies as claimed in claim 1, wherein: the method for correcting the fluid-solid coupling model of the virtual thermal test comprises the following steps:
(4.1) correcting the CFD model and adjusting the fluid grid scale;
(4.2) correcting the FEM heat transfer model of the test piece, and adjusting the scale of the structural grid;
(4.3) using the formulaThe heat flux density q at the fluid and structure interface is modified to ensure conservation of energy at the fluid and structure interface.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410601982.XA CN104461677B (en) | 2014-10-30 | 2014-10-30 | A kind of virtual thermal test method based on CFD and FEM technologies |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410601982.XA CN104461677B (en) | 2014-10-30 | 2014-10-30 | A kind of virtual thermal test method based on CFD and FEM technologies |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN104461677A CN104461677A (en) | 2015-03-25 |
| CN104461677B true CN104461677B (en) | 2017-09-29 |
Family
ID=52907781
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201410601982.XA Active CN104461677B (en) | 2014-10-30 | 2014-10-30 | A kind of virtual thermal test method based on CFD and FEM technologies |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104461677B (en) |
Families Citing this family (19)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106557597A (en) * | 2015-09-29 | 2017-04-05 | 中国飞机强度研究所 | A kind of hot virtual test method |
| CN105468846B (en) * | 2015-11-24 | 2019-01-11 | 北京宇航系统工程研究所 | A method of rocket bottom hot-fluid is determined using RADIATION ANGLE COEFFICIENT |
| CN106354970A (en) * | 2016-09-14 | 2017-01-25 | 郑州云海信息技术有限公司 | Analytical method and device for thermal environment |
| CN106570252B (en) * | 2016-10-26 | 2019-05-24 | 中国运载火箭技术研究院 | A kind of thermal protection system design method based on probabilistic technique |
| CN106841280B (en) * | 2016-11-14 | 2019-04-30 | 中国航天空气动力技术研究院 | Sharp leading edge determination method of heat-flow density under the conditions of a kind of arc tunnel |
| CN108120612B (en) * | 2016-11-29 | 2019-11-15 | 北京机电工程研究所 | A kind of processing method of heat test high temperature peak |
| CN108228921B (en) * | 2016-12-13 | 2021-07-13 | 北京空天技术研究所 | Wind tunnel test device and design method thereof |
| CN108304595B (en) * | 2017-05-04 | 2021-04-02 | 北京空天技术研究所 | Structural temperature analysis method for hypersonic aircraft semi-closed area |
| CN107563030A (en) * | 2017-08-22 | 2018-01-09 | 哈尔滨工程大学 | A kind of mesh free analogy method for being directed to two kinds of fluid heat transferrings and handing over mixed broken phase transition process |
| CN108332934B (en) * | 2017-11-15 | 2019-03-05 | 北京空天技术研究所 | A kind of arc tunnel test method of non-ablative thermally protective materials/structure |
| CN109670216B (en) * | 2018-11-30 | 2023-07-04 | 中国船舶重工集团公司第七一九研究所 | Passive residual heat removal condenser position optimization design method based on CFD technology |
| CN109856991B (en) * | 2019-01-22 | 2021-10-15 | 华东师范大学 | Dynamic virtual human simulation method based on kinetic energy and thermal distribution diagram |
| CN110261009B (en) * | 2019-07-09 | 2021-04-06 | 北京航空航天大学 | Wall surface heat flow density measuring method and device based on linear bottom layer hypothesis |
| CN110781629B (en) * | 2019-11-20 | 2023-04-21 | 桂林理工大学 | Method and system for determining convection heat dissipation coefficient |
| CN112035933B (en) * | 2020-09-03 | 2022-02-22 | 西北工业大学 | Solid rocket engine jet pipe thermal structure coupling analysis method considering structural clearance |
| CN112989485B (en) * | 2021-01-18 | 2022-04-12 | 中国空气动力研究与发展中心计算空气动力研究所 | Along trajectory heat flow interpolation method based on hot wall correction |
| CN114088253A (en) * | 2021-11-17 | 2022-02-25 | 华电国际电力股份有限公司十里泉发电厂 | Water-cooled wall backfire side heat flowmeter and online monitoring method |
| CN114840950B (en) * | 2022-07-04 | 2022-09-16 | 中国航空工业集团公司沈阳空气动力研究所 | Support layout design method for flexible plate system of wind tunnel flexible wall spray pipe |
| CN117963157B (en) * | 2024-03-28 | 2024-08-06 | 南京工业大学 | Thermal test method and system for multi-temperature-zone structure of full-size hypersonic aircraft |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5926399A (en) * | 1996-03-04 | 1999-07-20 | Beam Technologies, Inc. | Method of predicting change in shape of a solid structure |
| CN102254068A (en) * | 2010-12-01 | 2011-11-23 | 东南大学 | Multi-scale analyzing method for buffeting response of large-span bridge |
| CN103593518A (en) * | 2013-10-31 | 2014-02-19 | 中国运载火箭技术研究院 | Aircraft model modification system based on modal test data |
-
2014
- 2014-10-30 CN CN201410601982.XA patent/CN104461677B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5926399A (en) * | 1996-03-04 | 1999-07-20 | Beam Technologies, Inc. | Method of predicting change in shape of a solid structure |
| CN102254068A (en) * | 2010-12-01 | 2011-11-23 | 东南大学 | Multi-scale analyzing method for buffeting response of large-span bridge |
| CN103593518A (en) * | 2013-10-31 | 2014-02-19 | 中国运载火箭技术研究院 | Aircraft model modification system based on modal test data |
Also Published As
| Publication number | Publication date |
|---|---|
| CN104461677A (en) | 2015-03-25 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104461677B (en) | A kind of virtual thermal test method based on CFD and FEM technologies | |
| CN108304684B (en) | Rocket engine tail jet flow simulation method and system | |
| Marcantoni et al. | High speed flow simulation using openfoam | |
| CN107273593A (en) | A kind of turbulence model and its method for building up predicted for High Mach number intense shock wave flow field Aerodynamic Heating | |
| CN105095603A (en) | Multi-field coupling transient numerical method for hypersonic flow-heat transfer and structural response | |
| CN113792432A (en) | Flow field calculation method based on improved FVM-LBFS method | |
| CN114757070A (en) | New WENO format construction method under trigonometric function framework for numerical simulation | |
| Jiang et al. | On the use of thermally perfect gas model for heating prediction of laminar and turbulent SWBLI | |
| CN110826261A (en) | Buried gas pipeline leakage simulation method based on fluent | |
| CN111695307B (en) | Water hammer finite volume simulation method considering dynamic friction resistance explicitly | |
| CN113051846A (en) | Wall surface first layer grid thickness estimation method considering compressible and heat conduction effects | |
| Ding et al. | Thermal effect on mass flow-rate of sonic nozzle | |
| Guangyue et al. | Preliminary study on the influence of aerothermoelastic deformation on 2-D hypersonic inlet performance | |
| CN108228921B (en) | Wind tunnel test device and design method thereof | |
| CN111695219B (en) | Stress prediction method of skin plate covered with thermal protection coating under supersonic flight condition | |
| Lei et al. | Inverse determination of multiple heat sources’ release history in indoor environments | |
| CN110909511B (en) | Non-viscous low-speed streaming numerical simulation method without curved surface volume division | |
| CN114021499A (en) | Method for calculating heat conduction of aircraft thermal protection structure based on FVM-TLBFS method | |
| Didwania et al. | Analysis of turbulent flow over a 90 bend of duct using in centralized AC plant by CFD code | |
| Balam et al. | Numerical solution of natural-convection flow in enclosures: An implicit vorticity boundary condition type method | |
| Si et al. | Development of simulation program for temperature field of mass concrete structures | |
| Jang et al. | A scheme for improving computational efficiency of quasi-two-dimensional model | |
| Dreyer et al. | Towards characterization of relevant fidelity modeling of loads for maneuvering hypersonic vehicles | |
| Özkökdemir et al. | Investigation of the mechanical integrity of wings/fins under thermal loading | |
| Reinartz et al. | Computation of hypersonic double wedge shock/boundary layer interaction |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |