【原創】OpenFOAM 邊界條件系列解析—Slip邊界（1）

05-23

來自專欄研客知識原創文章

本文原創作者：國防科技大學葉帥、周后村博士

1 slip 的含義

1.1.1 slip 邊界的物理含義

slip boundary condition 指的是滑移邊界條件，是指在固壁上（邊界上），流體的速度與固壁的速度不相同，即產生滑移速度。一般來講，按照固壁與流體之間的滑移速度，可以將邊界條件分成三類，如圖1所示，即no slip、partial slip 和slip。

圖1: slip 邊界條件分類

在圖1中可以看出，no slip 指的是流體在固壁上（邊界上）的速度為0，即無滑移速度；slip 指的是，流體在固壁上速度與流體在流場中沿固壁方向上的速度相等，即邊界切向速度保持；partial slip 指的是，邊界切向速度成一定梯度。

1.1.2 slip 邊界產生的原因

在流體力學中有兩個假設，即連續介質假設：真實流體或固體所佔有的空間可以近似地看作連續地無空隙地充滿著「質點」，以及無滑移假設：貼附於固體表面的流體運動速度與固體保持一致。在氣體中，無滑移的假設是建立在連續介質的假設之上的，分子在邊界上與固壁發生碰撞（氣體中，能量的交換主要靠碰撞），當固壁表面絕對光滑時，分子在邊界處發生反射，從宏觀統計上看，「質點」切向速度不變（滑移），而當固壁表面不光滑時，分子在邊界處發生漫反射，由於質點中分子數量眾多，碰撞之後氣體分子在統計上向各個方向散射，從宏觀上看，切向的速度為0，即無滑移。而實際上在稀薄大氣中，由於空氣中分子較少，第一個假設不一定成立，在邊界處的漫反射在統計上也不一定均勻散射，實際上，Navier 指出，流體在固壁處的速度與流體速度沿邊界面法向的梯度成正比，即

如圖2所示，流體在固壁處的速度稱為滑移速度。

圖2: slip-length

而對於液體而言，流體與固壁之間的能量傳遞主要依靠分子間的作用力，一方面，固壁處的流體之間有內聚力（粘性力），流體與固壁之間也存在剪切應力，無滑移假設認為流體與固壁之間的剪切應力可以達到無窮大，因此要遠大於流體間的內聚力，因而不會發生滑移，有學者由此建立了極限剪切應力模型，即假設固液邊界面上存在一個極限剪切應力，當固壁處的流體受到的應力在極限剪切應力之下時，不發生滑移，而當在其之上時，發生滑移。在常規的尺度條件下，液體在固體表面流動的速度滑移由於很小而常常被忽略，隨著微納電子機械系統（MEMS/NEMS）的發展，納米尺度流動中的滑移現象研究也受到了廣泛的關注。如圖3所示，隨著knudsen number（克努森數）的變化，將流動分為連續區、滑移區、過渡區和自由分子區。

圖3: 流體流動隨knudsen 數的變化

knudsennumber 定義如下：

其中，λ 表示的是分子的平均自由程，L 表示的是物理的特徵尺寸，knudsen數一般被用來判斷流體是否適合連續假設以及無滑移假設，從上圖中可以看出，當knudsen 數趨近於0 時，可以用連續性假設中的歐拉方程來描述流體；當knudsen 數介於0.001-0.1 之間時，可以用有滑移邊界條件的Navier-Stokes 方程來描述；而當處於0.1-10 時，流體處於過渡區，連續介質假設本身不再合理；當大於10 時，採用分子假設，使用玻爾茲曼方程來描述流體。在1Pa 和273K 標準大氣條件下，空氣分子的平均自由程約為0.07μm，在亞微米的特徵尺度下，假設為0.5μm，knudsen 數的量級為0.1，因此滑移邊界條件經常在MEMS 等微納米尺度或者稀薄空氣中被使用。

1.2 slip 邊界的計算

1.2.1 slip 邊界上速度計算

設流體在流場中的速度為：

流體在邊界上的速度為：

邊界的法向量為：

如圖4所示：

圖4：固壁邊界幾何

對於slip 邊界（非滲透邊界）來講，流體在邊界法向上的速度為0，邊界切向上的速度為流體在速度在切向上的分量，則slip 邊界滿足的條件為：

（3）式也即：

1.2.2 有限體積的離散

在有限體積中，守恆公式可以寫成：

如表1其中，當? 為不同的值時，公式表示不同的守恆方程。

表1: ? 不同時所代表的守恆公式

穩態下的各守恆公式為：

在有限體積法中，控制體C 上的積分有：

由高斯散度定理可知，向量v的散度在體積V 的體積分，可以化為在圍成該體積的面上的面積分，即：

因此式（9）可化為：

1.2.3 對流項的離散

式（11）的左邊即散度項的積分，等於該Cell 的控制體周圍的面上的積分之和，可以寫成：

式（12）中右邊的積分項與插值的方式有關，在面上可以使用一個點來代表此面上的? 值，也可以使用多個點來計算，這就涉及到插值的精度問題，當我們使用面中心來計算時，式12就寫成：

