接触热阻对超声速飞行器结构响应分析的影响

2020-06-06张仡,成竹,秦强,陈宏

科学技术与工程 2020年12期

张仡,成竹,秦强,陈宏

(中国飞机强度研究所全尺寸飞机结构静力/疲劳航空科技重点实验室，西安 710065)

近年来，超声速飞行器日益受到各国的广泛关注[1]。与亚声速飞行器不同，超声速飞行器表面的气动加热现象更为显著，对飞行器结构温度与应力的影响更为突出。在气动加热过程中，热量通过蒙皮传递到内部梁、肋等结构中，界面间的接触热阻会直接影响结构的温度与应力。进行分析时，对接触热阻的影响估计过高会使结构过于笨重，而估计过低会带来安全隐患。因此，必须准确判断接触热阻对于结构响应的影响程度。

中外学者在接触热阻理论、仿真与试验方面已经进行了大量研究。在理论方面，通过描述接触面表面形貌与变形，并参照试验数据，建立了一系列预测接触热阻的理论模型[2-4]。Fletcher[5]、Sridhar等[6]和Lambert等[7]分别综述了部分理论模型。对接触热阻的理论研究拓展到微纳尺度，提出了基于热载子的理论模型[8]。在仿真与试验方面，Mantelli等[9-10]、Song等[11]通过仿真与试验研究了螺栓连接件中的接触热阻；Tirovic等[12]通过仿真与试验研究了汽车刹车片界面的接触热阻；Milanez等[13]对低压力下金属间接触热阻进行了试验研究并提出了相应的理论模型；Liu等[14-15]对C/C复合材料与GH600高温合金之间的接触热阻进行了试验研究，并通过数值仿真研究了接触热阻对疏导式热防护结构防热效果的影响。目前，对于接触热阻的研究主要集中在电子、卫星、核能等领域，而对于超声速飞行中的接触热阻问题关注较少[13]。但考虑到热问题在超声速飞行领域中的重要地位，有必要开展接触热阻对超声速飞行器结构响应分析影响的研究。

为此，以超声速飞行器外表面沉头螺栓局部模型为例，研究在数值仿真中，接触热阻对超声速气动载荷下结构响应的影响，为超声速飞行器设计及安全评估提供参考。

1 接触热阻理论

在Cooper等[2]提出的(Copper-Mikic-Yovanovich)CMY接触热阻模型中，假设金属表面粗糙峰的高度服从高斯分布，粗糙峰在金属表面上随机分布，粗糙峰变形均为塑性变形。基于以上假设，CMY模型中金属界面间接触热导表示为

(1)

(2)

ψ(ε)=(1-ε)1.5, 0<ε<0.3

(3)

(4)

式中：hc为接触热导；n为接触点密度；a为接触点平均半径；ks为等效热导率；ψ为热流约束参数；ε为接触点相对尺寸；k1与k2分别为界面两侧材料的热导率；Ar为接触面实际接触面积；Aa为接触面名义接触面积。

定义相对平均板间距λ=Y/σ，其中Y为实际平均板间距。对于塑性变形，Ar/Aa、n、a可用λ表示为

(5)

式(5)中：σ为等效均方根(RMS)表面粗糙度；m为等效粗糙峰斜率平均值。σ与m由界面两侧的表面形貌参数表示为

(6)

(7)

式中：σ1与σ2分别为界面两侧表面RMS粗糙度；m1与m2分别为界面两侧表面粗糙峰斜率平均值。

通过对金属表面形貌的试验研究，Lambert等[16]提出的m与σ的试验关系式为

(8)

由式(1)～式(5)得到无量纲形式的接触热导Cc为

(9)

式(9)中,无量纲接触热导只与相对平均板间距λ有关。在CMY模型中，λ只能通过试验获得。为解决这一问题，Yovanovich[3]提出λ与P/Hc(其中，P为接触面压力，Hc为材料微硬度)的关系为

