张 琳,冯妍卉,王智学

(1.中国航发四川燃气涡轮研究院， 成都 610500;2.北京科技大学 能源与环境工程学院， 北京 100083)

随着传热技术的不断发展和广泛应用，接触导热逐渐成为一个不可忽略的影响传热的因素，成为航空航天、机械制造、电子、低温超导等学科和动力工程领域中的研究热点。接触导热是一个受材料热物性、机械特性、表面形貌及接触环境(压力、温度和间隙填充介质等)等众多因素影响的非线性问题。在工程应用中，通常需要降低或提高接触热阻以控制接触传热，在设计过程中需要采用理论、实验或数值分析的方法对接触热导率或接触热阻进行预测。但在不同的应用中，材料种类繁多，应用环境也不同，用实验手段来获取接触热导率通常需要设计复杂的实验装置，具有较高的成本。近年来，由于计算机技术的发展，数值分析的方法越来越多地应用于接触导热的研究中[1-3]，但此类方法通常需要进行大量的前、后数据处理工作，时间成本较高。因此，工程领域通常采用理论模型和相应的计算公式对接触热导率进行预测。国内外学者对接触导热进行了大量的研究，产生了许多具有重要参考价值和工程应用价值的理论和方法。接触导热的理论模型和计算方法的研究包括对粗糙表面的形貌分析、接触形变分析和接触导热分析，研究者将不同的粗糙表面形貌模型、形变模型和导热模型进行组合，形成了不同的接触导热计算模型。王安良等[4]对接触导热的理论预测方法进行了较为全面的综述。

早期的接触导热计算模型主要以统计学参数描述固体粗糙表面形貌，认为固体粗糙表面是由无数个均匀分布的微凸体堆积而成，这些微凸体的高度符合高斯分布。以此为基础，许多学者建立了接触导热的计算公式，如Mikic[5]、Yovanovich等[6]、陈剑楠等[7]。此外，Kumar等[8]建立了接触热阻的蒙特卡洛模型，该模型根据粗糙表面粗糙峰分布符合高斯分布的特性，在粗糙表面上随机模拟N个符合高斯分布的粗糙峰，并在不同的接触界面间距下计算出发生接触的粗糙峰数目，每个发生接触的粗糙峰都形成一个单点接触热阻，界面间总的接触热阻即为这些单点接触热阻的并联。

需要指出的是，表面粗糙高度的分布并不完全符合高斯分布，以上所述接触导热计算模型都建立在以统计学参数描述的表面形貌模型的基础上，而研究表明，表面形貌的统计学参数很明显地受到测量仪器的分辨率和取样长度的影响[9]，不能唯一确定地表征一个粗糙表面，在此基础上建立的接触导热计算模型通常具有较大的不确定性，计算结果的准确性与工程应用需求差距较大。因此，研究者们开始寻求与仪器分辨率和取样长度无关的粗糙表面表征参数，而分形理论使得粗糙表面的唯一表征成为可能。分形理论采用反映表面不规则程度的分形维数D和反映表面轮廓幅值的比例参数G来确定一个表面轮廓，分形参数D和G都是与仪器分辨率和取样长度无关的参数。Majumdar等[9]采用修正的Weierstrass-Mandelbrot函数(W-M函数)描述粗糙表面轮廓，建立了固体接触界面的弹塑性分形接触模型和接触导热的M-T分形模型。Warren等[10]采用Cantor集分形形貌模型模拟粗糙表面，建立了弹性-塑性接触和完全塑性接触的Cantor集接触模型和接触导热模型。此外，赵兰萍[11]、徐瑞萍[12]、Zou等[13]、Ji等[14]、马丽娜[15]和李小彭等[16]也进行了基于表面分形理论的接触导热计算模型研究。

一系列研究表明，基于表面分形理论的接触导热计算模型，不受仪器分辨率和取样长度的影响，对接触导热的预测具有确定性，能够较为准确地预测接触导热。近年来将分形理论应用于接触导热预测的研究较少，现有的分形理论和模型较为复杂，与实际工程应用存在很大差距[17]。因此，本文提出了一种适合于工程应用的接触导热分形计算模型。模型基于固体粗糙表面形貌的随机特性和分形特性，采用蒙特卡洛法和W-M分形函数模拟表面形貌和接触，构建了固体界面间接触热阻的并联网络模型，简化了接触热阻的计算过程，模型计算结果与实验结果基本一致，能够较为准确地预测固体界面间的接触导热，具有较高的工程应用价值。

1 固体粗糙表面形貌的分形表征

固体粗糙表面形貌是一个非稳定随机过程，如图1所示[9]。表面轮廓曲线具有随机、多尺度和无序的特性以及统计自仿射和自相似的数学特性，分形几何学中的W-M函数可以满足表面轮廓曲线的所有数学特征，其表达式如下：

