二维两层介质界面接触热阻的反辨识研究

温晶晶, 李磊涛, 刘承骛, 吴斌

(西北工业大学 航天学院, 陕西 西安 710072)

两固体表面的接触实际只发生在一些离散点或微小面上。当热量从一个界面向另一个界面传导时，接触面附近的温度会发生突变。把这种由界面不一致接触引起的接触换热的附加阻力称为接触热阻[1]。接触热阻研究已经深入到机械制造、微电子、航空航天等诸多工程领域[2]。特别是航天领域：在真空、高低温交替的环境下，航天器部件之间的接触热阻对传热效率起着至关重要的作用，进而直接影响航天器的工作状态和寿命[3]。

目前，关于接触热阻计算方法的研究主要包括理论研究[4]、仿真计算[5-6]及试验三方面。考虑到接触传热机理复杂，影响因素多，个别物性参数定义不准确等原因，试验依然是求解接触热阻不可或缺的方法。而这种根据试验测得温度场分布，结合数据处理方法辨识出接触热阻的方法本质上属于热传导反问题范畴。其中，李建平等人[7]通过试验测得铝卷之间温度分布，再结合差分方法辨识出了铝卷间的接触热阻；刘冬欢等人[8]通过试验测得C/C复合材料和GH600高温合金之间的轴向温度值和热流密度，利用Fourier定律外推出接触界面温差，进而求得接触热阻；石友安等人[9]将灵敏度法和共轭梯度法应用到试验中，对相变材料之间的接触热阻进行了反辨识研究；El-Sabbagh、Gill等人[10-11]也分别利用共轭梯度法和灵敏度法对接触热阻进行了反辨识研究；王超[12]将顺序函数法应用到试验中，对高强度钢热成形过程中板料与模具之间的接触热阻进行了反辨识研究。总结来看，目前国内外关于接触热阻的反辨识研究都是将问题简化为两层介质之间的一维导热问题，并利用参数辨识的方法对接触热阻进行求解。该方法存在以下不足：①只能辨识出接触热阻的单一值，然而接触面附近的实际传热方式是三维的，真实的接触热阻是随接触面位置变化而变化的，特别在接触区域较大或接触面拓扑形貌变化较大的情况下必须考虑接触热阻的空间变化；②温度测点必须放置在介质内部轴线上，以最大限度模拟出一维温度场梯度变化，这就要求热电偶必须嵌入试件内部，而这会给温度测点的选择带来不便，同时也会破坏试件本身的温度场并带来测量误差，并且很难应用于尺寸较小的试件的测试[13]。

针对上述问题，本文采用边界元法对温度场进行离散求解，采用共轭梯度法进行参数辨识，二者相结合对接触热阻进行反辨识研究。本方法将问题抽象为两层介质之间的二维导热问题，可以求出随界面接触线位置变化的接触热阻；充分发挥边界元算法只需计算边界温度与热流不需计算内部温度与热流的优势，给温度测点的选取带来了很大的便利；算例分析表明：该方法可以有效求解出接触热阻，但作为接触热阻反辨识方法的一种，同样存在对温度测量误差敏感的不适定性问题，并且温度测点距离接触界面越远误差越大，采用最小二乘法对计算结果进行优化可以有效提高辨识的精度和稳定性。

1 二维热传导反问题的计算模型

1.1 建立二维热传导模型

建立如图1所示的二维热传导模型：导热系数分别为k1和k2的2块材料在y=l处接触;在接触界面附近选取一系列温度测点并测取温度值;大多数研究忽略了接触热阻在接触区域内沿接触面空间的变化,本文假设接触热阻在二维结构中沿接触线位置变化,记为R(x)。建立稳态热传导控制方程为:

图1 二维热传导模型

(1)

在不考虑热对流和热辐射的情况下,由第一类边界条件得:

(1)式～(5)式中，T=T(x,y)为温度变量;∂T/∂x=0为两侧(x=L或x=0处)的边界绝热条件;Tu,Td分别为上、下表面(y=L、y=0处)的给定温度值;Ti(x)为y=Li处的温度函数值,可由上述方程组解出。

1.2 建立目标函数和迭代停止标准

(6)

式中,Tl(x)为假定的接触界面的温度(以下简称界面温度);Ti(x)是根据T(x)由控制方程(1)计算得到的测点处的温度函数,Ti(xj)为其离散温度值。

