宋孟天，雷杰超，張建成

(1.廣西大學(xué) 土木建筑工程學(xué)院，廣西 南寧 530004；2.廣西大學(xué) 機械工程學(xué)院，廣西 南寧 530004)

0 引言

微流控是在微通道中操縱微量流體的運動以達到工程目的的技術(shù)，整個系統(tǒng)可集成在厘米級別甚至更小的芯片上，應(yīng)用領(lǐng)域有熱傳、微軸承潤滑、高通量藥物篩選等[1-3]。具有代表性的一種結(jié)構(gòu)是由外圓柱罩和黏性流體包裹的轉(zhuǎn)軸構(gòu)成的同心圓柱，微尺度下壁面效應(yīng)大大增強，工作過程中，由于器件之間的相對運動，并且流體存在黏性與腐蝕性，軸的表面會因為沉積、磨損和腐蝕原因變得粗糙[4-6]。表面粗糙度一方面會對微管內(nèi)流體的流動特性產(chǎn)生影響，改變其工作參數(shù)使其不再滿足原有的要求，另一方面對軸的機械性能以及作用于其上的力和轉(zhuǎn)矩有極大的影響，給設(shè)備的工作狀態(tài)帶來隱患。因此，發(fā)展表面粗糙度檢測技術(shù)極為重要，可為提前判斷設(shè)備是否需要更換或維修提供依據(jù)，避免設(shè)備失效破壞造成經(jīng)濟損失。傳統(tǒng)的粗糙度檢測方法有掃面隧道顯微鏡法、超聲波法和光學(xué)測量方法等[7-9],所需設(shè)備昂貴。同時，微流控設(shè)備中的元件往往被其他的裝配單元所包裹，檢測時需要進行拆卸，這就使得傳統(tǒng)的粗糙度檢測技術(shù)難以直接應(yīng)用。因此，需要發(fā)展基于流體的粗糙度檢測方法。

微流控領(lǐng)域的一個基本問題是粗糙微管內(nèi)流體的庫埃特流動和泊肅葉流動。PONJAVIC等[10]采用熒光漂白恢復(fù)法研究了彈流潤滑中流體的速度剖面和滑移長度，首次實現(xiàn)全厚度流速的直接測量。郝鵬飛等[11]發(fā)現(xiàn)雷諾數(shù)較小時，微觀粒子圖像測速技術(shù)本身帶有1 % ～ 2 %的誤差。姚華平等[12]采用有限差分法求解雷諾方程，研究了兩個方向的正弦/均方根隨機粗糙表面之間的微流道庫埃特流動，其結(jié)果表明油膜承載能力隨波長以二次曲線的規(guī)律變化。XIONG等[13]利用有限體積法求解三維隨機粗糙圓管中水流動的Navier-Stokes方程，結(jié)果表明主流區(qū)域受到粗糙度的強烈影響，泊肅葉數(shù)偏離高達11.9 %，劉趙淼等[14]也得到類似的結(jié)果。

在解析解的研究上，WANG等[2]利用邊界攝動法研究具有周向正弦表面粗糙度的兩同心圓柱間的庫埃特流，最早得到內(nèi)旋轉(zhuǎn)軸上的力和轉(zhuǎn)矩，此外，還利用邊界攝動法研究了具有兩個方向表面粗糙度的圓管泊肅葉流，表面粗糙度利用正弦函數(shù)的乘積來模擬[15]。SONG等[16]用余弦波、矩形波和三角波模擬圓管軸向壁面粗糙度，采用邊界攝動法推導(dǎo)泊肅葉流解析解，但是所研究的是單一方向的粗糙度。