(10)

(11)

(12)

(13)

(14)

2 考虑接触热阻的局部建模

在热力耦合接触问题中，接触热阻使得结构应力场与温度场之间存在复杂的耦合关系。根据式(9)、式(10)，面间压力对接触热阻有直接影响，而接触热阻直接影响温度场，从而使应力场间接对温度场造成影响；而温度场又直接影响应力场。为了避免间接耦合引起的计算误差，采用直接耦合法进行热力耦合分析。

为在计算中引入接触热阻，需在接触关系中，输入一系列不同材料、压力与温度下的hc数值。计算过程中，通过线性插值获得局部区域在给定材料、压力与温度下的hc，从而将接触热阻引入传热方程中，从而进行考虑接触热阻的热力直接耦合计算。

问题的计算流程图如图1所示。首先进行常温下应力分析，在接触面间形成初始接触状态，并通过插值获得常温下接触热阻。将常温下的应力场作为初始条件，开始直接热力耦合计算，对问题进行迭代求解。每次计算中，需判断计算前后的接触状态是否变化以及应力场与温度场是否收敛。如不收敛，则需根据计算得到的接触状态、应力场与温度场，重新计算材料的力热参数及界面间的接触热阻，进行新的直接热力耦合计算。如此反复迭代直至收敛，从而获得考虑接触热阻的稳态结构响应。

2.1 几何模型

选取超声速飞行器前缘中，以螺栓为中心，边长30 mm的方形局部模型进行研究。考虑到几何与载荷的对称性，取该方形模型的一半进行研究，如图2所示。为降低计算量，不考虑沉头螺栓与螺母间的接触。

图1 考虑接触的热力耦合计算流程图Fig.1 Flow chart of thermo-mechanical coupling calculation considering contact

图2 局部模型网格Fig.2 Grids of local model

图2所示的模型中，蒙皮与肋的厚度为3 mm，梁的厚度为2 mm。沉头螺栓公称长度l=8 mm，螺纹规格为d=M5，其余尺寸参照《MJ螺纹十字槽100°沉头螺栓》(GJB 3374/29)。螺母螺纹规格为D=M5，其余尺寸参照《抗剪型MJ螺纹高锁螺母》(GJB 3377—1998)。

采用考虑传热的六面体减缩积分实体单元进行有限元仿真。由于螺栓孔边、螺栓沉头处与倒角处可能出现应力集中，对这些区域的网格进行了加密处理。

2.2 材料及其物理属性

模型中，蒙皮、梁、肋的材料为TC4钛合金，沉头螺栓与螺母的材料为GH2132高温合金。两种材料的热导率、线膨胀系数、弹性模量与泊松比如表1～表3所示。TC4钛合金退火后的硬度为HBS 255～341，固溶时效硬度为HBS 293～361。GH 2132高温合金硬度统计值为HBS 293～306。计算时取TC4与GH 2132硬度均为HBS 300。

表1 不同温度下的材料热导率

表2 不同温度下的材料线膨胀系数

表3 不同温度下的材料弹性模量及泊松比

2.3 模型中的接触关系

模型中的接触面包括：蒙皮与梁之间的接触面；梁与肋之间的接触面；蒙皮、梁、肋与螺栓及螺母的接触面。根据接触面两侧的材料，将接触面分为TC4与TC4之间的接触面、TC4与GH2132之间的接触面。

给定所有接触面间切向摩擦系数均为0.15。

