0   引言

低温风洞通过降低试验气体温度可以提高试验雷诺数，能够满足飞行器气动性能准确模拟对雷诺数的需求，是十分重要的飞行器气动性能地面模拟设备。试验气体降温过程中，低温风洞结构与气流间存在着不均匀热交换，结构局部易出现较大温度梯度，引发热变形，且受各类边界约束以及气动载荷的作用，导致结构应力超出许用范围，从而威胁设备运行安全。另外，风洞结构变形会引起风洞回路轴线偏离、气动型面变形和内部气流表面出现阶差等诸多问题，影响风洞的气动性能。因此，分析和控制低温环境下风洞结构的热变形是低温风洞研制的关键技术问题。

复杂结构的热变形分析通常采用试验测量和数值仿真等方法。相对而言，数值仿真能够准确、全面地反映结构温度场和应力场分布情况，为结构热力学特性优化提供可靠、充分的参考依据，且成本较低，分析周期较短。所以，基于数值仿真方法，国内外学者对低温结构的热力学特性开展了广泛研究。朱鸿梅^[1-2]采用有限元方法，对冷却过程中非磁性低温杜瓦的螺纹粘接接头进行瞬态热应力仿真，分析了等效应力与冷却时间的变化关系。Kim、Kang^[3-4]等对复合材料/铝环状试件的热力学特性进行了研究，并将低温热应力与有限元仿真结果作出对比分析。朱立伟^[5]应用有限元分析软件分析了LNG船用超低温球阀在超低温条件下的应力特性，并提出改进应力集中的措施。上述研究均采用有限元分析方法，关注的是结构的热力学特性，而低温环境下考虑流固热耦合作用的热力学特性分析尚不多见。考虑流固热耦合作用的研究，也大多关注流体传热^[6-7]，或温度变化较小的结构传热^[8-9]。

本文以风洞典型结构部段--扩散段为研究对象，基于流固热耦合方法，建立低温风洞扩散段的流固热耦合仿真模型，重点分析低温下扩散段的换热特性，计算扩散段结构的温度及应力分布，并在风洞低温运行过程中对扩散段温度和应力进行测试研究。

1   试验装置

低温试验在中国空气动力研究与发展中心的0.3m跨声速低温风洞中进行，试验主要研究扩散段的热力学特性。低温风洞正常情况外壁面有保温层，风洞调试期间，壁面暂未加保温层，故本文计算按未加保温层考虑。风洞运行时，扩散段将气流动能恢复为压力能，以减少气流在扩散段下游各段的能量损失^[10]。扩散段结构关于垂直中心面左右对称，由内流道、驻室、驻室进气管道和支座组成。4个支座均为滑动支座，通过双头螺柱压紧于绝热块上，从而限制结构在竖直方向的位移，同时保证低温工况下在水平方向可自由收缩。内流道中通入快速流动的低温氮气，使流场和结构得以快速降温。为降低驻室与内流道之间的温度梯度，需从风洞回路向驻室引入冷却气体，其温度与内流道中的氮气接近，但流速相对较小。结构外壁面直接与空气接触。扩散段模型如图 1所示。

图  1  低温风洞扩散段模型

Fig.  1  Structure model of cryogenic wind tunnel diffuser section

下载: 全尺寸图片 幻灯片

2   计算方法

降温过程中，扩散段换热不均匀，导致局部存在温度梯度，引发热应力。本文采用流固热耦合方法，求解结构的温度分布；采用顺序热力耦合方法，求解结构的应力分布。

2.1   传热理论

传热有3种方式:导热、对流和热辐射。本文研究对象处于低温工况，可忽略辐射效应。

流体内部的传热涉及3个重要物理定律：质量守恒定律、动量守恒定律和能量守恒定律，此处不再赘述。

流体的管内流动状态取决于当地雷诺数Re，Re＞10⁴时为完全湍流^[11]。内流道入口雷诺数为2.96×10⁶，采用标准k-ε两方程模型进行模拟。

本文中流体通过耦合界面与固体发生热交换，壁温是时间的函数，为第一类边界条件。

结构内部传热方式为导热，直角坐标系中的导热方程通常写作：

式中：为单位体积产生的热能。无内热源，则=0。

边界条件为第三类边界条件：

式中：T_w为固体壁面温度，T_f为膜层温度，即对流边界层中流体的温度。考虑热边界层对换热的影响，用T_∞表示远场温度，T_f可用下式计算：