以上文獻報道的表面粗糙度模型多為一個方向，不能真實全面體現(xiàn)實際粗糙度的影響。兩個方向上的粗糙度雖有研究，但研究的是圓管內(nèi)的流動，對于幾何結(jié)構(gòu)較為復(fù)雜的兩圓柱之間的微管流，周向和軸向粗糙度同時存在的情況尚未考慮。另外，以上成果集中在粗糙度對流體流動特性的影響上，很少根據(jù)流體參數(shù)反推粗糙度。本文通過精確到二階的邊界微擾方法，研究微管內(nèi)轉(zhuǎn)軸表面兩個方向的粗糙度幅值和疏密程度對軸上的轉(zhuǎn)矩和流體流率的影響，以此提出一種新的表面粗糙度檢測模型。利用計算流體力學(xué)軟件COMSOL進行數(shù)值模擬，仿真結(jié)果和理論結(jié)果呈現(xiàn)良好的一致性。在此之前，本課題組已開展了表面粗糙度的相關(guān)微擾解析工作[17]，本文則進一步以數(shù)值計算承接解析的思路，驗證解析解，并拓展粗糙度檢測模型的應(yīng)用范圍。

1 物理模型與基本方程

圖1表示兩個同心圓柱，內(nèi)軸粗糙，外圓柱罩光滑。微流控設(shè)備在工作時，實際是內(nèi)軸轉(zhuǎn)動，例如：內(nèi)軸以角速度Ω繞z軸順時針旋轉(zhuǎn)，流體則因壁面摩擦力被驅(qū)動形成庫埃特流動，如圖1(b)所示。當內(nèi)軸上存在不同形式和大小的粗糙度時，其上的轉(zhuǎn)矩因此而變化，此問題是非穩(wěn)態(tài)的。本文的目的是為了檢測轉(zhuǎn)軸表面粗糙度，而軸上受到的轉(zhuǎn)矩與外圓柱罩上受到的轉(zhuǎn)矩相等，因此為了便于計算，將內(nèi)粗糙軸固定，使外光滑圓柱沿相反方向以同樣的角速度旋轉(zhuǎn)，則問題變?yōu)榉€(wěn)態(tài)，且外圓柱可以較低角速度Ω′旋轉(zhuǎn)，以便于進行層流分析，即為本文所求解的等價問題，如圖1(c)所示。此外，為了達到檢測軸表面粗糙度的目的，本文另一研究內(nèi)容是讓設(shè)備靜止，即內(nèi)軸和外圓柱都沒有旋轉(zhuǎn)，在微管間隙的兩端施加壓力差形成泊肅葉流動，研究內(nèi)軸表面粗糙度對流率的影響，從而得到更全面的信息以建立轉(zhuǎn)軸表面粗糙度模型。

外圓柱和內(nèi)粗糙軸半徑分別為R與r*[17]，

(1)

其中：Rb(b<1)為內(nèi)軸的平均半徑，ε為軸的粗糙度幅值與其平均半徑之比，λz為z*方向即軸向的波長。

以外圓柱半徑R為特征長度進行無量綱化，取外圓柱無量綱半徑為R=1，則內(nèi)軸無量綱半徑為

r=b+εcos(nθ)cos(αz)=b+εg(θ,z)

(2)

其中：n=2πR/λθ，α=2πR/λz，分別為周向θ和軸向z方向上的波數(shù)，λθ為θ方向上的波長。

(a) 物理模型立體圖

(b) 原始問題：Ω為工作時轉(zhuǎn)軸角速度

(c) 等價問題：Ω′為檢測時外圓柱罩角速度

對穩(wěn)態(tài)不可壓縮黏性流，連續(xù)方程為[18]

*·u*=0，

(3)

其中：u*=(u*,v*,w*)為流場速度向量，u*，v*和w*分別為徑向、周向和軸向的速度分量。動量方程為[18]

ρ(u*·*)u*=-*p*+μ*2u*，

(4)

其中：μ為流體動力黏度，p*為壓力，ρ為流體密度。

2 微擾分析

2.1 庫埃特流動

以U=Ω′R為特征速度，μU/R為特征壓力，則無量綱化的速度、壓力分別為：u=u*/U=u*/(Ω′R)，p=p*/(μU/R)。因此，連續(xù)方程式(3)和動量方程式(4)分別變?yōu)?/p>

·u=0，

(5)

(6)

小轉(zhuǎn)速下，旋轉(zhuǎn)雷諾數(shù)ReΩ≡ρΩ′R2/μ遠小于1，因此式(6)左側(cè)項可忽略，得到斯托克斯方程

2u=p。

(7)

將速度與壓力項進行精確到二階的微擾展開[17]，得

(8)

由于內(nèi)軸固定，外圓柱旋轉(zhuǎn)，r=b+εg處的無滑移速度邊界條件為

