载人航天器舱内流场与温度场松耦合计算方法研究

李志杰 果琳丽

（中国空间技术研究院载人航天总体部，北京 100094）

0 引言

与非载人航天器相比，载人航天器均含有密封舱，而密封舱内会有供人生存的大气存在[1]。因此在对舱内热环境进行预测分析计算时，不可避免的需要同时对流场与温度场进行耦合计算[2]。国内外学者针对此问题的紧耦合计算方法已经展开了很多研究工作[3]。文献[4]中Garimella等人完成了舱内流场与传热一体化耦合计算，得到了较为理想的结果。文献[5]中Hofacker等人利用PHOENICS软件模拟了欧空局赫尔墨斯（Hermes）航天飞机舱内对流换热温度场、流场和狭缝集中送风的通风效果。文献[6]对航天器密封舱内流动和传热问题进行了数值计算研究，得出了舱内对流换热系数与平均风速的经验关系式。文献[7]利用CFD数值计算方法研究了载人空间站舱内对流换热的问题。但上述工作中均采取的是CFD紧耦合数值计算方法[8-10]，即计算时采用流场的特征时间为计算步长，特点是流场与结构传热实时耦合，反映了物理实际，但计算需要耗费大量时间，在实际工程中的应用价值大打折扣，因此需要寻求一种快速便捷的仿真计算方法解决以上问题。

本文首先讨论了数值计算中的耦合传热问题，分析了流固界面耦合求解原理，指出了不同求解方法的优缺点。针对普遍使用的紧耦合方法计算时间过长的问题，结合一般松耦合求解思想提出了一种针对载人航天器舱内流场/温度场松耦合计算方法，并对其进行了详细的数值仿真验证。结果证明在流场随时间变化不大的情况下，本文提出的流场/温度场松耦合计算方法在求解瞬态耦合换热问题时不但具有很好的鲁棒性，而且在与紧耦合方法的计算结果进行对比时，精度上也只存在较小的误差。

1 耦合处理方法原理分析

紧耦合计算的特点是流场与固体计算域在接触面实时交互数据、耦合迭代，主要体现在流体与固体接触面上的边界条件处理，即各自的温度边界条件[11-12]。流/固接触面的耦合原理，如图1所示。

图1 流/固接触面的耦合原理Fig.1 Coupling principle of fluid/solid interface

图1 中，在直角坐标系下，流场区域和固体结构区域分别用FΩ和SΩ表示；Γ为流、固接触面；sq为接触面Γ传至结构内部的热流；xq和yq分别为沿x和y方向由流场传入至结构内部的热流。

热传导方程在边界Γ上的定解条件是[10]：

式中 k为固体热传导系数；T为温度；x为横向坐标；y为纵向坐标；nx为沿x方向的单位矢量；ny为沿y方向的单位矢量。

从能量守恒出发，在忽略辐射效应的假设前提下，流体区域传至固壁的热流应等于固壁传至结构体内部的热流，因此有[10]：

由上式建立起流场与结构传热计算的耦合关系，先给定时间为0、流场FΩ的初场和Γ上的温度分布，通过求解Navier-Stokes方程[13]得到tΔ时间后壁面Γ上的热流分布，此处tΔ为一个计算时间步长。然后以该热流分布为边界条件，用有限单元法计算tΔ时间后的结构温度分布，并得到壁面Γ上新的温度分布，至此完成耦合计算的一个周期。如此反复迭代，直到计算要求的截止时刻。

而在松耦合方法中，假定从t=0时刻开始，Γ上的热流分布为恒定值，即上述紧耦合计算的0时刻热流分布，结构区域SΩ从此作为边界条件进行热传导分析。利用松耦合方法求解流场与结构传热耦合的非定常问题，可按“准定常”方法处理。首先，把整个时间历程划分成一定的时间序列0t，1t，2t，…，1nt-，nt时间间隔tΔ的大小可视精度要求进行控制。在每一时间间隔tΔ内，认为固体表面接收到一个固定的外表面热流，具体步骤如下[14-15]：

1）给定初始t0时刻外表面壁温和结构内温度分布初场，算出稳态流场，得到外表面温度分布第一个时间步Δt1内的外表面壁温；

2）将外表面温度分布作为边界条件，求解得到t1时刻结构温度分布；