结构外壁面与大气传热方式为自然对流；内流道结构与流体之间的传热方式为强制对流。驻室腔体内氮气流速很低，与驻室壳体间的传热方式可视为自然对流。

等温水平长圆柱的大空间自然对流换热，可采用以下关联式^[11]估算：

式中:Ra为瑞利数，Pr为普朗特数，要求Ra≤10¹²。

2.2   流固热耦合计算方法

多物理场耦合计算通常分为直接耦合法和间接耦合法。直接耦合是对问题涉及的所有数学物理方程进行联合求解，一般适用于较简单的物理问题。本文采用间接迭代耦合法，分别对流场和结构进行分析，再将耦合界面的物理量通过映射与差值^[12-13]传递给外场进行计算，如此循环迭代。

采用多物理场代码耦合软件MpCCI，联合Abaqus与Fluent进行流固热耦合仿真。在结构与流场耦合界面上，交换的物理量为膜层温度T_f、换热系数h和壁温T_w。壁温由Abaqus计算，经MpCCI传入Fluent。耦合壁面的换热系数与膜层温度由Fluent计算，经MpCCI传入Abaqus。计算采用串行算法，即Abaqus收到Fluent传入的初始数据后开始迭代，经过固定耦合时间步长后，将计算得到的耦合量(节点温度)传入Fluent；Fluent收到数据后经相同时间步长的迭代，将计算得到的耦合量(换热系数与膜层温度)传入Abaqus，依此循环，典型串行算法的数据传递过程如图 2所示。

图  2  串行耦合算法

Fig.  2  Serial coupling algorithms

下载: 全尺寸图片 幻灯片

2.3   热应力计算

热应力是由于部件受热不均匀而存在温度差异，导致各处膨胀变形或收缩变形不一致，且相互制约而产生内应力。对各项同性材料，自由膨胀的变分量为:

式中：ε_x、ε_y、ε_z为沿X、Y、Z方向应变分量；α为线膨胀系数。

如果物体边界有约束条件，自由膨胀受到约束限制，微元体会产生热应力。根据线性应力原理和胡克定律，描述应力和应变关系的物理方程为：

式中：E为材料杨氏模量；μ为泊松比；σ_x、σ_y、σ_z为沿X、Y、Z方向应力分量；G为材料剪切模量，G=2E(1+μ)；ΔT为从无应力状态开始的温差。

扩散段中温差较大区域为驻室与内流道的过渡区域，这也是应力分析重点关心的区域。

3   建立流固热耦合仿真模型

本节分别建立流体与结构传热模型，在MpCCI中设置相关耦合参数，联合Abaqus与Fluent进行流固热耦合分析。将结构温度场作为载荷导入Abaqus中，进行顺序热力耦合分析。

3.1   流体传热模型

扩散段内的流场比较复杂，对其做出以下几点简化。首先，驻室内气流压差小，且大部分位置气流速度低于1m/s，对换热能力的促进作用较小，因此以自然对流作为后期结构热固耦合的边界条件。流体传热模型中仅考虑内流道的流场，内流道壁面为耦合界面。其次，扩散段内流道实际入口气流较紊乱，速度分布复杂，此处简化为简单的均匀来流。最后，不考虑介质的摩擦损失以及介质的物性参数随温度的变化。

流体传热模型入口马赫数为0.2；入口温度随时间变化，通过Matlab对实测入口温度进行多项式拟合，并写入UDF，如图 3所示。出口边界设置为outflow；壁面边界条件设置为不滑移，温度由Abaqus传入。

图  3  降温曲线与拟合降温曲线

Fig.  3  Cooling curve and fitting cooling curve

下载: 全尺寸图片 幻灯片

湍流模型采用k-ε两方程模型，壁面函数采用标准壁面函数，要求第一层网格厚度位于对数律区域，且Y^*>11.225。经Fluent计算后，调整第一层网格厚度，直到满足湍流模型对Y^*的要求。

网格为六面体结构网格，数量为293 436，如图 4所示。为验证网格无关性，对以下工况--入口马赫数为0.2，入口流体温度为270K，壁面温度为300K--进行验证，对比壁面平均对流换热系数，结果如图 5所示。

图  4  流体网格

Fig.  4  Fluid gridding

下载: 全尺寸图片 幻灯片

图  5  网格无关性验证

Fig.  5  Grid independence verification

下载: 全尺寸图片 幻灯片

3.2   结构传热模型