为获得界面间无量纲形式接触热导，需得知界面两侧材料的布氏硬度HB，界面两侧表面的RMS粗糙度σ1与σ2，以及界面平均压力。首先，由较软材料的布氏硬度HB，通过式(12)、式(13)算得维氏微硬度系数c1与维氏微硬度无量纲系数c2。之后，由σ1与σ2通过式(6)～式(8)分别算得等效RMS表面粗糙度σ与等效粗糙峰斜率平均值m。根据《GH2132 MJ螺栓螺纹和螺钉通用规范》(GJB 8618—2015)，参考机械设计中对于连接部位的相关要求，给定蒙皮、梁、肋及螺栓的表面粗糙度均为σ=3.2 μm。随后，由σ、m、c1、c2通过式(11)、式(12)算得相对接触压力P/Hc与相对平均板间距λ。最后，将λ代入式(9)获得无量纲形式的接触热导Cc。计算得到的Cc如表4所示。引入给定温度下的等效热导率ks，由式(9)最终求得对应温度下的接触热导hc。

表4 不同压力下的无量纲接触热导

2.4 边界条件与载荷

在超声速飞行中，图2所示的局部结构同时承受高温空气的气动加热以及飞行器其他部位传递来的剪力。因此，计算中对模型施加热学边界，以模拟气动热载荷；施加力学载荷与位移边界，以模拟结构中的剪力。

模型中，载荷与边界条件的位置如图3所示。图3中，A为模型外表面，在此处给定温度边界条件；B为蒙皮与肋的右侧端面，在这两个面上约束沿法向位移；C为肋内表面、肋右侧端面与模型中间截面上的交点，在该点上给定固支边界条件；D为模型内表面，在此处给定自然对流与辐射边界条件，其中对流换热系数为20 W/(m2·K)，环境温度为 20 ℃，表面发射率为0.5；E为梁左侧端面，在该面上定义沿法向指向面外的均布拉力；F为螺栓中垂直于轴向的截面，在该面上定义沿螺栓轴向大小为 1 932.25 N 的预紧力。

A为模型表面；B为蒙皮与肋的右侧端面；C为肋内表面、肋右侧端面与模型中间截面上的交点；D为模型内表面；E为梁左侧端面；F为螺栓截面图3 模型中边界与加载位置Fig.3 Boundary and loading position in model

3 接触热阻影响分析

3.1 接触热阻对温度场的影响

3.1.1 不考虑接触热阻时的温度场

不考虑接触热阻时，应力对传热过程没有显著影响，结构温度场只与热学边界条件有关。不同外表面温度下的结构相对温度分布如图4所示。相对温度为各点摄氏温度与外表面摄氏温度的比值。图中，等温线在所有交界面上连续，表明界面间不存在接触热阻。当外表面温度升高时，结构内表面附近的相对温度略有下降。这是由于温升引起的结构热阻下降不足以抵消内表面辐射换热量上升对传热过程的影响。

图4 不考虑接触热阻时的相对温度分布Fig.4 Relative temperature distribution without contact thermal resistance

相同外表面温度，拉力分别为5、50 MPa时的结构相对温度分布如图5所示。由图5可知，即使拉力增大10倍，结构的相对温度场分布仍未出现明显差异。表明不考虑接触热阻时，结构应力场对温度场没有明显影响。

3.1.2 考虑接触热阻时的温度场

拉力为50 MPa时，不同外表面温度下的结构相对温度分布如图6所示。由于接触热阻的影响，等温线在交界面上不连续。由表5可知，在蒙皮与梁的交界面上，蒙皮内表面的最低温度大于梁上表面的最高温度,并且两者的温差随着外表面温度的升高不断上升。梁与肋的交界面间也存在不低于 20 ℃ 的温差。界面两侧的温差表明接触热阻阻碍了蒙皮、梁、肋间的热传递。同时，外表面温度越高，界面两侧的温差越大。即气动加热越强，接触热阻对结构传热的影响越显著。

图5 不考虑接触热阻时不同拉力下的相对温度分布Fig.5 Relative temperature distribution under different tensile forces without contact thermal resistance

图6 考虑接触热阻时的相对温度分布Fig.6 Relative temperature distribution with contact thermal resistance

