基于2D 網格的軸對稱浸沒邊界法1)

蔡政剛 潘君華 倪明玖

(中國科學院大學工程科學學院,北京 100049)

引言

流體領域中的顆粒兩相流[1-2]一直是工程和學術研究上的重要內容.除了理論和實驗研究之外,數值模擬也是一種重要的研究手段.而在顆粒兩相流的相關研究中,復雜幾何和動邊界問題是十分普遍的.當我們采用數值模擬手段研究上述問題時,使用貼體網格是有一定挑戰的.例如固體幾何比較復雜,那么壁面邊界層附近的網格處理將會是麻煩的.同樣當固體邊界發生運動時,貼體網格會發生畸變,使得網格的正交性變差.但是浸沒邊界法(immersed boundary method,IBM)為相關問題提供了非常好的技術手段.

浸沒邊界法首先是由Peskin[3-4]為了模擬心臟力學以及相關的血液流動而發展出來一種方法.該方法對流體和固體采用兩套不同網格,流體一般采用正交的歐拉網格而固體一般采用拉格朗日網格.固體邊界上的信息通過直接在Navier-Stokes 方程中添加力源項反映到流場中,并引入狄拉克δ 函數來聯系固體域和流體域.隨后此方法也被后續的科研工作者推廣到了其他眾多領域,如文獻[5]將IBM方法推廣到了湍流領域、自然和強迫對流傳熱領域;文獻[6]將IBM 方法分別應用到顆粒液滴運動問題以及可壓縮黏性流動問題;Grigoriadis 等[7]則將其推廣到磁流體動力學(magnetohydrodynamics,MHD)領域.另外,為了對鈉冷快堆中的停堆組件落棒過程進行數值模擬研究,秦如冰等[8]開發了2D 的浸沒邊界法.同時,近年來由于格子玻爾茲曼方法興起,IBM 方法與格子法的結合使用也是熱門的領域.周睿等利用格子玻爾茲曼浸沒邊界法分別對雙層剛性植被明渠水流[9]以及變動水面水庫[10]進行了數值模擬和驗證.佟瑩等[11]發展了一種直接力格式的浸沒邊界格子模型用于動邊界流動數值計算.李橋忠等[12]提出了一種基于浸沒邊界-簡化熱格子玻爾茲曼方法的耦合模型,該模型簡化了邊界條件的實現形式,減小了計算成本.

通常根據浸沒固體邊界條件引入流體計算域的差異,后續發展的IBM 方法主要分為兩類:一類是連續力法如Peskin[3-4],另一類是離散力法如文獻[13].相比于連續力法,離散力法并不通過對控制方程改動來施加浸沒邊界條件,而是在浸沒邊界附近單元直接進行插值引入邊界條件.考慮到后者在滿足物理量守恒特性上要優于前者,本文所采用的IBM 方法也是屬于離散力法.但是使用離散力方式處理浸沒邊界時,固體界面將背景的流體網格切割成不規則形狀.如何處理切割后的網格單元的守恒性是一個難點.文獻[14]通過耦合IBM 以及切割網格技術有效地處理了由于動邊界引發的壓力震蕩,但其方法實現是比較困難的.隨后Pan 等[15]發展了一種相容守恒格式的浸沒邊界法用以研究MHD 流動以及動邊界問題.該方法既能夠抑制動邊界產生的壓力震蕩又相對前者更易實現.

同時經過本文調研發現,目前大部分關于IBM方法的研究[1-15]主要針對的是較為普遍的二維或三維問題,而對于軸對稱這種特定物理問題(比如雷諾數Re＜210 時的小球繞流[16],伽利略數G＜156 時的小球自由運動[17])的相關研究是比較少的.Lai 等[18]基于柱坐標發展了一種軸對稱IBM 用以研究伴隨不溶性表面活性物質的軸對稱界面流.Li 等[19]在研究細胞分裂增長的軸對稱過程中發展了一種軸對稱浸沒邊界模型.Li[20]采用渦量-流函數耦合控制方程發展了一種快速軸對稱浸沒邊界投影法用以處理小球顆粒與壁面正碰問題.上述三種模型均是基于連續力法來實現的,而目前基于離散力格式的軸對稱IBM 方法還未發展.