(9)

r=1處的邊界條件為

(10)

其中：速度u的下標0、1、2分別表示零階、一階和二階速度。

零階邊界條件(ε=0)：

u0|r=b=u0|r=1=0,v0|r=b=0,v0|r=1=1,w0|r=b=w0|r=1=0。

(11)

一階邊界條件：

(12)

二階邊界條件：

(13)

零階解：

(14)

壓力為常數(shù)，因此壓力梯度為零。為了簡便，我們?nèi)0=0。

一階解：

(15)

其中

超聲波測距模塊用來測量模塊距離地面的距離d;MPU-9250模塊用來測量拐杖運動的角速度w。通過大量的實驗來模擬老人摔倒時的狀況,發(fā)現(xiàn)當d>240 cm、w>5 rad/s且最終測量角度大于80°時,有99.2%的情況老人處于跌倒狀態(tài),將此作為判定老人跌倒的標志[9],同時控制GSM模塊將報警信息發(fā)送至遠程的手機監(jiān)測軟件,其軟件流程如10所示。

(16)

式(16)中：In和Kn為n階修正貝塞爾函數(shù)。

對于轉(zhuǎn)矩的計算，周向速度是主要貢獻，因此只給出V1c(r)的表達式，其余三項U1c(r)，W1c(r)和P1c(r)也是修正貝塞爾函數(shù)的組合，A1、B1、C1、D1、E1和F1為待定系數(shù)，根據(jù)邊界條件確定。

二階解：

(17)

其中：U2α，W2α和P2α均為零，V20(r)=A2I1(2αr)+B2K1(2αr)。

平均轉(zhuǎn)矩：

由于對稱性，有

(18)

對θ積分后平均速度為

(19)

(20)

外圓柱上單位面積的平均轉(zhuǎn)矩為

(21)

為了簡便，可以取αz1=π/2，αz2=3π/2，同時內(nèi)粗糙軸上的平均轉(zhuǎn)矩與外圓柱上的轉(zhuǎn)矩大小相等，方向相反，因此有內(nèi)軸上的轉(zhuǎn)矩為

(22)

M*=μΩ′R3M=μΩ′R3[M0+ε2η]。

(23)

在θ和z方向上均覆蓋一個波長的區(qū)域面積，可由下式給出

Sc=(2π)2/(nα)=(2π)2/Ac，

(24)

因此，Ac=nα可稱為粗糙度密度，通過定義粗糙度密度，可綜合分析波數(shù)n和α帶來的影響。以b=0.5，ε=0.05為例，平均轉(zhuǎn)矩M隨n和Ac的變化如圖2所示?？梢姰攏取1和2時，M隨粗糙度密度Ac的增加而增加，當M取到最小值，即Ac=20時，n和α的組合為(3, 20/3)；而當n大于等于3時，平均轉(zhuǎn)矩M基本不隨Ac的變化而變化。

圖2 b=0.5，ε=0.05時平均轉(zhuǎn)矩M隨n和Ac的變化情況Fig.2 Variation of the mean torque M versus n and Ac (b=0.5, ε=0.05)

2.2 泊肅葉流動

內(nèi)粗糙軸和外光滑圓柱靜止，微管間隙兩端施加壓力差G，問題為穩(wěn)態(tài)且流動為單方向，因此式(4)左側(cè)項為零。以G為特征壓力，GR2/(4μ)為特征速度對連續(xù)方程式(3)和動量方程式(4)進行無量綱化，得

(25)

無滑移邊界條件為

(26)

零階邊界條件(ε=0)：

u0|r=b=u0|r=1=0,v0|r=b=v0|r=1=0,w0|r=b=w0|r=1=0。

(27)

一階邊界條件：

(28)

二階邊界條件：

(29)

零階解：

(30)

一階解：

(31)

其中

(32)

對于總流率的計算，軸向速度是主要的貢獻，因此只給出W1c(r)的表達式，其余三項U1c(r)，V1c(r)和P1c(r)也是修正貝塞爾函數(shù)的組合。

二階解：

(33)

同樣地，只給出W20(r)的表達式：

(34)

總流率：

總流率可寫為

