多孔介質(zhì)粗糙表面分形-蒙特卡羅傳熱研究

徐靜磊，張賽，王昌進(jìn)，張杰，高偉業(yè)

（650500 云南省 昆明市 昆明理工大學(xué) 機(jī)電工程學(xué)院）

0 引言

多孔介質(zhì)是一種強(qiáng)化傳熱介質(zhì)[1-3]，被廣泛應(yīng)用于建筑、食品、化學(xué)品、煤炭、石油開采油氣貯藏、熱管虹吸材料、化工床、航空航天復(fù)合材料、人體和動物的組織器官、植物體內(nèi)的水分和營養(yǎng)運(yùn)輸系統(tǒng)等諸多領(lǐng)域，因此研究多孔介質(zhì)的傳熱在許多領(lǐng)域都具有重要的意義[4-7]。

由于多孔介質(zhì)具有自相似性，其任何表面，在微觀層面看來都具有大量的凸凹和幾何缺陷而顯示出粗糙的形態(tài)，即絕對光滑的表面是不可能存在的。對于固體表面，無論其表面是如何的形成方式，都不會是規(guī)則的，其高度在某一水平面來回波動。那些高度高出平均幾何平面的稱為微凸體，低于這個平面的稱為凹坑[8]。粗糙表面對于研究摩擦、輻射、導(dǎo)熱、導(dǎo)電、黏附及流體的流動傳熱都有相應(yīng)的影響，因此研究粗糙表面的傳輸機(jī)理，對于許多學(xué)科和領(lǐng)域都有很重要的意義[9]。

Shishkina O[10]等對粗糙形狀為矩形的邊界進(jìn)行了模擬研究，研究結(jié)果表明，冷熱板在粗糙形態(tài)下的邊界層厚度要明顯大于光滑表面，原因是粗糙單元間的間隙對流傳遞存在阻礙作用，所以傳熱不僅不會增強(qiáng)，甚至?xí)霈F(xiàn)減弱。Wei Ping[11]等對不同的粗糙壁面進(jìn)行了研究，結(jié)果表明，當(dāng)上下表面都為粗糙表面時，傳熱效果會現(xiàn)增強(qiáng)，但當(dāng)上下兩板依次為光滑表面與粗糙表面時，傳熱基本不發(fā)生增加。丁雪興[12]等針對干氣密封非穩(wěn)態(tài)下的摩擦特性對密封性能的影響進(jìn)行了研究，建立了三維粗糙實體與理想光滑剛體滑動摩擦熱力耦合模型。發(fā)現(xiàn)粗糙表面最高接觸溫度隨滑動時間增加呈逐步上升趨勢，并且溫度呈現(xiàn)一定的波動性。

近些年關(guān)于接觸熱導(dǎo)，很多學(xué)者進(jìn)行了大量的實驗研究[13-15]。理論上，人們主要研究接觸表面的微觀接觸狀態(tài)和形式，簡化成單點(diǎn)接觸模型，然后根據(jù)材料特性和熱阻定義等進(jìn)行求解。

雖然目前存在許多研究方法，但大部分學(xué)者研究的是粗糙表面與傳熱熱阻的關(guān)系，很少有人研究粗糙表面和傳熱系數(shù)之間的關(guān)系。本文基于分形-蒙特卡羅法研究傳熱系數(shù)與微凸體及孔隙率之間的關(guān)系。避開了傳統(tǒng)熱的模擬熱阻法，減少了經(jīng)驗常數(shù)的使用，運(yùn)用分形理論使得結(jié)果更接近真實情況。

1 關(guān)于蒙特卡羅法的模型建立

為了研究粗糙表面的具體形貌和內(nèi)在機(jī)制，就要刻畫出粗糙表面的分形特征，描述多尺度下的結(jié)構(gòu)。Majumdar 和Mandelbrot 提出的W-M 分形函數(shù)進(jìn)行修正得到了M-B 模型，粗糙表面的二維模型表達(dá)式如下[16]：

