CN107808065B - 三维复杂外形高速飞行器流-固-热快速计算方法 - Google Patents
三维复杂外形高速飞行器流-固-热快速计算方法 Download PDFInfo
- Publication number
- CN107808065B CN107808065B CN201711177571.2A CN201711177571A CN107808065B CN 107808065 B CN107808065 B CN 107808065B CN 201711177571 A CN201711177571 A CN 201711177571A CN 107808065 B CN107808065 B CN 107808065B
- Authority
- CN
- China
- Prior art keywords
- heat
- calculation
- flow
- temperature
- aircraft
- 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
- 238000004364 calculation method Methods 0.000 title claims abstract description 101
- 238000012546 transfer Methods 0.000 claims abstract description 46
- 238000000034 method Methods 0.000 claims abstract description 32
- 238000010438 heat treatment Methods 0.000 claims abstract description 28
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 16
- 230000008878 coupling Effects 0.000 claims abstract description 15
- 238000010168 coupling process Methods 0.000 claims abstract description 15
- 238000005859 coupling reaction Methods 0.000 claims abstract description 15
- 230000004907 flux Effects 0.000 claims abstract description 9
- 230000003068 static effect Effects 0.000 claims abstract 2
- 239000010410 layer Substances 0.000 claims description 54
- 239000000463 material Substances 0.000 claims description 20
- 230000000694 effects Effects 0.000 claims description 19
- 238000006243 chemical reaction Methods 0.000 claims description 16
- 239000000126 substance Substances 0.000 claims description 11
- 229910052757 nitrogen Inorganic materials 0.000 claims description 8
- 125000004433 nitrogen atom Chemical group N* 0.000 claims description 8
- 229910052760 oxygen Inorganic materials 0.000 claims description 8
- 125000004430 oxygen atom Chemical group O* 0.000 claims description 8
- 230000008859 change Effects 0.000 claims description 6
- 238000005516 engineering process Methods 0.000 claims description 5
- 230000008569 process Effects 0.000 claims description 5
- IJGRMHOSHXDMSA-UHFFFAOYSA-N Atomic nitrogen Chemical compound N#N IJGRMHOSHXDMSA-UHFFFAOYSA-N 0.000 claims description 4
- 238000012407 engineering method Methods 0.000 claims description 3
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 2
- QVGXLLKOCUKJST-UHFFFAOYSA-N atomic oxygen Chemical compound [O] QVGXLLKOCUKJST-UHFFFAOYSA-N 0.000 claims description 2
- 239000003054 catalyst Substances 0.000 claims description 2
- 230000003197 catalytic effect Effects 0.000 claims description 2
- 238000012512 characterization method Methods 0.000 claims description 2
- 239000001301 oxygen Substances 0.000 claims description 2
- 239000004576 sand Substances 0.000 claims description 2
- 239000002344 surface layer Substances 0.000 claims description 2
- 230000005855 radiation Effects 0.000 claims 1
- 238000004088 simulation Methods 0.000 abstract description 9
- 230000007547 defect Effects 0.000 abstract description 5
- 101000795074 Homo sapiens Tryptase alpha/beta-1 Proteins 0.000 description 6
- 102100029639 Tryptase alpha/beta-1 Human genes 0.000 description 6
- VYPSYNLAJGMNEJ-UHFFFAOYSA-N Silicium dioxide Chemical compound O=[Si]=O VYPSYNLAJGMNEJ-UHFFFAOYSA-N 0.000 description 5
- 238000004458 analytical method Methods 0.000 description 5
- 238000011160 research Methods 0.000 description 5
- 101000662819 Physarum polycephalum Terpene synthase 1 Proteins 0.000 description 4
- 239000011241 protective layer Substances 0.000 description 4
- 230000007123 defense Effects 0.000 description 3
- 238000013461 design Methods 0.000 description 3
- 238000012545 processing Methods 0.000 description 3
- 241000152447 Hades Species 0.000 description 2
- 101000830822 Physarum polycephalum Terpene synthase 2 Proteins 0.000 description 2
- 125000001797 benzyl group Chemical group [H]C1=C([H])C([H])=C(C([H])=C1[H])C([H])([H])* 0.000 description 2
- 238000005094 computer simulation Methods 0.000 description 2
- 238000000354 decomposition reaction Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 239000007789 gas Substances 0.000 description 2
- 230000007774 longterm Effects 0.000 description 2
- 239000002184 metal Substances 0.000 description 2
- 229910052751 metal Inorganic materials 0.000 description 2
- 235000012239 silicon dioxide Nutrition 0.000 description 2
- 239000000377 silicon dioxide Substances 0.000 description 2
- 230000007704 transition Effects 0.000 description 2
- 241000237524 Mytilus Species 0.000 description 1
- 241000321453 Paranthias colonus Species 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000001816 cooling Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- 238000005755 formation reaction Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000009413 insulation Methods 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 239000002245 particle Substances 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 230000001681 protective effect Effects 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 230000000284 resting effect Effects 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 230000035939 shock Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000002076 thermal analysis method Methods 0.000 description 1
- 238000012549 training Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/10—Geometric CAD
- G06F30/15—Vehicle, aircraft or watercraft design
-
- 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
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/08—Thermal analysis or thermal optimisation
-
- 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)
- Other Investigation Or Analysis Of Materials By Electrical Means (AREA)
Abstract
本发明公开了一种三维复杂外形高速飞行器流‑固‑热快速计算方法，属飞行器气动计算领域。本发明依据边界层理论，采用CFD数值计算与工程算法相结合的方法，静气动热的计算简化为绕飞行器的无粘外流数值解和边界层内热流求解两份部分，同时耦合防热系统结构传热计算模型。本发明结合了全长数值模拟与工程算法各自的优势，可以快速高效对三维复杂外形高速飞行器的多种飞行状态下气动加热与结构耦合传热特性的计算方法，给出飞行器全机表面热流密度与防热结构温度场时变特性，弥补了直接数值模拟方法计算代价高、效率低、周期长等缺陷，同时拓展了工程算法的应用范围。
Description
技术领域
本发明属飞行器气动计算技术领域，具体指代一种三维复杂外形高速飞行器流-固-热快速计算方法。
背景技术
高超声速飞行器在大气层中持续飞行时，飞行器表面将承受巨大的气动热(王江峰,伍贻兆.高超声速复杂气动问题数值方法研究进展[J].航空学报.2015,33(1):159-175)。气动热会导致飞行器结构温度升高，改变结构内温度场分布，从而改变结构材料的物理属性，而结构材料特性的改变会引起结构刚度和结构模态的变化，对飞行器的飞行安全造成极大隐患。因此热防护问题对于高超声速飞行器设计显得尤为重要(程克明,吕英伟.飞行器持续气动加热的耦合性分析[J].南京航空航天大学学报,2000,32(2):150-155)，越来越受到世界各国的高度重视。
气动加热问题(卞荫贵,徐立功.气动热力学[M].中国科技大学出版社.1997:15-20)的主要研究方法包括：数值模拟方法(Sinha K,K.Reddy D S.Effect ofChemicalReaction Rates on Aeroheating Predictions of Reentry Flows[J].Journal ofThermophysics and Heat Transfer,2011,25(1):21-33)、工程近似计算方法(Hamilton HH,Weilmuenster K J,DeJarnette F R.Ap-proximate Method for ComputingConvective Heating on Hypersonic Vehicles Using Unstructured Grids[J].Journalofspacecraft and rockets.2014,1288:1305)、地面风洞实验(田旭昂,王成鹏,程克明.Ma5斜激波串动态特性实验研究[J].推进技术.2014,35(8):1030-1039)及自由飞实验等，其中后两种方法因代价高昂等因素不适于工程设计的初期与选型、改型阶段。
数值模拟方法主要是对Euler方程、N-S(Navier-Stokes)方程及相关简化形式的控制方程进行求解，具有计算精度高、可处理复杂流动和全机外形等优点，但在气动热与结构传热耦合问题的求解方面对计算资源需求巨大且非常耗时。国内外在求解高超声速气动热相关问题的数值计算方法方面做了大量研究，主要工作集中在计算效率与精度等方面。
气动加热工程计算方法(Zhao J S,Gu L X,Ma H Z.A rapid approach toconvective aeroheating prediction ofhypersonic vehicles[J].Science ChinaTechnologicalSciences,2013:1-15)对简单外形的气动热求解具有高效、准确的特点，因此在实际工程应用中率先得到发展。但在复杂气动外形的气动热分析方面适应性较差，并且需要基于大量实验数据对计算方法与结果进行人工修正。在这方面比较著名的是NASA兰利研究中心研发的一套气动加热预测程序(MINIVER)(李会萍.高超声速飞行器气动加热特性及其计算方法研究[D].上海：上海交通大学,2010)，该程序中驻点区域使用经典的Fay-Riddle公式，采用参考焓方法来计算高速流动压缩性效应，另外可对过渡流区气动热进行计算。Hamilton(Hamilton H H,Weilmuenster K J,DeJarnette F R.Approximate methodfor computing laminar and turbulent convective heating on hypersonic vehiclesusing unstructured grids[C].The 41st Annual AIAA Thermophysics Conference,AIAA Paper.2009,4310:2009.)针对三维钝头体发展了一种适用于空气平衡气体的高超声速流动热流密度计算方法，该方法可计算不同高速流动状态(层流、转捩、湍流)下的热流密度。Zoby(Zoby E V,Simmonds A L.Engineering flowfield method with angle-of-attack applications[J].Journal ofspacecraft androckets,1985,22(4):398-404)等人开发了LATCH方法，该方法基于参考焓和修正雷诺比拟计算热流密度，可用于有化学反应参与的钝头体气动热计算。李建林(李建林,唐乾刚,霍霖,程兴华.复杂外形高超声速飞行器气动热快速工程估算[J].国防科技大学学报.2012,34(6):89-93)等人对气动热工程计算方法进行拓展应用，对乘波体构型高速飞行器进行气动热快速估算，得到较好结果。总体而言，工程计算方法存在需要大量人工干涉及以庞大实验数据为基础等缺陷，在处理复杂外形飞行器方面的通用性仍需完善。
发明内容
针对于上述现有技术的不足，本发明的目的在于提供一种三维复杂外形高速飞行器流-固-热快速计算方法，以解决直接数值模拟方法计算代价高、效率低、周期长等缺陷，同时拓展了工程算法的应用范围。
为达到上述目的，本发明采用的技术方案如下：
一种三维复杂外形高速飞行器流-固-热快速计算方法，其特征在于，包括步骤如下：
步骤1、无粘外流场数值解；
步骤2、工程表面热流率计算；
步骤3、工程结构传热耦合计算；
步骤4、弹道状态动态插值方法；
步骤5、高温化学非平衡效应。
步骤1具体包括：
控制方程采用三维欧拉方程对无粘外流场求解，采用块结构网格机翼基于分布式并行计算技术，取无粘流场解结果中的物面参数作为边界层外缘参数提供给工程方法中的边界层简化算法。输出的物面参数包括：坐标(m)、物面速度分量(m/s)、马赫数、压强(Pa)、密度(kg/m3)、静温(K)、总能(J/m3)。而边界层的工程算法可以得到用于防热系统中结构传热计算所需的飞行器表面热流及传热系数等参数。
步骤2具体包括：
本发明针对的是全机外形，根据热流密度工程计算方法将飞行器表面热流的计算划分为驻点区和非驻点区两个区域。驻点区热流计算采用目前广泛使用的Fay-Riddell公式(中国人民解放军总装备部军事训练教材编辑工作委员会.高超声速气动热和热防护[M].北京:国防工业出版社,2003:46-136.)：
式中：ρw、μw、hw分别表示物面密度、物面粘性系数和物面焓，ρs、μs分别为驻点密度及驻点粘性系数，hD为平均空气电离焓，计算中取普朗特数Pr＝0.71、路易斯数Le＝1.0。
非驻点区的热流密度计算公式(卞荫贵,钟家康.高温边界层传热[M].科学出版社.1986.)，本发明采用的是工程中常用的平板传热模型。
步骤3具体包括：
热防护结构内表面采用绝热壁边界条件，依据防热结构材料的毕奥数Bi的大小，采不同的传热计算模型(航天工业部部标准,“战术导弹气动加热工程计算方法”,QJ1734-89，1989.)：
式中α为传热系数，δ为材料结构厚度，λδ为当前材料热传导系数。
当Bi＜0.1时，采用如下式所示的热薄壁传热模型：
初始条件为：
Tw|t＝0＝T0 (4)
式中ρδ、cδ、ε分别为表示分别表示防热材料的密度、比热容和表面辐射系数。采用差分方法对式(3)进行加热时间t推进求解，即可得到热薄壁条件下热防护层温度随时间的变化。
当Bi＞0.1时，采用热厚壁传热模型。将热厚壁由内而外分为j层进行逐层推进求解：
表层：
最内层：
中间层：
初始条件为：
Tj|t＝0＝T1|t＝0＝Tn|t＝0＝T0 (8)
式中下标m为为材料类型，n表示层数，λm为材料m的热传导系数。利用差分离散可计算求得热厚壁条件下热防护结构各层材料的温度随时间变化的结果。
气动加热与结构传热过程的耦合特性(吴洁,阎超.气动热与热响应的耦合研究[J].导弹与航天运载技术,2009(004):35-39.)以物面温度Tw作为物面边界输入条件，然后通过简化边界层工程方法计算得到表面热流密度qn；而在计算结构传热时，又以热流密度qn作为热边界条件，计算得到物面温度。在时间步tn的求解中，以上一时间步tn-1的物面温度作为本时间步的物面温度边界条件，采用工程算法进行边界层的气动加热计算，得出tn时间步的热流密度然后以作为热边界输入条件提供给结构传热计算模块，计算得到tn时间步的物面温度如此循环以实现结构传热和气动热环境的耦合计算。
步骤4具体包括：
计算非定常弹道状态的气动加热与结构传热耦合计算需要足够多的无粘外流的数值解作为其输入参数。
本发明利用无粘外流解动态插值方法提高耦合算法计算效率：通过已经预先完成的有限个无粘外流解的流场结果，插值得到当前弹道时间点上飞行条件下的流场解，以供气动加热与结构传热计算所需，具体步骤如下：
1)根据飞行器物面上的当地流场参数与来流参数的比值在不同飞行高度下几乎保持不变(Bova S W,HowardMA.Coupling strategies forhigh-speedaeroheatingproblems[R].SandiaNationalLaboratories,2011.)这一规律，对同一飞行器，已经某高度H0下的边界层外缘参数P0以及该高度下的流动参数Pinf0，同时已知任意高度下的大气参数Pinfx，可以通过Px＝Pinfx×P0/Pinf0的简单换算获得该飞行器在任意高度上的边界层外缘参数。
2)由1)，对于整个弹道状态，无粘外流解只需计算设定参考飞行高度下的两个马赫数和两个迎角，即4个飞行状态，其他飞行状态则可根据插值方法得到。这四个计算状态分别为设定飞行高度下的两个马赫数和两个迎角的组合，即：(Ma1,α1),(Ma1,α2),(Ma2,α1),(Ma2,α2)。与这四个状态所对应的飞行器表面任一点处的流动参数的值分别记作P11、P12、P21、P22，则某个弹道时间点上飞行条件(Max,αx)状态下的飞行器表面同一位置处的流动参数Pxx可通过如下插值算法得到：
步骤5具体包括：
在实际高超声速飞行条件下，空气的高温非平衡效应对气动热影响显著(欧阳水吾,谢中强.高温非平衡空气扰流[M].国防工业出版社,2003)。高温空气非平衡边界层驻点传热与平衡边界层驻点传热表面热通量的关系式表示为：
式中，下标O、N分别表示氧原子和氮原子，φ为壁面催化因子；h0为生成焓；CO.s和CN.s分别表示氧原子和氮原子的质量浓度；Le为路易斯数，取值范围1～2。
对于组元参数特性，使用组元参数简化计算模型，在流场温度9000K以下时，只考虑O2、O、O+、N2、N、N+、NO、NO+、e-等9种组元及以下6个化学反应方程式：
其中Ki为摩尔密度平衡常数。本发明暂不考虑NO、NO+、e-、N+、O+。
氧原子和氮原子的壁面催化因子有如下关系成立：
其中，[O2]0，[N2]0分别表示为空气中氧分子和氮分子的摩尔密度。联立式(11)、(12)与反应方程①、②求解，得到氧原子和氮原子的摩尔密度：
取一级近似φO＝φN＝0。联立反应方程①②③，可得到[O2]、[N2]和[NO]的摩尔密度计算式为：
联立反应方程④⑤⑥和电荷守恒方程(唐国庆.高超声速多体分离动态气动加热计算技术[D].南京：南京航空航天大学,2013.)，可以得到电子的摩尔密度：
考虑高温化学非平衡效应下的气动加热计算方法可修正为：首先由Fay-Riddell驻点热流密度计算公式(1)得到平衡条件下热流密度qeq，然后根据化学反应计算模型得到各组元摩尔密度和质量分数，最后由驻点化学非平衡边界层与平衡边界层热流密度关系式(10)求得高温非平衡条件下的热流密度qne
本发明的有益效果：
本发明依据边界层理论，采用CFD(Computational Fluid Dynamics)数值计算与工程算法相结合的方法，将气动热的计算简化为绕飞行器的无粘外流数值解和边界层内热流求解两个部分，同时耦合防热系统结构传热计算模型，发展了可用于快速计算三维复杂外形高速飞行器流-固-热特性的计算方法，实现了复杂飞行条件下飞行器全机表面热流密度与防热结构温度场时变特性的计算。
该方法的优点在于：综合了全流场数值模拟与工程算法各自的优势，可以快速高效对三维复杂外形高速飞行器的多种飞行状态下气动加热与结构耦合传热特性进行计算与分析，给出热流密度与温度等参数分布，弥补了直接数值模拟方法计算代价高、效率低、周期长等缺陷，同时拓展了工程算法的应用范围。
附图说明
图1a为沿子午线的组元摩尔浓度分布(t＝1000s)。
图1b沿子午线表面温度分布(t＝1000s)。
图1c为表面温度：不考虑化学非平衡效应(t＝1000s)。
图1d为表面温度：考虑化学非平衡效应(t＝1000s)。
图2为计算模型及物面网格和TPS方案。
图3a为第1000s时得飞行器防热层外表面温度分布。
图3b为第1000s时得飞行器防热层内表面温度分布。
图3c为飞行器在计算初始时刻(t＝0s)的防护层外表面热流密度分布。
图3d为飞行器在计算结束时刻(t＝1000s)的防护层外表面热流密度分布。
图3e为驻点内/外层温度随时间的变化曲线。
图3f为驻点热流密度随时间的变化曲线。
图4为弹道状态动态插值方法计算流程图。
图5a为弹道状态t＝0s时的表面温度分布。
图5b为弹道状态t＝100s时的表面温度分布。
图5c为弹道状态t＝350s时的表面温度分布。
图5d为弹道状态t＝600s时的表面温度分布。
图5e为弹道状态t＝1000s时的表面温度分布。
图6为本发明计算方法总体流程示意图。
具体实施方式
为了便于本领域技术人员的理解，下面结合实施例与附图对本发明作进一步的说明，实施方式提及的内容并非对本发明的限定。
1、RAMC-II巡航状态长时加热(高温非平衡效应)
1.1技术参数
计算外形采用典型高速飞行器RAMC-II球锥体试验模型，其几何参数为：头部曲率半径0.152m，半顶角9°，模型长度1.3m。无粘外流场解的计算网格为总单元数约40万、物面单元约1.5万的块结构网格。设定的巡航状态为：Ma∞＝12，迎角α＝0°，飞行高度H＝33km，飞行时间t＝1000s。
长时加热中考虑结构传热，飞行器热防护结构材料为金属Ti，厚度2mm，表面发射率为0.8。结构传热计算中，时间步长取为0.05s(即总时间推进步数为2万步)，采用热薄壁传热模型.该算例在普通微机上计算耗时约17min CPU机时。
2.2数值模拟结果
图1给出了该算例数值计算结果。图1a为九组元模型O、N、O2、N2、NO、e-、N+、O+、NO+在长时加热第1000秒时刻沿子午线分布的组元摩尔浓度分布。由于在本算例计算状态下空气电离度很低，因此计算中未考虑电离反应对气动热计算的影响。因此，图中可以看出第1000秒时，带电粒子摩尔浓度都基本为零，沿子午线组元成分主要以O、N、N2为主，O2基本分解完毕，N2部分分解，计算中O2与N2的分解特性与空气高温非平衡流理论相符。图1b为子午线上考虑空气高温非平衡效应对比完全气体的表面温度分布，由图可见两者的温度差异明显，在驻点处化学非平衡效应导致了约300K的温降，这与理论分析趋势一致。图1c、1d分别给出了第1000秒时不考虑化学反应效应和考虑化学反应效应时的表面温度分布云图。在本算例条件下多组元的高温化学反应效应较明显，化学非平衡效应使得壁面热流密度降低，飞行器表面温度下降，这与理论分析相符。该算例表明本文所发展的计算方法可以用于高温化学非平衡流的气动加热问题。
2、X-37B巡航状态及弹道状态长时传热
2.1、技术参数
本发明的计算模型选取类X-37B外形，如图2所示，图2同时给出了计算所用的表面网格。计算状态如表1所示。
表1
用于无粘外流求解的块结构网格总单元数约512万，物面网格单元约25.1万。在热防护结构方面，该算例进行了更接近于工程设计的复杂方案，即飞行器头部区域设定为热防护区1，以TPS1表示，全机其他部位设为热防护区2，以TPS2表示(见图2)，不同的热防护区域使用表2给出的不同的热防护方案，其中热防护结构材料的最外层(飞行器表面)和最内层温度初值设为飞行高度上的大气环境温度245K。根据该算例的飞行条件，不考虑高温非平衡效应影响特性。该算例耗时约28min CPU机时。
表2
2.2、数值模拟结果
图3a和3b分别为第1000秒时飞行器防热层外表面(TWsurf)和内表面(TWin)温度分布，飞行器在计算初始时刻(第0s)和计算结束时(第1000s)的防护层外表面热流密度分布如图3c、3d所示。由计算结果可见，巡航1000秒时，类X-37B热防护层内表面温度在不同的热防护区域(TPS)出现明显差异(图3b)，而热防护层整个外表面温度分布均匀。这是由于TPS1区使用的是热薄壁与隔热层的组合防热结构，并且隔热层采用热沉性好的二氧化硅材料(SiO2)，故在该热防护区域内表面温度出现显著降低，大部分区域温度在410K左右，说明防热层隔热效果良好，与工程实际情况相符。而TPS2区仅采用了热薄壁结构，防热材料为0.002m的金属Ti，其热传导性很好，故防护层内外表面温度很快达到平衡，与外表面温度差别很小，且防护区域热流密度很小，该区域温度大部分在600K以上。图3e、3f分别为驻点内外层温度与热流密度随时间变化曲线，可以发现TPS1内外表面温度变化差别较大，计算结果显示驻点区域防热层外表面温度为最高温约为979K，内表面最高温度约为521K，造成此温度值差异的原因仍然是由于在这个区域使用了SiO2隔热层。
2.3、弹道状态长时加热
在以上算例的基础上，考虑弹道状态，即变高度、迎角、马赫数，弹道状态更加接近飞行器实际飞行包线状态。计算流程如图4所示。
本文根据相关文献设定的弹道参数如表3所示，弹道参数考虑了马赫数、飞行高速及迎角的变化，弹道飞行时间1000s，在此期间飞行高度由60km降至30km，飞行马赫数由Ma＝6降至Ma＝4，迎角变化范围为±5°。该算例全机表面采用同一种组合热防护方案模型，具体参数上述算例中TPS1区相同。该算例耗时约31min CPU机时(根据弹道飞行参数，不考虑高温非平衡效应)。
表3
图5a～图5e给出了该算例计算结果，从图中可以清晰的看到飞行器表面沿弹道飞行时逐渐加热到平衡状态的过程，并且随迎角从-5度至5度变化过程中飞行器头部与机翼前缘处的高温区往下表面移动，第1000秒时，最高温出现在飞行器头部下表面，约为1050K，这与高超声速飞行器防热结构材料特性理论分析相符。由于没有在公开发表文献中找到与此类似的算例进行对比，因此这里仅给出本文算例计算结果的分析。该算例表明本文所发展的计算方法可对复杂高超声速飞行器在弹道飞行状态下的气动加热与结构传热耦合进行快速高效计算与分析，可为高速飞行器的热防护设计与选型提高技术参考，其中，图6为本发明计算方法总体流程示意图。
本发明具体应用途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以作出若干改进，这些改进也应视为本发明的保护范围。
Claims (1)
1.一种三维复杂外形高速飞行器流-固-热快速计算方法，其特征在于，包括步骤如下：
步骤1、无粘外流场数值解；
步骤2、工程表面热流率计算；
步骤3、工程结构传热耦合计算；
步骤4、弹道状态动态插值方法；
步骤5、高温化学非平衡效应；
步骤1具体包括：
控制方程采用三维欧拉方程对无粘外流场求解，采用块结构网格机翼基于分布式并行计算技术，取无粘流场解结果中的物面参数作为边界层外缘参数提供给工程方法中的边界层简化算法，边界层简化算法用于快速计算飞行器表面热流及传热系数，输出的物面参数包括：坐标(m)、物面速度分量(m/s)、马赫数、压强(Pa)、密度(kg/m3)、静温(K)、总能(J/m3)；
步骤2具体包括：
三维复杂外形高速飞行器流-固-热快速计算方法针对的是全机外形，根据热流密度工程计算方法将飞行器表面热流的计算划分为驻点区和非驻点区两个区域，驻点区热流计算采用目前广泛使用的Fay-Riddell公式：
式中：ρw、μw、hw分别表示物面密度、物面粘性系数和物面焓，ρs、μs分别为驻点密度及驻点粘性系数，hD为平均空气电离焓，计算中取普朗特数
Pr＝0.71、路易斯数Le＝1.0；
非驻点区的热流密度计算公式采用的是平板传热模型；
步骤3具体包括：
热防护结构内表面采用绝热壁边界条件，依据防热结构材料的毕奥数Bi的大小，采不同的传热计算模型：
式中α为传热系数，δ为材料结构厚度，λδ为当前材料热传导系数；
当Bi＜0.1时，采用如下式所示的热薄壁传热模型：
初始条件为：
Tw|t＝0＝T0 (4)
式中ρδ、cδ、ε分别为表示分别表示防热材料的密度、比热容和表面辐射系数；采用差分方法对式(3)进行加热时间t推进求解，得到热薄壁条件下热防护层温度随时间的变化，当Bi＞0.1时，采用热厚壁传热模型；将热厚壁由内而外分为j层进行逐层推进求解：
表层：
最内层：
中间层：
初始条件为：
Tj|t＝0＝T1|t＝0＝Tn|t＝0＝T0 (8)
式中下标m为材料类型，n表示层数，λm为材料m的热传导系数，利用差分离散计算求得热厚壁条件下热防护结构各层材料的温度随时间变化的结果；
气动加热与结构传热过程的耦合特性：以物面温度Tw作为物面边界输入条件，然后通过公式(1)中热流计算公式计算得到表面热流密度qn；而在计算结构传热时，又以热流密度qn作为热边界条件，计算得到物面温度；在时间步tn的求解中，以上一时间步tn-1的物面温度作为本时间步的物面温度边界条件，采用公式(1)中热流计算公式进行边界层的气动加热计算，得出tn时间步的热流密度然后以作为热边界输入条件提供给结构传热计算模块，计算得到tn时间步的物面温度如此循环以实现结构传热和气动热环境的耦合计算；
步骤4具体包括：
三维复杂外形高速飞行器流-固-热快速计算方法利用无粘外流解动态插值方法提高耦合算法计算效率：通过已经预先完成的有限个无粘外流解的流场结果，插值得到当前弹道时间点上飞行条件下的流场解，以供气动加热与结构传热计算所需；
1)根据飞行器物面上的当地流场参数与来流参数的比值在不同飞行高度下几乎保持不变这一规律，对同一飞行器，已知某高度H0下的边界层外缘参数P0以及该高度下的流动参数Pinf0，同时已知任意高度下的大气参数Pinfx，通过Px＝Pinfx×P0/Pinf0的简单换算获得该飞行器在任意高度上的边界层外缘参数；
2)由1)，对于整个弹道状态，无粘外流解只需计算设定参考飞行高度下的两个马赫数和两个迎角，即4个计算状态，其他飞行状态则根据插值方法得到；这四个计算状态分别为设定飞行高度下的两个马赫数和两个迎角的组合，即：(Ma1,α1),(Ma1,α2),(Ma2,α1),(Ma2,α2)；与这四个状态所对应的飞行器表面任一点处的流动参数的值分别记作P11、P12、P21、P22，则某个弹道时间点上飞行条件(Max,αx)状态下的飞行器表面同一位置处的流动参数Pxx通过如下插值算法得到；
步骤5具体包括：
在实际高超声速飞行条件下，空气的高温非平衡效应对气动热影响显著，高温空气非平衡边界层驻点传热与平衡边界层驻点传热表面热通量的关系式表示为：
式中，下标O、N分别表示氧原子和氮原子，φ为壁面催化因子；h0为生成焓；CO,s和CN,s分别表示氧原子和氮原子的质量浓度；Le为路易斯数，取值范围1～2；对于组元参数特性，使用组元参数简化计算模型，在流场温度9000K以下时，只考虑O2、O、O+、N2、N、N+、NO、NO+、e-9种组元及以下6个化学反应方程式：
①②
③④
⑤⑥
其中Ki为摩尔密度平衡常数；暂不考虑NO、NO+、e-、N+、O+；
氧原子和氮原子的壁面催化因子有如下关系成立：
其中，[O2]0，[N2]0分别表示为空气中氧分子和氮分子的摩尔密度；联立式(11)、(12)与反应方程①、②求解，得到氧原子和氮原子的摩尔密度：
取一级近似φO＝φN＝0；联立反应方程①②③，得到[O2]、[N2]和[NO]的摩尔密度计算式为：
联立反应方程④⑤⑥和电荷守恒方程，得到电子的摩尔密度：
考虑高温化学非平衡效应下的气动加热计算方法修正为：首先由Fay-Riddell驻点热流密度计算公式(1)得到平衡条件下热流密度qeq，然后根据化学反应计算模型得到各组元摩尔密度和质量分数，最后由驻点化学非平衡边界层与平衡边界层热流密度关系式(10)求得高温非平衡条件下的热流密度qne。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711177571.2A CN107808065B (zh) | 2017-11-23 | 2017-11-23 | 三维复杂外形高速飞行器流-固-热快速计算方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711177571.2A CN107808065B (zh) | 2017-11-23 | 2017-11-23 | 三维复杂外形高速飞行器流-固-热快速计算方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107808065A CN107808065A (zh) | 2018-03-16 |
| CN107808065B true CN107808065B (zh) | 2019-12-31 |
Family
ID=61581312
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201711177571.2A Active CN107808065B (zh) | 2017-11-23 | 2017-11-23 | 三维复杂外形高速飞行器流-固-热快速计算方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107808065B (zh) |
Families Citing this family (18)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108629089B (zh) * | 2018-04-17 | 2021-08-03 | 西南交通大学 | 一种基于连续-离散单元耦合的强夯加固地基模拟方法 |
| CN109029907B (zh) * | 2018-07-18 | 2020-07-14 | 大连理工大学 | 一种气动热环境试验模拟条件的参数相似方法 |
| CN109583067B (zh) * | 2018-11-22 | 2020-04-07 | 北京空天技术研究所 | 基于温度平衡的高速飞行器转捩位置测量传感器设计方法 |
| CN109726432B (zh) * | 2018-11-26 | 2023-05-16 | 北京空天技术研究所 | 飞行器底部结构温度计算方法 |
| CN109918808B (zh) * | 2019-03-13 | 2023-10-20 | 北京宇航系统工程研究所 | 一种气热弹三场耦合仿真分析方法 |
| CN110119535B (zh) * | 2019-04-11 | 2023-05-26 | 上海交通大学 | 一种轴对称气固耦合传热模型、分析方法及应用系统 |
| CN111859532B (zh) * | 2020-06-16 | 2023-11-28 | 空气动力学国家重点实验室 | 考虑高超声速化学非平衡效应的改进热壁修正方法 |
| CN111780949B (zh) * | 2020-07-10 | 2021-04-30 | 南京航空航天大学 | 基于cfd分析的高速进气道前体风洞实验总压修正方法 |
| CN111931295B (zh) * | 2020-09-13 | 2024-08-20 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种全弹道整体迭代的气动热/传热耦合计算方法 |
| CN112329148A (zh) * | 2020-11-18 | 2021-02-05 | 西北工业大学 | 一种高超声速飞行器热防护结构优化方法及系统 |
| CN113065275B (zh) * | 2021-03-05 | 2022-05-24 | 浙江大学 | 一种飞行器飞行过程中驻点热流的在线预示方法 |
| CN114186330B (zh) * | 2021-11-04 | 2024-08-02 | 北京机电工程研究所 | 高速流场环境下轴类结构温度估算方法 |
| CN114021499B (zh) * | 2021-11-05 | 2024-06-11 | 沈阳飞机设计研究所扬州协同创新研究院有限公司 | 基于fvm-tlbfs方法的飞行器热防护结构热传导计算方法 |
| CN113867381B (zh) * | 2021-12-02 | 2022-02-22 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种飞行器姿态控制方法 |
| CN114626313B (zh) * | 2022-03-04 | 2023-05-16 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种可解析时变热响应的高速气动热cfd求解方法 |
| CN116341397B (zh) * | 2023-05-30 | 2023-09-05 | 湖南光华防务科技集团有限公司 | 一种基于深度神经网络的灭火弹温度场模拟方法 |
| CN118228521B (zh) * | 2024-05-24 | 2024-08-13 | 西安现代控制技术研究所 | 一种适应气动热耦合内热的固体发动机复合壳体设计方法 |
| CN118313069B (zh) * | 2024-06-11 | 2024-08-09 | 中国空气动力研究与发展中心计算空气动力研究所 | 飞行器流场与多壁面效应耦合模拟方法、装置、设备及介质 |
Family Cites Families (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103077259A (zh) * | 2011-10-26 | 2013-05-01 | 上海机电工程研究所 | 高超声速导弹多场耦合动力学一体化仿真分析方法 |
| WO2013078628A1 (zh) * | 2011-11-30 | 2013-06-06 | 天津空中代码工程应用软件开发有限公司 | 直升机旋翼飞行结冰的数值模拟方法 |
| CN103366052B (zh) * | 2013-06-27 | 2016-06-01 | 中国航天空气动力技术研究院 | 一种高超声速飞行器热静气弹分析方法 |
| CN106528990B (zh) * | 2016-10-27 | 2019-08-09 | 中国运载火箭技术研究院 | 一种基于泛函优化的高超声速尖锥外形热流密度建模方法 |
| CN106682392B (zh) * | 2016-11-24 | 2019-07-19 | 南京航空航天大学 | 复杂高超声速飞行器烧蚀效应快速计算方法 |
-
2017
- 2017-11-23 CN CN201711177571.2A patent/CN107808065B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN107808065A (zh) | 2018-03-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107808065B (zh) | 三维复杂外形高速飞行器流-固-热快速计算方法 | |
| CN106682392B (zh) | 复杂高超声速飞行器烧蚀效应快速计算方法 | |
| CN107944137B (zh) | 高超声速飞行器弹道状态多场耦合的热气动弹性计算技术 | |
| Xiao et al. | Numerical investigation of transpiration cooling for porous nose cone with liquid coolant | |
| Yao et al. | A heat transient model for the thermal behavior prediction of stratospheric airships | |
| Lu et al. | Investigation of thermal protection system by forward-facing cavity and opposing jet combinatorial configuration | |
| Lv et al. | A new design method of single expansion ramp nozzles under geometric constraints for scramjets | |
| CN109800488B (zh) | 关于液体火箭高空环境下底部热环境的数值计算方法 | |
| Gu et al. | Dynamic experimental investigations of a bypass dual throat nozzle | |
| Eggers | One-dimensional flows of an imperfect diatomic gas | |
| Wei | The development and application of CFD technology in mechanical engineering | |
| Khraibut et al. | Laminar hypersonic leading edge separation–a numerical study | |
| Lei et al. | Experimental model design and preliminary numerical verification of fluid–thermal–structural coupling problem | |
| CN106508020B (zh) | 一种可用于工程设计的复杂飞行器气动加热计算方法 | |
| Smith | Effect of gas radiation in the boundary layer on aerodynamic heat transfer | |
| CN107766620A (zh) | 一种基于降阶模型的气动‑热‑结构优化方法 | |
| Fan et al. | Prediction of hypersonic boundary layer transition with variable specific heat on plane flow | |
| Su et al. | Fluid–Thermal–Structure Interaction of Hypersonic Inlets Under Different Aspect Ratios | |
| Huang et al. | Coupled Fluid–Structural Thermal Numerical Methods for Thermal Protection System | |
| Hu et al. | Shock wave standoff distance of near space hypersonic vehicles | |
| Meng et al. | A hypersonic aeroheating calculation method based on inviscid outer edge of boundary layer parameters | |
| Li et al. | An engineering method of aerothermodynamic environments prediction for complex reentry configurations | |
| CN114021499A (zh) | 基于fvm-tlbfs方法的飞行器热防护结构热传导计算方法 | |
| Willems et al. | Coupled Fluid-Thermal Response in the Gap Region of a High-Speed Control Surface | |
| Kumar | Numerical simulation of chemically reactive hypersonic flows |
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 |