与不考虑接触热阻时的结构温度场相比，考虑接触热阻时结构各部分温度明显下降，如表6所示，其中肋的温度下降幅度最大。对比图4、图6发现越是接近内表面的位置，考虑接触热阻时温度下降幅度越大。这是因为越接近内表面，传热路径中的接触面越多，接触热阻的对传热过程的影响更为突出。

表5 不同外表面温度下交界面两侧温度极值Table 5 Extreme temperature on both sides of interface at different outer surface temperatures

表6 外表面温度500 ℃时结构各部分温度极值

相同外表面温度下，结构承受不同拉力时，肋中螺栓孔表面的温度分布如图7所示。由图7可知，拉力增大10倍后，螺栓孔左侧温度略有上升，右侧的温度降低了约20℃。这是因为拉力的增大使得肋与螺栓在左侧接触压力上升，右侧接触压力下降，从而引起界面上接触热阻的变化，进而对温度场产生影响。表明在考虑接触热阻时，应力场通过影响结构中的接触关系，间接影响了温度场分布。

图7 考虑接触热阻时肋中螺栓孔表面温度分布Fig.7 Surface temperature distribution of bolt hole in rib with contact thermal resistance

3.2 接触热阻对应力场的影响

考虑接触热阻时，接触热阻通过影响温度场从而影响结构热应力与热变形。根据表6，接触热阻对螺栓、梁、肋的温度场影响较大。因此，接触热阻对这些部位的应力场也有较大影响，主要表现为应力分布与幅值的变化。

3.2.1 结构应力分布

结构中，螺栓孔表面的Mises应力较大。外表面温度为500 ℃，结构承受5 MPa拉力时，螺栓孔表面应力分布如图8所示。由图8可知，不考虑接触热阻时，除沉头位置的出现了应力集中外，螺栓孔表面其他部位的应力梯度较小。考虑接触热阻时，螺栓孔表面应力梯度增大，主要表现为蒙皮、梁、肋的交界部位的局部应力变化。

考虑接触热阻时，螺栓孔表面的应力分布变化是由温度场变化引起的。外表面温度为500 ℃，结构承受5 MPa拉力时，结构的变形情况如图9所示。由图4、图6、表6可知，不考虑接触热阻时，蒙皮、梁、肋间温度差低于20 ℃，沿x方向变形基本一致；考虑接触热阻时，蒙皮与肋之间最大温差为129 ℃，肋沿X方向变形量不足蒙皮的2/3，导致螺栓转动角度增加，使得螺栓与蒙皮、梁、肋间的接触位置与压力产生变化，导致螺栓孔表面的应力梯度上升。

因此，考虑接触热阻时，接触热阻通过影响结构温度场，从而影响热应力与热变形，进而影响了结构应力分布。

图8 承受5 MPa拉力时螺栓孔表面应力分布Fig.8 Stress distribution on bolt hole surface under 5 MPa tension

图9 放大30倍的结构变形Fig.9 30-fold enlargement of structural deformation

3.2.2 结构应力水平

选取结构中各部分积分点上的最大Mises应力值代表其应力水平。不同的外表面温度与拉力下，各部分的Mises应力最大值如图10所示。图10中，含(termal contact resistance)(TCR)的图例代表考虑接触热阻时的情况。

图10 结构各部分应力水平随拉力变化情况Fig.10 Variation of stress level of structural parts with tensile force

由图10可知，外表面温度较低时，接触热阻对各部分应力水平的影响并不明显。外表面温度≥300 ℃时，接触热阻引起的结构温度变化较为明显，进而引起的应力变化也较为显著。同时，随着拉力的增大，接触热阻的影响也逐渐增大。其中，肋的应力水平变化最为显著。当外表面温度为500 ℃，拉力45 MPa时，由接触热阻引起的肋应力水平变化量达到89.6 MPa，相较于不考虑接触热阻时增长了23.0%。螺栓、蒙皮、梁中的应力水平变化率较小，均低于10%。不同温度下各部件应力水平的最大变化量如表7所示。