用迭代法搜索求解界面温度Tl(x),并规定迭代停止标准为:

(7)

若考虑温度测量误差,则ε也可以写成:

(8)

式中,σ为温度测量的方均根误差。

2 二维热传导反问题的参数辨识方法

采用共轭梯度法进行接触热阻的反辨识。共轭梯度法的求解策略分为:灵敏度算法求解迭代步长、伴随算法求解迭代方向和共轭系数、边界元法求解正问题。其迭代格式为:

(9)

2.1 灵敏度算法求解迭代步长

灵敏度是指当界面温度Tl(x)有一个增量ΔTl(x)时,区域内某点温度T(x,y)的变化量ΔT(x,y)如何变化。将T+ΔT,Tl(x)+ΔTl(x)代入(1)式～(5)式中得灵敏度控制方程及相应边界条件为:

(11)式表示两侧(x=L或x=0处)的热流变化量边界条件;(12)式表示上、下表面(y=L、y=0处)的温度变化量边界条件。可由上述方程组求解出y=Li处的温度变化量ΔTi(x)。

将(9)式代入(6)式得:

(14)

(15)

对(15)式中的βk求导使其导数为0,并将其沿x坐标离散得:

(16)

2.2 伴随算法求解迭代方向

(17)

(18)

式中，γk为共轭系数,规定γ0=0。

(19)

式中，λ为Lagrange算子;δ(x-xj)为Dirac函数。

用Ti(x)+ΔTi(x)代替Ti(x),用Tl(x)+ΔTl(x)代替Tl(x),经过一系列变换可得泛函增量为:

(20)

将方程(20)右边第二项分部积分,并利用灵敏度问题的边界条件,同时令ΔJ→0,可得伴随问题:

(22)式表示侧面(x=L或x=0处)的伴随问题边界条件;(23)式、(24)式表示上、下表面(y=L,y=0处)的伴随问题边界条件。同样,可由上述方程组求出伴随函数λ(x,y)的相应值。

根据Alifanov的定义[14],由泛函增量可以求得:

(25)

2.3 边界元法求解正问题

正问题是根据已知的边界条件求解微分方程,得到未知的热流或温度,方程组(1)～(5)、(10)～(13)、(21)～(24)作为偏微分方程组的边界值问题,其求解过程均可归结为正问题求解的范畴。本文采用边界元方法求解正问题,以方程组(1)～(5)为例,首先将控制方程化为边界积分形式:

(26)

式中,Γ为边界;Ti为场点(xi,yi)处的温度值;Tj为源点处的温度值;qj为源点处的热流值;q*=-k∂G*/∂n,G*为控制方程的基本解,在二维问题中,G*=ln(1/R)/(2πk),R为场点到源点的距离;Ci为自由项系数,当(xi,yi)点位于区域内时Ci=1,当(xi,yi)点位于边界时Ci=1/2。

如图2所示,将边界Γ划分成N个单元,每个单元记为Γj(j=1,2,…,N)。

图2 二维热传导问题的边界元模型

采用常单元插值,可将(26)式写成:

(27)

(28)

将(28)式中未知量移到等式左边,已知量移到等式的右边,写成通常的方程组形式为:

[A][X]=[F]

(29)

采用四点高斯积分公式，得到数值积分的系数。式中，[A]为数值积分系数Hij、Gij组成的矩阵；[X]为未知量组成的向量；[F]为已知量组成的向量。结合已知边界条件,按(29)式即可求解出边界未知温度和热流。

如果要求区域内的温度,可以将(27)式写成:

(30)

结合已知的边界处的热流值和温度值,由(30)式即可得到区域内部的温度值。

2.4 求解接触热阻

(31)

具体计算流程如图3所示。

图3 二维接触热阻计算流程图

3 算例分析

3.1 算法验证及反问题不适定性分析

如图4所示:2块1 m×1 m的钢板在y=0处接触;侧面(x=0 m或x=1 m处)为绝热;上、下表面恒温且Tu=500 K、Td=0 K;2个区域材料的导热系数为k1=k2=14.9 W/(m·K);在y=0.05 m处有一系列温度测点,设Ti(x)=282.87+0.25sin(6πx)(单位为K),其中0.25sin(6πx)为模拟的测量误差,可见测量误差不超过0.09%;假设界面接触热阻为:R(x)=0.001 5-0.004x+0.004x2+0.000 05·