因此本文發展了一種基于有限體積法(FVM)以及離散力格式的2D 軸對稱浸沒邊界法.此算法是以Pan 等[15]所發展的2D 和3D 相容守恒浸沒邊界法擴展而來.因此對于浸沒邊界條件的引入以及銳利階梯界面的劃分,本文沿襲其處理方式.但是Pan 等[15]所發展的原始IBM 方法并未考慮碰撞情況,當需考慮顆粒壁面碰撞,在小球離壁面只有一層網格時,由于插值數據信息不足導致數值計算崩潰.所以在其原始算法基礎之上,本文進一步完善了浸沒邊界離壁面邊界只有一層網格時的情況.

如圖1 所示為例,所有背景網格單元被分為三種:第一種是棕色表示的固體網格單元,其網格中心在浸沒邊界法內部;第二種是綠色表示的IB 網格單元,其為離浸沒邊界最近的第一層網格單元.第三種是剩下的藍色表示的流體單元(直接參與離散矩陣計算).除了上述網格單元之外,黃色的壁面面元代表了流體域的壁面邊界,面心的數據由邊界條件給出.圖中的紅色實線則是表示浸沒邊界,階梯狀虛線表示銳利界面(實際計算的邊界).

圖1 浸沒邊界和網格單元劃分示意圖Fig.1 Schematic diagram of immersed boundary and dividing cells type

由于IB 單元一般是不規則的切割單元(如圖1中所示,綠色IB 單元被紅色浸沒邊界切割),本文采用文獻[21]提出的加權最小二乘法(LSM)插值得到.而插值所用到的信息一般為浸沒固體邊界條件以及周圍流體網格中心的物理量.當浸沒邊界開始接觸壁面邊界時(相距只有幾層網格),流體側壁面邊界上的面心數據也會作為插值數據點.得到IB 網格單元中心數據之后,再通過引入階梯狀銳利界面作為近似邊界來封閉離散方程,從而實現固體邊界條件的引入.

一般根據物體運動速度的確定方式,本文將其分為兩種基本類型.第一種是運動物體的強迫運動,人為給出物體的速度隨時間的變化.二是通過耦合運動學方程和Navier-Stokes 方程,由流體-顆粒的相互作用來確定物體的速度.通常,運動學方程中的力系數在時間上是顯示離散的,以便將運動學方程從Navier-Stokes 方程中解耦[22].本文通過2D 軸對稱IBM 方法模擬Stokes 流小球近壁運動以及小球自由下落碰壁彈跳分別對上述兩種運動方式進行驗證.

1 數值模型

1.1 控制方程

本文算法考慮的是一種不可壓縮的、運動黏性為 ν、密度為 ρf的流體所產生的無環流軸對稱流動.在軸對稱假設之下,本文可以利用柱坐標系將三維流動的控制方程簡化成兩維形式.本文使用u=(ur,uz)以及p來分別表示子午面上的速度和壓力,其中ur和uz分別表示徑向和軸向的速度分量.因此無量綱的流體動力學控制方程可以寫成

這里無量綱參數雷諾數Re=UL/ν,其中U和L分別表示軸對稱流動的特征速度和特征長度.對于小球算例特征長度L=D,D為小球直徑.而對于圓盤特征長度L=D,D為圓盤徑向直徑.另外需要注意的是軸對稱柱坐標系下的拉普拉斯算符 Δ 以及梯度算符 ? 與直角坐標是不一樣的,分別為

另外,對比式(2)和式(3),我們能夠清楚地發現在徑向動量方程中由于坐標系的變化導致多出來的源項ur/r2.由于這一項在靠近對稱軸位置變得很大,因此對于該項的處理是十分重要的.

因為本文的數值模擬計算也涉及到小球顆粒自由下落的運動過程.所以除了上述流體的控制方程之外,顆粒運動的控制方程也需要計算.一般來說我們可以將顆粒的運動速度表示為UΓ=Up+ωp×rΓ,其中Up表示顆粒質心的平移速度、ωp表示顆粒的角速度、rΓ表示顆粒表面到質心的距離矢量.顆粒運動方程則分別是由牛頓動量方程(6)以及角動量方程(7)決定,即