3）以t1时刻固体边界温度为边界条件，求得t1时刻的稳态流场；

4）冻结t1时刻的流场，继续在t1～t2时间段内单独求解结构瞬态热传导，获得t2时刻的结构温度分布；

5）转向步骤3），继续计算直到整个时间推算完毕转向执行步骤2），直至整个时间推算完毕。

紧耦合方法中，虽然流场与结构传热实时耦合，反映了物理实际，但其计算效率较为低下，不满足实际应用中的高效性要求。而针对此类耦合问题，如果计算对象的流场随时间变化较为缓慢，即流场流动特征时间与结构传热特征时间相差较大，松耦合方法在结构体每一时间步长的时间间隔内可以认为流场是稳态，从而大大减少计算量。

2 流场/温度场松耦合求解思想

在前面所介绍的松耦合求解方法中，在计算固体导热的每一个时间步长内都将流场视为稳态，从而大大节省了整体的计算时间提高了效率。由此受到启发，提出一种新的流场/温度场松耦合求解方法，将其应用到载人航天器舱内对流换热的CFD数值计算当中，具体求解方法与步骤如下：

1）在一个时间步长1tΔ内同时计算流场+温度场，也就是说在1tΔ的时间步长内对能量方程和动量方程同时进行迭代求解；

2）固定流场，在下一个时间步长2tΔ内不再进行计算流场变化，即不再对流体动量方程进行计算，只对能量方程进行迭代求解一个步长2tΔ，得出温度场；

3）固定温度场不变，在1tΔ时间内继续同时求解流场；

4）重复步骤1）到步骤3），直到求解结束。

需要说明的是，一般的松耦合计算方法中是将1tΔ与2tΔ作为进行迭代计算的步长，在2tΔ的计算步长内将流场视为瞬态固定，而在本文中，1tΔ与2tΔ将作为进行计算的时间步长，而在每一个时间步长内会进行若干次相应的迭代计算，进一步提高整体计算效率。

在流场/温度场松耦合计算方法中，1tΔ与2tΔ时间步长的选取方式至关重要。根据非稳态导热问题数值解法原理[10]，时间计算步长越小，所需计算机内存及计算时间就大大增加；时间步长选取过大，计算精度则相应降低，综合考虑以上因素，本文选取9种1tΔ与2tΔ的组合，如表1所示。

表1 不同时间步长组合的选取Tab.1 Different combination of time steps

3 计算模型的选取

本文选取的计算对象为某升力体载人航天器设备舱的舱段模型，如图2所示。为方便进行数值计算，提高计算效率，对模型进行相应简化，简化后的模型及舱内7个设备的编号如图3所示，模型尺寸单位为m。

图2 计算模型示意Fig.2 Numerical simulation model

图3 简化后的计算模型及舱内设备编号Fig.3 Simplified model and serial numbers of devices in cabin

数值计算的边界条件设定如下：舱内设备材料为铝，舱蒙皮内壁面（弧面）给定温边界条件100℃，其它壁面设为绝热条件，CFD数值计算数学模型采取k-ε湍流计算模型[10]。

4 计算结果与分析

首先考察不同时间步长组合下基于本文松耦合求解思想的计算结果，与紧耦合计算结果对比，选取一组计算效率和计算精度均较优的时间步长组合，然后对此时间步长组合下不同时刻的计算结果进行对比分析。

4.1 不同时间步长下的计算结果

表2所示为采用本文松耦合方法在9种时间步长组合下进行瞬态计算1000s所消耗的计算时间，并与紧耦合方法的计算结果进行了对比。表中t为计算消耗时间，单位为h。

表2 不同时间步长组合下计算1000s所消耗的时间Tab.2 Time expended when calculating 1 000s at different time steps combination

从表2中可以看出，相比紧耦合计算方法，流场/温度场松耦合计算方法在进行瞬态计算时节省了大量的计算时间。结合表1中时间步长的选取可以看出，1tΔ越大，整体的计算时间就会越长，相反2tΔ的增大会减少计算时间，提高计算效率。也就是说，实际计算时间基本是与1tΔ/2tΔ的大小成正比关系。

很明显，不同时间步长的组合对数值计算所用时间的影响很大，但在实际计算中，只一味地选择计算耗时少的情况也是不可取的，因为也要同时考察不同计算方式下计算结果的精度问题。

