蔡政剛 潘君華 倪明玖

(中國(guó)科學(xué)院大學(xué)工程科學(xué)學(xué)院,北京 100049)

引言

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

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

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

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

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

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

圖1 浸沒邊界和網(wǎng)格單元?jiǎng)澐质疽鈭DFig.1 Schematic diagram of immersed boundary and dividing cells type

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

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

1 數(shù)值模型

1.1 控制方程

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

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

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

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

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

1.2 數(shù)值離散

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

圖2 柱坐標(biāo)系下的體積微元Fig.2 Volume element in cylindrical coordinates

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

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

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

第三步計(jì)算k+1時(shí)層壓力場(chǎng).首先根據(jù)上一步計(jì)算得到的網(wǎng)格中心的中間步速度場(chǎng),我們通過插值可以得到網(wǎng)格單元面上中間步質(zhì)量通量

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

第四步根據(jù)k+1 時(shí)層的壓力場(chǎng),修正中間步速度場(chǎng),得到k+1 時(shí)層的網(wǎng)格單元中心和面上的速度場(chǎng)以及面上的質(zhì)量通量分別為

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

另外對(duì)于在不可壓縮黏性流體中顆粒小球垂直靠近壁面運(yùn)動(dòng)時(shí),由于運(yùn)動(dòng)顆粒和靜止壁面之間的間隙很小,數(shù)值模擬所需的網(wǎng)格分辨率非常高.因此為了保證計(jì)算的穩(wěn)定性,需要在一個(gè)時(shí)間步長(zhǎng)內(nèi)重復(fù)多次執(zhí)行投影步計(jì)算,如圖3 所示.

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

根據(jù)本文測(cè)試的經(jīng)驗(yàn),這里的迭代次數(shù)一般設(shè)置為3～5 次即可.另外在前兩次迭代計(jì)算時(shí),一般會(huì)對(duì)壓力和速度做松弛處理來保證計(jì)算的穩(wěn)定性.

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

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

1.3 IB 邊界條件引入

因?yàn)楸疚牟捎秒x散力法,所以固體浸沒邊界的引入主要分為三個(gè)步驟.第一步通過最小二乘法插值得到IB 網(wǎng)格單元的值,第二步利用階梯狀銳利界面替代真實(shí)的固體浸沒邊界封閉控制方程求解流場(chǎng)信息,第三步利用新求解的流場(chǎng)信息修正IB 網(wǎng)格單元信息.接下來我們將分別說明這三個(gè)步驟的實(shí)現(xiàn).

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

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

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

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

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

這里IB 單元被當(dāng)做源項(xiàng)處理放在公式右側(cè).

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

1.4 2D 自適應(yīng)網(wǎng)格技術(shù)實(shí)現(xiàn)

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

圖4 給出的流場(chǎng)子午面上的加密網(wǎng)格.粗網(wǎng)格和細(xì)網(wǎng)格加密界面上的法向和網(wǎng)格中心連線存在一定傾角.針對(duì)這一現(xiàn)象,本文采用傾斜校正處理使得計(jì)算更加穩(wěn)定.

圖4 2D 自適應(yīng)網(wǎng)格Fig.4 An 2D AMR grid around immersed boundary

2 數(shù)值算法驗(yàn)證

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

2.1 靜止小球繞流驗(yàn)證

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

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

這里考慮到本文的工況最高雷諾數(shù)為200 以及Pan[26]的數(shù)值經(jīng)驗(yàn),一個(gè)網(wǎng)格無關(guān)性驗(yàn)證如表1所示.

表1 雷諾數(shù)Re=200 網(wǎng)格無關(guān)性驗(yàn)證Table 1 Grid independence test at Re=200

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

圖6 不同雷諾數(shù)流線對(duì)比Fig.6 Streamline at different Re with IBM

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

圖7 小球回流區(qū)長(zhǎng)度和阻力系數(shù)對(duì)比Fig.7 Comparisons of the length of recirculation zone Lre and drag coefficient Cd of sphere

2.2 靜止圓盤繞流驗(yàn)證

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

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

這里關(guān)于邊界條件的設(shè)置與上一小節(jié)的小球繞流中的設(shè)置是一樣,便不再?gòu)?fù)述.我們同樣采用兩種不同分辨率的網(wǎng)格進(jìn)行計(jì)算,用以驗(yàn)證網(wǎng)格無關(guān)性,如表2 所示.

表2 雷諾數(shù)Re=100 網(wǎng)格無關(guān)性驗(yàn)證Table 2 Grid independence test at Re=100

根據(jù)表2 中的阻力系數(shù)以及回流區(qū)長(zhǎng)度對(duì)比,我們發(fā)現(xiàn)粗網(wǎng)格的數(shù)值計(jì)算結(jié)果已經(jīng)足夠求解該問題,所以后續(xù)的低雷諾數(shù)下的工況我們均采用粗網(wǎng)格進(jìn)行計(jì)算.

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

圖9 不同雷諾數(shù)流線對(duì)比Fig.9 Comparison of streamline at different Re

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

圖10 圓盤回流區(qū)長(zhǎng)度和阻力系數(shù)對(duì)比Fig.10 Comparisons of the length of recirculation zone Lre and drag coefficient Cd of circular disk

2.3 Stokes 流小球近壁運(yùn)動(dòng)驗(yàn)證

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

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

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

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

圖13 數(shù)值模擬阻力和理論結(jié)果對(duì)比Fig.13 Comparison of numerical and theoretical results

這是因?yàn)殡S著間隙越來越小,壓力場(chǎng)梯度劇烈變化,想要精確地解析流場(chǎng),所需要的網(wǎng)格分辨率是非常高的.雖然這種極限小間隙情況很難解析正確,但是在網(wǎng)格分辨率足夠的情況下,我們的結(jié)果依然是有效的.這可以說明本文的軸對(duì)稱算法對(duì)于給定速度的強(qiáng)迫運(yùn)動(dòng)產(chǎn)生的動(dòng)邊界問題也是可以正確求解的.

2.4 小球自由下落碰壁彈跳驗(yàn)證

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

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

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

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

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

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

表3 物理計(jì)算參數(shù)Table 3 Physical and computational parameters

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

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

3 結(jié)論

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