這里 τ=-Ip/ρf+ν(?u+?uT)表示應力張量,其中I和p分別為單位張量和流體動壓(不含靜水壓力).而下標p和f分別對應顆粒和流體.Vp,mp和np則分別表示固體顆粒的體積、顆粒質量以及顆粒表面的單位外法向矢量.fp為顆粒所受外力之和,包含流體作用力、重力以及浮力.Ip和rp分別表示固體顆粒的慣性張量以及顆粒所受到的外力作用點到顆粒重心的距離,而Tp則表示顆粒所受到的合力矩.因為本文只考慮無環流軸對稱運動,顆粒在自由下落運動過程中是沒有轉動發生的即固體顆粒的角速度 ωp=0 .這也意味著我們只需要單獨求解顆粒的牛頓運動方程 (6)即可.

1.2 數值離散

本文的算法是基于柱坐標系以及有限體積法來開展的.因此對于無環流軸對稱流動,子午面上生成兩維的平面網格后,將網格沿軸線(z軸)分別向正反環向方向旋轉 dθ/2 角度后便形成了一層楔形三維網格(夾角為 dθ).如圖2 所示:其中坐標r軸位于網格中界面,S f表示網格面面積,nf表示網格面單位法向量.

圖2 柱坐標系下的體積微元Fig.2 Volume element in cylindrical coordinates

網格體積微元表示為 Ωc=drdzrcdθ,z方向網格表面積計算表示為Sz=rfdθdr,r方向網格表面積計算表示為Sr=rfdθdz,θ 方向網格表面積計算表示為Sr=dθdz.這里rc和rf分別為網格微元中心和網格面心徑向坐標.雖然在進行空間離散時,由于有限體積法的使用引入了一層三維網格單元,但在計算過程中本質上仍然是2D 計算.然后對流體控制方程進行的完整離散,本文采用一個經典的具有二階時間和空間精度的分步式投影法[23].主要計算過程可以分為以下四個步驟.

第一步分別計算預測步徑向速度ur和軸向速度uz.對于對流項以及黏性項中的拉普拉斯項均采用二階中心差分格式離散,而對于徑向方程中多出來的ur/r2項采用隱式處理.則上述徑向和軸向動量方程的離散形式為

這里需要注意的是式(8)和式(9)中的上標k和 * 分別表示當前時層(已知)和預測步時層(未知),下標c和f則分別表示網格單元中心和網格單元面.表示網格單元面上的質量通量.?pr和?pz則分別表示壓力梯度的徑向和軸向分量.Ωc表示網格單元的體積(在環向 θ 方向單位長度取為1,只有一層網格).

第三步計算k+1時層壓力場.首先根據上一步計算得到的網格中心的中間步速度場,我們通過插值可以得到網格單元面上中間步質量通量

這里為了防止速度壓力求解失耦,發生棋盤壓力震蕩現象,我們需要將網格單元中心的速度插值到網格單元面心上().其次通過上述的網格單元面上的中間步質量通量來求解壓力泊松方程得到新時層的壓力場.這里壓力泊松方程為

第四步根據k+1 時層的壓力場,修正中間步速度場,得到k+1 時層的網格單元中心和面上的速度場以及面上的質量通量分別為

至此完整的一次投影步迭代計算結束.需要說明的是對于網格中心處的梯度計算是采用高斯定理方式 ?pc=,而對于網格單元界面上的法向壓力梯度是采用二階中心差分形式,下標A和nb分別表示網格A與網格A的相鄰網格nb,dnb-A表示網格A與相鄰網格nb中心距離.

另外對于在不可壓縮黏性流體中顆粒小球垂直靠近壁面運動時,由于運動顆粒和靜止壁面之間的間隙很小,數值模擬所需的網格分辨率非常高.因此為了保證計算的穩定性,需要在一個時間步長內重復多次執行投影步計算,如圖3 所示.

圖3 流體方程計算流程圖Fig.3 Flow chart of fluid equation calculation

根據本文測試的經驗,這里的迭代次數一般設置為3～5 次即可.另外在前兩次迭代計算時,一般會對壓力和速度做松弛處理來保證計算的穩定性.

除了求解流體的Navier-Stokes 方程之外,本文也要求解浸沒固體顆粒的牛頓運動方程.這里我們直接采用Euler 的顯式離散格式處理,如下式所示

