0   引言

低温风洞利用低温气体作为试验介质，提高风洞雷诺数范围，是先进飞行器创新研制必不可少的基础试验设施^[1]。低温环境下，风洞洞体结构存在一定的温度梯度，产生温度附加应力和热变形，严重时可能会引起结构局部损伤，引发安全事故。同时，低温风洞运行消耗大量液氮和电力，探索合理的降温策略有助于提高风洞运行经济性。因此，建立可靠的低温风洞热力学模型对研究风洞运行安全性和经济性是必不可少的。对于热力学模型修正问题，工程上普遍采用基于设计经验的试凑法，即直接依据经验对热力学模型偏差处的参数进行感性调整和试算。近年来，热力学模型修正技术取得一定发展，其中以热网络方程模型修正、基于蒙特卡洛原理的模型修正和基于响应面法模型修正的发展最为引人注目。

李鹏采用试验稳态工况的试验数据对对接机构热力学模型进行修正, 研究了热网络方程修正的实用化方法^[2]。因为系统复杂性以及某些关键参数的不确定性, 航天器系统的热设计和热分析一直是一个难以很好解决的关键问题，国内外均有采用基于蒙特卡洛原理的随机近似方法进行模型修正研究^[3-5]。响应面法是近年快速发展的基于统计分析的模型修正技术，邓小雷使用响应面模型和多目标遗传算法，根据热平衡试验所获得的数控机床主轴系统热态特性数据，利用多目标遗传算法对二阶响应面模型进行循环逼近优化，提高原有限元热力学模型的准确度^[6]。

试凑法每次试算需要调用有限元软件进行解算，因而需要很大的人力和时间投入，效率不高，并且要调整哪些参数全凭经验和感性认识，并没有理论依据^[7]。热网络方程待修正的未知数过多, 修正模型为不定方程组，试验测点与模型节点不统一, 部分模型节点无温度测点等, 导致数据处理和工程实现困难^[8-10]。蒙特卡罗法是一种采用统计抽样理论近似地求解数学或物理问题的方法，主要用于多目标的敏感度分析，且多用于传热分析。响应面方法是一种数理统计学技术，目的在于改善和优化输入与响应过程的数据。响应面方法实施起来相对容易，相对于其他修正方法提高了计算效率，并且避免了灵敏度分析过程，逐渐成为模型修正研究的热点^[11-12]。

本文以确保低温风洞降温安全性和提高风洞运行经济性为工程背景，研究建立可靠性高的有限元热力学模型的方法。以低温风洞扩散段为方法研究对象，在降温试验的基础上，进行有限元热应力分析，对比分析温度、应力监测点试验数据和仿真数据的差别，确定修正参数；通过中心复合试验设计生成有限元样本空间，以监测点试验数据和仿真数据的残差均方和为考核指标，在样本空间内对残差均方和进行非线性回归分析；以所有监测点残差均方和总和为目标函数，在样本空间内进行优化分析，得到最优解；使用修正后的模型进行验证分析，检验修正可靠性。

1   工程背景及降温试验简介

1.1   工程背景

低温风洞的运行成本主要包括3个方面：A是能源消耗(液氮和电力)，B是人员成本，C是维护、保险和其他成本，根据欧洲跨声速低温风洞(European Transonic Wind Tunnel，简称ETW)运行经验，3者比例约为50%: 31%: 19%，可以看到，能源成本占据主要，其中又以液氮成本为主。ETW 1次实验平均需消耗液氮620吨，更换模型车另外需要每次增加30吨，按现在液氮市场价880元/吨，仅液氮一项就需要57万元。低温风洞不仅需要类似传统风洞的常温静态和动态调试(简称为静调和动调)，还需要低温下的静调和动调，其中最重要的是低温下设备安全静调和流动参数精确控制动调。由已有的低温引导风洞调试经验可知，低温下的静调和动调难度大，成本高。根据ETW调试经验，采用有限单元法建立洞体精确仿真模型，有助于提高设备调试的经济性。另外，控制结构件在低温下的热应力和热变形也是低温风洞建设过程中的难题之一。美国国家跨声速低温设备(National Transonic Facility，简称NTF)和ETW建设过程中，均采用计算机程序解决大型、强迫对流和热应力问题。对选定的内部部件进行试验研究和数值分析，验证设计方法的可行性。