(1)

式中：D是分形维数，它反映的是轮廓z(x)在所有空间尺度上的不规则性；G是决定轮廓高度幅值的比例系数，决定了z(x)的具体尺寸；γn为空间频率，曲线的最低频率取决于取样长度L，并有γn1=1/L。对于任意取样长度，决定轮廓高度z(x)的是分形参数D和G。

图1 表面形貌的统计自仿射和自相似特性

由式(1)可以得出

z(γx)=γ(2-D)z(x)

(2)

若位置坐标x放大γ倍，则轮廓高度放大为原来的γ(2-D)倍。粗糙表面的分形参数D和G是尺度独立的，可以提供存在于分形表面上的所有尺度范围内的全部粗糙度信息[10]。

图2为某粗糙表面实测形貌与采用分形分析后得到分形参数并采用W-M分形函数模拟的形貌对比，可以看出二者表面形貌曲线波动幅度和频率具有较高的符合性，W-M函数对固体粗糙表面的形貌具有较好的模拟效果。

图2 W-M函数模拟形貌与真实表面形貌对比

2 基于蒙特卡洛法和分形理论的接触热阻计算

2.1 基本假设

本文研究在一定压力下相接触的固体表面间的接触热阻，基本假设如下：

1) 粗糙表面轮廓具有随机性和分形特性，固体表面之间的接触可等效为一个当量粗糙表面和刚性光滑平面之间的接触。

2) 固体表面间的接触和传热仅仅发生在离散分布的一系列大小不同的接触点上，接触点发生塑性变形，忽略间隙介质的导热和热对流以及非接触部分的辐射换热。

3) 表面间各接触点的传热互不影响，每个接触点相当于一个热流通道，并形成一个单点接触热阻，总的接触热阻是所有单点接触热阻的并联。

2.2 接触界面形貌的模拟

接触界面形貌模拟的目的是确定粗糙表面接触点的数量、尺寸及分布等参数。粗糙表面轮廓本质上是一个随机过程，因此本文将分形理论与蒙特卡洛法相结合，对表面粗糙峰的分布进行模拟。

首先将2个粗糙表面等效成一个当量粗糙表面，将粗糙表面的接触等效为当量粗糙表面与刚性光滑平面的接触，从而求出当量表面分形参数。一条数字化轮廓曲线可以视为一个时间序列，在z(x)曲线上取时间间隔为Δτ的N个采样点，令τ=nΔτ，则该曲线的结构函数为

E(τ)=〈[z(x+nΔτ)-z(x)]2〉=

(3)

等价于

E(τ)=CG2(D-1)τ4-2D

(4)

其中

(5)

式中函数Γ为Gamma函数，当1

在双对数坐标上，结构函数E(τ)与τ呈线性关系，即

lgE(τ)=B+klgτ

(6)

式中：

k=4-2D

(7)

B=lgCG2(D-1)

(8)

对于当量粗糙表面，其结构函数为2个粗糙表面结构函数之和，即

E(τ)=E1(τ)+E2(τ)

(9)

当量硬度H、当量弹性模量E和当量热导率k通过以下各式求得：

H=min(H1,H2)

(10)

(11)

(12)

式中ν为材料的泊松比，下标1和2分别代表相接触的2个表面。

研究表明，表面粗糙峰均方根高度σ与分形参数具有如下关系[9]：

(13)