這里的應力張量 τk+1=,其中浸沒邊界上的壓力以及速度梯度分別為,而dIB-ib表示IB 網格單元到浸沒邊界上的映射ib點的距離.需要注意的是可由當前浸沒顆粒速度即給出,而需要通過由求解流體離散方程(8)～(14)得到的流體單元數據經過最小二乘法插值得到.關于IB 單元和映射ib點的位置可以參考圖1.在下一小節中我們會說明如何插值得到IB 單元數據.因此通過使用式(15),我們很容易求解得到k+1 時層的顆粒速度.進而更新k+1 時層的浸沒邊界條件,繼續求解下一時層流場.

1.3 IB 邊界條件引入

因為本文采用離散力法,所以固體浸沒邊界的引入主要分為三個步驟.第一步通過最小二乘法插值得到IB 網格單元的值,第二步利用階梯狀銳利界面替代真實的固體浸沒邊界封閉控制方程求解流場信息,第三步利用新求解的流場信息修正IB 網格單元信息.接下來我們將分別說明這三個步驟的實現.

第一步:我們對IB 單元中心物理量 φ 在浸沒邊界投影點(圖1ib點)在子午面上做二維平面泰勒展開,即

其中投影點坐標用 (xib,yib),Δx=x-xib和Δy=y-yib則分別表示IB 單元中心到投影點的各個方向上的距離.這里需要說明一下,盡管本文流場的方程離散計算是在軸對稱柱坐標系上進行的,但是IB 網格單元中心的插值還是在二維直角坐標基礎上完成的.那么只要確定了式(16)中的泰勒展開系數,等,就能得到IB 單元中心的值.針對這些未知系數,本文采用IB 單元周圍的已知流體網格單元以及壁面邊界面元信息構建一個代數方程組(未知系數作為待求解變量).通過最小二乘法求解方程系數,以此確定式(16)中的泰勒展開系數.當插值搜索的數據點多于插值方程未知系數的個數時,我們通過帶權重的最小二乘法[21]求解:γ=.這里M和m分別表示數據點的總數和第m個數據點,而表示距離權函數.其中表示數據點與投影點之間的間距.R表示這些間距中的最大值.

針對紐曼邊界條件也是類似處理,這里不再贅述,詳情可以參考文獻[15].考慮到插值的精度和效率,本文的插值精度采用三階形式.

第二步:階梯近似界面引入,如圖1 所示.首先忽略浸沒邊界的存在,則離散方程的最終離散形式在k+1 時層的一般表達為

其中下標nb表示網格單元p周圍的網格單元,nfl/(nIB)和nIB則分別表示相鄰網格單元中不包括IB 單元的總數以及IB 網格單元的數量.這里方程(17)因為缺少邊界條件是無法求解的.但是我們可以利用第一步得到的IB 單元中心的插值信息(k時層)以及階梯近似邊界來封閉流體離散方程.此時新的離散方程最終離散形式為

這里IB 單元被當做源項處理放在公式右側.

第三步:修正IB 單元的信息.通過第二步計算,我們得到了k+1 時層的流場信息.但是IB 單元還是k時層的,我們需要對其再做一次最小二乘法插值,將其更新到k+1 時層.至此IB 邊界條件引入完成.

1.4 2D 自適應網格技術實現

對于動邊界問題,由于顆粒在靠近壁面或與壁面發生碰撞彈跳之前需要運動較長距離,而顆粒附近的黏性邊界層和顆粒碰壁間隙處相對于遠處流場需要較密的網格才能解析.因此如果在顆粒運動路徑都采用較密網格顯然會造成很高的計算成本,這顯然是不能接受的.而自適應網格技術很好地替我們解決這一困難,因此本文開發了一套基于浸沒邊界法的2D 自適應網格技術.對于浸沒邊界分層加密策略,主要通過設置三個加密參數來執行.首先是最大加密次數k1,其決定最小網格尺寸為,其中Δh為背景基礎網格尺寸.另外為了保證浸沒邊界附近的網格加密層過度平滑(有助于計算穩定),我們針對在浸沒邊界附近的加密層,設置前k2個加密層中每一個加密層具有N層網格.一個k1=3,k2=2,N=8的加密結果如圖4 所示.

圖4 給出的流場子午面上的加密網格.粗網格和細網格加密界面上的法向和網格中心連線存在一定傾角.針對這一現象,本文采用傾斜校正處理使得計算更加穩定.

