引言

当两个粗糙表面相互接触时,其真实的接触仅发生在一些离散的点或微小面积上,其余的间隙部分是真空或介质,即使在10 MPa量级的接 触压强下,真实的接触面积仅占 名义接触面积的1%~2%^[1,2]. 这种接触状况引起热流收缩,使两个相互接触的表面间的温度产生一个阶跃变化,从而产生接触热阻,因此,将接触热阻R定义为接触 面温度变化值和热流的比值^[1],即

$R\left( {x,y} \right) = \dfrac{\Delta T\left( {x,y} \right)}{q\left( {x,y}\right)} (1)$

式中, ΔTx,y为接触面的温度变化值, qx,y为界面热流值.

接触热阻的研究包括正问题和反问题两个方面,其中正问题研究包括理论研究和试验研究^{[3,4,5,6,7,8,9]}. 目前,国内外学者对接触热阻的理论研究主要集中在正问题研究,而对反问题的研究比较少. 正问题中的理论研究是建立接触热阻的解析模型,具有代表性的接触热阻模型有CMY模型、GW模型、BGT模型、WA模型、MT模型等^[2]. 之后一些学者对接触面形状进行模拟,将数值求解与粗糙表面形貌描述相结合,发展了接触热阻的数值模型^{[10,11,12,13,14,15,16,17,18,19,20,21]}. 尽管接触热阻可以通过理论进行预测,但是试验研究仍是一种更快速准确的方法^[7,9,22]. 然而,目前的试验着重研究接触热阻随接触压强、界面温度、表面粗糙度等因素的变化规律^{[23,24,25,26,27,28]}.

接触热阻反问题研究是指依据部分边界条件和部分测量点的温度值,通过优化算法反演计算得到接触面的温度和热流,进而得到接触热阻. 相对于接触热阻的正问题研究,它可作为实际工程中预测接触热阻的一个简洁有效的方法.

目前求解反问题主要依靠数值方法^[29],其数值方法主要有有限差分法、有限体积法、有限元法和边界元法4种^{[30,31,32,33]}. 相比较于其他数值方法,边界元法具有简便且计算量小的特点,不需要计算区域内的温度和热流,而只需求解边界温度和热流,因此,对于实际问题求解具有较大优势^[32,33]. 同时,反问题的计算还需要结合优化迭代算法来反演计算得到所求信息^[29-30,32],在众多优化算法中,由于共轭梯度法具有简单,计算量小等优势^[30,31,32]. 因此,本文将边界元法和共轭梯度法结合起来,研究二维空间随坐标变化的接触热阻的反问题,通过算例验证方法准确性和可行性. 同时,由于求解反问题过程中不可避免其不适定性,即对误差敏感性^[30,31,32],为处理该问题,本文采用最小二乘法(LSM),以避开反问题结果中一些偏离实际值较大的测量点,使得问题的求解结果与真实值更接近.

1 接触热阻反问题求解

1.1 正问题求解

计算模型域 Ω内

$\dfrac{\partial ^2T}{\partial x^2} + \dfrac{\partial ^2T}{\partial y^2} = 0 (2)$

边界条件为

$\dfrac{\partial T}{\partial x} = 0$, 在 x=0和 x=1.0处

T=T1, 在 y=-1.0处

T=T2, 在 y=1.0处

T1(x,y)=Ti, 在 y=0.05

T2(x,y)=Tj, 在 y=-0.05

本文采用边界元法求解上述偏微分方程,因此上式可将其转化为边界点的积分方程^[30,32-33],即

CiTi+∫Γj‍Tjq*dΓ=∫Γj‍qjG*dΓ(3)

式中, Γ为域 Ω的边界; T为温度; q为热流密度; G*为基本解,二维情况下,为$\dfrac{1}{2\pi }\ln \dfrac{1}{r}$; q*为 G*的法向导数, i为 Ω内或 Γ上的点,当 i为 Ω区域内的点时, Ci=1,当 i为光滑边界 Γ上的点时, Ci=0.5. 将边界离散后得到

${ C}_i { T}_i + \sum_{j = 1}^n { T}_j \int_{\varGamma _j } { q}^\ast d\varGamma = \sum_{j = 1}^n { q}_j \int_{\varGamma _j } { G}^\ast d\varGamma $(4)

记

Ĥij=∫Γj‍q*dΓ(5)

Gij=∫Γj‍G*dΓ(6)

则式(4)可以写成

${ C}_i { T}_i + \sum_{j = 1}^n \hat { H}_{ij} { T}_j = \sum_{j = 1}^n { G}_{ij} { q}_j (7)$

式(7)最终可化为如下形式的标准边界元离散方程

${\pmb H}\left( {x,y} \right){\pmb T}\left( {x,y} \right) = {\pmb G}\left( {x,y} \right){\pmb q}\left( {x,y} \right) (8)$