1.2   降温试验简介

低温风洞扩散段(简称扩散段)主要作用是把气体的动能恢复为压力能，从而减少气流在扩散段下游各段的能量损失。扩散段主要由驻室、内流道、驻室进气管道和支座组成，如图 1所示。冷却气流在内流道中快速流动，使内流道结构快速降温。为了降低驻室和内流道之间的温度梯度，在驻室进气管道中通入冷却气体，但相对流量较小。驻室内的流动较弱，可视为自然对流，内流道气流流速较快，气流与结构之间是强迫对流换热。

图  1  扩散段三维模型

Fig.  1  Diffuser section model

下载: 全尺寸图片 幻灯片

依托中国空气动力研究与发展中心0.3m低温高雷诺数风洞，开展降温试验，如图 2(a)所示，白色方框为温度测点位置，白色五星为应变测点位置。在驻室与内流道连接处布置一个应变测点，监测应力变化；在驻室锥形体外表面水平母线上，均布6个温度测点，图中由左至右分别为1~6号测点，监测温度梯度。在内流道、驻室靠近内表面处布置温度探针，监测与结构发生对流换热的流体温度，降温时间约2.5h，气流温度从室温降至107K。

图  2  降温试验

Fig.  2  Cool-down test

下载: 全尺寸图片 幻灯片

降温试验测试结果如图 2(b)所示，图中，T00、T0表示内流道和驻室温度探针测试得到的温度曲线，T1~T6分别表示图 2(a)中所示的6个温度测点测试曲线。温度探针直接测试流体温度，内流道中为冷却气体，且流速较快，驻室中冷却气体是由驻室进气阀门自内流道中引导而来，流速相对较慢，故T00比T0降温要快。T1~T6为结构外表面温度测点，T1靠近驻室进气口，T6靠近内流道，故T6温度最低，T1至T6温度曲线应该为先升后降，试验结果与预测一致。S1为图 2(a)中所示应变片测得的应力曲线，随着温度降低，驻室结构外表面轴线方向温度梯度增大，应力上升，实测结果与预测一致。

2   有限元分析及修正参数选择

使用有限元软件，在降温试验的基础上，对低温风洞扩散段进行有限元热应力分析。根据仿真结果，对比分析温度、应力监测点试验数据和仿真数据的差别，确定修正参数；

2.1   有限元热应力分析

扩散段有限元热力学模型如图 3所示，采用半模分析，内流道截面积变化不大，内流道内表面对流换热系数设置成统一的，对于圆管内湍流的对流换热系数估计，常用的较为精确的经验公式是格列林斯基(Gnielinski)公式，表述为：

图  3  扩散段有限元热力学模型

Fig.  3  Finite element model of diffuser section

下载: 全尺寸图片 幻灯片

式中:h为对流换热系数；L为特征长度；k_f为热导率；Re为雷诺数；Pr为普朗特数。

根据格列林斯基公式可以得到随流体变化的圆管内湍流强迫对流换热系数^[13-14]。驻室内的流动缓慢，流场和结构之间传热视为自然对流，对流换热系数取为5，降温曲线为试验实测曲线，如图 4所示。为验证此模型的准确性，在驻室锥形体外表面母线上均匀布置温度监测点，在驻室和内流道交界处布置应变监测点。

图  4  流场降温曲线

Fig.  4  Flow field cool-down curve

下载: 全尺寸图片 幻灯片

计算单元采用三维连续体20节点温度耦合缩减积分单元，即C3D20RT，局部网格采用单向偏置细化保证应力计算准确，计算结果对比总应变能和伪应变能的比值，验证计算准确性。

分析时长为10 000s，在9300s左右到达温度最低点，约108K。图 5为有限元计算结果，左侧为应力云图，右侧为温度云图。可以看出：在内流道和驻室连接处存在明显的温度梯度和应力梯度。

图  5  有限元计算结果

Fig.  5  Results of FEA

下载: 全尺寸图片 幻灯片

2.2   修正参数选择

对比试验与仿真数据，图 6为对应点的应力随时间变化曲线，仿真首先进行静力学分析，故起始阶段不为零，总体看来，计算所得应力偏大。图 7为各测点温度随时间变化曲线，实线为仿真值，曲线为试验值，可以看出第6点仿真与试验吻合较好，其余均为较大误差。

图  6  试验与仿真应力曲线