圖4 2D 自適應網格Fig.4 An 2D AMR grid around immersed boundary

2 數值算法驗證

為了驗證本文所提出的軸對稱IBM 算法的有效性與可靠性.本章節將分別給出小球繞流、圓盤繞流、Stokes 流小球近壁運動以及小球自由下落碰壁彈跳等四種軸對稱流動的數值模擬驗證結果.

2.1 靜止小球繞流驗證

根據Wu[24]的試驗研究以及Johnson[16]和Tomboulides[25]數值模擬結果,當雷諾數Re=U∞D/νf小于210 時,均勻來流通過一個固定小球所形成的流動是保持軸對稱的.這里U∞,D,νf分別表示遠場均勻來流速度、靜止小球直徑以及流體的運動黏度.因此本文分別測試了雷諾數為1,10,50,100,150,200 等六種不同工況,并將本文的數值模擬結果與前人結果進行比較.首先我們給出計算模型如圖5所示.模型流向總長為 40D,小球中心距離入口和出口截面距離均為 20D,小球中心距離上側壁為 20D.其次本文對于邊界條件設置:進口采用均勻來流邊界條件U∞,出口設置為對流出口邊界條件,對稱軸設置為對稱邊界條件,側壁設置為遠場邊界條件,浸沒小球表面設置為無滑移邊界條件.模擬過程中本文采用自適應時間步長,保證最大的庫朗數(CFL)不超過0.5.

圖5 小球繞流幾何參數Fig.5 Configuration for flow past a fixed sphere

這里考慮到本文的工況最高雷諾數為200 以及Pan[26]的數值經驗,一個網格無關性驗證如表1所示.

表1 雷諾數Re=200 網格無關性驗證Table 1 Grid independence test at Re=200

表中 Δx為網格最小尺寸,Cd=為阻力系數(Fz為流向阻力),Lre表示小球尾流回流區長度.通過粗細網格阻力系數和回流區長度結果對比,我們發現粗網格對于本問題是足夠解析的.因此后續工況的數值結果均是基于粗網格實現.在圖6中,從上到下依次給出雷諾數Re為50,100,150,200 時采用軸對稱IBM 計算得到的流線結果,白色半圓表示浸沒小球.從圖中我們可以看到四種雷諾數下,隨著雷諾數的增加回流區長度也是增加的.

圖6 不同雷諾數流線對比Fig.6 Streamline at different Re with IBM

圖7(a)和圖7(b)中,將采用本文算法計算得到的回流區長度Lre以及阻力系數Cd與Johnson[16]和Magnaudet[27]的結果進行定量對比.根據圖7(a)可以發現:我們計算的回流區長度在雷諾數為50 和100 時與前人結果均符合較好,在雷諾數為150 和200 時與Johnson[16]結果也符合很好.而根據圖7(b)可以發現:我們通過軸對稱IBM 算法計算的阻力系數在低雷諾數和高雷諾數與Johnson[16]和Magnaudet[27]的均符合很好.綜上所述,本文的軸對稱浸沒邊界算法對于求解小球繞流中軸對稱流動是可靠的、有效的.

圖7 小球回流區長度和阻力系數對比Fig.7 Comparisons of the length of recirculation zone Lre and drag coefficient Cd of sphere

2.2 靜止圓盤繞流驗證

除了小球這種特殊幾何結構之外,圓盤也是一種典型的軸對稱結構.為了更好地說明本文的算法對于固體幾何的適應性,我們對于低雷諾數下的軸對稱圓盤繞流也進行了相應的測試.根據文獻[28-29]的數值研究:對于縱橫比χ=D/Hz=10 的圓盤,當雷諾數Re=U∞D/νf小于135 時,遠場均勻來流通過固定圓盤后形成穩態軸對稱流動.這里U∞,D,Hz,νf分別表示遠場均勻來流速度、靜止圓盤徑向直徑、靜止圓盤軸向厚度以及流體的運動黏度.本文分別計算了雷諾數為10,20,40,60,80,100 等六種工況.為了更好地與文獻[28]研究對比,我們在計算模型的設置上采取與之相一致的參數.計算模型如圖8 所示:進口距離圓盤 2.5D,出口距離圓盤 15D,側壁距離對稱軸 6D.