式中,H和G是边界元离散矩阵,T和q分别是温度和热流的列向量. 根据已知的边界条件,求解上式可以得到边界温度和热流.

1.2 反问题目标函数

在接触热阻反问题研究中,接触界面温度与热流是未知的. 但其他边界温度与热流是已知的,可选择在接近接触界面适当的位置布置温度测量点. 假设在N个位置布置测量点,记测量点的温度为$\bar{\pmb T}_i (x)$ $ \left( {i = 1,2, \cdots, N} \right)$,则反问题研究就可以表示为：用测点温度和部分已知的边界条件来反推出接触界面上未知的温度和热流密度,进而获得接触热阻值.

反问题求解可以转化为以下目标函数泛函变分的最优化问题^{[29,30,31,32]},即

$J({ T} (x)) = \left\| {{ T}_i - \bar{ T}_i } \right\|^2 = \sum_{i = 1}^N {\left[ {{ T}_i (x) - \bar{T}_i \left( x \right)} \right]} ^2 (9)$

式中, T̅i(x)是在模型测量点上实际测量的温度值, Ti(x)是根据估计的初始边界温度 T0反演计算得到的测量点处的温度值.

用共轭梯度法迭代来搜索边界温度 Ti(x,y)涉及迭代停止的标准问题,并同温度测量的误差有关,在没有测量误差的情况下,停止标准为

J[Tk+1(x)]<ε(10)

式中, ε是一个较小的正数. 存在测量误差的情况下,假设温度测量的标准差为 σ,在测量点加上测量误差后,温度差值近似可表示为

TiR-T̅i≈σ(11)

则 ε可以写成

$\varepsilon = \sum_{i = 1}^N {\sigma ^2} = N\sigma ^2 (12)$

根据计算值 Ti(x)和测量值 T̅i(x),计算反问题目标函数式 (3),判断目标函数结果是否满足迭代停止标准,进而判断是否还需要进行迭代.

1.3 共轭梯度法

求解反问题所涉及的方法较多,鉴于共轭梯度法的简单性和高效性^[29,30],本文采用共轭梯度法来优化目标函数. 未知 的边界温度 Ti(x)是初始的估计值或迭代过程中通过搜索得到的,即

$\hat{\pmb T}^{n + 1}\left( x \right) = \hat{\pmb T}^n\left( x \right) - \beta ^n {\pmb p}^n\left( x \right) (13)$

式中, βn为 n到 n+1次迭代的搜索步长; pnx为搜索方向,即

${\pmb p}^n\left( x \right) = {\pmb J}'^n\left( x \right) + \gamma ^n{\pmb p}^{n -1}\left( x \right) (14)$

式中, ${\pmb J}'^n\left( x \right)$为梯度方向,因此第 n次的搜索方向 pnx是梯度方向${\pmb J}'^n\left( x \right)$和第 n-1次的搜索方向 pn-1x的线性组合,其中 γn为共轭系数,可由下式计算得到

$\gamma ^n = \dfrac{\int_{ 0}^L {\left( {{\pmb J}'^n} \right)^2 x} }{\int_{ 0}^L {\left( {{\pmb J}'^{n -1}} \right)^2 x} } (15)$

式中, γ0=0. 从式(14)可知,当 γn=0时,无论 n为何值, pnx变成梯度方向.

灵敏度方程可由控制方程(2)求其增量得到

$\dfrac{\partial ^2\Delta T}{\partial x^2} + \dfrac{\partial ^2\Delta {T}}{\partial y^2} = 0 (16)$

采用边界元法求解灵敏度方程,即可得到式(13)中的搜索步长 βn

$ \beta ^n = \dfrac{\sum_{i = 1}^n {\left( {{ T}_i - \hat{ T}_i } \right)\Delta {T}_i } }{\sum_{i = 1}^n {\Delta { T}_i ^2} } (17)$

伴随方程可由控制方程(2)的泛函求导得到

$\dfrac{\partial ^2\lambda }{\partial x^2} + \dfrac{\partial ^2\lambda}{\partial y^2} = 0 (18)$

同样采用边界元法求解伴随方程,结合正问题的计算结果,得到泛函 J(T(x))的导数为

${J}'(T(x)) = - \dfrac{\partial \lambda }{\partial y}\dfrac{\partial {\pmb T}}{\partial y}\left| {_{y = 0} }\right. (19)$

其中, λx,y为拉格朗日算子. 结合正问题得到的温度在 y方向上的导数 ∂T/∂y,可计算出来泛函 J(T(x))的导数 J'(T(x)). 根据式(15),可计算出共轭系数 γn,结合 n-1次的搜索方向 pn-1x,就可求解新的搜索方向 pnx. 由式(13)可知,依据灵敏度方程求解得到的搜索步长 βn,就可计算出新的边界温度函数 T̂n+1x.