取舱内设备3～4的温度分布结果作为分析对象，在1000s时，表1中所示的9种不同时间步长组合下的计算结果如图4～6所示。图中T为温度；横坐标为沿舱段横向的距离，参照图3所示。

图4 计算时间1000s时算例1～3的温度计算结果Fig.4 Temperature results of case 1～3 in 1 000s

图5 计算时间1000s时算例4～6的温度计算结果Fig.5 Temperature results of case 4～6 in 1 000s

图6 计算时间1000s时算例7～9的温度计算结果Fig.6 Temperature results of case 7～9 in 1 000s

从图4中不难发现，当流场与温度场同时求解的时间步长1tΔ固定时，计算结果与紧耦合结果之间的误差会随着求解温度场的时间步长2tΔ的增大而增大，当1tΔ较大时，计算结果与紧耦合结果较为符合，这意味着相应的计算时间也会较大，对计算效率没有太多提高。

从图5和图6中可以看出，减小1tΔ的大小，同时选取合适的2tΔ的值，计算结果与紧耦合结果也会符合得很好。另一方面，在保证精度的同时也应将计算时间尽量减少，结合表2可以发现在算例5的情况下计算1000s耗时为3h，比算例4情况下的4.3h和算例7情况下的3.5h要小，而三者的误差均在1%～1.5%（见表3）。结合上述分析，最终选取算例5即1tΔ=10s和2tΔ=10s作为时间步长应用到本文的流场/温度场松耦合计算方法中，作进一步计算分析。

上述分析结果与1000s时舱内设备的平均温度计算结果见表3，从中进一步可以看出，与紧耦合计算方法相比，虽然算例1、算例4和算例7的计算误差较算例5要小，但计算效率较差，因此算例5的时间步长组合选择为最优。设备表面的在流场中的平均对流换热系数[13]见表4，从表4中发现不同时间步长组合对对流换热系数的计算结果影响不大，当算例9中引起的温度误差为7.51%时换热系数的误差却仅为2.87%。

表3 瞬态1000s时刻，设备平均温度在不同时间步长情况下的计算结果Tab.3 Tab. 3 Average temperature at t=1 000s at different time steps

表4 设备表面平均对流换热系数的计算结果Tab.4 Average heat convection coefficient on the surface of devices

4.2 同时间步长、不同时刻计算结果

根据上一节中得出的结论，选取算例5中的Δt1=10s，Δt2=10s作为流场/松耦合计算方法中的时间步长。在瞬态计算时长为200、500、1000和1500s时，对比分析温度场的（设备1到设备2之间）计算结果，如图7～10所示。

图7 t=200s时刻温度场计算结果Fig.7 Temperature results at t = 200s

图8 t=500s时刻温度场计算结果Fig.8 Temperature results at t = 500s

图9 t=1000s时刻温度场计算结果Fig.9 Temperature results at t = 1 000s

图10 t=1500s时刻温度场计算结果Fig.10 Temperature results at t = 1 500s

可以看出本文流场/温度场松耦合计算方法的结果与紧耦合方法符合得较好，但随着计算时间的增大，二者计算结果的误差也相应增大，4个图中设备的计算温度出现最大误差的时刻在1500s，但也仅为2%左右。舱内设备在不同时刻下的温度计算结果见表5所示。从表中可看出，在选取算例5进行计算同时间步长下不同时刻的设备温度时，与紧耦合方法相比，计算结果均能够满足精算精度要求，进一步说明了本文流场/温度场松耦合计算方法的合理性。

舱内设备表面平均对流换热系数的计算值见表6所示。与表4中的结果类似，舱内设备表面的瞬态对流换热系数基本不受计算方法和时间的影响。究其原因，一方面从对流换热系数的定义分析，物体与流体接触面的对流换热系数h的大小主要取决于表面的几何形状、流体的运动特性及流体的热力学与输运性质，而受其它因素的影响较小[10,13]。

表5 不同时刻下设备平均温度的计算结果Tab.5 Average temperature at different time

表6 不同时刻下设备表面平均对流换热系数的计算结果Tab.6 Average heat convection coefficient at different time

5 结束语