圖8 圓盤繞流幾何參數Fig.8 Configuration for flow past a circular disk

這里關于邊界條件的設置與上一小節的小球繞流中的設置是一樣,便不再復述.我們同樣采用兩種不同分辨率的網格進行計算,用以驗證網格無關性,如表2 所示.

表2 雷諾數Re=100 網格無關性驗證Table 2 Grid independence test at Re=100

根據表2 中的阻力系數以及回流區長度對比,我們發現粗網格的數值計算結果已經足夠求解該問題,所以后續的低雷諾數下的工況我們均采用粗網格進行計算.

圖9 中,我們分別給出了雷諾數為10 和100 兩種工況下的流場流線圖,白色矩形表示浸沒圓盤.我們可以直觀地看出回流區長度也是隨著雷諾數增加而增加的.

圖9 不同雷諾數流線對比Fig.9 Comparison of streamline at different Re

除此之外,我們同樣對比了圓盤尾跡回流區長度Lre以及阻力系數Cd,對比結果如圖10(a)和圖10(b)所示.我們發現二者在不同雷諾數時均與文獻[28]符合較好,并且阻力系數Cd隨著雷諾數的增加是減小的,這與小球繞流結果是相似的.綜上所述,本文的軸對稱IBM 算法對于圓盤幾何的軸對稱繞流問題求解結果也是有效和可靠的,同時也說明本文算法對復雜的軸對稱幾何也是適用的.

圖10 圓盤回流區長度和阻力系數對比Fig.10 Comparisons of the length of recirculation zone Lre and drag coefficient Cd of circular disk

2.3 Stokes 流小球近壁運動驗證

在前面兩個小節中,本文對于靜止浸沒邊界問題已經做了詳細的驗證和討論,充分說明了本文的軸對稱IBM 算法對于各種復雜幾何靜邊界軸對稱流動問題是能夠正確求解的.但是工業工程更多的是動邊界問題,本小節將通過數值模擬的方法計算勻速小球靠近壁面運動所受阻力,并與經典的Brenner[30]小球碰壁潤滑力解析解式(19)和式(20)進行對比

圖11 小球勻速靠近壁面幾何參數Fig.11 Configuration for a sphere with a constant velocity approaching a wall

圖12 給出了不同間隙時的流場壓力云圖.從圖中我們可以看出由于小球不斷靠近壁面,間隙流體被擠壓向徑向外側流動的同時,間隙處的壓力會急劇升高導致小球阻力上升.當間隙特別小時,壓力會變得非常大,相應所需要的網格分辨率也就會很高.本文阻力數值模擬結果與Brenner[30]理論結果對比如圖13 所示:采用了800,1600 和3200 三種不同分辨率的網格進行計算,很明顯三種網格的數值模擬結果在 ε＞0.1 時和理論結果都能符合很好,但是在小球無限接近壁面時,即 ε→0 時結果開始出現偏差,而且網格分辨率越小偏差越大.

圖12 壓力云圖間隙為(a)ε=1.1和(b)ε=0.1Fig.12 Contour of pressure around sphere with (a)ε=1.1 and(b)ε=0.1

圖13 數值模擬阻力和理論結果對比Fig.13 Comparison of numerical and theoretical results

這是因為隨著間隙越來越小,壓力場梯度劇烈變化,想要精確地解析流場,所需要的網格分辨率是非常高的.雖然這種極限小間隙情況很難解析正確,但是在網格分辨率足夠的情況下,我們的結果依然是有效的.這可以說明本文的軸對稱算法對于給定速度的強迫運動產生的動邊界問題也是可以正確求解的.

2.4 小球自由下落碰壁彈跳驗證