其中，

代表f 面上的面積，實際上，

為量? 在面f上的通量，設

表示面f上的質量通量（單位時間內通過單位面積的質量），即

則有：

至此，我們已經將對流項（散度項）化成了各個面上的通量和，在計算時，只需要對面上的值

進行插值，這部分內容將在後面舉例說明。

1.2.4 擴散項的離散

同樣由高斯定理，式（11）中右邊的擴散項可以寫成：

式中，

這一項與該面相鄰的cell 的值插值而得。

1.2.5 一維穩態無源問題離散過程

無源項的一維穩態對流擴散問題對一般場變數? 應滿足：

式中，u 為流體在x 方向的流動速度。對於一維穩態問題，應該滿足連續方程，即:

設其空間離散如圖5所示，

圖5：一維空間麗離散

則設場內共1 和2 兩個體中心點，因此未知量為?1 和?2，w,e 等虛線圍成的矩形為結點1 的控制體積，設截面積為S，流場方向為x 正向，則對於結點1 而言，由式15 可知，散度項的離散為：

其中，

這是因為面上的法向量方向指向體積外側，又因為連續定理，流過w,e 面的質量通量守恆，因此有：

因此式（20）可化為：

對w 和e 面對?e 和?w 進行插值，插值方法有很多種，這裡我們使用簡單的線性插值，可得：

因而在結點1 處，對流項的離散為：

同理，在結點2 處，有：

且

則

由式17知，結點1 處的擴散項的離散為：

其中

是流場的物理參數，一般是已知量，為方便表述，我們假設流場是各向均質的，流場中的擴散係數相同，因此上式可化為：

同樣採用一階插值來近似求解梯度，則

其中，d12 表示結點1 與2 之間的距離，因此

同樣地，在結點2 處，其擴散項的離散為：

在結點1 處，聯立(18)、(21) 以及(24) 可得：

即：

同樣地，在結點2 處有：

假設在該問題中，邊界上的值固定為0，則上式即：

在邊界處，一般給定邊界條件，對?W 及?E 進行插值，即利用與邊界相鄰的cell 的體中心的值（第二類邊界條件：梯度值）以及常數值（第一類邊界條件：固定值）進行插值，一般地

對於對流項（隱式散度運算元）來說，由式(20) 可知，需對?b 插值；而對於擴散項（隱式拉普拉斯運算元）來說，由式(23) 可知，需對r?b 插值。在OpenFOAM 中，對流場中的? 場進行求解的方程最終可以寫成：

對流項和擴散項中對邊界的插值，將會對其相鄰的cell 的變數在矩陣A中的係數（一般在對角線上）以及常數向量b產生影響。

以對流項為例，在OpenFOAM 中，fvm:div 運算元（隱式散度運算元）在gaussConvectionScheme.C 中有實現，實現的代碼如圖6所示，它需要對邊界面上的值進行插值。

圖6：隱式散度運算元

假設

則由式(20) 知，該項的面積分值為

在實現中，用inernalCoeffs() 來表示插值所產生的相鄰Cell 的變數值前的係數，用boundaryCoeffs() 來表示插值所產生的常數項，internalCoeffs()將加在矩陣A 中對應變數的係數中，boundaryCoeffs() 則將加在b 中（注意移項）。因此有：

1.2.6 slip 邊界條件在OpenFOAM 中的實現

要實現邊界條件，主要是利用邊界上的物理條件，對邊界面上的值進行插值，以計算散度運算元和拉普拉斯運算元，在OpenFOAM 中，散度運算元和拉普拉斯運算元在邊界面上主要分別利用valueInternalCoeffs()、valueBoundaryCoeffs() 以及gradientInternalCoeffs()、gradientBoundaryCoeffs()來表示，例如散度運算元對面上?b 的插值為：

在slip 的邊界條件下，對於速度場

，在面上對?b 進行計算（注意維度是3），有：

亦即

上式中分別對各方向上的?b 值進行插值，可以看出，右邊第一項與該方向上的?c 值相關，第二項與?c 值無關。以x 方向為例，有：

由於valueBoundaryCoeffs()_x 是可計算的值（需要加入常數矩陣b 中），而?cy 和?cz 都是本次帶計算的值，在此處可以用上一個時刻的值(以0 來表示上一時刻) 來進行計算，即：

如圖7所示，slipFvPatchField 中這兩個參數的繼承了在transformFvPatch-Field 類中的函數，在transformFvPatchField 中，valueInternalCoeffs() 和valueBoundaryCoeffs() 兩個函數的聲明如圖8所示。

圖7：slipFvPatchField 類的繼承關係

圖8：transformFvPatchField 類中兩個函數的聲明

其中，nGradTransformDiag()函數在slip 邊界條件時，由basicSymmetryFvPatchField 實現，其代碼如圖9所示，Type 為vector，（圖中對snGradTransformDiag 的實現中，最後一個transformFieldMask <Type>pow<vector,pTraits<Type>::rank>diag)) 尚未深入了解，待補充）。