文中首先讨论了数值计算中的耦合对流换热问题，分析了流固界面耦合传热求解原理，指出了不同求解方法的优缺点。针对普遍使用的紧耦合方法计算时间过长的问题，结合一般松耦合求解思想提出了一种针对舱内流场/温度场的松耦合计算方法，并对其进行了详细的数值仿真验证。结果证明此种方法在求解瞬态耦合换热问题时不但具有很好的鲁棒性，与紧耦合方法的计算结果进行对比，精度上也只存在较小的误差。尤其是在计算对流换热系数的时候发现，其计算结果不受到计算时间步长和求解时间的影响，这样进一步验证了本文提出的流场/温度场松耦合计算方法的合理性。

References)

[1]于莹潇,田佳林.美国新型载人火星探测技术方案[J].航天返回与遥感,2009,30(3):1-7,15.YU Yingxiao, TIAN Jialin. USA

[2]黄家荣,范含林.载人航天器生活舱内热湿环境的数值模拟[J].中国空间科学技术,2004(6):7-13.HUANG Jiarong, FAN Hanlin. Numerical Simulation of Thermal/Humidity in Manned Spacecraft Cabin[J]. Chinese Space Science and Technology, 2004(6): 7-13. (in Chinese)

[3]郑忠海.空间站舱内气流分布特性数值模拟与热舒适性评价[D].哈尔滨:哈尔滨工业大学,2005.ZHENG Zhonghai. Numerical Simulation of Air Distribution in Space Station and Thermal Evaluation[D]. Harbin: Harbin Institute of Technology, 2005. (in Chinese)

[4]Garimella S. Flow and Heat Transfer in Space Vehicle Tile Gaps[R]. AIAA Paper No. 92-4061, Houston: American Institute of Aeronautics and Astronautics, 1992.

[5]Hofacker W, Lucas J. Thermal Analysis of Spacecraft by Combining Different Analysis Tools[R]. SAE911584, Detroit:Society of Automotive Engineers, 2009.

[6]钟奇,刘强.航天器密封舱流动和传热的数值研究[J].宇航学报,2002,23(5):44-48.ZHONG Qi, LIU Qiang. A Numerical Investigation on Heat Transfer and Flow in a Pressured Cabin of Spacecraft[J]. Journal of Astronautics, 2002, 23(5): 44-48. (in Chinese)

[7]张吉礼,梁珍,郑忠海,等.载人航天空间站舱内通风对流换热数值研究进展[J].暖通空调,2006,36(1):28-34.ZHANG Jili, LIANG Zhen, ZHENG Zhonghai, etal. Numerical Simulation Progresses of Convection Heat Transfer in Manned Spacecraft Cabin[J]. Heating, Ventilating and Air Conditioning, 2006, 36(1): 28-34. (in Chinese)

[8]Patankar S V. Numerical Heat Transfer and Fluid Flow[M]. New York: McGraw-Hill, 1980: 125-129.

[9]张政.传热与流体流动的数值计算[M].北京:科学出版社,1984:25-28.ZHANG Zheng. Numerical Simulation of Heat Transfer and Fluid Flow[M]. Beijing: Science Press, 1984: 25-28. (in Chinese)

[10]陶文铨.数值传热学[M].第2版.西安:西安交通大学出版社,2001:112-174.TAO Wenquan. Numerical Heat Transfer[M]. Second Edition. Xi’an: Xi’an Jiaotong University Press, 2001: 112-174. (inChinese)

[11]Zabrodskiy S S. Several Specific Terms and Concepts in Russian Heat Transfer Notation[J]. Heat Transfer-soviet Research,1972, 4(2): 1-4.

[12]Versteeg H K, Malasekera W. An Introduction to Computational Fluid Dynamics: The Fine Volume Method[M]. Essex:Longman Scientific Technical, 1995: 35-68.

[13]Bejan A. Convection Heat Transfer [M]. New York: Wiley, 2004: 334-384.

[14]Sun Y S, Emery A F. Multigrid Computation of Natural Convection in Enclosures with a Conductive Baffle[J]. Numerical Heat Transfer, Part A, 1994(25): 575-592.

[15]Anand N K, Chin C D, Mcmath J G. Heat Transfer in Rectangular Channels with a Series of Normally in-line Positioned Plates [J]. Numerical Heat Transfer, Part A, 1995(27): 19-34.