sin(4πx)(单位为K·m2/W),称为准确接触热阻。

在不考虑温度测点误差的情况下,利用本文第2节的参数辨识方法可以计算出上界面温度,再通过界面接触热阻反推出下界面温度;同理,在考虑温度测点误差的情况下,也可以计算出上、下界面温度。利用温度测点无误差情况下计算出的上界面温度和温度测点有误差情况下计算出的下界面温度可以计算出考虑温度测点误差的接触热阻,称为计算接触热阻,与准确接触热阻的对比如图5所示。

图4 计算模型

图5 准确接触热阻和计算接触热阻的对比图

在离散化的准确接触热阻和计算接触热阻中各取n个值,求取平均相对误差,计算公式为:

(32)

式中,R*(xi)为计算接触热阻离散值;R(xi)为准确接触热阻离散值。

再由Bessel公式可得计算接触热阻的方均根误差计算公式为:

(33)

计算可得:v=6.15%,σ=5.31×10-4Km2/W。该算例表明:本文提出的边界元法和共轭梯度法相结合的参数辨识方法可以解决接触热阻反问题;但反问题的求解结果对温度测点处的误差比较敏感,很小的输入误差(不到0.09%)会导致较大的输出误差(6.15%),并且输出结果的稳定性也不高(方均根误差较大),即反问题的求解存在不适定性。

3.2 测点位置对求解精度的影响

结合上述模型,分别在距离上界面0.05 m,0.10 m,0.15 m,0.20 m,0.25 m,0.30 m处布置温度测点,并计算接触热阻及其相应的平均相对误差,如图6所示。

图6 接触热阻平均相对误差随测点位置y变化的曲线

由图6结果可知:测点位置距离接触界面越远,接触热阻计算误差越大,并且相对误差增加幅度也越大,该结论与文献[13]中结论一致。因此,在靠近接触界面处布置温度测点并读取温度测量值可以提高接触热阻计算精度。

3.3 基于最小二乘法的不适定性处理方法

利用区域Ⅰ中测点温度的准确值计算出上界面温度,利用区域Ⅱ中测点温度的准确值计算出下界面温度,结合计算出的热流,即可计算出准确接触热阻值R(x);同理,可以计算出考虑测量误差的计算接触热阻值R*(x)。下面采用最小二乘法对计算接触热阻进行优化,得到修正接触热阻值Rcr(x):

在不失一般性的情况下,不妨假设:

(34)

式中,αn为未知系数。定义:

(35)

分别求(35)式对α1,…,αn的偏导,并令结果为0,可得到如下方程组:

(36)

求解(36)式即可解出未知的系数αn,进而代入(34)式求出Rcr(x)。本文取n=6,求得:

(37)

对比R(x)、R*(x)、Rcr(x),如图7所示。

图7 R(x),R*(x),Rcr(x)对比图

利用(32)式、(33)式计算可得:R*(x),Rcr(x)的平均相对误差分别为5.06%,2.96%;二者的方均根误差分别为6.11×10-5Km2/W,5.69×10-5Km2/W。可见经过最小二乘法处理后的接触热阻的精度得到了提高,并且计算结果的稳定性也有所提高。

4 结 论

1) 考虑到大多数研究忽略了接触热阻随接触面空间的变化,本文采用边界元法和共轭梯度法相结合的方法反辨识出二维热传导问题中沿接触线位置变化的接触热阻。对于接触区域较大或接触面拓扑形貌变化较大等不宜将模型简化为一维热传导问题的情况,本方法可以计算出更准确的接触热阻;

2) 考虑到边界元算法不需要对内部区域进行离散的特点,本方法可以全域选择温度测点;

3) 算例分析表明,本方法可以有效辨识出接触热阻,但计算结果对温度测量误差非常敏感,并且温度测点离接触界面越远敏感程度越高;

4) 利用最小二乘法对接触热阻的辨识结果进行优化处理,提高了接触热阻的辨识精度和稳定性;

5) 本方法基于二维两层介质模型建立,具有一定的代表性,可以用于类似的接触热阻测试工作,为计算接触热阻提供了新的思路;

6) 可以更进一步将本方法推广至三维接触热阻的分析情况,计算出更精确的随接触界面位置变化的接触热阻。

[1] Madhusudana C V, Fletcher L S. Contact Heat Transfer——The Last Decade[J]. AIAA Journal, 1986, 24(3):510-523

