CN107832494B - Hypersonic aircraft leading edge flow-heat-solid integrated calculation method - Google Patents

Hypersonic aircraft leading edge flow-heat-solid integrated calculation method Download PDF

Info

Publication number
CN107832494B
CN107832494B CN201710951385.3A CN201710951385A CN107832494B CN 107832494 B CN107832494 B CN 107832494B CN 201710951385 A CN201710951385 A CN 201710951385A CN 107832494 B CN107832494 B CN 107832494B
Authority
CN
China
Prior art keywords
calculation
flow
heat
temperature
flow field
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.)
Expired - Fee Related
Application number
CN201710951385.3A
Other languages
Chinese (zh)
Other versions
CN107832494A (en
Inventor
李佳伟
王钰涵
杨天鹏
王江峰
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing University of Aeronautics and Astronautics
Original Assignee
Nanjing University of Aeronautics and Astronautics
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Nanjing University of Aeronautics and Astronautics filed Critical Nanjing University of Aeronautics and Astronautics
Priority to CN201710951385.3A priority Critical patent/CN107832494B/en
Publication of CN107832494A publication Critical patent/CN107832494A/en
Application granted granted Critical
Publication of CN107832494B publication Critical patent/CN107832494B/en
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Automation & Control Theory (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Investigating Or Analyzing Materials Using Thermal Means (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

The invention discloses a hypersonic aircraft leading edge flow-heat-solid integrated calculation method, and belongs to the field of aircraft pneumatic calculation. Aiming at the problem of complex coupling of hypersonic flow pneumatic heating and structural heat transfer, the method avoids the complex data exchange and calculation amount caused by the fact that the traditional pneumatic heating/structural heat transfer coupling solving method carries out flow field and structural coupling alternative iterative calculation in a time domain, takes the flow field and the structure as a physical field, and adopts a unified control equation. Redefining physical parameters of the flow-solid interface, carrying out space dispersion of a finite volume method in a full physical field, and carrying out time advance by adopting implicit time iteration. Compared with a coupling algorithm, the method does not need additional data exchange and coupling strategies, the calculation result is closer to the experimental value, the calculation amount and the grid dependency are relatively small, and the method has better stability and calculation accuracy.

Description

Hypersonic aircraft leading edge flow-heat-solid integrated calculation method
Technical Field
The invention belongs to the technical field of aerodynamic calculation of aircrafts, and particularly relates to a hypersonic speed aircraft leading edge flow-heat-solid integrated calculation method.
Background
Hypersonic flow generally refers to flow with a mach number greater than 5. When the aircraft goes in and out of the atmosphere at hypersonic speed or continuously flies in space, due to the compression effect and the violent friction between the surface of the aircraft and the air (Wangjiang peak, five-mytilus, Jiweiwan, Xiaoxian peak, Zhao Faming, Lu dete 2015 Ainsliao 36(1):159 and 175' numerical method research progress of hypersonic speed complicated aerodynamic problem), the head part, the front edge of an air inlet channel and other key parts of the aircraft can bear huge aerodynamic heating, and strong aerodynamic force, aerodynamic heat and structural coupling problems can be generated, thereby bringing great hidden danger to flight safety. Therefore, the physical process of accurate prediction, pneumatic heating and structural heat transfer plays an important role in the light weight design of the heat protection system of the high-speed aircraft. Because the ground experiment difficulty of the problems is high and the period is long, the numerical simulation technology is mainly adopted for analyzing the problems at present.
At present, numerical simulation of the aerodynamic heat/structure heat transfer coupling problem of the hypersonic aircraft mainly comprises two methods of partition coupling calculation and integrated solution. A traditional multi-field partition coupling method (Xia Gang, Liu Xin, Cheng Wen science and technology university's bulletin 25(1): 35-39' numerical calculation of blunt body hypersonic pneumatic heating and structural heat transfer coupling '; Yao Xiaohu, Han Qiang 2008' Physics bulletin 57(8):5056 and 5062 'twisting and buckling of double-layer carbon nano tubes under the action of thermal coupling') divides a flow field and a structure into two independent parts, and performs data exchange of two parameters of the flow field heat flow density and the structure surface temperature on a coupling interface in a coupling alternating iteration mode on a time domain. Deshaumphi (Wieting A R, Deshaumphai P, Bey K S1991 Thin wall.struct.11112) of the NASA Lanli research center considers that the mathematical model of the partition coupling method needs an additional data transmission strategy, and artificially divides the originally continuous physical process, thereby generating calculation errors and influencing the accuracy of the calculation results.
The idea of the integrated solving method is proposed in the last 70 th century, and researchers at home and abroad also have some research results. Professor Thornton (Dechaumaphi P, Thornton E A, Wieting A R1989J. Spacecraft 26201209) is equal to that in 1988, a finite element method is adopted to carry out integrated solution research on a flow field and a solid structure, and a two-dimensional circular tube pneumatic heating test (Allen R W1987 NASA TM-100484) is carried out in an 8-foot NASA Lanli high-enthalpy wind tunnel, so that the effectiveness of the method is verified, but because the discrete process of the finite element method still only adopts a simple interpolation method, obvious non-physical shock occurs at shock wave equal discontinuities during calculation, the shock wave resolution is low, so that shock waves need to be captured by encrypting a grid at the shock wave, the calculation method is poor in flexibility, and the application is limited. Related researches are carried out on the two-dimensional flow field, thermal and structure integrated numerical simulation by Huangtang and the like (Huangtang, Mao nationality, Jiang Guiqing, ZhouWeijiang 2000, aero-dynamics newspaper 18(1):115 + 119 'two-dimensional flow field, thermal and structure integrated numerical simulation'), wherein the flow field adopts a finite difference method based on a TVD format to carry out numerical dispersion, the structure heat transfer adopts a mature finite element method, the two satisfy an energy balance equation at a coupling interface, but the coupling calculated amount and the calculation error are large and the like due to the difference in the order of magnitude of the flow field and the structure time step length. The method considers that in order to really realize the integrated calculation of the flow field, the heat and the structure in engineering, the integrated numerical calculation concept of the flow field, the heat and the structure must be expanded, the continuous and integral physical change conditions of the flow field, the heat and the structure are met, and the calculation precision of the thermal structure can be ensured. In addition, Jiang Gui Qing et al (Jiang Gui Qing, Tong's class, Cao Tree Sound 1992 mechanics and practice 14(3):1-8 "computational pneumatic thermodynamics mainly based on finite element method") think that along with the development and deepening of aerospace technology and pneumatic thermodynamics, modern pneumatic thermodynamics must solve the problem of solving the three-dimensional integration of gas-heat-structure, analyze the advantages of the finite element method for solving the problem, and simultaneously show that the strong interruption problem is not solved well in the finite element method. Gunn in China aerodynamic research and development center and the like (Gunn, Zhang Xun, Shenqing 2002, journal of aerodynamic science 20(4): 422-. The method omits the time change details of the smaller characteristic time, has the problem of solving rigidity, and has to be improved in the calculation adaptability to the complex aerodynamic shape.
Disclosure of Invention
In view of the defects of the prior art, the invention aims to provide a hypersonic aircraft leading edge flow-heat-solid integrated calculation method to solve the problem of complex data exchange and calculation amount caused by the fact that the traditional pneumatic heating/structural heat transfer coupling solving method carries out flow field and structural coupling alternative iterative calculation in a time domain.
In order to achieve the purpose, the technical scheme adopted by the invention is as follows:
the invention discloses a hypersonic aircraft leading edge flow-heat-solid integrated calculation method, which comprises the following steps of:
the flow field and the structure are used as the same physical field, physical parameters and thermodynamic properties of the flow field and the structure are calculated simultaneously, an interface of the flow field and the structure is used as an internal boundary of the whole physical field, a heat transfer control equation of the flow field and the structure is established simultaneously, and a unified numerical calculation method is adopted for solving simultaneously.
Preferably, for structural heat transfer, the integral of the structural heat transfer control equation over the control volume Ω s is formed as follows, regardless of the heat source:
Figure BDA0001432928850000021
wherein dS is a control body surface unit, CsIs the specific heat capacity of a solid material, T is the structure temperature,
Figure BDA0001432928850000022
is the temperature gradient, rho is the material density, k is the thermal conductivity coefficient,
Figure BDA0001432928850000023
is a boundary; the flow field calculation adopts a compressible Reynolds average N-S (RANS) equation, and on a control body omega s, a flow field control equation and a structural heat transfer control equation are unified into a control equation in the same integral form:
Figure BDA0001432928850000024
wherein W is a constant, FcFor convection flux, FvIs the viscous flux; it is defined as follows:
Figure BDA0001432928850000025
for a streamAnd field calculation is carried out, the density rho, the pressure p and the temperature T meet an ideal gas state equation p which is rho RT, u, v and w are speeds of the control body in three directions respectively, E is total energy of the control body of the fluid unit, H is total enthalpy H which is E + p/rho, k is heat conductivity coefficient, and tau isijFor viscous stress tensor, the structure is undeformed for structure heat transfer calculations with u-v-w-0 and zero F for the above formula convection flux c0, wherein E ═ CsT is in the solid unit control body;
an SST k-omega two-equation model is selected as the turbulence model, a damping function is not needed in the model, the development of near-wall turbulence can be well simulated, the model is applied to a region with strong boundary layer viscous interference, and meanwhile, a free shear layer can be better simulated and is applied to turbulence far away from the wall surface; and during space dispersion, the convection flux is dispersed by adopting an AUSM + format based on a lattice format, and the viscous flux dispersion adopts a second-order central difference format. In order to accelerate the calculation, in the time dispersion, the non-constant calculation adopts double-time implicit time iteration, and the constant calculation adopts LU-SGS implicit time iteration.
Preferably, the calculation method comprises the following boundary conditions:
(1) a flow field far field boundary condition which adopts a Riemann boundary condition;
(2) an object plane boundary condition comprising a solid surface flow boundary condition and a thermodynamic boundary condition;
a. the boundary condition of the solid surface flow interface meets the non-slip boundary condition, and the pressure gradient is zero;
b. structural heat conduction boundary conditions are thermodynamic boundary conditions, which include Dirichlet temperature boundary conditions and Neuman heat flow density boundary conditions.
Preferably, the calculation method further comprises: the interface temperature calculation was calculated using the center-mean method as follows:
T＝(Tl+Tr)/2
in the formula, Tl、TrThe temperatures of the interface left and right control units are respectively; the calculation of temperature gradient ∑ T requires correction, the calculation method being as follows:
Figure BDA0001432928850000031
Figure BDA0001432928850000032
in the formula (I), the compound is shown in the specification,
Figure BDA0001432928850000033
temperature gradients, L, of left and right control units of the interface, respectivelylrIs the distance between cell centers, rlrThe unit vector is from the center point of the left control unit to the center point of the right control unit; in addition, the calculation method of the temperature gradient employs a gaussian green method, which calculates by the value of the boundary of the cell and the normal vector:
Figure BDA0001432928850000034
for a lattice format, this is:
Figure BDA0001432928850000035
wherein Ω is the volume of the control body, i is the number of the cell, j is the number of the adjacent cell, and nijIs a unit normal vector, Δ SijThe surface area of the control body is shown, and N is the unit surface number of the control body; in order to improve the accuracy of the calculation, the concept of thermal impedance in solid heat transfer is introduced, and the thermal impedance is defined as:
Figure BDA0001432928850000041
in the formula, d is the thickness in the heat flow direction, A is the cross-sectional area perpendicular to the heat flow direction, and k is the heat conductivity coefficient; the thermal impedance relation of the interface is obtained by introducing the concept of thermal impedance:
Rt,bnd＝Rt,l+Rt,r
the calculated equation of the interface thermal conductivity k is obtained as follows:
Figure BDA0001432928850000042
wherein k isl、krRespectively, the heat transfer coefficients of the left and right control units, Ll、LrThe distances from the center of the left and right control units to the center of the boundary, respectively.
The invention has the beneficial effects that:
the invention avoids the complex data exchange and calculation amount caused by the flow field and structure coupling alternative iterative calculation in the time domain by the traditional pneumatic heating/structure heat transfer coupling solving method, takes the flow field and the structure as a physical field and adopts a unified control equation. Redefining physical parameters of the flow-solid interface, carrying out space dispersion of a finite volume method in a full physical field, and carrying out time advance by adopting implicit time iteration. The integrated calculation method is verified by adopting a typical hypersonic velocity circular tube steady-state/unsteady-state flow-heat-solid coupling calculation example, the maximum stagnation temperature of the circular tube in a steady state reaches 648K, and the heat flow density and the structure temperature in an unsteady state are well matched with reference documents and experimental values, so that the reliability and the correctness of the method are proved. Meanwhile, the comparison and analysis results of the coupling calculation method show that the calculation result obtained by the integrated calculation method is closer to an experimental value, the calculation amount and the grid dependency are relatively small, and the stability and the calculation precision are better.
Drawings
FIG. 1a is a schematic diagram of a computational grid.
FIG. 1b is a schematic diagram of boundary conditions.
FIG. 2a is a schematic diagram of the temperature distribution of a round tube structure.
Fig. 2b is a schematic view of the temperature distribution of the flow field.
FIG. 3a is a schematic comparison of a density contour plot and a test schlieren plot.
FIG. 3b is a schematic diagram showing the comparison of the pressure distribution on the surface of the round tube with the test value.
FIG. 4a is a temperature cloud graph of the round tube structure at 0.1 s.
FIG. 4b is a temperature cloud of the tube structure at 0.5 s.
FIG. 4c is a temperature cloud of the tube structure at 1.0 s.
FIG. 4d is a temperature cloud of the tube structure at 2.0 s.
FIG. 5 is a schematic representation of the temperature at the stagnation point of a round tube as a function of time.
FIG. 6 is a comparison of the temperature distribution of the round tube structure at the 2s time with that of the literature.
Fig. 7a is a cloud graph of flow field temperature at 0.1 s.
Fig. 7b is a cloud graph of the flow field temperature at 0.5 s.
Fig. 7c is a cloud graph of flow field temperature at 1.0 s.
Fig. 7d is a cloud graph of flow field temperature at 2.0 s.
Fig. 8 is a schematic diagram of the temperature change of the flow field along the symmetry line at different moments.
FIG. 9 is a comparison of the heat flow distribution on the surface of the round tube at time 0 s.
FIG. 10 is a graph illustrating the change in stagnation heat flow with time.
Fig. 11 is a schematic diagram of a serial coupling iterative method.
FIG. 12 is a flowchart of a coupling solution method calculation.
FIG. 13a is a schematic of the change in stagnation temperature over time.
FIG. 13b is a schematic of the temperature difference as a function of time.
FIG. 14 is a diagram of a physical model of the flow-heat-solid coupling problem.
Detailed Description
In order to facilitate understanding of those skilled in the art, the present invention will be further described with reference to the following examples and drawings, which are not intended to limit the present invention.
1. Technical parameters
Wieting. A.R (Allen R W1987 NASA TM-100484) A pneumatic heating test of the leading edge of stainless steel round tubes was performed in 1987 in an 8-foot high enthalpy wind tunnel at the NASA Lanli research center, USA (Wieting A R, Holden M S1987 AIAA 22nd Thermophysics Conference Honolulu Hawaii June 8-10,198)7) This test has been used many times to verify the accuracy of the flow field structure heat transfer coupling calculations. The invention also selects infinite-length stainless steel round tube models with the same conditions, and the size of the round tube is the outer diameter Rout0.0381m, inner diameter Rin0.0254m, the structural material is AISI321 series (1Crl8Ni9Ti) stainless steel. The thermodynamic parameters are shown in table 1, and the incoming flow condition calculation parameters are given in table 2 as follows:
TABLE 1
Figure BDA0001432928850000051
TABLE 2
Figure BDA0001432928850000061
2. Model mesh
The flow field calculation and the structure calculation are the same set of grids, fig. 1a, the flow field grid amount is about 37100, the structure grid amount is about 3800, and the grid height of the first layer of the near object plane is about 1 multiplied by 10-6And m is selected. Fig. 1b is a boundary condition, the two ends of the flow field are pressure outlet boundary conditions, the two ends of the circular tube are heat insulation walls, the inner wall is an isothermal wall, and the initial temperature of the circular tube is 294.4K.
3. Numerical simulation results
And (3) performing hypersonic flow-solid-heat integrated steady and steady numerical calculation on the stainless steel circular tube, wherein in an unsteady state, the actual calculation physical time is 2s, and the actual time step length is equal to delta t which is 0.001 s.
3.1 Steady State results analysis:
the structural temperature distribution of the round tube in the steady state is shown in fig. 2a, and it can be seen from the graph that the high-temperature area of the round tube structure in the steady state is distributed in the stagnation point area, the temperature of the stagnation point area is 648K at the maximum, the temperature rise reaches 353.6K, the temperature of the inner wall is raised from 294.4K initially to 305K, and the temperature rise is 10.6K less than that of the inner wall. Fig. 2b shows the temperature distribution of the flow field outside the circular tube in a steady state, and the maximum temperature of the flow field reaches 2263K after the flow field is heated by bow shock wave. Fig. 2a and fig. 2b well show that one of the advantages of the method of the present invention is that the temperature distribution of the steady-state structure and the flow field can be calculated quickly, and the hypersonic flow-solid-thermal steady-state solving problem is well solved. The current coupling calculation method is to perform parameter exchange of heat flow and temperature on the interface of a flow field and a structure, if the temperature distribution of the flow field and the structure in a steady state is solved, iteration is required for multiple times until the temperature of the structure is converged, and repeated iteration of the flow field and the structure brings great calculated amount.
FIG. 3a is a comparison of the calculated density contour plot with a test schlieren plot, the upper half being a photograph of the test schlieren plot and the lower half being a cloud plot (Dechaumaphai P, Thornton E A, WietingA R1989J. Spacecraft 26201209), showing that the shock positions are substantially coincident. Meanwhile, the comparison between the pressure distribution (normalization) on the surface of the circular tube and the test value is shown in fig. 3b, and the coincidence degree between the pressure distribution along the circular tube and the test value is better, so that the accuracy of the calculation of the steady-state flow field is illustrated.
3.2 analysis of temperature field characteristics of the unstable structure:
fig. 4 shows the temperature distribution cloud charts of the circular tube structure at different times (t ═ 0.1s, 0.5s, 1.0s, 2.0s) in the unsteady state calculation. The change of the temperature distribution in the circular tube structure within 2 seconds can be obviously seen from the graph, the heat generated by pneumatic heating is conducted in the structure along with the advance of time, the whole temperature of the structure is increased, the high-temperature area is gradually increased from the stagnation point area, and the stagnation point temperature is always highest. Within 2 seconds, the stagnation point temperature is increased from 294.4K to 390.2K, which is slightly higher than 388.8K calculated by Dechaumpai, the error is about 0.4%, and meanwhile, the temperature of the inner wall of the circular tube basically keeps the initial temperature of 294.4K, which indicates that the structural thermal response caused by aerodynamic heating within 2 seconds does not affect the inner wall of the structure temporarily.
Fig. 5 shows the change of the stagnation point temperature of the circular tube with the calculation time, and compared with the reference, the calculation result is in good agreement with the reference (huangjie 2013 scholar thesis hypersonic aircraft flow thermal solid multi-physical field coupling calculation research, harbin industrial university), and the maximum error of the stagnation point temperature at 2s is 3.2K, which is detailed in table 3. It can be seen from the figure that the initial time stagnation point temperature rises sharply, and the temperature rise degree gradually tends to be gentle as time advances, because the initial time stagnation point heat flow is the largest, and the structure temperature rises, the stagnation point heat flow rate gradually decreases, and the temperature rise trend gradually decreases as heat is transferred in the structure. Table 3 is as follows:
TABLE 3
Figure BDA0001432928850000071
In table 3, document 1: gunn, zhangxian, shenqing 2002 aerodynamic reports 20422; document 2 huangjie 2013 thesis of hypersonic aircraft flow thermal solid multi-physics coupling calculation research, harbin university of industry; document 3: dechaumaphai P, Thornton E A, Wieting A R1989 J.Spacecraft 26201.
In order to better analyze the accuracy of the calculation method, fig. 6 shows a comparison graph of the temperature distribution contour line of the round tube when t is 2s, the upper half part is a reference result, and the lower half part is a calculation result of the method, so that the calculation result is well matched with the reference result.
3.3 analysis of unsteady flow field characteristics:
fig. 7 a-7 d are contour graphs of the flow field temperature in the stagnation area at different times, and fig. 8 is a curve of the flow field temperature variation along a symmetrical line (y is 0) at different times. As can be seen from the figure, the temperature rises sharply from 241.5K to 2163K after the incoming flow is heated by the bow shock, which is not much different from 2166.7 of the above-mentioned document 2. It can be seen from the figure that the position of the bow shock wave is about-54.5 mm, which is better matched with-54.5 mm calculated by an empirical formula (Billig F S1967J. Spacecraft 4822) in a complete gas state, and the condition that the chemical unbalance effect is not considered in the calculation of the invention is reasonable. A temperature boundary layer area exists near the stagnation point, the thickness of the temperature boundary layer area is about 3% of the thickness of the bow shock wave, a large temperature gradient exists in the temperature boundary layer, the temperature drops sharply from 2163K to 294.4K, and therefore a large heat flow density is generated in the stagnation point area.
FIG. 9 is a comparison of the heat flow distribution on the surface of the round tube at the initial time (t ═ 0) and the test valueThe heat flow calculation distribution result (normalization processing) of the invention is consistent with the experimental value. The heat flow value at the initial moment is maximum at the stagnation point, and the maximum heat flow obtained by calculation is 49.71 multiplied by 104W/m2Slightly higher than Fay-Riddell (Fay J A, Riddell F R1958 Journal of the Aeronoutic Sciences 257385) formula, 48.27 × 104W/m2And the viscous shock layer formula (Holcomb J E, Curtis J T, scope F L1985 TN AEDC-TMR-85-V7) to obtain 47.02X 104W/m2But much smaller than the experimental value of 70.07X 104W/m2Since the heat flow calculation is sensitive to the incoming flow turbulence, the difference may be caused by the difference between the incoming flow turbulence and the calculated incoming flow state, which is not considered at present, and the specific numerical comparison is detailed in table 1. Table 3 compares the 2s time stagnation temperature and initial time heat flow values to the reference. Fig. 10 is a curve of the stationary point heat flow changing with time, and it can be seen that the stationary point heat flow gradually decreases and the decreasing trend gradually becomes gentle as time advances, because the heat is conducted in the circular tube structure, the stationary point temperature gradually increases, the thickness in the temperature boundary layer increases, the temperature gradient gradually decreases, and the heat flow density accordingly gradually decreases, and the heat flow density decreases by about 6.3% at 2s compared with the initial time, which is close to 8% of the heat flow decrease in 2s in the literature.
The mutual comparison of the calculation method and the coupling solution of the invention
The pneumatic heating/structure heat transfer coupling solving method (Neire, Liuwei Qiang 2012 'physical science report' 61184401 'hypersonic aircraft leading edge fluid-solid coupling calculation method research') is a calculation method for alternately iterating a flow field and a structure in a time domain. The method mainly comprises the following two processes: firstly, calculating a flow field by taking the object surface temperature of a boundary interface as a flow field boundary condition to obtain the heat flux density of the boundary interface; and secondly, calculating the temperature field of the structure by taking the heat flow density as a boundary condition of the structure.
The time iterative computation flow of the coupling solution method is shown in fig. 11, the coupling scheme is a serial iterative coupling scheme, and the scheme is based on that the characteristic time of the structural heat transfer is far longer than the characteristic time of the flow field and the heat transfer is a slow process, so that compared with the structural heat transfer computation, the flow field can be assumed to be transient stable, and the structure and the flow field only perform data exchange on a real time node.
For the above coupling iteration method, fig. 12 shows a calculation flow chart of the coupling solution method, in which data exchange is performed on the flow field and the structure once in each real physical time step, and the calculation flow is as follows:
(1) will tiTransferring the boundary temperature distribution Ti of the physical moment to a flow field grid;
(2) performing pseudo-time iteration in the flow field until convergence;
(3) calculating tiPhysical time flow field boundary heat flow interpolation Qi；
(4) Interpolating Q the Heat flowiTo the structural grid;
(5)ticarrying out pseudo-time iteration on transient heat transfer of the moment structure until convergence;
(6) judging whether the calculation requirement is met, if so, ending the program, otherwise, ti＝ti+ Δ t, return (1).
The calculation accuracy and grid sensitivity of the coupling algorithm and the calculation method are compared through an unsteady pneumatic heating calculation example with the heating time of 2s (a plurality of problems of the format and the grid effect in the heat flow CFD calculation are researched in Hayao, Yujiajun, Li Jun philosophy 2006 "aerodynamic reports" 24(1):125- > 130-. The time step is calculated as 0.001s, and the total time is calculated as 2 s. Two sets of different grids are respectively taken, and the flow field grid scale and the structural grid scale at the interface of the first set of grids are both 1 multiplied by 10-6m, the flow field grid dimension and the structural grid dimension at the second set of grid interface are both 5 multiplied by 10-5m。
FIG. 13a is a plot of stagnation temperature versus time at different grid scales, and FIG. 13b is a plot of temperature difference calculated over time at both algorithms. As can be seen from fig. 13a and 13b, the calculation results of the two algorithms are substantially coincident when the network is subdivided, and the stagnation temperature at 2s is about 388K, which indicates that the calculation results of the two algorithms are reasonable under the grid scale. When the network rough time division has large difference between the calculation results of the two algorithms, the stagnation temperature of the method for 2s is about 388K, and the calculation result is correct; the stagnation temperature in the coupling algorithm 2s is less than 360 degrees, and the calculation result is unreasonable. The calculation results of the two examples show that the dependence of the calculation method on the grid scale is smaller than that of a coupling algorithm, and the calculation method has higher grid adaptability. In fig. 13a and 13b, the stagnation temperature difference between the coupling algorithm and the calculation method of the present invention changes sharply at the initial stage, decreases rapidly with the passage of time, and finally tends to be constant. The stagnation point temperature value calculated by the calculation method of the invention is smaller than the stagnation point temperature value calculated by the coupling algorithm at the initial stage of heat transfer, and the value calculated by the calculation method of the invention is larger than the value calculated by the coupling algorithm at the later stage. The reason is that an interface stagnation point temperature value is assumed during the calculation of the coupling algorithm in the initial iteration, then the flow field is calculated to obtain interface heat flow, and the stagnation point heat flow value obtained at the moment is larger than the real heat flow value. In the calculation method of the invention, the stagnation point temperature value of the interface is calculated in real time, so that the obtained stagnation point heat flow value is closer to the real heat flow. Thus, the stagnation temperature value calculated by the calculation method of the invention at the initial stage of heat transfer is lower than the stagnation temperature value obtained by the coupling algorithm. In the following iteration, it is inevitable that the stagnation point heat flow value obtained by the coupling algorithm is smaller than the stagnation point heat flow value obtained by the calculation method of the present invention. Therefore, the stagnation temperature calculated by the calculation method is higher than the stagnation temperature value obtained by the coupling algorithm.
The main reason why the dependence of the calculation method on the grid scale is less than that of the coupling algorithm is that the temperature, the temperature gradient and the heat transfer coefficient on the interface of the flow field and the structure in the calculation method are obtained by the parameter interpolation of the flow field and the inner near object plane of the structure, and the boundary temperature and the heat flow calculation in the coupling algorithm are only related to the flow field or the structure.
While the invention has been described in terms of its preferred embodiments, such as that shown in fig. 14, it will be understood by those skilled in the art that various changes in form and detail may be made without departing from the spirit and scope of the invention.

Claims (2)

1. A hypersonic aircraft leading edge flow-heat-solid integrated calculation method is characterized by comprising the following steps:
taking the flow field and the structure as the same physical field, simultaneously calculating physical parameters and thermodynamic properties of the flow field and the structure, taking an interface of the flow field and the structure as an internal boundary of the whole physical field, simultaneously establishing a heat transfer control equation of the flow field and the structure, and simultaneously solving by adopting a unified numerical calculation method;
for structural heat transfer, the integral of the structural heat transfer control equation over the control volume Ω s is formed as follows, regardless of the heat source:
Figure FDA0002618510570000011
wherein dS is a control body surface unit, CsIs the specific heat capacity of a solid material, T is the structure temperature,
Figure FDA0002618510570000012
is the temperature gradient, rho is the material density, k is the thermal conductivity coefficient,
Figure FDA0002618510570000013
is a boundary; the flow field calculation adopts a compressible Reynolds average N-S equation, so that a flow field control equation and a structural heat transfer control equation are unified into a control equation in the same integral form on a control body omega:
Figure FDA0002618510570000014
wherein W is a constant, FcFor convection flux, FvIs the viscous flux; selecting an SST k-omega two-equation model from the turbulence model; during space dispersion, the convection flux is dispersed by adopting an AUSM + format based on a lattice format, and the viscous flux dispersion adopts a second-order central difference format;
the interface temperature calculation was calculated using the center-mean method as follows:
T＝(Tl+Tr)/2
in the formula, Tl、TrThe temperatures of the interface left and right control units are respectively; temperature gradient
Figure FDA0002618510570000015
The calculation of (2) needs to be corrected, and the calculation method is as follows:
Figure FDA0002618510570000016
Figure FDA0002618510570000017
in the formula (I), the compound is shown in the specification,
Figure FDA0002618510570000018
temperature gradients, L, of left and right control units of the interface, respectivelylrIs the distance between cell centers, rlrThe unit vector is from the center point of the left control unit to the center point of the right control unit; in addition, the calculation method of the temperature gradient employs a gaussian green method, which calculates by the value of the boundary of the cell and the normal vector:
Figure FDA0002618510570000019
for a lattice format, this is:
Figure FDA00026185105700000110
wherein Ω is the volume of the control body, i is the number of the cell, j is the number of the adjacent cell, and nijIs a unit normal vector, Δ SijFor control body surfaceAccumulating, wherein N is the unit surface number of the control body; the concept of thermal impedance in solid heat transfer is introduced, and the thermal impedance is defined as follows:
Figure FDA0002618510570000021
in the formula, d is the thickness in the heat flow direction, and A is the cross-sectional area perpendicular to the heat flow direction; the thermal impedance relation of the interface is obtained by introducing the concept of thermal impedance:
Rt,bnd＝Rt,l+Rt,r
the calculated equation of the interface thermal conductivity k is obtained as follows:
Figure FDA0002618510570000022
wherein k isl、krRespectively, the heat transfer coefficients of the left and right control units, Ll、LrThe distances from the center of the left and right control units to the center of the boundary, respectively.
2. The hypersonic aircraft leading edge flow-heat-solid integrated calculation method of claim 1, wherein the calculation method comprises the following boundary conditions:
(1) a flow field far field boundary condition which adopts a Riemann boundary condition;
(2) an object plane boundary condition comprising a solid surface flow boundary condition and a thermodynamic boundary condition;
a. the boundary condition of the solid surface flow interface meets the non-slip boundary condition, and the pressure gradient is zero;
b. structural heat conduction boundary conditions are thermodynamic boundary conditions, which include Dirichlet temperature boundary conditions and Neuman heat flow density boundary conditions.
CN201710951385.3A 2017-10-13 2017-10-13 Hypersonic aircraft leading edge flow-heat-solid integrated calculation method Expired - Fee Related CN107832494B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710951385.3A CN107832494B (en) 2017-10-13 2017-10-13 Hypersonic aircraft leading edge flow-heat-solid integrated calculation method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710951385.3A CN107832494B (en) 2017-10-13 2017-10-13 Hypersonic aircraft leading edge flow-heat-solid integrated calculation method

Publications (2)

Publication Number Publication Date
CN107832494A CN107832494A (en) 2018-03-23
CN107832494B true CN107832494B (en) 2021-02-19

Family

ID=61648117

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710951385.3A Expired - Fee Related CN107832494B (en) 2017-10-13 2017-10-13 Hypersonic aircraft leading edge flow-heat-solid integrated calculation method

Country Status (1)

Country Link
CN (1) CN107832494B (en)

Families Citing this family (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109885885B (en) * 2019-01-22 2023-01-06 南京航空航天大学 Method for estimating wall temperature of nozzle rod based on gas-solid-liquid three-phase coupling heat transfer
CN109918808B (en) * 2019-03-13 2023-10-20 北京宇航系统工程研究所 Three-field coupling simulation analysis method for gas-thermal bomb
CN110057536B (en) * 2019-04-12 2021-04-02 北京空天技术研究所 Method for simulating inner and outer flow coupling of air-breathing aircraft under engine combustion condition
CN110309541B (en) * 2019-05-31 2023-04-07 中国航天空气动力技术研究院 Method for constructing different-medium multi-component flow field interface conditions of variable-specific-heat gas
CN110321602B (en) * 2019-06-17 2020-12-11 大连理工大学 Full-field temperature calculation method for mine magnetic coupler
CN111339672B (en) * 2020-03-02 2021-06-08 上海索辰信息科技股份有限公司 Method for analyzing aerodynamic thermal simulation of shock wave at front edge of air inlet channel
CN111460578A (en) * 2020-03-23 2020-07-28 南京航空航天大学 High-precision flow-solid coupling calculation method for hypersonic aircraft nose cone thermal environment
CN111695219B (en) * 2020-06-12 2023-05-12 燕山大学 Stress prediction method of skin plate covered with thermal protection coating under supersonic flight condition
CN111859532B (en) * 2020-06-16 2023-11-28 空气动力学国家重点实验室 Improved hot wall correction method considering hypersonic chemical unbalance effect
CN112989485B (en) * 2021-01-18 2022-04-12 中国空气动力研究与发展中心计算空气动力研究所 Along trajectory heat flow interpolation method based on hot wall correction
CN113850008B (en) * 2021-12-02 2022-03-11 北京航空航天大学 Self-adaptive grid disturbance domain updating acceleration method for aircraft aerodynamic characteristic prediction
CN113899527B (en) * 2021-12-06 2022-03-01 中国空气动力研究与发展中心低速空气动力研究所 Method for correcting surface temperature of test model
CN114218833B (en) * 2021-12-16 2023-11-10 西北工业大学太仓长三角研究院 Method and system for predicting performance of flow field in secondary light gas gun
CN114329319B (en) * 2021-12-27 2024-05-28 北京航空航天大学 Stream thermosetting coupling calculation method based on physical neural network
CN115238397B (en) * 2022-09-15 2022-12-02 中国人民解放军国防科技大学 Method and device for calculating thermal environment of hypersonic aircraft and computer equipment
CN115544675B (en) * 2022-12-01 2023-04-07 北京航空航天大学 Multi-scale prediction method for surface catalytic properties of heat-proof material of hypersonic aircraft
CN116070554B (en) * 2023-04-06 2023-06-09 中国人民解放军国防科技大学 Hypersonic aircraft aerodynamic heat load calculation method, device and equipment
CN117933144B (en) * 2024-03-22 2024-05-28 中国空气动力研究与发展中心超高速空气动力研究所 Multiple grid method for solving grid flow field of complex topological structure

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7197442B2 (en) * 2002-08-30 2007-03-27 Fujitsu Limited Mesh creating device, mesh creating method and mesh creating program
CN104298826A (en) * 2014-10-10 2015-01-21 南京航空航天大学 Aerodynamic stability predicting and estimating method of aerial engine under counter thrust state
CN105667811A (en) * 2016-01-27 2016-06-15 南京航空航天大学 Design method for multi-stage coupling integrated structure of front body and air inflow channel of hypersonic aircraft
CN106682392A (en) * 2016-11-24 2017-05-17 南京航空航天大学 Technology for rapidly calculating ablation effect of complex hypersonic flight vehicle

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7197442B2 (en) * 2002-08-30 2007-03-27 Fujitsu Limited Mesh creating device, mesh creating method and mesh creating program
CN104298826A (en) * 2014-10-10 2015-01-21 南京航空航天大学 Aerodynamic stability predicting and estimating method of aerial engine under counter thrust state
CN105667811A (en) * 2016-01-27 2016-06-15 南京航空航天大学 Design method for multi-stage coupling integrated structure of front body and air inflow channel of hypersonic aircraft
CN106682392A (en) * 2016-11-24 2017-05-17 南京航空航天大学 Technology for rapidly calculating ablation effect of complex hypersonic flight vehicle

Also Published As

Publication number Publication date
CN107832494A (en) 2018-03-23

Similar Documents

Publication Publication Date Title
CN107832494B (en) Hypersonic aircraft leading edge flow-heat-solid integrated calculation method
Sun et al. Numerical investigation on the unsteady cavitation shedding dynamics over a hydrofoil in thermo-sensitive fluid
Yu et al. Improved treatment of the open boundary in the method of lattice boltzmann equation: general description of the method
CN107832260B (en) Numerical simulation method for flat plate impact jet flow heat transfer problem
Gao et al. A high-order lifting collocation penalty formulation for the Navier-Stokes equations on 2-D mixed grids
Rahbarimanesh et al. Development and validation of a homogeneous flow model for simulating cavitation in cryogenic fluids
Xu et al. A parallelized hybrid NS/DSMC-IP approach based on adaptive structured/unstructured overlapping grids for hypersonic transitional flows
Zhao et al. Numerical study of total temperature effect on hypersonic boundary layer transition
Liu et al. Investigation of heat transfer characteristics of high-altitude intercooler for piston aero-engine based on multi-scale coupling method
CN107766620A (en) A kind of Aerodynamic Heating structural optimization method based on reduced-order model
Yang et al. Numerical study of flow and heat transfer in a three-dimensional metal foam considering different direction micropores in skeleton structure
Kamali et al. Development and validation of a High-Fidelity Aero-Thermo-Elastic analysis capability
Tongqing et al. CFD/CSD-based flutter prediction method for experimental models in a transonic wind tunnel with porous wall
Bourgault-Cote et al. Two-dimensional/infinite swept wing ice accretion model
Jiang et al. Numerical simulation of three-dimensional high-speed flows using a second-order nonlinear model
König et al. Influence of non-uniform injection into a transpiration-cooled turbulent channel flow: a numerical study
Yang et al. Behaviors of hypersonic wing under aerodynamic heating
Nima et al. Numerical study of heat transfer enhancement for a flat plate solar collector by adding metal foam blocks
Vasconcellos et al. Flow past a circular cylinder: a comparison between commercial finite volume and finite element codes
Zope et al. Towards More Efficient Fluid-Thermal Interaction Analysis for Hypersonic Trajectory Flights
Remaki et al. Hermite-based mesh adaptation for functional outputs improvement in fluid flow simulation
Zhao et al. A fast natural convection algorithm based on dividing fluid development stages
Levy et al. The EZNSS CFD Code—Theoretical and User’s Manual
Qiu et al. Numerical investigation of geometric parameters effects on vortexes and aerodynamic heating environment in transverse gaps on hypersonic vehicle
Yuan et al. An efficient topology optimization method for steady gas flows in all flow regimes

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20210219

Termination date: 20211013