式中：z（x）——粗糙表面輪廓幅值函數(shù)；x——測量坐標(biāo)。此處，D——表面輪廓分形維數(shù)，（1

式（1）具有連續(xù)性、處處不可微和自仿射性，能滿足粗糙表面輪廓曲線的所有特性。

通過式（1）得到的粗糙表面垂直剖面的二維模型仿真圖如圖1 所示。通過該圖可以看出粗糙表面的二維剖面輪廓，該輪廓線能夠很好地符合粗糙表面的特性，即處處連續(xù)、不可微、自仿射性。通過放大局部輪廓，可以看出粗糙表面在任意一個長度尺寸下都是不平滑的。

圖1 二維粗糙表面剖面模擬圖Fig.1 2D rough surface profile simulation

對于粗糙表面的三維模型圖也存在相關(guān)的蒙特卡羅表達(dá)式[17]，如式（2）：

式中:M——曲面褶皺重疊數(shù)；φm，n——隨機(jī)相位（0 ≤φm，n≤2π）；n1——最低空間頻率序數(shù)，n1=0；nmax——最高頻率序數(shù)，nmax=int[lg（L/Ls）/lgγ]；Df——面積分形維數(shù)（三維空間2

圖2、圖3 是通過MATLAB 仿真得到的三維模型圖。圖2 為關(guān)于x-y 軸方向的正視圖，分別在x，y 軸方向設(shè)置了800 個點(diǎn)，每兩個點(diǎn)之間間距為1 nm；圖3 是一般視角看到的三維圖，此處在z 軸方向設(shè)置的粗糙表面上微凸體高度取值在3.2×10-8m 左右，圖3 中顏色深淺代表微凸體在x-y 面內(nèi)的面積大小，占據(jù)面積越大顏色就越重，微凸體尖端的x-y 截面的面積較小，所以顏色較淺。由圖3 可知三維蒙特卡羅法能夠顯示出粗糙表面的表面輪廓，與真實粗糙表面很接近。

圖2 三維粗糙表面關(guān)于X-Y 方向的正視圖Fig.2 Elevation view of 3D rough surface in X-Y direction

圖3 三維粗糙表面仿真Fig.3 3D rough surface simulation

2 考慮多孔介質(zhì)接觸粗糙表面分形—蒙特卡羅方法

多孔介質(zhì)表面的微凸體是隨機(jī)分布的，可以把分形理論與蒙特卡羅方法結(jié)合起來，研究多孔介質(zhì)粗糙表面的傳熱特性。

由關(guān)于微凸體分布的概率密度函數(shù)[18-19]

式中：λ——微凸體直徑。

其中f（λ）≥0，λmin≤λ≤λmax。根據(jù)概率理論，在λmin～λ區(qū)間微凸體的累積概率為

式 中：當(dāng)λ→λmin時，R=0；當(dāng)λ→λmax時，R ≈1。因為最大微凸體與最小微凸體直徑相差很大，所以式（4）中R 的取值在0～1。

式（4）可變形為

因為微凸體在多孔介質(zhì)壁面上是隨機(jī)分布的，且大小服從分形標(biāo)度，所以式中的R 是一組0～1 的隨機(jī)數(shù)，通過隨機(jī)抽樣，可以得到微凸體直徑分布的概率模型。

當(dāng)微凸體直徑從小到大排序后，可根據(jù)式（7）求得微凸體的個數(shù)。

單個微凸體的名義面積可表示為π（Lc/2）2，其中Lc為微凸體表面的名義特征長度，因此可得總的微凸體名義面積Aa

φ由式（9）決定[18]。

式中：φ——有效孔隙率；r ——孔隙直徑，這里把孔隙表面上的微凸體分布類比于分形多孔介質(zhì)的孔隙分布；dE——?dú)W幾里得維數(shù)，二維空間內(nèi)為2，三維空間取3。由于實際計算情況的需要，孔隙率φ常直接在0～1 之間進(jìn)行取值。