上一小節中,本文驗證了強迫運動的動邊界問題,本節我們通過耦合顆粒運動學方程以及流體控制方程來數值模擬小球顆粒在不可壓縮黏性流體中自由下落碰撞壁面過程.針對這樣的一個球壁碰撞運動,Gondret 等[31]通過實驗研究,給出了小球碰撞后的運動軌跡以及相應的恢復系數.文獻[32-33]通過數值模擬的手段再現了Gondret 等[31]的實驗結果.為了驗證本文軸對稱IBM 算法的正確性,我們選取與之相同的工況參數設置并進行結果對比.由于本節問題和上一小節類似,上一小節的計算模型在本節仍然繼續使用.在進行碰撞數值模擬之前,我們需要引入兩個理論模型:一個是潤滑力模型,另一個是碰撞模型.首先對于潤滑力模型,結合上一小節,我們發現間隙越來越小所需網格分辨率越來越高.直接模擬碰撞將會耗費巨大的資源,而且還不能求解準確,因此在小間隙時一個潤滑力修正模型[34-35]的引入是有必要的.再加入潤滑力模型后,顆粒運動方程(6)修正為式(21),其中潤滑力由式(22)[35]和式(23)[36]給出(式(23)為式(20)的近似形式).其中Un為顆粒垂直壁面法向(z軸方向)速度大小,εal為激活潤滑力模型間隙條件.由于本文的模擬參數Re＜210,流場保持為軸對稱狀態,因此小球速度Up只有垂直壁面法向分量Un

這里因為 1/ε 的存在,導致當間隙趨近于零時校正潤滑力趨近于無窮大.這顯然不符合物理實際,事實上由于顆粒表面粗糙度的存在,一般當 ε 接近于表面粗糙度時碰撞已經發生(ε=εw).此時潤滑力不會繼續增大,顆粒進入固體碰撞階段.這一階段流體對于顆粒小球的作用力相比于固體的碰撞力可以忽略.為了解析這樣的碰撞過程,使用一個經典的軟球模型[34],即

其中kn和 ηn分別為剛度系數和阻尼系數且由式(25)給出,en,d表示空氣中的干碰撞恢復系數與顆粒材料有關,一般直接參考實驗數據.根據Gondret 等[31]文獻,本文碰撞算例均采用en,d=0.98 的設置.δs-w表示顆粒與壁面重疊位移.碰撞時間步長 Δtc=[33],Δtf表示流體時間步長.NΔtf為人為設定的碰撞時間,本文中N=10[32-34].實際上碰撞發生的時間是比較短的,這里人為地放大了碰撞接觸時間,避免速度改變過于劇烈導致計算誤差太大

完整小球與壁面正碰撞數值模擬過程如圖14所示:當間隙 ε＞εal時,顆粒只受到來自流體的作用力與重力;當 εw＜ε＜εal時,顆粒除了受到來自流體的作用力及重力之外,還有一項潤滑修正力;當ε＜εw時,顆粒和壁面已經發生碰撞接觸,顆粒運動由方程(24)控制.

圖14 小球與壁面正碰過程Fig.14 Process of a sphere impacting normally on a wall

這里我們對比了沖擊壁面Stokes 數為27 和152兩種工況結果,具體的物理參數如表3 所示.其中Stokes 數Stc=ρpUinD/(9ρfν),雷諾數Re=UinD/νf,Uin表示小球自由下落后在壁面影響前達到的穩定速度.

表3 物理計算參數Table 3 Physical and computational parameters

圖15 (a)和15(b)給出了兩種工況下的數值模擬結果與Gondret 等[31]實驗結果對比.圖中和T*=表示無量綱高度和時間,其中tc為首次碰撞發生時間.我們可以發現小球的軌跡均符合較好.因此本文的軸對稱算法對于耦合顆粒運動學方程的動邊界問題也是準確的.

圖15 顆粒碰撞軌跡Fig.15 Trajectory of the sphere after colliding with a wall

3 結論

本文針對顆粒兩相流中的無環流軸對稱流動問題,發展了基于2D 笛卡爾網格的軸對稱IBM 算法.該IBM 算法通過一個階梯銳利狀界面實現固體浸沒邊界引入,對于界面上的IB 網格單元通過最小二乘法插值得到.本文通過柱坐標系建立流體控制方程,并采用一個有限體積格式的經典投影算法進行方程離散.同時對方程中由于柱坐標系的使用而產生的源項ur/r2采用隱式格式處理.針對動邊界,本文也發展了一個2D 自適應網格技術以提升計算效率.通過小球繞流、圓盤繞流、Stokes 流小球近壁運動以及小球自由下落碰壁彈跳等數值模擬算例,驗證該軸對稱IBM 算法對于復雜幾何邊界以及動邊界軸對稱流動的求解是高效的和準確的.