Fig.  6  Stress curves of test and simulation

下载: 全尺寸图片 幻灯片

图  7  试验与仿真降温曲线

Fig.  7  Temperature curves of test and simulation

下载: 全尺寸图片 幻灯片

试验值表明，驻室锥形体外表面存在温度梯度，但是在仿真中表现不明显，尤其是1~3号测点。原因在于仿真中，驻室采用统一的降温曲线，实际上驻室内流体温度存在一定温差。同时考虑到驻室内的内流道对驻室内气体的冷却作用，以及驻室冷气进气管道对局部的冷却作用，如图 8所示，即区域① 和区域② 对驻室局部气流的冷却作用，在驻室锥面不同地方，冷却作用效果不一。因此，驻室不同区域单位时间降温量存在一定的误差，故需要对驻室不同区域对流换热量进行修正，提高热力学模型准确性。

图  8  误差分析

Fig.  8  Error analysis

下载: 全尺寸图片 幻灯片

由牛顿冷却公式可知，单位时间的对流换热量等于温差与对流换热系数的乘积，表述为：

式中：q″为热流密度；T_∞为流体温度；T_s为固体表面温度；h为对流换热系数。若选择对温差修正，即修正图 4中驻室降温曲线，难度较大。通过修正对流换热系数，同样可以达到修正对流换热量的目的。因此，本文选择对驻室内对流换热系数进行修正，原来热力学模型中使用统一的对流换热系数，修正为不同区域使用不同对流换热系数。

3   热力学模型修正

通过对比分析试验与仿真数据，结合试验点布置方式，确定驻室局部对流换热系数为待修正参数，进行扩散段热力学模型修正。图 9为待修正对流换热系数对应的驻室区域，在温度梯度较大的区域设置5个对流换热系数。修正步骤为：

图  9  待修正参数说明

Fig.  9  Parameters updated

下载: 全尺寸图片 幻灯片

(1) 使用中心复合试验设计构造采样点，在每个采样点调用有限元模型计算输出参数；

(2) 利用输入和输出参数构造响应面，并对响应面进行回归分析；

(3) 使用分析和试验结果构造目标函数；

(4) 在响应面模型的范围内对目标函数进行迭代优化分析，得到优化后的修正参数；

(5) 模型验证，使用修正参数检验新模型与试验数据误差。

3.1   中心复合试验设计

基于回归分析的响应面拟合是对样本数据进行操作, 利用试验设计方法，可以用较少的样本点数(降低有限元分析的计算量)保证较高的响应面模型的精度^[15-17]。在众多的样本空间设计方法中，中心复合设计应用最为广泛^[18-19]。本次修正参数为5个，使用中心复合试验设计，各目标有5个水平，即5目标5水平中心复合试验设计，共需要27个样本，修正目标水平如表 1所示，表中h为对流换热系数。

表  1  修正目标水平表(单位：W/(m²·K))

Table  1  Modal updating parameters and level

Level	-2	-1	0	1	2
h	2	7.75	13.5	19.25	25

下载: 导出CSV 

| 显示表格

在试验时间内每隔50s选择一个时间点，得到每个时间点上试验值和仿真值的残差，对同一模型每一个测点的所有残差求均方和，即：

式中:y_test、y_sim分别为对应的试验和仿真数据。R_Error²越小，则表明试验值与仿真值误差越小，两者更吻合，热力学模型更准确。

3.2   建立响应面模型

响应面法就是根据研究对象的特点，在试验设计的基础上，用多项式或其它响应面模型近似描述设计变量和响应特征之间的复杂关系，得到响应特征的响应面模型，利用该模型来预测非试验点的响应值。实际中根据工程经验，通常选取二次多项式形式的响应面模型^[20]。完全二次多项式模型公式为：

式中:x_i是预测变量，对应热力学模型中的待修正的参数；ε为误差项。假设待修正参数个数为m，对应的二次多项式响应面展开为：

假设试验次数为n，用矩阵形式，令：

则：Y-Xa=ε

用最小二乘法拟合估计a，假设

令，即，可得：

求解式(8)，即可得到，即响应面函数系数。

针对本次研究，根据公式(5) 建立5目标值响应面模型，采用完全二次多项式为非线性回归分析拟合函数，生成的响应面模型为：