根據(jù)多孔介質(zhì)中的孔隙直徑和孔隙數(shù)目滿足的標(biāo)度關(guān)系[19]可知，接觸點(diǎn)直徑分布也遵循同樣的標(biāo)度關(guān)系，具體如式（10）：

可以看出，接觸點(diǎn)直徑的概率的密度函數(shù)n（λ）可由式（11）求得

可知，所有接觸點(diǎn)面積之和Ar近似為

此處Ar為所有接觸點(diǎn)的面積，即流體流過粗糙面時流體與粗糙面能夠接觸到的面積，本文也是用它來形容粗糙表面的粗糙程度的。當(dāng)流體流經(jīng)相同截面的光滑表面和粗糙表面時，與流體的接觸面積粗糙表面要比光滑表面大；反之，當(dāng)截面一定時，與流體的接觸面積越大的表面也就越粗糙。

3 多孔介質(zhì)粗糙表面?zhèn)鳠嵯禂?shù)的求解

多孔介質(zhì)粗糙表面的傳熱，可近似看作多孔介質(zhì)某一截面上的過程傳熱，熱流仍需經(jīng)過孔相和固相，只不過此時的固相與經(jīng)過截面時的固相不太一樣?？紤]截面?zhèn)鳠釙r，熱流經(jīng)過的固相就是固體骨架的截面；考慮粗糙表面時，流體除經(jīng)過孔隙外，還要考慮流體與固相接觸的面積，而非全部固相的截面積。流體與粗糙固相接觸面積由式（12）表達(dá)。

對于孔隙的面積，可由微元法得到，多孔介質(zhì)橫截面上孔隙總面積為

由傅里葉定律知傳熱系數(shù)k 如式（14）：

A 為多孔介質(zhì)傳熱截面積，Ls為表面的名義特征長度，表達(dá)式為

熱量Q 為多孔介質(zhì)粗糙面吸收的熱量。

式中：Q1——多孔介質(zhì)孔隙吸收熱空氣的熱量；Q2——多孔介質(zhì)粗糙固體面吸收的熱量。

式中:Ax——孔道的平均截面積；ca——空氣的比熱容；ρa(bǔ)——空氣密度；v ——熱空氣流動的速度；t——加熱時間；ΔT——加熱前后的溫差。

式中：cl——固體的比熱容；ρl——固體的密度。

總熱量Q 為

因為多孔介質(zhì)結(jié)構(gòu)復(fù)雜，所以熱量通過多孔介質(zhì)的不同相時，導(dǎo)熱路徑不同，得出改進(jìn)后的傅里葉導(dǎo)熱定律[17]。

聯(lián)立式（12）—式（15）、式（19）、式（20）可得傳熱系數(shù)k 如下：

4 結(jié)果分析與討論

通過對聯(lián)立后的傳熱系數(shù)計算，以及對變量的調(diào)整對比，可得分析圖如圖4 所示。

圖4 有效傳熱系數(shù)模擬值與實驗值比較Fig.4 Comparison of simulated and experimental values of effective heat transfer coefficient

圖4 是以氧化鋁為實驗對象，其傳熱系數(shù)隨孔隙率變化的實驗數(shù)據(jù)[20]與模擬值的對比圖。從圖中可以看出，有效傳熱系數(shù)隨孔隙率的增大而減小。實驗數(shù)據(jù)較模擬數(shù)據(jù)稍稍偏小，考慮到實驗存在一定的誤差，所以該偏差是在實驗誤差允許范圍內(nèi)的，該模擬值能夠很好地反映出多孔介質(zhì)粗糙表面?zhèn)鳠嵯禂?shù)隨孔隙率變化的關(guān)系，可用于計算多孔材料粗糙表面的等效導(dǎo)熱系數(shù)。本文是基于粗糙表面的分形蒙特卡羅法進(jìn)行的推導(dǎo)驗算，并且與實驗數(shù)值有很好的驗證，所以可用于多種多孔介質(zhì)粗糙表面的傳熱分析，如:粗糙墻體保溫、地下石油開采、機(jī)械零部件的傳熱等。