对结构作出以下几点简化。风洞扩散段是薄壁壳体，有限元采用实体单元对实际模型简化少，能较好反映实际结构，但厚度方向节点不易控制，网格质量不好，计算量大。因此本文采用壳体单元建模，可提高计算效率。此外，降温过程中，空气中的水蒸气会在扩散段外壁面凝华结霜，对传热会产生一定影响。建模中通过设置较小的对流换热系数来近似模拟，不考虑结霜过程。最后，不考虑支座处摩擦阻力对位移的约束。

扩散段结构关于垂直中心面左右对称，对于有限元计算，只需建立半模即可，如图 6所示。采用Abaqus前处理器对结构划分网格，并细化驻室与内流道交界处网格，保证壳单元厚度方向至少有5个积分点^[14-15]，总网格数为19 714。传热分析采用传热单元DS4。结构传热分析中只需考虑模型的换热问题。初始温度设置为300K；外壁面自然对流换热系数可按公式(4)进行估算，驻室外表面不同截面处的直径、温度均变化较大，计算时采用算数平均值。降温末期，结构外表面结冰结霜会减弱传热效果，因此换热系数的取值应适当减小。结构外壁面换热系数取为5W/(m²·K)；驻室内冷却气体引自风洞回路，温度与内流道氮气温度接近，因此将内流道入口测量温度作为驻室气体温度。驻室内气流速度小，在大部分位置速度低于1m/s，因此简化为自然对流边界处理。自然对流换热系数通常在10W/(m²·K)以内，参照封闭腔内自然对流换热工况^[16]下的换热水平，将驻室内换热系数取为3.5W/(m²·K)。

图  6  结构网格

Fig.  6  Structure gridding

下载: 全尺寸图片 幻灯片

3.3   流固热耦合模型

在流固热耦合模型中，耦合界面为内流道壁面。总的降温时间为6400s，耦合物理量为流体温度、换热系数和固壁温度。

温度从300K降至110K，每秒钟温度下降0.03K。软件中，流固热耦合时间步长设置为1s。在该时间步长内，耦合量的变化都足够小，从而可准确反映瞬态过程中的换热特性。

3.4   结构应力模型

热应力分析中采用S4R壳体单元，支座处限制Y方向自由度，出口端面限制Z方向自由度，将传热分析得到的温度场作为温度载荷加载到应力分析步。

4   试验测试

测试系统包括低温电阻温度传感器、低温应变传感器、数据采集器和显示系统。

温度和应变传感器布置在驻室与内流道外壁面的水平位置，实时监测结构的温度与应力变化。粘贴传感器之前，先要对粘贴表面进行处理，保证传感器与表面的紧密贴合。然后在传感器与钢结构之间涂上刚性较好的低温胶，在传感器背面涂上弹性较好的防护胶。通过在测点附近粘贴补偿块的方式，消除温度变化造成的传感器误差。

数据采集仪型号为吉时利2701，将其与工控机箱串联，通过电脑显示器即可设置相关参数，实时监测结构的温度与应力变化。

试验中测点分布如图 7所示，其中T1~T13为温度测点。测点T1~T9位于驻室椎形体外表面水平母线上，测点T10位于冷气入口管道与驻室连接处附近，测点T11~T13位于冷气入口管道外壁；示意图中三角形代表应力测点S1，位于驻室与内流道连接处。在内流道、驻室靠近内表面处布置温度探针，监测与结构发生换热的气体温度。

图  7  测点分布图

Fig.  7  Diagram of measuring points distribution

下载: 全尺寸图片 幻灯片

降温时间为6400s，气流温度由300K降至130K。此间，结构表面会结霜。随后，风洞恢复常温，降温过程如图 8所示。

图  8  降温过程

Fig.  8  Cooling process

下载: 全尺寸图片 幻灯片

5   试验与仿真结果分析

本节首先将低温风洞试验温度测量结果与流固热耦合仿真结果进行对比，分析扩散段结构的传热特性。然后将低温风洞试验应力测量结果与顺序热力耦合仿真结果进行对比，分析结构应力状态。

5.1   耦合传热分析

由降温曲线的斜率可判断降温速率：0~2500s，温度均匀下降；2500~3400s，降温速率逐渐减慢；3400~4500s温度均匀下降；4500~4800s，降温速率突然变大，随后恢复之前的降温速率。

流固热耦合仿真计算得到结构在1000、2500、5000和6400s的温度场，如图 9所示。由温度分布可知，驻室与内流道过渡区域的温度梯度较大。驻室内气流速度低，换热效果差，且热量传导至温度较低的内流道需经过结构导热，以上原因导致驻室外壁温度明显高于内流道外壁温度。