2 接触热阻反问题方法验证

为了验证接触热阻反问题方法的准确性和可行性,假定真实接触热阻的函数形式已知,并且还获得部分测点温度,采用边界元法和共轭梯度法反演计算得到界面的温度和热流,进而得到接触热阻值.

2.1 计算模型

考虑一个二维的稳态情况的接触热阻问题. 计算模型如图1所示,是一个 1 m × 1 m和1 m × 1 m的模型,在 y=0处两模型接触^[9]. 模型的左右边

显示原图| 下载原图ZIP| 生成PPT

图1   接触热阻模型

Fig. 1   Mode of thermal contact resistance

界为绝热,上下表面保持温度不变, Tu=500K,下边界 Tb=0K,两个区域材料为导热系数相同 λ=14.9W/(m ⋅K)不锈钢,环境温度 T=293K,在接触界面上方 y=0.05m处设置一系列温度测点,而在接触面位置假定接触热阻为一个已知函数.

2.2 边界条件

计算模型域 Ω内

$\dfrac{\partial ^2T}{\partial x^2} + \dfrac{\partial^2T}{\partial y^2} = 0$

边界条件为

$\dfrac{\partial T}{\partial x} = 0$, 在 x=0和 x=1.0处

T=0K, 在 y=-1.0处

T=500K, 在 y=1.0处

内部测量点的温度值为

T1(x,y)=285+0.25sin6πx,y=0.05

式中, T=285K为准确值, 0.25sin6πx为测量误差,以模拟实际测量等造成的误差. 为刻画其分布情况,可采用通常的函数形式,如多项式或三角函数,这里,采用三角函数形式. 其相对于准确值,测量误差很小,约 ξ=0.09%.

假定真实的接触热阻随空间的变化函数为

Rx=0.0015-0.004x+0.004x2+0.00005sin4πx

2.3 计算结果

首先考虑没有测量误差的情况,反演获得界面的接触热阻如图2所示.

从图2可以看出,在没有测量误差的情况下,边界元法和共轭梯度法可以准确计算出接触热阻值,与真实值没有差别,说明了所采用 的方法在没有测量误差的情况下,计算结果的准确性.

显示原图| 下载原图ZIP| 生成PPT

图2   计算结果比较

Fig. 2   Comparison of calculation results

现在考虑有 ξ=0.25sin6πx的误差情况,计算得到接触热阻值,为了反映出结果的准确性,定义以下各类误差.

将计算结果的相对误差定义为

$\varepsilon = \left| {\dfrac{R_i ^\ast - R_i }{R_i }} \right|\times 100\% (20)$

将计算结果的平均相对误差定义为

$\bar {\varepsilon } = \dfrac{1}{N}\sum_{i = 1}^N {\left|{\dfrac{R_i^\ast - R_i }{R_i }} \right|} \times 100\% (21)$

将计算结果的标准差定义为

$\sigma = \sqrt {\dfrac{1}{N - 1}\sum_{i = 1}^N {\left({R_i^\ast - R_i } \right)^2} } (22)$

式中, Ri*为计算值, Ri真实值.

接触热阻的真实值和考虑有测量误差情况下的计算值对比如图3所示. 相对误差沿着 x轴变化曲线见图4.

显示原图| 下载原图ZIP| 生成PPT

图3   计算结果比较

Fig. 3   Comparison of calculation results

显示原图| 下载原图ZIP| 生成PPT

图4   相对误差曲线

Fig. 4   Relative error curve

由图3和图4可以看出,正问题的测量误差会随着反问题的求解而被放大. 在 x=0.325m,相对误差最大量值为20.23%;在 x=0.225m,相对误差最小为0.43%,平均的相对误差 ε̅=6.15%,标准差为 σ=5.31×10-4K ⋅m 2/W,其结果的相对误差的最大值和最小值与测量误差相比,分别被放大了约225倍和5倍. 本节算例说明本文所采用的边界元法和共轭梯度法可以解决接触热阻反问题,但反问题的求解对测量误差比较敏感,结果会因测量输入误差而产生较大偏差.

3 接触热阻反问题算例

3.1 基于最小二乘法的结果校正

由2.3节可知,由于测量误差的存在,若采用上述方法直接求解,可能导致计算结果 不准确^[30,32]. 这时就需要采用特殊 的校正方法,来校正因测量等造成的误差对计算结果的影响,本文采用最小二乘法来优化计算结果.

为了简化问题,假设接触热阻R是x的连续函数R(x). 为了确定未知的接触热阻R(x),又不失一般性,可假设

$ R(x)=\sum^n_{n=1} \alpha_n x^n $(23)

其中, αn是未知的系数.