式中:y为同一测点仿真和试验残差均方和R_Error²; x_i为待修正的对流换热系数; a_i为响应面系数，共21个。对所求得的响应面模型进行拟合度检验，采用相关指数R²作为衡量拟合后曲线效果好坏的指标^[21]，即：

使用式(10) 对式(9) 所求响应面进行评估，求得相关指数分别为0.9286、0.9717、0.9726、0.9474、0.9865、0.9756和0.9621，均较接近1，拟合度较高。

3.3   优化分析

构造优化目标函数为所有残差均方和的总和，约束条件为响应面区间内，即：

函数优化过程如图 10所示，经过52次运算后，函数收敛在容差范围内。优化后的参数水平值及对应值如表 2所示，h1~h5优化后水平为1.88，-1.95，-1.55，-1.66和0.35，对应的实际对流换热系数为24.34，2.32，4.59，3.95和15.48，单位为W/(m²·K)。

图  10  优化过程

Fig.  10  Optimization process

下载: 全尺寸图片 幻灯片

表  2  参数修正结果(单位：W/(m²·K))

Table  2  Model updating result

Parameter	h1	h2	h3	h4	h5
Updated level	1.88	-1.95	-1.55	-1.66	0.35
Updated h	24.34	2.32	4.59	3.95	15.48

下载: 导出CSV 

| 显示表格

3.4   模型验证

采用表 2修正后参数，更新有限元热力学模型，使用修正后的热力学模型，进行热应力计算，与试验值和初始热力学模型数据作对比，图 11为应力对比，图 12为1~6号测点温度对比，图中Test表示试验曲线，Update表示修正后的仿真曲线，Sim表示初始热力学模型仿真曲线。由图可以看出，修正后的热力学模型应力和各监测点温度均与试验值更接近，说明修正后有限元热力学模型与实际更吻合，提高了有限元热力学模型的可靠度。

图  11  修正后应力对比

Fig.  11  Stress comparison after updating

下载: 全尺寸图片 幻灯片

图  12  修正后温度对比

Fig.  12  Temperature comparison after updating

下载: 全尺寸图片 幻灯片

4   模型应用

使用修正后的热力学模型分析另外一次降温试验，降温时间为2h，对比应力监测点数据，如图 13所示。由图可见，仿真应力数据与试验应力数据较为吻合，说明修正后的模型可以用于其它降温试验分析，达到预期目的。但可以看出，试验和仿真仍存在一定的误差，表 3为每隔1000s，试验应力值和修正应力值的误差分析，误差随时间逐渐减小，而应力值随时间是增大的，即高应力时，仿真与试验的误差是较小的，说明模型预测应力可靠性较高。

图  13  应用验证

Fig.  13  Applying validate

下载: 全尺寸图片 幻灯片

表  3  试验应力值和修正应力值误差分析

Table  3  Stress error analysis of test and updated model

Time/s	Update value/MPa	Test value/MPa	Error/%
1000	1.90	15.35	107.8
2000	47.74	37.07	28.8
3000	53.88	43.16	24.8
4000	59.20	48.27	22.7
5000	73.83	68.44	7.9
6000	90.79	88.56	2.5

下载: 导出CSV 

| 显示表格

经过分析，可以确定误差主要来源有2点：(1) 稳定流场的建立需要一段时间。仿真假设流场是均匀稳定的，实际上稳定流场的建立是需要一段时间的，随着时间的增加，稳定流场逐步建立，误差也随之减小；(2) 仿真没有考虑驻室内表面对流换热系数在径向上的变化。实际上驻室内表面对流换热系数是空间坐标的函数，初始仿真模型采用的是常量，修正模型假设对流换热系数是驻室轴向坐标的分段函数，简化了模型修正的难度和工作量，同时也带来了一定的误差。综上所述，修正后的模型能有效预测实际结构的热应力，预测误差随试验时间减小，在合理范围之内。

5   总结

(1) 联合中心复合试验设计和响应面法来提高低温风洞扩散段热力学模型的可靠性是可行的，通过验证，修正后的模型与实际模型误差更小，相关性更高。

(2) 修正后的模型可以用作低温风洞其他降温试验，说明在一定试验范围内，修正模型都是可靠的。

(3) 扩散段热力学模型修正研究，为低温风洞其他部段和整个洞体的模型修正提供了可行性方法研究，为下一步开展更大部段，更复杂模型修正提供基础。