图  9  降温过程结构温度分布

Fig.  9  Distribution of structural temperature when cooling down

下载: 全尺寸图片 幻灯片

测点T2位于内流道外壁，T4位于驻室与内流道过渡区域。通过对比发现，仿真曲线趋势与测量曲线保持一致，最大偏差为20K左右。2500~3400s，降温速率减慢，温度测量与仿真曲线的斜率也随时间减小，符合实际情况。4500~4800s，降温速率突然变大，在仿真曲线中可明显看到斜率增大，而实际测量曲线并未反映出该降温速率的变化，如图 10和11所示。

图  10  测点T2温度变化曲线

Fig.  10  Temperature change curve of T2

下载: 全尺寸图片 幻灯片

图  11  测点T4温度变化曲线

Fig.  11  Temperature change curve of T4

下载: 全尺寸图片 幻灯片

测点T7、T9位于驻室外壁面，由本节前面的分析可知该处降温速率较慢。实际测量发现，温度在250K以上，计算结果与试验数据吻合较好，如图 12和13所示。

图  12  测点T7温度变化曲线

Fig.  12  Temperature change curve of T7

下载: 全尺寸图片 幻灯片

图  13  测点T9温度变化曲线

Fig.  13  Temperature change curve of T9

下载: 全尺寸图片 幻灯片

测点T10位于驻室上，靠近冷气入口管道底部。测点T12位于冷气入口管道上。由于管道内气流速度较低，换热不充分，所以该处结构温度较高。试验测量与仿真计算偏差较小，在10K以内，如图 14和15所示。

图  14  测点T10温度变化曲线

Fig.  14  Temperature change curve of T10

下载: 全尺寸图片 幻灯片

图  15  测点T12温度变化曲线

Fig.  15  Temperature change curve of T12

下载: 全尺寸图片 幻灯片

总的来说，温度测量与仿真曲线吻合较好。最大偏差在20K以内，发生在降温末期。

5.2   结构应力分析

结构应力的大小与温度梯度成正比。由5.1节传热分析可知，驻室与内流道过渡区域应力较大。图 16为结构在1000、2500、5000和6400s时刻的应力分布。由应力分布可知，支座、法兰以及驻室与内流道过渡区域应力较大。由于采用壳体单元建模，支座、法兰处过于简化，可能导致应力计算不准确。

图  16  降温过程结构应力分布

Fig.  16  Distribution of structural stress when cooling down

下载: 全尺寸图片 幻灯片

对比应力仿真曲线与结构降温曲线可知，应力变化与温度变化基本对应。在0~2500s，均匀降温，应力随时间增大；2500~3400s，降温速率逐渐减小，由于结构内部导热，局部温度梯度减小，热应力减小；3400~4500s，均匀降温，应力随时间增大；4500~4800s，降温速率突然变大后恢复，S1处应力曲线也出现一个凸起。

由图 17可知，应力仿真值与实际测量值吻合较好。最大偏差为10MPa左右，发生在降温末期。温度应力最大值121MPa出现在降温末期，在驻室与支座之间。使仿真与计算产生偏差的原因，可能有以下3点：(1)风洞回路中，低温氮气经过试验段后均匀性较差，导致结构传热不均匀；(2)降温末期，扩散段外结霜结冰，影响应力测量；(3)防护胶经过低温后粘合度下降，甚至出现脆裂。

图  17  降温过程测量与仿真应力变化曲线

Fig.  17  Simulation and measuring curve of stress when cooling down

下载: 全尺寸图片 幻灯片

通过仿真分析和试验测试，扩散段温度应力低于材料的许用应力(115MPa)，在低温和给定的试验温度变化条件下，扩散段结构是安全的。

图  18  最大应力位置

Fig.  18  The position of maximum stress

下载: 全尺寸图片 幻灯片

6   结论

通过对低温风洞扩散段的数值计算与试验得出如下2点结论：

(1)采用壳体单元对风洞扩散段的结构进行有限元建模是可行的，可以减小计算规模。尽管局部位置可能过于简化，但应力计算值与仿真值基本在一个量级。

(2)利用第三方软件MpCCI联合Abaqus与Fluent进行流固热耦合仿真，得到结构的温度场分布；再将温度场作为载荷，利用Abaqus进行顺序热力耦合分析，求解结构的应力场。通过风洞低温试验测试结果表明，该方法能够准确反映低温风洞结构的热力学特性，可为低温风洞的结构安全性能优化提供可靠的仿真分析方法。