表面粗糙峰高度具有随机性，并符合W-M函数，因此对于一个边长为L的方形接触表面，可随机选取N个x值(0

N=L2η

(14)

式中η为表面粗糙峰密度。根据Hsieh[20]的结果，表面粗糙峰密度

(15)

式中m为表面粗糙峰平均斜率。

2.3 表面接触参数的计算

本文采用的接触点形变模型为塑性圆锥模型，如图3所示。假设当量粗糙表面与刚性光滑平面的距离为d，则高度为zi的粗糙峰的峰顶与刚性光滑平面的距离δi为

δi=zi-d

(16)

当δi>0时，粗糙峰与刚性光滑平面发生接触，其接触半径ai、接触面积Aci和接触载荷Fi分别为：

ai=δi/m

(17)

(18)

(19)

通过对所有粗糙峰的计算，可以得到发生接触的粗糙峰的个数n，进而计算出表面距离为d时的实际接触面积Ar和总接触压力Fc：

(20)

(21)

上述计算过程是在已知d的情况下求界面间接触载荷Fc等参数，然而在进行实验研究或实际工程应用中很难知道接触界面的实际距离，但很容易知道接触载荷Fc，因此在接触热阻的理论计算中可迭代求得d，进而求出其他参数。

图3 塑性圆锥形变模型

2.4 单点接触热阻和总热阻的计算

研究表明，半无限长流管换热模型更适用于接触热阻的计算，Yovanovich[6]采用的单点接触热阻计算公式为

(22)

式中

ψ=(1-ε)1.5

(23)

并且有

(24)

式(24)中Aa为名义接触面积，

Aa=L2

(25)

根据本文的基本假设可知，总热阻为各单点接触热阻的并联，因此可以得出总接触热阻Rc和接触热导率hc：

(26)

(27)

3 结果与分析

3.1 计算方法的验证

本文采用文献[11]中的2组实验数据对接触热阻的计算模型进行了验证，并与Mikic模型、Yovanovich模型和MT分形模型进行了对比。实验材料为铝合金Al5052，实验温度为155 K，表1和表2分别为材料Al5052的机械性能和热力学性能参数以及4个Al5052试样的表面分形参数。

图4和图5分别为理论计算结果与文献[11]中的实验结果的对比，可以看出：Yovanovich模型的接触热导率计算值高于实验值，并有很大偏差，而Mikic模型的接触热导率计算值小于实验值。本文的分形模型和MT分形模型的接触热导率计算值与实验值比较接近，这表明分形模型能够较为准确地预测接触热阻。在Yovanovich模型中，认为粗糙表面微凸体的形变是塑形形变，但并没有考虑变形后体积的损失，所以在计算过程中，每对微凸体接触点的接触半径都要小于实际变形后的接触半径，接触点的对数就会相应地增加，因此得出的接触热导率计算值偏大。Mikic模型认为微凸体的变形为弹性形变，而在实际接触过程中，随着压力的增大，必然会引起塑形形变，所以接触热导率求解的结果会偏低。

表1 铝试样Al5052的热力学和机械性能

表2 铝试样Al5052的表面分形参数

图4 接触热导率理论计算值与实验值对比(试样1/试样2)

图5 接触热导率理论计算值与实验值对比(试样3/试样4)

表面形貌是接触导热的主要影响因素，传统的粗糙度参数受到仪器分辨率和取样长度的影响，使用尺度独立的分形参数模拟粗糙表面形貌建立的接触热阻计算模型更有利于对接触传热的预测。

3.2 分形参数对接触热阻的影响分析

计算模拟了2个具有相同性能的尺寸为2.5 mm×2.5 mm的材料表面的接触(材料参数见表3)，计算中设定的分形参数为当量分形参数。

表3 模拟材料的热力学和机械性能

图6为G=1.0×10-5m时，不同压力下的接触热导率与分形维数D的关系。可以看出：压力越大界面间的接触热导率越大，这反映了压力增大使界面间的接触面积增大。在定压下看，接触热导率随着分形维数D的增大而增大，这是因为当D增大时，表面轮廓的精细结构增多，精细结构增加导致整个轮廓的整体高度降低，减低了轮廓粗糙度，导致接触导热的增强，压力越大接触导热的增强越明显。

图7为在接触压力p=2 MPa时、不同分形参数G下，接触热导率与分形维数D的关系。可以看出：接触热导率随分形维数D的增大而增大，分形比例参数G越小，增大越明显。图8为接触压力p=2 MPa时、不同分形维数D下接触热导率与分形参数G的关系，可以看出：分形维数D一定时，分形参数G越大，接触热导率越小。这是因为分形参数G越大，表面粗糙度越大，增加了界面间的传热阻力。分形维数D越大，G的变化对接触热导率的影响越明显。

图6 不同接触压力下分形维数D对接触热导率的影响(G=1.0×10-5 m)

图7 不同分形比例参数G条件下分形维数D对接触热导率的影响(P=2 MPa)

图8 不同分形维数D的条件下分形比例参数G对接触热导率的影响(P=2 MPa)

综上所述，在一定的接触压力下，接触界面间的热导率随着分形维数D的增大而增大，随着分形参数G的增大而减小。这是因为在相同的分形参数G下，D越大表面粗糙度越小，在相同的分形参数D下，G越小表面粗糙度越小。表面粗糙度越小，材料表面越光滑，表面之间的接触越充分，接触热导率越大。

4 结束语

本文采用蒙特卡洛法和W-M分形函数模拟表面形貌，建立了一种能够反映固体粗糙表面形貌的随机特性和分形特性的固体接触界面间接触导热的计算模型，其计算结果与实验结果基本一致，能够较为准确地预测固体界面间的接触导热。分析了分形参数对固体界面间接触导热的影响。分析结果表明：在一定的接触压力下，分形参数G一定时，接触热导率随着分形维数D的增大而增大，G越小，增大越明显；分形参数D一定时，接触热导率随着分形参数G的增大而减小，分形维数D越大，G的变化对接触热导率的影响越明显。由于分形参数能够唯一确定地表征固体粗糙表面的形貌，本文提出的模型对接触热阻的预测具有确定性，为接触传热的工程设计提供了一种有效的方法。