二维两层介质界面接触热阻的反辨识研究
温晶晶, 李磊涛, 刘承骛, 吴斌     
西北工业大学 航天学院, 陕西 西安 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所示的二维热传导模型:导热系数分别为k1k2的2块材料在y=l处接触; 在接触界面附近选取一系列温度测点并测取温度值; 大多数研究忽略了接触热阻在接触区域内沿接触面空间的变化, 本文假设接触热阻在二维结构中沿接触线位置变化, 记为R(x)。建立稳态热传导控制方程为:

图 1 二维热传导模型
(1)

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

(2)
(3)
(4)
(5)

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

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

假设在y=Li处布置n个温度测点, 测得的温度记为Ti(xj)(j=1, 2, …, n), 则反问题可以表述为:用测点温度和其他边界条件反推出未知边界温度。反问题可以转化为泛函最优控制问题来求解, 优化目标函数为:

(6)

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

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

(7)

式中, Tlk+1(x)为第k+1次迭代所得的界面温度; ε为一个较小的数, 可根据具体收敛情况确定。

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

(8)

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

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

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

(9)

式中,Tlk+1(x)与(7)式中含义相同; βk, pk分别为第kk+1次迭代搜索的步长和方向。

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

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

(10)
(11)
(12)
(13)

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

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

(14)

式中, Ti[Tlk(x)-βkpk(x)]表示由第k+1次迭代所得的界面温度计算出的温度测点处的温度。

Ti[Tlk(x)-βkpk(x)]Taylor展开并取线性项代入(14)式中得:

(15)

式中,Ti[Tlk(x)]表示由第k次迭代所得的界面温度计算出的温度测点处的温度; ΔTi[pk(x)]表示由第k次迭代所得的界面温度增量计算出的温度测点处的温变化量, 它和pk(x)相关。

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

(16)

根据计算出的Ti(x)和ΔTi(x), 再结合测量值Ti(xj)可由(16)式计算出βk

2.2 伴随算法求解迭代方向

k次迭代搜索方向pk是梯度J′[Tlk(x)]和第k-1次迭代搜索方向pk-1的共轭, 计算方程为:

(17)
(18)

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

梯度J′[Tlk(x)]的求解称为伴随问题。为导出伴随问题的方程, 构造如下泛函形式:

(19)

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

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

(20)

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

(21)
(22)
(23)
(24)

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

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

(25)

由伴随函数在y方向上的导数和温度函数在y方向上的导数可以计算出J′[Tlk(x)], 代入(2)式可以求得共轭系数γk, 再结合第k-1次迭代搜索方向pk-1, 可由(17)式求解出pk

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)

并代入方程(27)得:

(28)

式中,

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

(29)

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

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

(30)

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

2.4 求解接触热阻

有了以上3个问题的解决, 就可以实现反问题的求解。按照上述方法求出上、下接触界面的温度Tlu(x)、Tld(x)和热流ql(x), 则界面接触热阻的计算公式为:

(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), 称为准确接触热阻。

图 4 计算模型

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

图 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 基于最小二乘法的不适定性处理方法

建立模型:在图 4所示模型的y=-0.05 m处再添加一系列温度测点, 并设y=0.05 m处温度测点的测量值为Tiu(x)=285+0.25sin(6πx)(单位为K, 0.25sin(6πx)为模拟测量误差), y=-0.05 m处温度测点的测量值为Tid(x)=257+0.25sin(6πx)(单位为K, 0.25sin(6πx)为模拟测量误差), 其余条件与本文3.2小节中的算例相同。

利用区域Ⅰ中测点温度的准确值计算出上界面温度, 利用区域Ⅱ中测点温度的准确值计算出下界面温度, 结合计算出的热流, 即可计算出准确接触热阻值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-5 Km2/W, 5.69×10-5 Km2/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. DOI:10.2514/3.9298
[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. DOI:10.2514/1.12294
[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. DOI:10.1016/j.ijheatmasstransfer.2005.02.033
[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. DOI:10.6052/j.issn.1000-4750.2010.12.0916 (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. DOI:10.1002/(ISSN)1521-4052
[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. DOI:10.1016/j.enganabound.2008.04.001
[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) http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=D608549
[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.
Study on Thermal Contact Resistance Estimation in Planar Medium
Wen Jingjing, Li Leitao, Liu Chengwu, Wu Bin     
School of Astronautics, Northwestern Polytechnical University, Xi'an 710072, China
Abstract: Thermal contact resistance(TCR) is one of the important parameters in heat transfer problems of engineering, and it is necessary to estimate the value of TCR effectively in many engineering fields. Considering the limitation of current estimation methods of TCR such as only focusing on one-dimensional thermal conduction, getting a single value of TCR merely, and the temperature measuring points only being placed in temperature gradient direction of mediums, boundary element method(BEM) and conjugate gradient method are combined to estimate the TCR in planar mediums. The value of TCR in relation to the position of contact interface line is estimated with this method, and the positions of temperature measuring points can be selected randomly because of the characteristic of BEM that there is no necessity to discrete the inner area and it is sufficient to discrete the boundary. The analysis of calculation examples base on heat transfer model of planar medium demonstrates that:this method can estimate the TCR effectively, but the ill-posedness is also existed in this method which is one of the inverse problems, and the calculation error of TCR is increased with the distance from temperature measuring points to contact interface, the estimation precision and stability can be improved after optimization with least square method.
Key words: thermal contact resistance     two-dimensional conduction     boundary element method     conjugate gradient method     least square method    
西北工业大学主办。
0

文章信息

温晶晶, 李磊涛, 刘承骛, 吴斌
Wen Jingjing, Li Leitao, Liu Chengwu, Wu Bin
二维两层介质界面接触热阻的反辨识研究
Study on Thermal Contact Resistance Estimation in Planar Medium
西北工业大学学报, 2018, 36(1): 28-34.
Journal of Northwestern Polytechnical University, 2018, 36(1): 28-34.

文章历史

收稿日期: 2017-04-25

相关文章

工作空间