[2] 王安良, 赵剑锋. 接触热阻预测的研究综述[J]. 中国科学: 技术科学, 2011, 41(5):545-556

Wang Anliang, Zhao Jianfeng. The Research Overview of Prediction of Thermal Contact Resistance[J]. Scientia Sinica Techologica, 2011, 41(5):545-556 (in Chinese)

[3] Ramamurthi K, Kumar S S, Abilash P M. Thermal Contact Conductance of Molybdenum-Sulphide-Coated Joints at Low Temperature[J]. Thermophys Heat Transf, 2007, 21(4):811-813

[4] Bahrami M, Yovanovich M M, Culham J R. Thermal Contact Resistance at Low Contact Pressure: Effect of Elastic Deformation[J]. International Journal of Heat and Mass Transfer, 2005, 48(16):3284-3293

[5] 刘冬欢, 王飞, 曾凡文,等. 高温接触热阻的有限元模拟方法[J]. 工程力学, 2012, 29(9):375-379

Liu Donghuan, Wang Fei, Zeng Fanwen, et al. Finite Element Simulation Method of High Temperature Thermal Contact Resistance[J]. Engineering Mechanics, 2012, 29(9):375-379 (in Chinese)

[6] 李磊涛, 郑晓亚. 一种瞬态接触热导数值模拟方法[J]. 机械工程学报, 2016, 52(6):174-180

Li Leitao, Zheng Xiaoya. A Numerical Simulation Method of Transient Thermal Contact Conductance[J]. Journal of Mechanical Engineering, 2016, 52(6):174-180 (in Chinese)

[7] 李建平, 刘涛, 王伯长,等. 铝卷径向接触热阻的仿真研究[J]. 金属热处理, 2009, 34(4):91-95

Li Jianping, Liu Tao, Wang Bochang, et al. Simulation of the Radial Contact Thermal Resistance for the Aluminum Coil[J]. Heat Treatment of Metals, 2009, 34(4):91-95 (in Chinese)

[8] 刘冬欢, 郑小平, 黄拳章,等. C/C复合材料与高温合金GH600之间高温接触热阻的试验研究[J]. 航空学报, 2010, 31(11):2189-2194

Liu Donghuan, Zheng Xiaoping, Huang Quanzhang, et al. Experimental Investigation of High-Temperature Thermal Contact Resistance between C/C Composite Material and Superalloy GH600[J]. ACTA Aeronautica et Astronautica Sinica, 2010, 31(11): 2189-2194 (in Chinese)

[9] 石友安, 桂业伟, 杜雁霞,等. 相变材料热控系统内部接触热阻的辨识方法研究[J]. 实验流体力学, 2012, 26(4):54-58

Shi You′an, Gui Yewei, Du Yanxia, et al. Study on Thermal Contact Resistance Estimation in the Phase-Change Material Thermal Control Device[J]. Journal of Experiments in Fluid Mechanics, 2012, 26(4):54-58 (in Chinese)

[10] El-Sabbagh A M, Fieberg C, El-Marg E, et al. Modeling of Transient Thermal Contact Resistance out of Conjugate Gradient Method[J]. Materialwissenschaft und Werkstofftechnik, 2007, 38(4):288-293

[11] Gill J, Divo E, Kassab A J. Estimating Thermal Contact Resistance Using Sensitivity Analysis and Regularization[J]. Engineering Analysis with Boundary Elements, 2009, 33(1):54-62

[12] 王超. 高强钢热成形接触导热和零件力学性能及工艺优化研究[D]. 武汉:华中科技大学, 2014

Wang Chao. Investigation on Thermal Contact Behavior and Prediction of Mechanical Properties and Process Optimization in Hot-Stamping[D]. Wuhan, Huazhong University of Science and Technology, 2014 (in Chinese)

[13] 张平, 宣益民, 李强. 界面接触热阻的研究进展[J]. 化工学报, 2012, 63(2):335-349

Zhang Ping, Xuan Yimin, Li Qiang. Development on Thermal Contact Resistance[J]. CIESC Journal, 2012, 63(2):335-349 (in Chinese)

[14] Alifanov O M. Solution of an Inverse Problem of Heat Conduction by Iteration Methods[J]. Journal of Engineering Physics and Thermophysics, 1974, 26(4): 471-476 (in Chinese)