圖5 是面積分形維數(shù)對粗糙面的影響，這里的粗糙面的粗糙程度并不是通過粗糙度體現(xiàn)出來的，而是由流體與粗糙面的接觸面積的大小而體現(xiàn)的，當(dāng)接觸面積越大時，粗糙面就相對越粗糙，當(dāng)接觸面積越小時粗糙面就相對沒有那么粗糙。通過圖5 可以看出，所有接觸點(diǎn)的面積隨面積分形維數(shù)的增大而增加，也就是說面積分形維數(shù)越大時，相對接觸面越大粗糙表面結(jié)構(gòu)也就越復(fù)雜。通過圖5 還能看出，所有接觸面積的增長速度剛開始增長較慢，隨面積分形維數(shù)的增大而增加得越來越快。

圖5 面積分形維數(shù)對粗糙表面的影響Fig.5 Influence of surface integral shape dimension on rough surface

圖6 是面積分形維數(shù)對有效傳熱系數(shù)的影響變化關(guān)系圖。從圖6 可以看出，當(dāng)孔隙率一定時，有效傳熱系數(shù)隨面積分形維數(shù)的增大而減小，減小趨勢越來越大。因為當(dāng)面積分形維數(shù)較小時，對有效傳熱系數(shù)影響較小，孔隙率對傳熱占主導(dǎo)因素；當(dāng)面積分形維數(shù)較大時，其影響程度相對于孔隙率的影響程度較大。

圖6 面積分形維數(shù)對有效傳熱系數(shù)的影響Fig.6 Influence of surface integral shape dimension on effective heat transfer coefficient

在面積分形維數(shù)較大時，不同孔隙率下，有效導(dǎo)熱最后趨近于同一值；當(dāng)分形維數(shù)一定時，有效傳熱系數(shù)隨孔隙率的增加而減小。因為孔隙越大，其中的氣相占比越大，氣相的傳熱系數(shù)比固相的傳熱系數(shù)小，所以孔隙越大，有效傳熱系數(shù)越小。

5 結(jié)論

（1）根據(jù)分形原理，考慮粗糙表面?zhèn)鳠崽厥庑?，結(jié)合蒙特卡羅法，建立粗糙表面的模型，并推導(dǎo)出關(guān)于粗糙表面的傳熱系數(shù)的表達(dá)式。該模型既考慮了傳統(tǒng)的多孔介質(zhì)的傳熱方式，又結(jié)合了蒙特卡羅法，能夠更詳細(xì)地顧及到細(xì)微結(jié)構(gòu)對傳熱系數(shù)的影響。

（2）通過與實驗數(shù)據(jù)對比發(fā)現(xiàn)，當(dāng)面積分形維數(shù)一定時，有效傳熱系數(shù)隨孔隙率的增大而減小，在孔隙率較小時減小的速率較快，在孔隙率較大時，減小的速率較小并逐漸趨于平穩(wěn)，因為當(dāng)孔隙率較大時，空氣與粗糙表面接觸的面積相對減小，而空氣的傳熱系數(shù)較小，所以變化也就較小并趨于平穩(wěn)。

（3）通過對多孔介質(zhì)粗糙表面的蒙特卡羅-分形傳熱研究發(fā)現(xiàn)，面積分形維數(shù)越大時，粗糙表面越不平整，傳熱過程也就越復(fù)雜，所以會導(dǎo)致傳熱系數(shù)減小，同理，當(dāng)面積分形維數(shù)較小時，表面區(qū)域比較平整，傳熱路徑較短，耗時也就較少，傳熱系數(shù)較大。