(35)

微管中任意界面流率相同，為了簡便，令z=π/(4α)，則W2α項積分為零，式(35)變?yōu)?/p>

(36)

其中

(37)

Q0為粗糙度為零時的流率，χ為粗糙度對總流率產(chǎn)生的凈效應(yīng)?？梢姡⒐苤械目偭髀时沪哦A修正。

有量綱的流率為

(38)

與庫埃特流動一樣，以b=0.5，ε=0.05為例，總流率Q隨n和Ac的變化如圖3所示?？梢妼τ谕瑯拥膎值，當粗糙度密度Ac增加，即波數(shù)α增加，總流率Q隨之減小；n變大時，總流率Q隨Ac變化的減小趨勢放緩。當Q取到最小值，即Ac=3.5時，n和α的組合為(2, 1.75)。

圖3 b=0.5，ε=0.05時總流率Q隨n和Ac的變化情況Fig.3 Variation of the total flow rate Q versus n and Ac (b=0.5, ε=0.05)

3 數(shù)值驗證與拓展

為了驗證理論結(jié)果的有效性，本文使用CFD軟件COMSOL進行數(shù)值模擬，采用方法為有限元法。外圓柱和內(nèi)粗糙軸為壁面，均為無滑移邊界條件。外圓柱半徑為1×10-4m，計算域的軸向長度為一個粗糙度波長，即2πR/α，流動介質(zhì)為水。邊界設(shè)定和網(wǎng)格如圖4所示。

(a) 庫埃特流動設(shè)定

(b) 泊肅葉流動設(shè)定

以b=0.5，ε=0的情況為例，進行網(wǎng)格無關(guān)性驗證，結(jié)果見表1和表2?？梢姰斁W(wǎng)格加密到一定程度時，所得結(jié)果很接近。因此，對于庫埃特流動，選擇normal網(wǎng)格類型；對于泊肅葉流動，選擇finer網(wǎng)格類型。庫埃特流的轉(zhuǎn)矩計算只需在外圓柱表面積分，而泊肅葉流的流率計算需要在微管截面進行積分，因此泊肅葉流內(nèi)部網(wǎng)格需要劃分得更加密集，才能得到更準確的結(jié)果。

數(shù)值模擬的參數(shù)設(shè)定為：粗糙軸半徑b=0.2，0.4，0.5，0.6，0.7，0.8，粗糙度幅值ε=0，0.01，0.025，0.050，0.1，0.150，0.200。對于庫埃特流動，對應(yīng)的(n,α)值的組合分別為(2, 10)，(2, 10)，(3, 20/3)，(3, 20/3)，(3, 20/3)和(2, 10)。對于泊肅葉流動，(n,α)值的組合分別為(4, 2)，(4, 2)，(2, 1.75)，(3, 7/6)，(4, 2)和(8, 1.5)。微擾分析和數(shù)值模擬的結(jié)果對比如圖5所示，可見，庫埃特流動中粗糙軸受到的轉(zhuǎn)矩M隨著軸半徑b和粗糙度幅值ε的增加而增大，而泊肅葉流動中的總流率Q則相反。這是因為當b和ε逐漸增大時，粗糙軸與外圓柱之間間隙變小，壁面效應(yīng)變得更加顯著。而對于泊肅葉流動，b和ε越小，間隙越大，因此流體流率越大。當ε=0，0.010，0.025、0.050和0.100時，兩者之間誤差較小。當ε> 0.1時，誤差逐漸增大，說明此時邊界微擾法已不再適用。這是因為ε過大，一階和二階解的絕對值更加接近零階解，使得微擾解失去了有效性。

因此，邊界微擾法的適用范圍為粗糙度幅值ε在0到0.100之間。ε=0.150，0.200和0.250時，庫埃特流動和泊肅葉流動的數(shù)值延伸結(jié)果如表3和表4所示，以拓展粗糙度檢測模型的應(yīng)用范圍。b=0.8時，由于尺寸限制，只計算到ε=0.200。

表1 庫埃特流動網(wǎng)格無關(guān)性驗證Tab.1 Grid independence verification of Couette flow

表2 泊肅葉流動網(wǎng)格無關(guān)性驗證Tab.2 Grid independence verification of Poiseuille flow