将最小二乘法准则 Sαn定义为

$S(\alpha_n)=\sum^N_{i=1} ( \bar R_i(x,y)-R)^2 $(24)

对于其中的未知系数 αn,可通过最小化 Sαn得到,对 Sαn进行微分,使其结果等于零,即

$ \dfrac{\partial S}{\partial \alpha_n} =\sum^n_{i=1} \dfrac{\partial R}{\partial \alpha_n}(\bar R_i -R) =0 (25)$

再结合反问题得到的结果 R̅i,由式(25)就能确定未知的系数 αn,根据式(23)计算接触热阻 R(x).

3.2 计算模型

实际工程中,通常在接触界面两边布置温度测点,反演出接触热阻. 仍然采用2.1节的计算模型,在距离接触界面 ±0.05 m处布置温度测点,具体形式如图5所示.

显示原图| 下载原图ZIP| 生成PPT

图5   接触热阻模型

Fig. 5   Mode of thermal contact resistance

3.3 边界条件

计算模型域 Ω内

$\dfrac{\partial ^2T}{\partial x^2} + \dfrac{\partial ^2T}{\partial y^2} = 0$

边界条件为

$\dfrac{\partial T}{\partial x} = 0$, 在 x=0和 x=1.0处

T=0K, 在 y=-1.0处

T=500K, 在 y=1.0处

内部测量点的温度值为

T1(x,y)=285+0.25sin6πx, y=0.05

T2(x,y)=257+0.25sin6πx, y=-0.05

式中, 0.25sin6πx同样为测量误差,其相对于准确值测量误差很小,约 ξ=0.09%.

3.4 计算结果

由2.3节可知,在没有测量误差情况下,采用边界元法和共轭梯度法反演获得界面的真实接触热阻,见图6中标记符为五角星的 曲线所示. 有测量误差情况下,采用边界元法和共轭梯度法反演获得的接触热阻见图6中标记符为三角形的曲线所示.

显示原图| 下载原图ZIP| 生成PPT

图6   计算结果比较

Fig. 6   Comparison of calculation results

基于最小二乘法对考虑测量误差反演出的接触热阻值进行校正,取 n=6. 由最小二乘法计算出来的接触热阻函数为

$ R\left( x \right) = 0.001 - 0.000 3x - 0.010 9x^2+ \\ 0.078 4x^3 - 0.197 8x^4 + 0.213 5x^5 - 0.083 3x^6 (26)$

校正后的接触热阻见图6中标记符为圆圈的曲线所示.

校正处理前后部分点的相对误差的减小量如表1所示,相对误差沿着 x轴变化曲线如图7所示.

显示原图| 下载原图ZIP| 生成PPT

图7   相对误差曲线

Fig. 7   Relative error curve

从表1可以看出,在一些位置上,经过最小二乘法优化后,其计算结果的相对误差下降明显,与真实值更接近. 由图7可以看出,采 用最小二乘法优化后,整体的相对误差有所下降. 在处理之前平均相对误差为 ε̅=5.06%, σ=6.11×10-5K ⋅m 2/W,处理后平均相对误差为 ε̅=2.96%, σ=5.69×10-5K ⋅m 2/W,经过最小二乘法处理后,平均相对误差减小41.5%,标准差减小2.45%. 最小二乘法拟合的曲线既反映接触热阻的变化趋势,还消除了其局部波动,避开了计算结果中偏离真实值较大的测量点,使其总体误差进一步减小,计算结果更准确. 本节算例说明了最小二乘法可以作为处理不适定性问题的一种有效方法,最小二乘法优化算法可以使误差减小,使计算结果较不适定处理前的计算值更接近真实值,显著提高了计算结果的准确性.

4 结 论

本文将边界元法和共轭梯度法结合起来,研究了二维空间随坐标变化的接触热阻反问题,具有一定的实际参考价值,所得结论如下：

(1)在接触热阻的反问题识别中,采用边界元法和共轭梯度法可求解反问题. 由于边界元法不需要知道和计算内部区域的温度和热流,因而大大减少了计算量.

(2)在无测量误差的情况下,本文采用的方法准确可行,可以真实地反演计算得到接触面上的接触热阻值.

(3)反问题对测量误差较敏感,存在误差的情况下,计算结果会偏离真实值. 且正问题的误差会随反问题的求解而被放大,使得反问题求解的数值结果误差较大.

(4)在正问题中的误差在反问题中被放大的情况下,为处理不适定性问题,采用了最小二乘法,使曲线避开反问题计算结果中一些偏离实际值较大的点,显著提高了反演计算结果的准确性.

The authors have declared that no competing interests exist.

参考文献

文献选项

• 原文顺序
• 文献年度倒序
• 文中引用次数倒序
• 被引期刊影响因子