(a) 庫埃特流動轉(zhuǎn)矩M對比圖

(b) 泊肅葉總流率Q對比圖

表3 不同b和ε下庫埃特流動無量綱轉(zhuǎn)矩數(shù)值結(jié)果Tab.3 Non-dimensional numerical simulations torque results of Couette flow under different b and ε

表4 不同b和ε下泊肅葉流動無量綱流率數(shù)值結(jié)果Tab.4 Non-dimensional numerical simulations flow rate results of Poiseuille flow under different b and ε

4 粗糙度檢測模型

首先必須指出，轉(zhuǎn)軸和外圓柱之間的間隙固定，即b是已知的，因此M0和Q0易于計算。由此，基于平均轉(zhuǎn)矩M和總流率Q的導(dǎo)出表達式(22)和(36)，在實踐中，結(jié)合測量結(jié)果，建立了表面粗糙度預(yù)測模型

(39)

對于給定的b值，粗糙軸不同的ε，n和α值得到不同的平均轉(zhuǎn)矩M和總流率Q。以b=0.5為例，不同粗糙度下(ε=0.01, 0.05, 0.08, 0.1)n和α變化時的M和Q值如圖6所示。由于平均轉(zhuǎn)矩M和總流率Q由實際測量得到，因此可以反推出ε和n的值，并得到α的值。例如，圖6中的星號處M=1.111 0，Q=0.118 2，此時ε=0.1，n=12，α=27。此時，如果外圓柱半徑為R=1×10-6m，則表面粗糙度幅值為ε×bR=0.1×0.5×10-6m=5×10-8m。不同的n和α的值組合下，M或Q值可能為相同的值，但這種情況對于M和Q不會同時發(fā)生，因此n和α的值可以唯一確定。又如，對于數(shù)值拓展部分，若M=0.754 7，Q=0.177 9，此時ε=0.2，n=3，α=20/3，如果外圓柱半徑為R=1×10-6m，則表面粗糙度幅值為ε×bR=0.2×0.5×10-6m=1×10-7m。

此外，可單獨使用庫埃特流進行表面粗糙度檢測。對于同一種流體，有量綱轉(zhuǎn)矩隨著轉(zhuǎn)速Ω′線性增大。由圖6(a)及式(39)可知，若測試兩種不同轉(zhuǎn)速下的軸上的轉(zhuǎn)矩，必定對應(yīng)于相同的n和α，然后根據(jù)轉(zhuǎn)矩的大小確定相對粗糙度ε。類似地，對于泊肅葉流動，對于同一種流體，有量綱流率隨著壓力差G線性增大。由圖6(b)及式(39)可知，若測試兩種不同壓力差下的流率，必定對應(yīng)于相同的n和α，然后根據(jù)流率的大小確定相對粗糙度ε。

(a) 平均轉(zhuǎn)矩M隨n和α的變化圖

(b) 總流率Q隨n和α的變化圖

5 結(jié)語

本文針對粗糙度沿周向和軸向分布的轉(zhuǎn)軸和外光滑圓柱罩之間微管的液體微流動，以粗糙度幅值為擾動項，求解了Stobes方程的解析解，得到了不同半徑下，庫埃特流動中粗糙軸受到的轉(zhuǎn)矩M和泊肅葉流動中的總流率Q隨粗糙度幅值和波數(shù)的變化規(guī)律。結(jié)果表明，轉(zhuǎn)矩M隨著軸半徑和粗糙度幅值的增加而增大，而總流率Q則相反。通過COMSOL軟件進行數(shù)值模擬，驗證了解析解的準確性，并拓展了應(yīng)用范圍。基于以上結(jié)果提出了利用流場參數(shù)反推表面粗糙度的方法。

本文的創(chuàng)新點在于，引入周向和軸向兩個方向的粗糙度模型，將兩圓柱間微管流邊界微擾理論推廣到二階，提出了基于流體方法的粗糙度檢測模型與檢測方法，所檢出的粗糙度可視為整體表面的平均，同時采用數(shù)值方法拓展了模型的應(yīng)用范圍。本文的結(jié)論可用于微流控設(shè)備表面粗糙度的不拆卸快速檢測。