非球形顆粒兩相流的數(shù)值模擬研究進展

崔智文 王 澤 蔣新宇 趙立豪

清華大學航天航空學院工程力學系, 北京 100084

1 引 言

顆粒兩相流常見于工業(yè)生產、自然環(huán)境和生物醫(yī)學等領域. 20 世紀以來, 顆粒兩相流相關的研究主要集中在球形顆粒問題, 且現(xiàn)在依然是兩相流領域中最活躍的方向之一. 在這些研究中顆粒的形狀因素通常被忽略, 而簡化成球形, 并且研究關注的重點主要是顆粒在流動中的聚集和輸運現(xiàn)象. 然而, 在實際流動問題中顆粒形狀往往不規(guī)則, 顆粒的非球形特性對顆粒取向與轉動行為的影響不可忽略, 比如在紙張制造過程中紙纖維的取向(Lundell et al. 2011a, H?kansson et al. 2014)以及纖維增強材料制備過程中纖維的排列(Cheung et al. 2009, Soutis 2005)均會對最終成品的力學性能產生重要的影響; 大量的實驗與數(shù)值模擬結果表明, 纖維以及聚合物對壁湍流的減阻和調制作用與其取向行為密切相關(Paschkewitz et al. 2004, 2005a, 2005b); 不同環(huán)境溫度和濕度下, 形成的冰晶具有不同的形狀, 會影響冰晶的碰撞融合以及雨雪的形成(Heymsfield 1977, Jucha et al. 2018); 海洋和湖泊中浮游生物往往具有形狀多樣性等(Pedley & Kessler 1992,Saintillan 2018, Qiu et al. 2019, Ardekani et al. 2017a).

近代有關非球形顆粒兩相流的理論研究最早可以追溯到Jeffery(1922)在Stokes 流中推導的橢球顆粒受剪切作用的力矩方程, 隨后在Batchelor(1970)和 Brenner(1961, 1963b)的努力下, 逐漸完成了Stokes 流動中橢球顆粒的阻力及顆粒對流體產生的額外應力(Batchelor 1970, Brenner 1974)的理論表達式. 由于當時計算機計算能力的限制, 這些研究主要集中在橢球形顆粒在極低雷諾數(shù)流動中運動的理論或實驗研究. 隨著計算機硬件水平的飛速提升以及高精度數(shù)值模擬方法的快速發(fā)展, 復雜流場中大規(guī)模非球形顆粒運動的高精度數(shù)值模擬得以實現(xiàn). 研究的顆粒形狀、尺寸以及流場類型日益多樣化. 然而, 相較于球形顆粒兩相流的研究, 非球形顆粒兩相流的數(shù)值模擬研究依然處于起步階段, 而非球形顆粒的各向異性的形狀特性也引入了更多的運動行為 (比如轉動、取向等)以及顆粒參數(shù)范圍(比如顆粒形狀等). 這使得非球形顆粒兩相流問題變得更加復雜且難以分析, 進而給非球形顆粒兩相流的研究帶來了較大的挑戰(zhàn). 目前, 有關非球形顆粒兩相流的數(shù)值模擬研究主要集中以下三個方面: (1) 單個顆粒行為的研究, 建立對非球形顆粒在不同流場中的運動學和動力學特征以及顆粒與流場之間的相互作用的基本認識, 輔助點顆粒模型的建立; (2) 多顆粒行為研究, 以了解顆粒與顆粒、顆粒與流場間的相互作用, 以及建立相關模型; (3) 研究大規(guī)模顆粒群在復雜流場 (如湍流) 中的運動行為, 通過統(tǒng)計學方法描述顆粒在流場中的行為特征, 以建立與完善非球形顆粒的歐拉模型.

目前, 非球形顆粒兩相流的數(shù)值模擬主要基于歐拉-拉格朗日的求解框架. 在該框架下, 流體相的Navier-Stokes 方程在歐拉網(wǎng)格上進行求解, 而顆粒相則通過拉格朗日追蹤方法去描述每一個顆粒的行為. 顆粒在拉格朗日觀點下的求解能準確地描述顆粒的平動、取向與轉動行為, 進而可以幫助研究者更加直觀地認識非球形顆粒復雜的運動行為與流場的關系. 基于該框架, 根據(jù)不同的顆粒計算方法, 顆粒相的模擬可劃分為點顆粒方法與全分辨顆粒方法. 前者用質點替代具體形狀的顆粒, 流體與顆粒之間的相互作用均通過相應的理論模型實現(xiàn). 而后者需要完全分辨顆粒與流體的界面, 流體與顆粒之間的相互作用通過全分辨數(shù)值模擬直接實現(xiàn). 近些年來, 國內外學者使用這兩類方法進行了大量直接數(shù)值模擬的研究, 使得人們對非球形顆粒在復雜流動中的行為規(guī)律有了初步認識.

本文針對非球形顆粒兩相流數(shù)值模擬的研究進展進行綜述. 文章的主體分為三個部分: 第一部分介紹非球形顆粒兩相流研究的基本概念、數(shù)值方法和理論模型; 第二部分介紹非球形顆粒在簡單基本流場(均勻來流、靜止流場和剪切流)和復雜湍流場(均勻各向同性湍流和壁湍流)中的運動學行為; 第三部分針對湍流中非球形顆粒取向與轉動行為、非球形顆粒對湍流的減阻調制作用的機理展開討論. 最后, 對未來值得研究的工作提出一些建議.

2 數(shù)值模擬方法

在本綜述討論的顆粒兩相流中, 顆粒均為離散相. 本節(jié)將主要介紹歐拉-拉格朗日方法中非球形顆粒的描述方法、理論模型以及相關的數(shù)值模擬手段. 通常, 流體的求解可以采用有限差分方法、有限體積方法、譜方法等常見的數(shù)值方法. 而顆粒的求解根據(jù)建模的方法可以大致分為點顆粒方法與全分辨顆粒方法. 點顆粒方法用質點替代具體的顆粒, 流體與顆粒之間的相互作用均通過相應的力學模型實現(xiàn). 而全分辨顆粒方法則完全分辨顆粒的形狀, 流體與顆粒之間的相互作用也通過數(shù)值模擬直接實現(xiàn). 兩種方法的使用條件主要由研究的主要問題確定. 通常點顆粒方法適合顆粒尺寸小于流體的最小特征尺度, 比如湍流中的Kolmogorov 長度尺度, 此時顆粒雷諾數(shù)較小, 顆粒尾部不出現(xiàn)明顯的渦脫落, 顆粒的受力以及顆粒對流場的影響可以通過點顆粒模型比較準確地給出. 而全分辨顆粒方法適合顆粒尺寸大于流體最小特征尺寸, 以捕捉點顆粒模型無法表征的流體慣性特征和對流場的反饋影響. 點顆粒方法通常使得單顆粒的計算量較低, 因此,非常適合大量顆粒 (百萬量級) 的計算, 但點顆粒的模型適用范圍有限, 例如顆粒尺寸足夠小、顆粒雷諾數(shù)須滿足模型要求等. 而全分辨方法通?？梢暂^為精確地模擬單個顆粒的行為, 但每個顆粒的計算成本高, 尤其是對于尺寸遠小于流場尺度的顆粒, 難以實現(xiàn)大規(guī)模顆粒的計算. 目前這兩種方法在研究中是相輔相成的關系, 一方面可以通過全分辨方法提煉點顆粒的模型, 另一方面可利用完善后的點顆粒模型進行大規(guī)模顆粒兩相流的模擬.

2.1 非球形顆粒的基本方程

2.1.1 坐標系的定義

為了準確描述非球形顆粒的行為, 一般需要描述顆粒空間位置以及空間的取向(或姿態(tài)).在本文使用三種參考系來描述單個顆粒的運動 (如圖1 所示) , 它們分別是:

(1) 慣性參考系 : 用來描述連續(xù)流場以及顆粒的平動行為.

(2) 顆粒參考系 : 以顆粒質心位置為原點, 軸分別對齊顆粒的主軸. 顆粒的轉動行為通常在該參考系下求解.

(3) 隨動參考系 : 以顆粒質心位置為原點,x′′,y′′,z′′軸分別與慣性參考系中的x,y,z軸平行. 與慣性參考系之間相差一個平動位移.

其中, 顆粒參考系與隨動參考系之間通過一個旋轉張量R進行轉換, 即

圖1描述顆粒運動的參考系: 慣性參考系 〈 x,y,z〉 , 顆粒參考系 〈 x′,y′,z′〉 以 及隨動參考系〈x′′,y′′,z′′〉

同時, 顆粒的取向、轉動與速度均由矢量進行描述. 對于矢量而言, 其在慣性參考系與隨動參考系下的表達相同. 為了簡化標記, 在不引入歧義的情況下, 通常用慣性參考系的描述方法來取代隨動參考系, 進而對顆粒的運動行為以及相關變量進行描述. 因此, 對于慣性參考系描述的矢量p滿足

其中‘(·)′’代表顆粒參考系下的變量, 下同. 對于慣性參考系中的二階張量H滿足

2.1.2 平動與轉動方程

根據(jù)牛頓第二定律, 顆粒的平動與轉動方程分別表示為

式中,mp為顆粒質量,I為顆粒的轉動慣量,v為顆粒的平動速度,為顆粒的轉動角速度,F與T分別為作用在顆粒上的合外力與合外力矩. 顆粒的轉動行為一般在顆粒坐標系中進行求解, 即

如果顆粒的合外力與合外力矩可以準確得到, 顆粒的運動行為可以通過式(1)與式(2)進行求解得到.

2.1.3 顆粒取向求解

從顆粒參考系到慣性參考系的旋轉矩陣可由歐拉四元數(shù)E=[e0,e1,e2,e3] 來表示, 即

歐拉四元數(shù)隨時間推進的方程為

式中

2.1.4 顆粒形狀

非球形顆粒通常具有復雜的外形. 在非球形顆粒的研究中, 主要使用橢球顆粒模型對非球形顆粒進行近似. 主要的原因包括:

(1) 橢球顆粒是最簡單的非球形顆粒模型. 近一個世紀以來, 科學家們通過橢球坐標系變換以及相關橢球諧函數(shù)等數(shù)學方法得到一些橢球顆粒的作用力的解析表達式. 因此, 橢球顆粒作為非球形顆粒動力學的理論研究的基礎模型, 擁有比較強大的數(shù)學工具. 盡管如此, 非球形顆粒兩相流的研究道路依然非常困難, 目前得到的理論模型依然有限.

(2) 大量的直接數(shù)值模擬以及實驗研究表明, 橢球顆粒模型能比較好地捕捉到部分非球形顆粒(非橢球形狀)的運動特性, 比如圓盤、纖維、圓柱等. 研究橢球顆粒的運動行為能夠幫助人們認識形狀的各向異性所帶來的復雜動力學行為.

在本文中, 橢球顆粒表示為

2.1.5 無量綱參數(shù)

在非球形顆粒兩相流中常用的無量綱參數(shù)如下.

(1) 顆粒雷諾數(shù)Rep: 基于顆粒與流體之間相對滑移的速度與顆粒直徑定義的無量綱數(shù), 即Rep=2|uf-v|rp/ν. 其中rp為顆粒的參考半徑(通常選取等效半徑或者顆粒的最大半徑).

(3) 平動Stokes 數(shù)St: 用來反映顆粒的平動慣性,St=τp/τf.τp與τf分別表示顆粒的平動弛豫時間與流體特征時間尺度.

(4) 轉動Stokes 數(shù)Str: 主要在剪切流中使用, 用來反映顆粒的轉動慣性,即Str=ρpsrp2/μ.

圖2橢球顆粒模型示意圖.(a)三軸不等顆粒, (b)桿狀顆粒, (c)碟狀顆粒

(5) 密度比D: 表示顆粒與流體密度之比, 即D=ρp/ρf.

需要注意的是在不同的文獻中雷諾數(shù)的定義稍有差別. 在絕大部分問題中特征長度選為顆粒的直徑, 但在有些理論推導的模型中往往選取顆粒的半徑作為特征長度. 本文已經(jīng)將該差別進行了統(tǒng)一, 但為了使形式與原文獻保持一致以方便讀者查閱, 一些地方保留了原文獻中的定義,相關的差別將會在需要的時候提及.

2.2 點顆粒方法

2.2.1 基本思路

點顆粒, 顧名思義, 就是用質點替代復雜形狀顆粒且同時考慮顆粒取向與轉動行為. 如果研究的顆粒的尺寸遠小于流動的特征尺度, 點顆粒是一種非常有效的方法(Maxey 2017). 該方法通常需要適當?shù)睦碚摶蚪?jīng)驗模型來描述流體與顆粒之間的相互作用, 顆粒與流體的相互作用主要受到顆粒相的體積分數(shù)和質量分數(shù)的影響(Balachandar & Eaton 2010). 常見的耦合方式如下.

(1) 單向耦合: 當體積分數(shù)和質量分數(shù)很小時, 即顆粒稀疏懸浮, 顆粒相對流體的反饋作用可以忽略, 只需考慮流體對顆粒的作用, 即流體作用在顆粒上的力與力矩;

(2) 雙向耦合: 當體積分數(shù)或質量分數(shù)較大時, 顆粒相對流場的反饋作用不能忽略, 此時需要考慮顆粒對流體的作用并進行耦合, 比如力耦合(Crowe et al. 1977)、力矩耦合(Andersson et al.2012)以及應力耦合(Batchelor 1970, Brenner 1974);

(3) 四向耦合: 當顆粒相的體積分數(shù)進一步增大, 顆粒間的相互作用 (如碰撞, 團聚和破碎等) 也需要考慮. 此外, 還需要考慮顆粒與顆粒之間的相互作用力(比如靜電力等).

點顆粒求解的基本流程與思路如下.

(1) 獲取顆粒處未擾動流場. 當一步流場更新后, 需要通過適當?shù)牟逯捣椒?比如線性插值、二階拉格朗日插值等)從流場中獲取顆粒位置處的未擾動流場的速度與速度梯度. 對于非球形顆粒而言, 顆粒位置處的速度梯度信息是非常必要的.

(2) 代入合適的顆粒受力模型, 計算顆粒的動力學方程.

(3) 如果進行雙向耦合, 則需要計算顆粒對流場的反饋, 并將反饋力、力矩與應力耦合到流場中.

(4) 如果進行四向耦合, 則需要進行顆粒與顆粒的接觸檢測, 并施加相應的碰撞模型或者相互作用模型.

非球形顆粒兩相流中常見的流體對顆粒的作用力以及顆粒對流體的反饋力模型將在第2.2.2 和2.2.3 小節(jié)進行介紹. 因為非球形顆粒之間的相互作用模型及機制相對復雜, 目前相關研究工作較少, 本文暫不討論顆粒與顆粒之間的相互作用, 而有關顆粒與顆粒碰撞的模型可以參見第2.3.4 小節(jié).

2.2.2 流體對顆粒的作用

流體對顆粒的作用通常體現(xiàn)為力與力矩的作用.

(1)力的作用

式中,Ks為無量綱Stokes 阻力系數(shù)矩陣,μ為流體的動力黏性系數(shù),uf為顆粒處流體速度,為顆粒速度,rp為顆粒的參考半徑(通常采用最大半徑或者等效半徑). 對于任意形狀的微小橢球顆粒, 其Stokes 阻力系數(shù)矩陣是對稱正定的(Brenner 1963b). 對于橢球顆粒模型, 在顆粒坐標系下, Stokes 阻力系數(shù)矩陣對角線上的元素非零, 且表示為

對于具有回轉軸的橢球顆粒而言, 令形狀參數(shù)ka=kb=1 ,kc=λ, 其中λ定義為顆?；剞D軸長度與赤道軸長度之比. 那么與形狀相關的參數(shù)α,β,γ和χ可以解析地用λ表達.

由于非球形顆粒的阻力系數(shù)矩陣的各向異性, Stokes 的阻力方向通常不與滑移速度共線. 與球形顆粒相比, 非球形顆粒受到的阻力不僅依賴顆粒與流體之間的速度滑移, 而且還取決于顆粒的取向.

同時, 為了反映顆粒趨近流體速度所需要的特征時間, 通常使用顆粒的弛豫時間進行描述,記為τp. 一般定義τp/τf為顆粒的Stokes 數(shù), 記為St.當St趨近于零時, 意味著顆粒幾乎完全跟隨流體運動. 為了評估非球形顆粒的弛豫時間, 在湍流問題中, 假設顆粒在取向隨機時(Brenner 1963b,Shapiro & Goldenberg 1993), 評估顆粒運動對流場變化的響應. 此時, 顆粒的弛豫時間τp可以定義為

(c) Saffman 升力

在顆粒雷諾數(shù)有限小的情況下, 強剪切會對運動的顆粒產生一個側向的升力, Saffman(1965, 1968)最早通過奇異攝動理論得到了球形顆粒在簡單剪切流中的升力模型, 由此稱為Saffman 升力. 該模型預測的運動與實驗中觀察到的結果符合較好(Taylor 1923), 因此該模型逐漸得到了推廣與應用. McLaughlin (1991), Bagchi 和Balachandar (2002), Magnaudet 與他的合作者們(Legendre & Magnaudet 1998, Magnaudet 2003, Magnaudet et al. 2003)對球形顆粒的Saffman 升力進行了改進以適應更加廣泛的剪切流場. 而Harper 和Chang (1968)最早將Saffman 升力拓展到非球形顆粒, 并發(fā)現(xiàn)相關的升力系數(shù)矩陣與顆粒形狀無關. 隨后Fan 和Ahmadi (1995),Feng 和Kleinstreuer (2003), Cui Y 等 (2018, 2019, 2020) 在Harper 和Chang (1968)的基礎上繼續(xù)發(fā)展非球形顆粒的剪切誘導升力模型. 其中Cui Y 等(2020) 最近提出了一個比較高效求解桿狀顆粒剪切升力的近似模型

然而, 目前的剪切誘導的升力模型與全分辨模擬的結果依然存在比較明顯的差距(Cui Y et al. 2019). 至于該模型式(14)是否適合碟狀顆粒依然有待研究. 目前在絕大部分點顆粒非球形顆粒兩相流的理論研究問題中均未考慮剪切升力模型. 升力對非球形顆粒的動力學行為的影響依然有待進一步討論.

(d) 重力與浮力

當考慮顆粒沉降等問題, 或顆粒Froude 數(shù)小于1, 即Frp=U/gτp ＜1 , 其中U為流體的特征速度,而g為重力加速度, 此時需要進一步考慮重力和浮力的影響(Sardina et al. 2012, Milici et al. 2014)

式中Fg為顆粒受到的重力與浮力的合力, 與mf分別為顆粒質量和顆粒排開流體的質量,g為重力加速度矢量.

(e) 經(jīng)驗公式

然而, 在一些實際工程問題中, 顆粒的雷諾數(shù)通常遠大于1, 此時顆粒的受到的流體作用力無法通過理論的辦法獲得. 為了能夠繼續(xù)采用點顆粒方法, 就需要通過實驗或者直接數(shù)值模擬的數(shù)據(jù)分析建立相應的經(jīng)驗模型. 而關于非球形顆粒受流體力和力矩的經(jīng)驗模型非常多, 在本文中不進行詳細展開,具體可見Ouchene et al. (2015, 2016). 本文主要介紹一種基于橢球顆粒的經(jīng)驗模型(Zastawny et al. 2012). 不同的經(jīng)驗模型在形式上都具有相似性, 不同點在于系數(shù)的組合形式. 關于力的模型, 主要被分為阻力FD與升力FL. 在經(jīng)驗模型中, 顆粒阻力與滑移速度共線, 而升力作用在與滑移速度垂直的方向. 其對應的阻力系數(shù)與升力系數(shù)定義為

在式(19)與式(20)中的擬合系數(shù)ai與bi可參見(Zastawny et al. 2012). 其他相關的模型形式類似, 具體可參見(Mand? & Roséndahl 2010, H?lzer & Sommerfeld 2008).

(2) 力矩作用

(a) 黏性力矩

對于中心對稱(即具有三個對稱面)的非球形顆粒, 在Stokes 流的情況下, 其受到的黏性力矩為(Brenner 1974)

上式中,KR為旋轉阻力系數(shù)張量,Θ為剪切阻力系數(shù)張量, 其中Θ為三階張量,KR與Θ均與顆粒的形狀有關. 對于任意三軸不等橢球顆粒而言, 在顆粒參考系中受到的黏性力矩為(Jeffery 1922)

對于軸對稱橢球顆粒, 形狀有關的參數(shù)α,β與γ由表1 的表達式給出.

表1 軸對稱橢球顆粒相關的形狀參數(shù)

(b) 流體慣性力矩

(c) 經(jīng)驗公式

對于高雷諾數(shù)問題, 可將力矩分解成俯仰力矩Tp(pitch torque)與旋轉力矩TR(rotation torque). 俯仰力矩主要是當顆?；剞D軸的方向與流場速度方向不共線時產生, 而旋轉力矩為顆粒相對于流體旋轉時流場對顆粒產生的力矩(Zastawny et al. 2012). 在經(jīng)驗公式模型中, 通過定義俯仰力矩系數(shù)為

在Zastawny 等(2012)的工作中, 該系數(shù)與雷諾數(shù)以及攻角的關系為

上式中擬合系數(shù)c1～c10可參見Zastawny 等(2012)的工作. 旋轉力矩系數(shù)定義為

式中Ωs為流體與顆粒之間的相對滑移角速度, 即Ωs=Ω-ωp. 該系數(shù)與顆粒旋轉雷諾數(shù)的關系為

式中顆粒旋轉雷諾數(shù)為ReR=4rp2|Ωs|/ν, 系數(shù)r1～r4參見Zastawny 等(2012)的工作. 其它相關的擬合形式的力矩形式可參見文獻(Ouchene et al. 2015, 2016; Yin et al. 2003).

(3) 其他受力模型

根據(jù)經(jīng)典的Maxey-Riley 小球受力模型(Maxey & Riley 1983), 非定常流動中點顆粒模型還需要進一步考慮壓力梯度力、附加質量力、Basset 歷史力效應項等. 當顆粒密度遠大于流體密度時, 壓力梯度力與附加質量力的作用可以忽略不計. 對于顆粒密度接近流體密度或者小于流體密度時, 壓力梯度力以及附加質量力的作用不可忽略. 在當前綜述的工作中描述的點顆粒均是密度比遠大于1 的情況, 暫不考慮壓力梯度力與附加質量力的影響. Basset 力反映了顆粒在加速過程中邊界層發(fā)展滯后所導致的“歷史”效應的作用, 但往往因為其在計算中實現(xiàn)困難以及該項僅在具有較大滑移速度變化時比較重要, Basset 力通常在計算中被忽略. 近期的研究文章(Olivieri et al. 2014, Daitche 2015, Prasath et al. 2019)表明Basset 力在顆粒輸運中的重要影響. 然而,Basset 歷史力效應項是基于球形顆粒得到的, 對于非球形顆粒而言, 目前還沒有針對一般橢球顆粒的模型. Lawrence 和Weinbaum (1986)通過攝動理論對近球形顆粒的小偏心率進行展開, 推導了近球形橢球顆粒的歷史力模型, 并發(fā)現(xiàn)對于非球形顆粒而言, 除了Basset 歷史效應力外還存在新的歷史力項. 相關的推導與解釋可參見論文Lawrence 和Weinbaum (1986, 1988). 目前,在大部分非球形點顆粒模擬中, 依然不考慮歷史力效應, 而歷史力對非球形顆粒的影響還有待進一步分析.

2.2.3 顆粒對流體的作用

顆粒對流體的反饋作用主要通過力耦合、力矩耦合和應力耦合等方法實現(xiàn), 對應的實現(xiàn)方法如下.

(1) 力耦合

力耦合是基于牛頓第三定律的傳統(tǒng)雙向耦合方法, 廣泛應用于慣性顆粒與流體相互作用的數(shù)值模擬中. 該方法的基本思想是將顆粒所受到的流體作用力作為顆粒對流體的反饋力. 反饋力的實現(xiàn)是將某個網(wǎng)格內所有顆粒的反作用力 -F進行體積平均 (或者加權平均) 以得到作用于流體相的歐拉網(wǎng)格點的體積力( Crowe et al. 1977, Maxey et al. 1997, Boivin et al. 1998, Vreman 2015 )

其中,Δ為網(wǎng)格體積,Fil為該網(wǎng)格體積內第l個顆粒所受流體作用力,np為該網(wǎng)格體積內的顆粒數(shù). 將體積力fip添加到流體動量方程的右端項, 便實現(xiàn)了基于力耦合方法的雙向耦合.

(2) 力矩耦合

Andersson 等 (2012)在力耦合的基礎上提出了力矩耦合方法, 認為除了將顆粒的點力作為反饋力外, 還應該考慮顆粒所受力矩對流體的反作用. 和力耦合方法相似, 單個顆粒受到流體施加的力矩N, 則該顆粒對周圍流體產生反饋力矩 -N, 其可表示為作用于流體微團的表面力. 體積為Δ的網(wǎng)格內, 顆粒產生的反對稱應力為

(3) 應力耦合

不同于大密度比的慣性顆粒, 小尺寸的近中性懸浮顆?？梢暈闊o慣性顆粒, 因此顆粒對流體的反饋力和力矩可以忽略. 但剛性顆粒由于不可變形, 還會通過額外的顆粒應力影響附近流體.稀疏懸浮液中的這種應力可以通過對應力子(stresslet)進行局部體積平均得到(Batchelor 1970).盡管應力子不像力矩和點力那樣與顆粒的動力學直接相關, 但應力子實質上代表剛性顆粒對流體應變運動的抵抗 (Guazzelli & Morris 2011). Brenner (1974)推導了一般軸對稱顆粒的應力子形式, 這里只簡單介紹軸對稱橢球顆粒的對應表達式

2.2.4 無慣性橢球顆粒模型

式中,vs為顆粒的Stokes 沉降速度. 對于小慣性顆粒的沉降問題, 顆粒的速度可以通過顆粒當?shù)氐牧鲌鏊俣燃由铣两邓俣全@得(Qiu et al. 2019). 若不考慮重力作用, 顆粒近似跟隨流體軌跡線.

此時顆粒的轉動, 也直接由當?shù)氐乃俣忍荻扰c顆粒取向決定

類似慣性顆粒的求解, 顆粒的取向相關的歐拉四元數(shù)由式(5)推進.

對于軸對稱的橢球顆粒, 顆粒的取向行為也可以直接通過求解顆粒回轉軸方向方程獲得, 即

2.3 全分辨顆粒方法

對于尺寸大于流動最小尺度 (kolmogorov 尺度) 的顆粒, 以及一些形狀不規(guī)則的顆粒, 它們與流體的相互作用就難以再用點顆粒的方法描述了, 這時就需要借助全分辨的方法來對顆粒兩相流進行模擬. 這種模擬方法在顆粒表面施加無滑移邊界條件, 完全分辨出顆粒的形狀、尺寸,進而分辨出顆粒周圍的流場以及顆粒-流體間的相互作用. 基于完全分辨的方法研究顆粒兩相流能夠最準確地捕捉顆粒-流體之間的相互作用(Tenneti & Subramaniam 2014), 同時也可以用于對基于點顆粒的方法進行建模、驗證, 具有重要的學術價值和研究意義.

完全分辨的顆粒兩相流問題中顆粒通常被視為剛體, 其與流體的相互作用就可以描述為流體中自由運動剛體的流固耦合問題, 如圖3 所示, 其中Si代表顆粒表面, 小箭頭代表顆粒運動方向. 對于該系統(tǒng), 流體相用基于歐拉觀點的不可壓Naiver-Stokes(NS)方程描述, 顆粒相由基于拉格朗日觀點的動力學方程(1)和(2)描述, 在不考慮其他外力 (如電磁力、電場力等) 的情況下,方程中的右端項合外力、合外力矩的形式為F=FH+Fint,T=TH+Tint, 其中FH,TH表示流體對顆粒的力和力矩,Fint,Tint表示顆粒之間相互碰撞的力和力矩. 而顆粒表面為流體提供了無滑移邊界條件

式中uF代表顆粒表面的速度, 和ω分別為顆粒中心點平動速度和轉動角速度,X為顆粒表面點的坐標, 而Xc為顆粒中心點坐標. 所以全分辨方式模擬顆粒兩相流自然是雙向 (不考慮顆粒碰撞) 或四向 (考慮顆粒碰撞) 耦合系統(tǒng).

對于全分辨顆粒兩相流問題的求解主要有兩類方法, 即貼體網(wǎng)格類方法和非貼體網(wǎng)格類方法(Haeri & Shrimpton 2012, Kim & Choi 2019). 貼體類方法的思想很簡單, 即在每一時刻都根據(jù)顆粒的位置劃分流體網(wǎng)格, 并且根據(jù)顆粒的位置和當前運動設置流場的邊界條件, 然后在流體區(qū)域內求解NS 方程. 到下一個時間步, 顆粒的位置發(fā)生變化, 則需要對空間網(wǎng)格進行重新劃分以對流動進行求解, 如圖4(a)所示. 這類方法的典型代表就是任意-拉格朗日-歐拉法 (ALE) .而后一類方法則通常需要兩套網(wǎng)格, 即描述流場的背景笛卡爾網(wǎng)格和描述顆粒表面的拉格朗日網(wǎng)格, 如圖4(b)所示的這類方法中, 背景的笛卡爾網(wǎng)格并不會因為顆粒的運動而發(fā)生改變, 因此不需要流體網(wǎng)格隨顆粒運動的不斷重新生成. 但是需要通過一些特殊的手段來處理流體-固體界面的邊界條件. 這類方法的典型代表包括浸沒邊界方法 (immersed boundary method, IBM), 虛擬區(qū)域方法 (fictitious domain method, FDM) 和格子-玻爾茲曼方法 (lattice Boltzmann method,LBM). 貼體類方法在追蹤運動的顆粒的過程中要不斷更新求解流體方程的網(wǎng)格, 這個過程需要很大的計算量, 所以這類方法對顆粒兩相流的數(shù)值模擬應用目前還僅限于二維(Haeri & Shrimpton 2012), 這里不再展開. 而非貼體類方法省去了網(wǎng)格更新過程的計算量, 同時背景的笛卡爾網(wǎng)格的簡潔性使得流體相NS 方程求解過程中的一些快速算法得以應用. 這兩點使得該方法的計算效率得到顯著提升, 因而這類方法是目前模擬運動物體流固耦合 (如顆粒兩相流) 問題的主流方法. 以下對非貼體類的浸沒邊界方法和格子-玻爾茲曼方法進行展開介紹. 由于虛擬邊界法與浸沒邊界法思路相似, 這里不再展開描述, 相關詳細過程可參考(Yu & Shao 2007b, 2010).

2.3.1 浸沒邊界方法

浸沒邊界方法最早是由Peskin 提出用來模擬心臟瓣膜流固耦合問題(Peskin 1972, 1977). 如圖5 所示, 浸沒邊界法使用背景的歐拉網(wǎng)格和物體表面的拉格朗日網(wǎng)格分別描述流體和物體(Griffith & Patankar 2020). 對于流體-固體界面的邊界條件, 該方法在NS 方程動量方程中施加一個額外的體積力項來實現(xiàn), 即將原本NS 方程的動量方程修改為

圖3全分辨顆粒兩相流問題示意圖

圖4網(wǎng)格方法示意圖. (a)貼體網(wǎng)格類方法, (b)非貼體類網(wǎng)格方法

圖5浸沒邊界方法的歐拉網(wǎng)格和拉格朗日網(wǎng)格

式中f稱為浸沒邊界力 (IB 力) . 適當?shù)厥┘覫B 力后就可以使得流固無滑移邊界條件即式(39)得到滿足. 根據(jù)IB 力的施加方法, 浸沒邊界法可以分為兩類, 即連續(xù)力法和離散力法(Mittal &Iaccarino 2005). 連續(xù)力法中的IB 力先在物體表面的拉格朗日網(wǎng)格上計算, 然后再通過離散的Delta 函數(shù)“傳播”(spread)到背景歐拉網(wǎng)格上. 根據(jù)IB 力的具體計算方法, 這類方法包括Peskin 最早提出浸沒邊界法時的經(jīng)典方法(Peskin 2002)、基于彈簧模型的方法(Lai & Peskin 2000)、罰函數(shù)法(Huang et al. 2011)、多孔介質模型法(Angot et al. 1999)、反饋力法(Goldstein et al. 1993)、直接力法(Uhlmann 2005)和投影法(Taira & Colonius 2007)等. 離散力法則是對沒有IB 力的原始NS 方程先進行離散, 然后再根據(jù)流固界面邊界條件在背景歐拉網(wǎng)格上直接施加IB 力. 這類方法又包含直接力法(Fadlun et al. 2000)、鬼點法(Majumdar et al. 2001)和切割網(wǎng)格法(Udaykumar et al. 1996)等. 連續(xù)力法中物體邊界的識別是模糊的, 離散力法中邊界的識別是清晰的, 但由于離散力法在處理運動物體時會由于邊界位置改變導致插值過程不連續(xù),進而導致虛擬的IB 力振蕩(Mittal & Iaccarino 2005), 所以在處理顆粒兩相流問題中, 通常使用連續(xù)力法. 這里以顆粒兩相流全分辨模擬中廣泛使用的直接力法(Uhlmann 2005, Breugem 2012)為例介紹浸沒邊界法的具體實施步驟.

首先對不施加IB 力的NS 方程動量方程進行離散得到預測速度u*

隨后, 對顆粒的運動 (包括平動和轉動) 利用顆粒動力學方程進行推進

通常, 對于非球形顆粒, 轉動方程在顆粒參考系下求解, 詳見2.1 節(jié).

2.3.2 格子玻爾茲曼方法

格子玻爾茲曼方法 (lattice boltzmann method, LBM) 也是全分辨顆粒兩相流問題數(shù)值模擬的有效手段. 格子玻爾茲曼方法本質上是一種流體求解器. 但是與傳統(tǒng)的有限差分、有限體積等基于對描述宏觀流體運動的NS 方程進行離散的流體求解器不同, 格子玻爾茲曼方法是從介觀角度的氣體動理學Boltzmann 方程出發(fā)來描述流體動力學(Chen & Doolen 1998), 即

方程(47)實際描述的是流體系統(tǒng)中“流體粒子”的運動, 式中的fi是流體粒子的在第i 個方向的速度分布函數(shù),x代表“晶格” (lattice) 的位置,Ωi代表碰撞算子,i為晶格的方向指標 (通常對于二維問題可取9 個方向, 對于三維問題可取15,19 或27 個方向, 如圖6 所示) , |ei|=dx/dt. 求解方程 (47) 涉及到“對流”(streaming)和“碰撞”(collision)兩步, 常用Bhatnagar-Gross-Krook(BGK)單弛豫時間的碰撞算子(Chen et al. 1992, Qian et al. 1992)

式中, 上標“eq”代表平衡態(tài)的分布函數(shù), 具體形式可見參考文獻(Chen & Doolen 1998). 得到方程 (47) 的解后, 對于不可壓流動問題, 基于Chapman-Enskog 展開(Chen & Doolen 1998), 就可以得到流體力學中的宏觀量 (如密度、壓強、速度等) . 關于LBM 的具體求解步驟和推導過程,可以參見相關綜述文章(Chen & Doolen 1998, Aidun & Clausen 2010)或相關書籍可參見 (何雅玲等 2009) .

LBM 方法較之于傳統(tǒng)的方法有以下三個優(yōu)點: 其一, 由于其求解過程局部性強, 求解過程簡單 (無需求解微分方程) , 因此求解程序并行化容易, 并行效率高(Aidun & Clausen 2010). 其二,由于LBM 求解器是基于介觀尺度的方法, 因此它不僅能夠處理連續(xù)介質問題, 對于稀薄流動等連續(xù)介質假設失效的問題也可以進行模擬. 其三, 對于邊界條件的處理也較容易, 例如可以使用回彈方法處理無滑移壁面(Ladd 1994a, 1994b), 使用Lees-Edwards 邊界條件方法處理無界域內的剪切流動邊界條件(Lees & Edwards 1972)等. 此外, 對于復雜流動問題 (如湍流、帶界面的流動) 都可以使用基于LBM 的求解器進行處理(Aidun & Clausen 2010).

對于全分辨顆粒兩相流問題, LBM 的回彈方法就可以直接地處理顆粒-流體邊界處的無滑移條件(Ladd 1994a, 1994b), 并且可以通過計算碰撞過程中的動量交換求出顆粒受到的流體力和力矩(即FH,TH) , 從而用以推進顆粒運動的求解. 但是使用這種方法處理的顆粒邊界呈臺階狀(Tenneti & Subramaniam 2014), 并存在質量守恒無法滿足、動邊界力振蕩以及伽利略不變性無法保證的問題(Aidun & Clausen 2010). 關于這些問題的討論以及改進可以參看相關文獻(Aidun &Clausen 2010, Peng et al. 2016).

圖6格子結構與速度的示意圖. (a) D2Q9, (b) D3Q19

2.3.3 浸沒邊界與格子玻爾茲曼混合方法

除了使用回彈模式處理顆粒-流體的邊界之外, 也可以將LBM 流體求解器和浸沒邊界法相結合 (即LBM-IBM 方法) , 通過在顆粒表面處流體施加體積力的方法來滿足無滑移邊界條件.這種方法的好處在于既利用了LBM 求解器的易編程、易并行化的優(yōu)點, 同時能夠避免使用回彈方法帶來的質量不守恒以及邊界動量力振蕩的問題(Aidun & Clausen 2010). 因此這種方法也被廣泛應用于顆粒兩相流全分辨數(shù)值模擬之中(Derksen 2011, Eshghinejadfard et al. 2016).

2.3.4 顆粒碰撞

對于全分辨的顆粒兩相流模擬, 由于顆粒的尺寸不可忽略, 通常還要考慮顆粒之間的碰撞,即四向耦合. 當不考慮流體相時, 已經(jīng)有相關綜述文章詳細討論了非球形顆粒的碰撞問題 (稱為干碰撞) , 并且對其中涉及到的接觸檢測和碰撞模型等重點問題進行了闡述(Zhong et al. 2016).但是當有流體存在時, 除了碰撞力, 顆粒間距很小時還存在潤滑力作用. 即便是全分辨的數(shù)值模擬也無法分辨顆粒距離很近時的潤滑力作用, 因此有流體存在時的顆粒碰撞 (稱為濕碰撞) 更加復雜. 一般有三種思路對濕碰撞過程進行數(shù)值模擬. 最簡單的是用簡化的排斥力來模擬顆粒間的碰撞和相互作用Fip,j(Glowinski et al. 1999), 對于球形顆粒表達式為

式中?是碰撞力產生的閾值,?p ?1 決定了碰撞力的大小, 這兩個量均需要人為給定.di,j為兩顆粒中心位置的距離,Ri與Rj分別為兩顆粒半徑,Xi與Xj分別為兩顆粒的中心位置. 這種碰撞模型形式簡單, 但是只能考慮法向碰撞力, 而且碰撞過程沒有能量損失. 第二種思路是潤滑力模型加硬球碰撞模型(Derksen 2011, Jain et al. 2019). 其中的潤滑力模型是基于Jeffrey 對球形顆粒間潤滑力的理論推導(Jeffrey 1982), 在顆粒間距小于一定范圍 (通常是1 個歐拉網(wǎng)格) 時施加.而當顆粒相互接觸之后, 再通過碰撞前兩顆粒的速度、碰撞過程的動量守恒關系、和先驗給定的恢復系數(shù)確定兩顆粒碰后的速度大小. 這種方法不需要分辨顆粒接觸后的碰撞過程, 理論上可以考慮法向和切向的動量交換, 但是對于多體碰撞 (三個或三個以上顆粒間的碰撞) 難以處理.最后一種思路是潤滑力加軟球碰撞模型(Xia et al. 2020). 這種思路同樣在顆粒間距很近時首先施加潤滑力作用, 而當顆粒產生相互接觸后使用基于彈性力學理論推導得到的碰撞力 (例如經(jīng)典的Hertz 模型) 來分辨碰撞過程. 這種方法與物理真實情況最為接近, 可以考慮碰撞過程的切向動量交換, 也可以考慮多體碰撞問題, 但是會由于碰撞過程剛度系數(shù)過大影響數(shù)值穩(wěn)定性, 需要很小的時間步或用多時間步的方法來處理碰撞過程. 所以有些學者通過一種“人工延長碰撞時間”的方法減小碰撞過程的剛度, 從而改善碰撞模擬的數(shù)值穩(wěn)定性(Kempe & Froehlich 2012,Costa et al. 2015).

3 非球形顆粒在靜止流場沉降或均勻來流中轉動

非球形顆粒在均勻來流中受到的力矩與其方位角的關系, 將決定其在均勻來流中的轉動方式, 進而影響其在流體中的取向和受力. 影響非球形顆粒在均勻來流下的轉動的關鍵無量綱數(shù)是顆粒雷諾數(shù) . 本節(jié)將基于考慮不同Re的情況對該問題進行討論, 并且以單個顆粒在靜止流場中的沉降問題為例解釋該問題的物理意義.

Re=U∞D/ν

3.1 Stokes 流動情況

當忽略流體慣性效應時, 根據(jù)Jeffery 的理論(Jeffery 1922), 在均勻來流中無論顆粒的方位角如何, 都不會受到流體力矩的作用, 也就是說顆粒會保持方位角不變. 根據(jù)這樣的結論, 對于一個非球形顆粒, 以任意的初始方位角釋放, 都會保持該角度不變. 而再根據(jù)非球形顆粒的受力公式(7), 就可以推導出顆粒以不同方位角沉降時的沉降速度(Ardekani et al. 2017a)

式中eg,n分別是重力方向和顆粒的回轉軸方向單位矢量, 而v1,v3分別是回轉軸與重力垂直和回轉軸與重力平行時顆粒的沉降速度, 即對于細長顆粒v1對應最大阻力沉降速度,v3對應最小阻力沉降速度; 對于扁平顆粒則v1對應最小阻力沉降速度,v3對應最大阻力沉降速度. 也就是說,顆粒的沉降方位角、沉降速度矢量均可以不和重力方向重合, 顆粒以不同的取向沉降時也會有不同的沉降速度. Candelier 和Mehlig (2016)就利用細長體理論得到的顆粒受力表達式, 推導出一根細針 (如圖7) 受重力作用沉降時, 速度矢量方向β和細針的取向角α的關系

3.2 有限小Re 情況

圖7細針以任意取向在流體中沉降示意圖

在考慮流體慣性效應后, 結果會有顯著的不同. Khyat 和Cox(1989)基于細長體近似理論,利用小參數(shù)漸進展開方法, 推導出了細長體 (κ=a/l ?1 ) 在均勻來流中, 當κRe ?1 時, 受到的力矩與取向角的關系. 隨后Dabade 等(2015)給出了軸對稱橢球顆粒的流體慣性力矩在有限小Re下 (Re ?1 ) 的解析表達式式(23). 根據(jù)流體慣性力矩的形式與方位角的關系, 利用顆粒轉動的穩(wěn)定性分析, 可以推導放置在均勻來流中的顆粒, 在流體慣性力矩的作用下, 最終到達的平衡狀態(tài)一定是阻力最大的取向方向.

對于沉降問題來說, 上述結論的推論就是, 以任意初始方位角釋放的非球形顆粒, 即便在Re ?1的情況下, 也不會像完全忽略流體慣性效應的Stokes 情況一樣保持其方位角不變沉降,而是會在流體慣性的作用下逐漸改變姿態(tài)直至穩(wěn)定到阻力最大的取向 (即細長顆?；剞D軸和重力方向垂直, 扁平顆?；剞D軸和重力方向平行) , 并且沉降速度矢量也會和重力方向平行, 沉降速度大小即前面所述的最大阻力沉降速度v1.

3.3 有限大Re 情況

盡管力矩系數(shù)的擬合公式定量上誤差較大, 但是非球形顆粒在均勻來流中受到的力矩與其方位角的定性關系是明確的, 即隨著回轉軸和來流方向的夾角從0°增大到90°, 顆粒受到的力矩會先增大后減小. 并且如果顆?？梢宰杂赊D動, 最終穩(wěn)定的取向將是阻力最大的方向(Ardekani et al. 2016). 當然如果雷諾數(shù)足夠高, 使得顆粒后方出現(xiàn)周期性的渦脫落, 顆粒方位角將在其穩(wěn)定取向方向附近振蕩.

4 非球形顆粒在線性剪切流中的運動

4.1 顆粒中心處滑移速度為零

該部分的研究最早可以追溯到1922 年, Jeffery (1922)首先在Stokes 流中研究了橢球顆粒在簡單線性剪切作用下的轉動行為. Jeffery 認為忽略流體與顆粒的慣性時, 顆粒是處于不受外力的狀態(tài). 此時, 軸對稱橢球顆粒在剪切流中的轉動軌跡為穩(wěn)定的封閉曲線, 又稱Jeffery 軌跡. 該軌跡的形態(tài)與初始的顆粒取向有關. 如果橢球顆粒的三個主軸的長度均不相等, 橢球顆粒在線性剪切流中會呈現(xiàn)出復雜的運動狀態(tài), 而且在特定的形狀參數(shù)下會呈現(xiàn)混沌的轉動的行為(Hinch &Leal, 1979, Yarin et al. 1997), 如圖8(b)～(d). Einarsson 等(2016)通過微槽道實驗觀察到顆粒因為軸對稱的缺失引起的復雜的轉動行為. 當考慮顆粒慣性且依然忽略流體慣性作用時, Lundell 和Carlsson(2010)與Challabotla 等(2015a)將Jeffery 力矩作為顆粒的合外力矩. 他們通過求解橢球顆粒的剛體轉動的歐拉方程獲得顆粒的轉動行為, 他們發(fā)現(xiàn)桿狀顆粒的回轉軸傾向性地垂直于渦軸并在剪切平面內轉動, 而碟狀顆粒的回轉軸的方向傾向朝著渦軸方向且繞著回轉軸在剪切平面內轉動, 且顆粒的轉動形態(tài)隨顆粒慣性的不同而呈現(xiàn)不同的取向軌跡線 (Lundell &Carlsson 2010, Challabotla et al. 2015a) 如圖8(e)和圖8(f). 對于三軸不等顆粒, Lundell(2011b)發(fā)現(xiàn)顆粒慣性會讓顆粒從混沌的轉動的狀態(tài)轉變?yōu)轭w粒繞自身最短軸在剪切平面內周期轉動的狀態(tài), 如圖8(g). 為了進一步解釋顆粒慣性與形狀對顆粒行為的作用, Yarin 等(1997)最早使用了線性穩(wěn)定性分析中的Floquet 理論研究了顆粒形狀對轉動行為的影響, 并發(fā)現(xiàn)趨近桿狀的三軸不等顆粒的運動行為更容易出現(xiàn)不穩(wěn)定轉動的行為. Lundell (2011b)也通過Floquet 分析理論給出顆粒慣性對三軸不等顆粒的轉動穩(wěn)定性的影響. 而Einarsson 等(2014)對顆粒方程基于弱顆粒慣性進行攝動展開, 并利用Floquet 理論解析地得到軸對稱橢球顆粒在弱顆粒慣性下自旋與翻轉模態(tài)運動的穩(wěn)定性. Cui Z 等(2019)系統(tǒng)地研究了軸對稱與三軸不等橢球顆粒繞各個顆粒主軸方向轉動的穩(wěn)定性, 并從理論上提供了更加廣泛的顆粒慣性與形狀的參數(shù)下顆粒不同轉動模態(tài)的穩(wěn)定性. 與此同時, 當顆粒在簡單線性剪切平面中進行準二維轉動時, 不管顆粒形狀與慣性如何變化, 顆粒均是周期轉動. Lundell 和Carlsson (2010)發(fā)現(xiàn)顆粒的轉動周期在顆粒慣性的臨界值附近存在較大的階躍, 當顆粒慣性遠小于臨界值時, 顆粒的轉動周期只與形狀有關但與慣性無關, 而在大于臨界值時, 不同形狀的顆粒轉動周期均趨于相同值. 當流體的剪切形式隨時間進行周期變化, 在特定的顆粒形狀與慣性范圍內, 軸對稱橢球顆粒的轉動也會呈現(xiàn)出混沌的運動現(xiàn)象(Nilsen et al. 2013).

圖8(a)顆粒在簡單剪切流示意圖; (b)～(g)橢球顆粒在剪切流中的轉動時代表橢球顆粒方向的軌跡圖. 其中(b)～(d)無慣性橢球顆粒, (e)～(g)以最長軸定義的顆粒轉動慣性分布為 S tr =669, 646, 100 ;形 狀 參 數(shù)(b)(e) λ =5 , (c)(f) λ =1/5 , (d)(g)ka =1,kb =10-0.65,kc =0.1

如果考慮流體慣性的影響, 非球形顆粒的轉動將會變得更加復雜. Qi 和Luo (2002, 2003)較早利用格子玻爾茲曼方法研究了有限尺寸橢球顆粒在剪切流中的轉動行為, 并發(fā)現(xiàn)隨著顆粒雷諾數(shù)的增大以及顆粒形狀的變化, 顆粒呈現(xiàn)出多種轉動與取向模態(tài). 隨后, Yu 等(2007a), Huang等(2012) 以及Rosén 等 (2014) 進一步研究了軸對稱橢球顆粒在剪切流中的行為, 發(fā)現(xiàn)并總結出在不同顆粒剪切雷諾數(shù)下橢球顆粒的不同運動模態(tài) (如圖9) , 隨著顆粒慣性的增強, 桿狀與碟狀顆粒的運動模態(tài)發(fā)生變化. 表2 總結了長寬比為λ=2 和4 時桿狀顆粒隨顆粒剪切雷諾數(shù)變化時不同的運動模態(tài). 表3 總結了長寬比λ=1/2 的碟狀顆粒在不同雷諾數(shù)下的轉動模態(tài). Rosén 與其合作者(2015a, 2015b, 2016, 2017a, 2017b)進一步從定量角度深入探討有限尺寸橢球顆粒在剪切流中的復雜行為 (包括不同模態(tài)以及模態(tài)之間的轉換條件) , 并給出了定量研究的網(wǎng)格要求.他們發(fā)現(xiàn)顆粒慣性與流體慣性對顆粒轉動行為的影響呈現(xiàn)出競爭的關系, 進而導致相同形狀顆粒在不同初始位置下可能存在不同的轉動行為. 同時, Rosén (2017a)發(fā)現(xiàn)了具有較大長寬比的桿狀顆粒在流體慣性與顆粒慣性的相互作用下, 其自旋模態(tài)(回轉軸與渦軸方向平行)可能會出現(xiàn)Shilnikov 分岔甚至混沌的行為. 當顆粒為三軸不等橢球顆粒時, 在極低顆粒雷諾數(shù)下三軸不等顆?？赡軙霈F(xiàn)混沌行為(Hinch & Leal 1979, Yarin et al. 1997, Lundell 2011, Cui Z et al.2019), Rosén 等 (2017b)較系統(tǒng)地研究了流體慣性對三軸不等橢球顆粒在剪切流中的運動行為.研究表明在流體慣性作用下三軸不等橢球顆粒的行為與極低雷諾數(shù)的情況下類似, 不過流體慣性會對顆粒繞中間軸的轉動有一定穩(wěn)定作用. 在理論研究方面, Meibohm 等(2016)通過攝動理論得到了軸對稱橢球顆粒在自旋時的角速度隨顆粒剪切雷諾數(shù)的關系. 而Dabade 等(2016)與Einarsson 等(2015)利用互等原理系統(tǒng)地研究在小剪切雷諾數(shù)下流體慣性對顆粒轉動行為的影響.

圖9顆粒在剪切流中轉動的不同模態(tài)(Rosén et al. 2014)

表2 有限尺寸桿狀顆粒在剪切流中不同 R es 對應的轉動模態(tài)

表3 有限尺寸碟狀顆粒在剪切流中不同 R es 對應的轉動模態(tài)

在Jeffery (1922)推導出橢球顆粒在Stokes 流中的力矩方程的同時, 也提出一個關于顆粒最終穩(wěn)定轉動狀態(tài)應該處于系統(tǒng)最小能量耗散的猜想. Jeffery 認為桿狀顆粒的回轉軸平行于渦軸進行自旋和碟狀顆粒的回轉軸垂直于渦軸在速度梯度平面翻轉為最終的穩(wěn)定狀態(tài). 盡管Taylor(1923)通過實驗觀察顆粒轉動的軌跡驗證了Jeffery 提出的最小能量耗散的猜想, 但在近些年的直接數(shù)值模擬的結果中均不支持該猜想(Qi & Luo 2003, Rosén et al. 2014, Huang et al. 2012,Mao & Alexeev 2014). Huang 等 (2012) 的直接數(shù)值模擬的結果表明, 桿狀顆粒的回轉軸在速度梯度平面的翻轉和碟狀顆粒的回轉軸平行于渦軸進行自旋是最終穩(wěn)定的狀態(tài), 他們對應的是最大能量耗散的狀態(tài). Mao 和Alexeev (2014) 則發(fā)現(xiàn)對于長寬比為2 的桿狀顆粒, 當顆粒剪切雷諾數(shù)大于60, 顆粒的轉動會從翻轉變?yōu)樽孕B(tài). Einarsson 等 (2015) 通過攝動理論對顆粒運動方程基于小剪切雷諾數(shù)進行展開, 理論表明非常扁平的碟狀顆粒的翻轉與自旋在小剪切雷諾數(shù)下都可能是穩(wěn)定的模態(tài). 這些研究成果意味著Jeffery 提出的最小能量耗散的猜想可能存在一定的適用范圍, 比如顆粒形狀范圍與顆粒剪切雷諾數(shù)范圍等.

4.2 顆粒中心處滑移速度非零

對于大多數(shù)實際情況, 顆粒在流場中輸運過程中一定會存在平動的滑移速度. 僅具有滑移速度的顆粒問題已經(jīng)被大量的研究, 相關內容可以參見第3 節(jié). 但同時考慮剪切與滑移的研究目前還缺少系統(tǒng)研究. Li 等(2019)研究了將顆粒中心固定在僅有一個壁面滑動的Couette 槽道流動的中心位置的轉動情況. 其結果展示了橢球顆粒在平動滑移產生的流體慣性力矩與剪切產生的力矩作用下產生的轉動行為, 但尚未對其運動行為展開系統(tǒng)與細致的深入研究工作(圖10).

5 非球形顆粒在均勻各向同性湍流中的運動

均勻各向同性湍流是一種最簡單的湍流場, 其各階湍流統(tǒng)計特性在空間處處相等且不隨坐標系的剛性轉動而改變(張兆順等 2017). 為了研究非球形顆粒與湍流的相互作用, 均勻各向同性湍流是一個比較合適的選擇. 近些年, 有關非球形顆粒在均勻各向同性湍流中運動的研究比較活躍, 接下來將主要針對點顆粒與全分辨顆粒方法, 對目前的研究現(xiàn)狀進行總結.

5.1 點顆粒橢球顆粒

5.1.1 無慣性橢球顆粒

圖10顆粒中心處滑移速度非零時的剪切流模型示意圖

對于顆粒尺寸遠小于湍流的Kolmogorov 尺寸, 可忽略流體慣性以及顆粒慣性的作用. 此時顆粒軌跡可近似跟隨流體質點, 顆粒的轉動可以根據(jù)當?shù)氐牧黧w速度梯度和顆粒取向信息通過Jeffery 方程求得, 而顆粒取向需要沿顆粒軌跡進行積分求得. 無慣性橢球顆粒因為其模型簡單并且具有很好的跟隨性而被廣泛應用在非球形顆粒湍流的理論研究中. 盡管實際問題中顆粒慣性和流體慣性并不能完美與無慣性顆粒的假設符合, 但大量的實驗觀測結果表明, 無慣性橢球顆粒模型能較好預測微小顆粒(流體慣性和顆粒慣性影響非常小)的復雜行為 (Parsa et al. 2012,Marcus et al. 2014, Fries et al. 2017, Einarsson et al. 2016). Pumir 和Wilkinson (2011)和Ni 等(2014)報道了桿狀顆粒的取向與流體的渦量方向具有強相關. 同時, 他們還發(fā)現(xiàn)桿狀顆粒沿變形率張量的最大拉伸方向相關較弱, 但與變形率張量的第二特征方向相關反而較強. 但是顆粒取向與渦量的相關性主要在顆粒尺寸小于耗散尺度時較強. 如果顆粒尺寸超過流體的耗散尺度, 且在流體的慣性尺度的范疇, 那么桿狀顆粒的取向是與變形率張量的最大拉伸方向相關最強(Pujara et al. 2019). 而對于碟狀顆粒, 其回轉軸傾向垂直于渦量方向, 且與變形率的最大壓縮方向具有較強相關 (Gustavsson et al. 2014, Byron et al. 2015) . Ni 等(2014)發(fā)現(xiàn)桿狀顆粒的取向與流體拉格朗日拉伸方向具有非常強的相關性. 他們認為長寬比大于10 的桿狀顆粒與流體拉格朗日拉伸方向具有非常好的一致性. 與此同時, Zhao 等(2019b)發(fā)現(xiàn)顆粒的取向在空間的分布并不連續(xù), 這意味著兩個距離很近(小于Kolmogorov 尺寸)的顆粒依然存在大的夾角. 在顆粒轉動方面, Parsa 等(2012)、Marcus 等(2014)和Byron 等(2015)通過實驗和數(shù)值模擬得到了均勻各向同性湍流中顆粒翻轉率(tumbling)與自旋率(spinning)隨顆粒形狀變化的函數(shù)關系 (如圖11(a)所示) . 該圖表明碟狀顆粒的翻轉率強于桿狀顆粒, 但碟狀顆粒的自旋率弱于桿狀顆粒. 而當顆粒長寬比大于3 和小于1/3 時顆粒轉動變化不明顯. 然而, 在關于顆粒轉動的前期的統(tǒng)計理論模型的研究中, 顆粒的轉動行為對于碟狀與桿狀顆粒沒有區(qū)別, 其轉動率隨顆粒形狀的變化應該是基于λ=1 對稱的(Shin & Koch 2005, Byron et al. 2015). 在該模型中, 湍流脈動隨機且顆粒取向隨機, 顆粒取向與流動相互獨立. 而Parsa 等(2012)的工作從實驗與數(shù)值上均說明了顆粒取向與流動結構的相關的重要性. Gustavsson 等 (2014) 從理論上解釋在傳統(tǒng)基于隨機湍流的模型預測的顆粒行為與真實湍流中顆粒轉動行為差異的原因. 他們認為碟狀顆粒的持續(xù)翻轉與顆粒軌跡經(jīng)過高渦量區(qū)有關, 而傳統(tǒng)的模型還不能捕捉. Pujara 等(2021)研究不同形狀的顆粒的轉動隨著對直接數(shù)值模擬的湍流場過濾尺寸變化的關系. 研究結果表明無論用Kolmogorov 尺寸還是用積分尺度過濾流場, 顆粒的轉動行為與顆粒形狀的關系與過濾尺寸影響不大, 這表明非球形顆粒的轉動可能體現(xiàn)了一種拉格朗日尺度的不變性. Chen 等(2016)發(fā)現(xiàn)目前的大渦模擬方法中流體相的亞格子模式對非球形顆粒的轉動行為的影響較大, 使得其結構與直接數(shù)值模擬的結果具有較大偏差. 近日, Pujara 等(2021)提出的不同尺寸過濾后的拉格朗日尺度不變性可能會對非球形顆粒的亞格子模型的建立具有參考意義. 如果顆粒尺寸大于Kolmogorov 尺寸, Pujara 等(2019)發(fā)現(xiàn)桿狀顆粒的翻轉率與顆粒自身的長度滿足-4/3 的冪律關系. 對于三軸不等顆粒, Chevillard 和Meneveau (2013)研究了三軸不等顆粒在各向同性湍流中的取向行為并與隨機模型進行對比. Pujara & Variano (2017)研究了三軸不等橢球顆粒在各向同性湍流中的轉動行為. 研究結果表明顆粒的擬渦能主要被分配給最長軸, 并且顆粒繞著其最長軸的轉動是最持久的. 同時, 研究發(fā)現(xiàn)湍流中渦量-變形率的相關使得顆粒的總擬渦能與顆粒形狀的關系較弱.

圖11無慣性橢球顆粒的(a)轉動率隨形狀變化的關系(數(shù)值模擬數(shù)據(jù)來自Byron 等 (2015), 兩處實驗數(shù)據(jù)分別來自Marcus 等 (2014)與Parsa 等(2012))與(b)在湍流中運動示意圖

5.1.2 慣性橢球顆粒

當考慮顆粒慣性的作用時, 近期的實驗發(fā)現(xiàn), 纖維自身的慣性會降低纖維自身的轉動(Bounoua et al. 2018). 然而, 非球形慣性顆粒在各向同性湍流的數(shù)值模擬研究卻比較少. Roy 等(2018)僅研究了顆粒平動慣性對顆粒運動的影響, 卻忽略了顆粒轉動慣性的作用. 他們認為在各向同性湍流中非球形顆粒的轉動弛豫時間要明顯小于平動的弛豫時間. Zhao 等(2015)在槽道湍流的中部(流動趨近于各向同性湍流)研究了慣性非球形顆粒的轉動行為隨顆粒慣性的變化. 其結果也展示了顆粒慣性會削弱顆粒自身的轉動行為. Gustavsson 等(2014)則是通過建立隨機湍流模型研究了慣性非球形顆粒的轉動行為. 關于顆粒慣性對顆粒在各向同性湍流中的取向與轉動行為的影響依然缺少系統(tǒng)與深入的研究.

5.1.3 橢球顆粒沉降

在早期關于球形顆粒的沉降問題中, 湍流會加速顆粒的沉降(Maxey 1987, Wang & Maxey 1993, Bec et al. 2014). 這是因為重顆粒會傾向性地聚集在流動向下的區(qū)域(Wang & Maxey 1993, Bec et al. 2014). 對于非球形顆粒而言, 顆粒沉降的過程比球形顆粒的沉降更加復雜. 一方面是形狀的各向異性導致沉降速度不僅與湍流強度有關還與顆粒取向有關(Siewert et al.2014a); 另一方面, 非球形顆粒的碰撞頻率要比球形顆粒更加強烈(Siewert et al. 2014b). 為忽略顆粒慣性的影響, Ardekani 等(2017a)直接在無慣性顆粒的模型上添加了重力沉降速度. 其研究結果表明顆粒形狀的增加會影響顆粒聚集以及平均的沉降速度. 其背后可能的原因有兩個: 一個是形狀對顆粒轉動的影響, 另一個是提高了阻力的各向異性, 而前者的作用更加明顯. Jucha 等(2018) 通過研究碟狀顆粒的沉降過程時發(fā)現(xiàn)在湍流能量耗散率較小時顆粒的碰撞主要由相鄰顆粒的取向差異引起的重力沉降導致. 隨著湍流能量耗散率的增加, 湍流的脈動主要誘導了顆粒的碰撞行為. Gustavsson 等(2017)利用簡單的高斯隨機模型就準確預測了非球形重顆粒在湍流中的取向行為.

隨著Dabade 等(2015)推導出軸對稱橢球顆粒流體慣性力矩的解析表達式, Gustavsson 等(2019)、Sheikh 等(2020) 和Anand 等(2020)意識到流體慣性力矩在非球形顆粒在湍流中沉降的重要性. 流體慣性力矩主要由顆粒與流場的滑移速度誘導, Gustavsson 等(2019)發(fā)現(xiàn)當Rep ?1時, 流體慣性力矩的作用依然不可忽略. Sheikh 等(2020)引入無量綱參數(shù) ? 衡量流體慣性力矩與黏性力矩. 當 ??1 時, 顆粒傾向朝著阻力最大的方向進行沉降.

5.2 有限尺寸顆粒

而在各向同性湍流中進行全分辨顆粒的模擬的研究到目前為止依然比較少. Schneiders 等(2017, 2019)在各向湍流中研究了尺寸接近流體Kolmogorov 尺寸的非球形顆粒運動行為以及顆粒與流體之間的相互作用, 并與傳統(tǒng)的拉格朗日雙向耦合點顆粒模型進行比較(圖12). 結果表明非球形點顆粒模型在預測取向行為上的準確性, 但在顆粒速度、顆粒與流場的相互作用等方面依然還存在明顯差距. 如何基于全分辨顆粒的研究結果去改進拉格朗日點顆粒模型在未來依然是一個值得研究的方向.

6 非球形顆粒在壁湍流中的運動

在實際工程問題中, 湍流具有各向異性且在空間上分布不均勻. 壁湍流就是一種典型的強剪切和強各向異性的湍流模型. 由于壁面的存在, 在壁面附近不僅存在強的速度剪切, 而且還存在非常豐富的湍流結構, 比如流向渦、條帶、猝發(fā)等(許春曉 2015). 隨著雷諾數(shù)的升高, 壁湍流的外區(qū)存在大尺度結構, 這些大尺度結構的流向長度通常為2～3 倍、13～15 倍甚至20 倍外區(qū)尺度(許春曉 2015). Zhang 等(2001)最早在壁湍流中進行非球形顆粒兩相流的直接數(shù)值模擬研究, 隨后研究學者們陸續(xù)在壁湍流中展開了全方面的非球形顆粒兩相流研究. 目前對非球形顆粒在壁湍流的運動行為特征有了基本的認識. 本節(jié)圍繞近些年關于非球形顆粒在壁湍流運動的研究進展展開討論.

6.1 點顆粒橢球顆粒

6.1.1 無慣性橢球顆粒

圖12有限尺寸桿狀顆粒在湍流中分布示意圖(a)與(b)自適應網(wǎng)格加密(Schneiders et al. 2019)

由于壁面的存在, 壁面附近的流動呈現(xiàn)出明顯的各向異性與不均勻性的特點. Challabotla等(2015b) 最早研究了無慣性橢球顆粒在槽道湍流中的取向行為. 在近壁面附近, 桿狀顆粒的回轉軸方向傾向朝著流向方向, 而碟狀顆粒的回轉軸傾向指向壁面的法向方向(如圖13(a)(b)).但在遠離壁面的槽道中部附近, 顆粒的取向行為趨于同顆粒在各向同性湍流中的行為類似, 這是因為在槽道中部的湍流具有局部各向同性的特性(Andersson et al. 2015). 同時, 近壁面附近存在復雜的湍流相干結構, Cui Z 等(2021)首次觀察到非球形顆粒取向與流向渦結構、上拋和下掃事件以及速度條帶之間的相關, 并提出顆粒的取向行為在近壁面附近可以被定性地劃分為剪切主導區(qū)、結構主導區(qū)以及各向同性區(qū)(如圖14). 在剪切主導區(qū), 顆粒在強剪切的作用下其旋轉會受到抑制, 桿狀顆粒與碟狀顆粒會在特定的朝向停留較長時間 (Challabotla et al. 2015b, Zhao et al. 2015) . 逐漸遠離壁面, 流場的剪切減弱, 湍流結構增加, 進入結構主導區(qū)域. 此時, 桿狀顆粒與碟狀顆粒均會受到結構的影響, 進而使得顆粒在流向渦周圍、上拋與下掃事件中以及高低速條帶中的取向行為產生明顯差異(Cui Z et al. 2021). 該類似現(xiàn)象在早期的實驗中也得到了觀察, 即桿狀顆粒在高速條帶中取向相對隨機而在低速條帶取向更趨近于流向方向(Abbasi Hoseini et al. 2015). 但是該工作將此現(xiàn)象簡單地歸結為顆粒與壁面碰撞的有限尺寸的影響, 并未深入分析和探討背后的原因. 而Cui Z 等(2021)的工作已經(jīng)消除了潛在的顆粒與壁面之間由于顆粒有限尺寸導致的碰撞效應的影響. 當顆粒位置進一步遠離壁面, 剪切趨近于零, 而湍流也趨于各向同性, 此時顆粒的轉動與取向行為逐漸與其在各向同性湍流中的行為類似(Zhao et al.2015). 不過, 在中高等雷諾數(shù)槽道湍流中, 由于流場中存在大尺度結構(比如等動量區(qū)與非等動量區(qū)), Jie 等(2019)發(fā)現(xiàn)等動量區(qū)的內外顆粒的轉動統(tǒng)計行為存在明顯差異, 圖15 展示了Reτ=1000槽道湍流中桿狀顆粒的取向分布瞬時圖.

圖13橢球顆粒在槽道湍流中的取向分布. (a)(b)無慣性橢球顆粒(Challabotla et al. 2015b); (c)(d)St =30 橢球顆粒

同時, 由于壁面剪切的存在, 桿狀顆粒與碟狀顆粒在近壁面附近的取向均與流體的渦量方向接近垂直 (Zhao et al. 2015) . Zhao 和Andersson (2016)發(fā)現(xiàn)桿狀顆粒與碟狀顆粒分別與流體拉格朗日的拉伸與壓縮方向有關, 但相關程度并沒有其在各向同性湍流中強. Cui Z 等(2020)進一步研究了桿狀顆粒的取向與流體拉格朗日拉伸方向的相關, 并發(fā)現(xiàn)Ni 等(2014)提出的長寬比大于10 的桿狀顆粒與流體拉格朗日拉伸方向一致的結論在槽道中并不是處處成立. Cui Z 等(2020)發(fā)現(xiàn)長寬比大于10 的桿狀顆粒行為與流體拉格朗日拉伸方向在黏性底層存在明顯差異.以圖16 為例, 其展示了長寬比為23.3 的桿狀顆粒與流體拉格朗日拉伸方向沿同一條顆粒軌跡線的取向行為. 圖16 的結果顯示流體拉格朗日拉伸方向在靠近上壁面時主要朝著流向方向而桿狀顆粒在壁面附近持續(xù)翻轉. Cui Z 等(2020)認為這一現(xiàn)象主要與近壁面的速度剪切和速度梯度脈動有關. 當考慮三軸不等顆粒時, Challabotla 等(2016a)發(fā)現(xiàn)顆粒同時具有碟狀與桿狀顆粒的特性, 即顆粒像桿狀顆粒一樣繞著最長軸轉動也像碟狀顆粒一樣繞著最短軸轉動.

如果壁面剪切趨于零, Yang 等(2018)通過構造Couette-Poiseuille 流動發(fā)現(xiàn)無慣性的桿狀顆粒在零剪切壁面附近取向趨向隨機而碟狀顆粒依然指向壁面法向. 同時, 他們的研究結果還表明顆粒的各向異性的轉動模態(tài)主要與顆粒取向的各向異性有關, 而壁面平均剪切對轉動模態(tài)的影響起次要作用.

6.1.2 慣性橢球顆粒

圖14顆粒在近壁流向渦結構影響下的取向行為及區(qū)域劃分(Cui Z et al. 2021). (a)(b)分別為細長桿狀與扁平顆粒在瞬時流向渦附近的取向分布; (c)(d)條件系綜平均后的顆粒取向分布, 其中(c)細長桿狀顆粒與流向夾角余弦值 | cosθx| ; (d) 扁平顆粒與展向夾角余弦值 | cosθz| ; (e)依據(jù)顆粒取向行為特點進行的區(qū)域劃分示意圖

圖15桿狀顆粒在雷諾數(shù) R eτ =1000 槽道湍流分布圖. 其中云圖為流向速度, 白色等值線代表0.95 倍槽流中心平均速度

圖16無慣性桿狀顆粒在壁面取向行為與流體拉格朗日拉伸方向的差異(Cui Z et al. 2020)

在壁湍流中, 考慮顆粒慣性的研究要早于無慣性顆粒的研究, Zhang 等(2001)最早實現(xiàn)了壁湍流中慣性桿狀顆粒的輸運與沉積特性. 隨后, Mortensen 等(2008a, 2008b)、Marchioli 等(2010,2013, 2016)、Zhao 等(2013, 2014, 2015)、Challabotla 等(2015c)和Yuan 等(2018)的研究表明慣性橢球顆粒傾向性地聚集在低速度條帶, 且與顆粒形狀的影響較小, 這與慣性球形顆粒的輸運特性類似. 然而顆粒的轉動與取向行為依賴顆粒形狀. 對于桿狀顆粒, 隨著顆粒長寬比的增大,桿狀顆粒在近壁面傾向朝著流向方向, 但隨著顆粒慣性的增大, 傾向性取向行為減弱, 進而桿狀顆粒傾向在強剪切平面內翻滾地向下游輸運. 對于碟狀顆粒, 隨著顆粒非球形度的增大, 弱慣性的碟狀顆粒在近壁面傾向朝著壁面法向, 但隨著顆粒慣性的增大, 碟狀顆粒像車輪一樣在壁面滾動著向下游輸運. Zhao 等(2014)與Marchioli 等(2016)分別研究了桿狀顆粒與流場之間平動滑移速度與轉動滑移角速度隨顆粒慣性與形狀的影響. 其結果表明隨著慣性逐漸趨近于零, 顆粒與流場之間的轉動滑移角速度并不像平動滑移速度那樣趨近于零, 依然存在明顯的滑移角速度.Zhao 等(2015)對槽道湍流中慣性橢球顆粒的轉動進行分析, 發(fā)現(xiàn)壁湍流中的近壁轉動模態(tài)與遠離壁面的轉動模態(tài)是完全不同的. 在遠離壁面的轉動模態(tài), 碟狀顆粒傾向繞非回轉軸垂直于渦量方向進行翻轉, 而桿狀顆粒繞回轉軸沿渦量方向自旋. 隨著顆粒慣性的增大, 這種趨勢被削弱.然而在近壁區(qū)域, 顆粒的慣性反而增加了碟狀顆粒繞回轉軸在剪切平面轉動而桿狀顆粒繞非回轉軸轉動. 進一步, Zhao 等(2019a)系統(tǒng)研究軸對稱橢球顆粒在槽道不同高度時的轉動模態(tài)隨剪切、湍流強度與顆粒慣性的影響, 并研究了從壁面到槽道中部兩種轉動模態(tài)的轉變. 最近,Michel 和Arcen (2021a)的研究表明非球形顆粒的統(tǒng)計量需要較長的時間才能統(tǒng)計穩(wěn)定, 同時,Marchioli 等(2016)的報道表明顆粒的轉動和取向行為比平動的統(tǒng)計更快統(tǒng)計穩(wěn)定. Michel 和Arcen (2021b)研究了不同槽道湍流的雷諾數(shù)(從Reτ=180 到550)對慣性桿狀顆粒的行為進行研究. 其結果表明隨著雷諾數(shù)的增大, 顆粒在壁面的傾向性取向與轉動行為有輕微的減弱.

6.1.3 重力沉降考慮重力沉降時, 目前大部分相關研究中主要考慮重力方向與槽道流向平行的情況, 即豎直槽道. 當槽道流動方向與重力方向一致, 稱為向下流動的豎直槽道, 當流動方向與重力相反, 則稱為向上流動的豎直槽道, 如圖17. Challabotla 等(2016b, 2016c)和Yuan 等(2017)的研究表明重力會影響到顆粒的聚集形態(tài), 使得顆粒有顯著的法向輸運. 當流動方向向上時, 非球形顆粒的分布會更加均勻, 而流動方向向下時, 顆粒聚集十分明顯, 而該影響隨著顆粒慣性的增大而減弱,且與顆粒形狀的關系不大. 除了非?？拷诿娴奈恢? 重力并不會對顆粒的取向與轉動行為有明顯的影響. Yuan 等(2018)進一步發(fā)現(xiàn)重力的作用會影響湍泳聚集(Turbophoresis)的機制且使得遠離壁面的運動顆粒比靠近壁面的運動顆粒更聚集在流體的低速條帶. 如果忽略顆粒慣性的作用, Qiu 等(2019)發(fā)現(xiàn)的非球形顆粒在向上流動的豎直槽道中傾向聚集在壁面, 而在向下流動的豎直槽道中傾向聚集在槽道中部. 這說明, 如果重力作用與基于湍泳機制聚集的作用相反, 那么顆粒的分布將會趨于均勻 (即向下流動情況) ; 如果兩種作用相同, 那么會增強顆粒在壁面的聚集(即向上流動情況). 而重力方向導致不同的法向輸運特點可能與壁面附近的流動上拋與下掃事件以及流動結構相關. 其具體的影響值得進一步深入研究.

6.2 有限尺寸顆粒

Do-Quang 等(2014) 對慣性較小的纖維使用全分辨方法研究其在槽道湍流中的行為, 并發(fā)現(xiàn)纖維會因為其與壁面的碰撞效應而傾向性的聚集在高速條帶, 且較長的纖維在槽道中部有朝著展向的趨勢, 隨著逐漸靠近壁面, 較長的纖維在逐漸在剪切平面內轉動. 而在壁面附近纖維主要朝著流向方向. Eshghinejadfard 等(2017, 2018)以及Zhu 等(2018)進一步研究了中性懸浮桿狀與碟狀顆粒在槽道中的行為以及其與流體的相互作用. 他們發(fā)現(xiàn)在緩沖區(qū)時, 顆粒的運動速度要快于流體速度, 且桿狀顆粒傾向性朝著流向而碟狀顆粒傾向性地朝著壁面法向. 而在遠離壁面的位置, 顆粒的取向又趨于各向同性, 這與點顆粒方法得到的結論類似. Eshghinejadfard 等(2019)通過改變顆粒與流體的密度比來改變顆粒的慣性, 其研究結果表明顆粒與流體間的密度比增大會影響流體與顆粒的統(tǒng)計行為. 其中顆粒的取向統(tǒng)計行為與點顆粒類似, 隨著顆粒慣性的增大(密度比增加), 桿狀顆粒與流向、碟狀顆粒與法向的取向行為被削弱. Zhu 等(2020)研究軸對稱桿狀顆粒與碟狀顆粒在豎直向上流動的槽道湍流中的沉降問題. 其結果顯示沉降的顆粒傾向于朝槽道中部遷移, 這種效應對于球形顆粒最明顯. 而非球形顆粒在近壁附近依然使其最長軸沿流向方向, 但在槽道中部由于沉降作用呈現(xiàn)出垂直流向的關系. 同時, Zhu 等(2020)還發(fā)現(xiàn)顆粒傾向性會聚集在高速條帶.

7 顆粒取向與轉動行為的力學機理

非球形顆粒盡管外形各向異性, 但是對于大部分問題中顆粒的聚集行為與顆粒形狀關系不明顯, 聚集的機理基本可以由湍流泳動的機制以及近壁相干結構進行解釋. 本文將主要討論由于形狀各向異性導致的復雜顆粒取向與轉動行為背后的機理.

圖17豎直槽道示意圖.(a)向下流動, (b)向上流動

7.1 渦量與變形率張量的影響

目前, 一些工作嘗試解釋各向同性湍流和壁湍流中的顆粒復雜行為的力學機理, 其中有關顆粒行為與渦量以及變形率張量的第二主軸的研究較為普遍(Pumir & Wilkinson 2011, Parsa et al. 2012, Byron et al. 2015). 比如, 在均勻各向同性湍流中, 無慣性的桿狀顆粒與流體的渦量方向具有較強的一致, 而碟狀顆粒垂直于渦量方向(Pumir & Wilkinson 2011, Byron et al. 2015).Pumir 和Wilkinson(2011)認為桿狀顆粒之所以與渦量方向一樣, 是因為渦量方向的演化方程與桿狀顆粒的運動方程具有類似性, 其中渦量方向的方程只比桿狀顆粒取向方程多了額外的黏性擴散項. 因為桿狀顆粒與渦量方向的關系, 在統(tǒng)計意義上桿狀顆粒的繞回轉軸轉動的速率(Spinning)要強于碟狀顆粒, 但桿狀顆粒翻轉速率(Tumbling)要弱于碟狀顆粒(Parsa et al. 2012,Marcus et al. 2014, Gustavsson et al. 2014, Byron et al. 2015). 不過, 直覺認為桿狀顆粒應該更加傾向性地朝著當?shù)氐淖冃温蕪埩康淖畲罄旆较? 然而, 學者們(Pumir & Wilkinson 2011, Gustavsson et al. 2014, Guala et al. 2005, Chevillard et al. 2013, Ni et al. 2014)卻發(fā)現(xiàn)桿狀顆粒和流體渦量并不會朝著變形率張量的最大主軸方向. 相反, 桿狀顆粒與流體渦量與變形率張量的第二特征向量相關較強. Xu 等(2011)發(fā)現(xiàn)并提出渦量會跟隨變形率張量的最大拉伸方向變化, 但是當渦量轉動到相應位置時, 原來的流體變形率張量的主軸方向已經(jīng)隨流體轉到其他方向. 在壁湍流中, 壁面附近的強剪切的使得壁面處的渦量方向一直朝著壁面展向方向, 但是無慣性桿狀顆粒與碟狀顆粒的回轉軸方向卻分別傾向性地朝著流向與壁面法向方向, 與渦量方向的夾角均接近直角, 如圖18(a)所示. 如果考慮顆粒慣性, 壁面的強渦量給顆粒提供了持續(xù)不斷的能量, 顆粒為了維持穩(wěn)定, 在自身轉動慣量的作用下, 桿狀顆粒與碟狀顆粒的轉動均傾向性地在流向與法向平面轉動(Zhao et al. 2015, 2019). 此時, 在顆粒慣性的作用下, 桿狀顆粒垂直于渦量但碟狀顆?；剞D軸卻與渦量對齊(Zhao et al. 2015, 2019), 如圖18(b). 但是, 因為在槽道中部槽道中部的湍流趨近于均勻各向同性的湍流的流動特征(Andersson et al. 2015), 顆粒的取向規(guī)律基本與各向同性湍流中的相同(Zhao et al. 2015, 2019).

圖18不同形狀和顆粒慣性的非球形顆?；剞D軸方向與渦量的夾角分布(Zhao et al. 2015). (a)槽道中部, (b)近壁區(qū)

由此得知, 顆粒取向行為與渦量和變形率的統(tǒng)計關系并不具備普適性, 比如, 桿狀顆粒傾向性朝著渦量方向僅在槽道中部或者均勻各向同性湍流中成立, 但在各向異性的近壁湍流附近桿狀顆粒卻與渦量方向垂直(Zhao et al. 2015). 而且, 顆粒在二維湍流或準二維流動中均不存在上述關系, 因為在二維流動中渦量始終垂直于二維平面.

7.2 拉格朗日結構

無論是分析顆粒取向與流體渦量以及變形率張量的關系, 還是顆粒取向與標量梯度存在的聯(lián)系, 這些研究均考慮顆粒取向與歐拉觀點下的流場物理量進行的比較分析. 實際上, 對于顆粒而言, 其方程均在拉格朗日觀點下進行點顆粒追蹤求解. 對于無慣性的微小橢球顆粒, 顆粒軌跡就是流體跡線, 而顆粒的轉動可以由當?shù)氐乃俣忍荻刃畔⑴cJeffery 方程直接得到. 如果僅采用顆粒與歐拉物理量(比如渦量與變形率張量) 進行比較, 將得到的是顆粒與瞬時的歐拉流場的關系, 其并不能真實反映顆粒的拉格朗日行為特性.

如何從拉格朗日觀點研究顆粒與流體拉格朗日量進行比較？在拉格朗日觀點下, 一個初始是球形的流體微團沿著流體軌跡線會發(fā)生拉伸變形, 其中最強拉伸與壓縮方向可以被認為是流體拉格朗日的拉伸與壓縮方向(Ni et al. 2014, Zhao & Andersson 2016). Parsa 等(2011)則通過實驗數(shù)據(jù)獲得桿狀顆粒與流體拉格朗日拉伸方向之間存在明顯關系, 并展示了流體拉格朗日拉伸結構對桿狀顆粒取向行為的影響 (如圖19(a)(b)).

圖19纖維與拉格朗日拉伸結構的關系. (a)周期流動(Parsa et al. 2011), (b)非周期流動(Parsa et al.2011), (c)拉格朗日拉伸(黑色箭頭)與壓縮(紅色箭頭)與拉格朗日結構的關系(Cui Z & Zhao 2021)

流體拉格朗日拉伸與桿狀顆粒取向存在怎樣的關系呢? 如果考慮一個細長桿狀顆粒, 其可以被近似地認為是一個物質線段, 其方向與流體拉格朗日拉伸方向漸近一致(Pumir et al. 2011,Parsa et al. 2012, Batchelor 1952, Johnson et al. 2017). 在早期研究工作就發(fā)現(xiàn)物質線段方向以及流體拉格朗日的拉伸方向就與渦量以及變形率張量的第二特征值方向有較強的關系(Girimaji& Pope 1990, Guala et al. 2005). 桿狀顆粒與拉格朗日拉伸方向以及物質線段存在一些差別, 比如桿狀顆粒始終是有限長寬比, 且不會在流體作用下發(fā)生拉伸與壓縮的變形. Ni 等(2014)驗證了桿狀顆粒(長寬比大于10)的取向可以用流體拉格朗日拉伸方向近似替代. 與此同時, Zhao 和Andersson(2016)的研究結果也發(fā)現(xiàn)在槽道湍流中顆粒的朝向與Ni 等(2014)在各向同性湍流中的結論類似, 即細長桿狀顆粒的取向與流體拉格朗日最大拉伸方向有關, 而扁平碟狀顆粒取向與流體拉格朗日最大壓縮方向相關. 但是, Cui Z 等(2020)的研究表明顆粒形狀(長寬比大于10)依然會對槽道的近壁附近的桿狀顆粒取向行為產生顯著影響. 而產生這一現(xiàn)象的原因可以歸結為槽道不同位置處流場的平均剪切與速度梯度脈動強度的比值較大. 如果考慮桿狀顆粒在流體拉格朗日拉伸與壓縮方向組成的坐標系中的取向分布, Cui Z 等(2020)的工作表明無論在壁面附近還是在槽道中部不同形狀的桿狀顆粒的取向分布均存在一個平臺區(qū)和冪律接近-2 的尾部區(qū).平臺區(qū)的存在就意味著顆粒在該區(qū)間內取向分布相對均勻, 對應著顆粒取向與流體拉格朗日拉伸方向的差別. 越靠近壁面, 平臺區(qū)的范圍更大, 顆粒取向與流體拉格朗日拉伸方向的差別就更明顯.類似的平臺與冪律區(qū)的分布也意味著桿狀顆粒與流體拉格朗日拉伸方向的關系在全場是類似的. 不過, 越靠近壁面, 流場的平均剪切與速度梯度脈動的比值變大, 分布函數(shù)中的平臺區(qū)變寬, 桿狀顆粒與流體拉格朗日拉伸方向的差別就變得更加顯著. Ni 等(2014)、Zhao 和Andersson (2016)以及Cui Z 等(2020)的研究均說明了一個共同的規(guī)律, 即隨著長寬比的增加/趨于零, 桿狀/碟狀顆粒逐漸與流體拉格朗日拉伸/壓縮方向一致. 只不過在剪切湍流中, 若想讓顆粒與流體拉格朗日方向完全一致, 則需要更小的形狀差異δΛ=1-|Λ| . 有關這一部分的闡述可見崔智文和趙立豪(2021)關于該內容的詳細綜述.

隨著人們對無慣性橢球顆粒取向行為研究的深入, 流體的拉格朗日拉伸與壓縮方向與非球形顆粒取向之間的數(shù)學關系也逐漸明晰. 在Pumir 等(2011)、Ni 等(2014)以及Cui Z 等(2020,2021)的工作中均將Jeffery 方程中的Λ=1 或 - 1 去近似替代流體拉格朗日拉伸或壓縮以及細長的桿狀或扁平的碟狀顆粒. 其背后的數(shù)學關系可以由Batchelor (1952)無限小物質線與物質面的演化得到. 不過 Balkovsky 和Fouxon (1999)的工作發(fā)現(xiàn)當積分時間足夠長后左柯西格林張量的最大與最小的特征向量(即流體拉格朗日拉伸與壓縮方向)的方程形式與Jeffery 方程中令Λ=1或 - 1 一樣. Cui Z 和Zhao (2021)則在該工作的基礎上, 進一步推導并完善了左柯西格林張量特征方向的時間演化方程, 并發(fā)現(xiàn)其在數(shù)學形式上與無慣性三軸不等橢球顆粒主軸的時間演化方程的高度相似性, 見式(49)和式(50). 在公式中, 隨著積分時間的延長, 柯西格林張量方程中的參數(shù)?i(i=1,2,3) 將 由無窮趨于 1,-1 和 1 . 這恰好與長軸比中間軸無限大、中間軸比短軸無限大的三軸不等橢球顆粒的取向方程一致, 即Λ1=1,Λ2=-1,Λ3=1 . 此時, 式(49)與式(50)完全一致. 隨后Cui Z 和Zhao (2021)利用該特性結合三軸不等橢球顆粒求解過程提出了一種新穎的左柯西格林張量的方法, 詳細內容可參見(Cui Z & Zhao 2021).

左柯西格林張量三個特征向量〈e1,e2,e3〉的時間演化方程

對應的三軸不等橢球顆粒的主軸〈p1,p2,p3〉的時間演化方程

8 顆粒在湍流中的減阻控制作用

湍流的減阻控制在工業(yè)生產和交通輸運等工程領域具有重要意義. 大量實驗發(fā)現(xiàn), 湍流里放入少量的添加劑, 如柔性聚合物、表面活化劑、氣泡等, 可能產生一定的減阻效應. 其中柔性聚合物減阻明顯, 且使用較為方便, 在許多領域得到了廣泛應用. 但其不受強剪切, 不適用于剪切較強的流動中(Reddy & Singh 1985). 部分研究者通過實驗發(fā)現(xiàn)剛性纖維可經(jīng)受強剪切, 且同樣具有減阻作用(Radin et al. 1975, Gyr & Bewersdorff 1995). 在數(shù)值模擬和理論分析中, 細小的剛性纖維通常被近似為無慣性的細長橢球顆粒 (如Brenner 1974, Gillissen et al. 2008) . 無慣性橢球顆粒通過額外的應力影響附近流體, 而反饋力和力矩則可忽略(Guazzelli & Morris 2011). Den Toonder 等(1997)首先利用直接數(shù)值模擬方法研究含細長顆粒的圓管湍流, 盡管他們簡單地假定顆粒的取向沿著當?shù)氐牧黧w速度向量的方向, 但仍然發(fā)現(xiàn)細長的剛性顆粒具有湍流減阻效應,且湍流統(tǒng)計結果與柔性聚合物的實驗結果在定性上基本一致. Paschkewitz 等(2004)通過歐拉-歐拉雙向耦合方法研究了最小尺寸槽道湍流中細長顆粒的減阻作用, 探討了細長顆粒的布朗運動、體積濃度、長寬比和顆粒取向的封閉模型對減阻強弱的影響, 實現(xiàn)最大減阻率可達26%, 還提出了一種可能的湍流減阻機理. 他們認為細長顆粒在準流向渦之間的拉伸區(qū)會產生強烈的應力脈動, 從而減弱了近壁渦旋結構和改變了近壁湍流的自維持過程, 最終導致了湍流減阻現(xiàn)象,如圖20 所示. Gillissen 等(2008)進一步發(fā)現(xiàn)細長顆粒在不同流動雷諾數(shù)下都具有湍流減阻效應.在這些數(shù)值模擬中, 顆粒相是在歐拉觀點下結合封閉模型進行求解, 而非使用拉格朗日點顆粒法直接求解單個顆粒的動力學方程. Moosaie 和Manhart (2013)首次通過歐拉-拉格朗日雙向耦合方法數(shù)值模擬了無慣性細長橢球顆粒與槽道湍流的相互作用, 同樣發(fā)現(xiàn)了壁湍流減阻的現(xiàn)象和流動特征. 不過, 相比于細長顆粒, 扁平顆粒對壁湍流的影響的相關研究則非常少. 單向耦合的數(shù)值模擬結果表明, 長寬比較小的扁平顆粒也會產生強烈的顆粒應力脈動, 這意味著一定濃度的扁平顆粒或許能夠調制湍流(Wang & Zhao 2020). Wang 等(2021)最近利用雙向耦合的直接數(shù)值模擬, 研究了體積分數(shù)為0.75%的無慣性扁平顆粒和細長顆粒對槽道湍流的影響, 發(fā)現(xiàn)長寬比為100 的細長顆粒減阻14.93%, 長寬比為0.01 和0.002 的扁平顆粒分別減阻1.92%和7.15%. 他們提出了無慣性軸對稱橢球顆粒的湍流減阻機理, 認為顆粒通過在展向和壁面法向對流體做的功抑制了準流向渦結構, 從而減弱了上拋下掃事件的強度, 導致了雷諾切應力的減小, 同時顆粒切應力會部分補償雷諾切應力的減小, 兩者之和的加權積分確定了阻力系數(shù)的大小. 對于扁平顆粒而言, 雖然雷諾切應力出現(xiàn)了明顯減小, 但近壁區(qū)處較大的顆粒切應力減弱了湍流減阻的效果.

另一方面, 一些學者利用全分辨方法數(shù)值模擬含有有限尺寸顆粒的壁湍流, 發(fā)現(xiàn)有限尺寸的橢球顆粒也可能減小流動阻力. Ardekani 等(2017b)發(fā)現(xiàn)中性懸浮的扁平顆粒在近壁區(qū)轉動較慢且回轉軸傾向沿壁面法向, 這導致了流體雷諾應力和湍流產生項的減小, 從而產生湍流減阻.Eshghinejadfard 等(2018)觀察到體積分數(shù)為10%的有限尺寸的扁平顆粒使得流體平均速度增加了1.3%. 他們提出, 展向速度脈動的減小和條帶結構的展向間距增大是扁平顆粒減阻的主要原因, 而顆粒對渦結構的增強作用則抑制了減阻效果. Zhu 等(2018)數(shù)值模擬了含有限尺寸的細長和扁平顆粒的槽道湍流, 發(fā)現(xiàn)顆粒偏離球形的程度越大減阻越明顯, 他們認為近壁區(qū)較低的顆粒體積分數(shù)和較小的流體雷諾應力導致了湍流減阻.

9 總結與展望

本文主要綜述了非球形顆粒兩相流的數(shù)值模擬方法、非球形顆粒在不同類型的簡單流場與復雜湍流場中的運動行為的研究進展, 并對非球形顆粒取向與轉動行為背后的機理進行了闡述與討論. 根據(jù)目前的研究現(xiàn)狀, 未來非球形顆粒兩相流可能的研究方向包括:

圖20纖維減阻機制示意圖(Paschkewitz et al. 2004)

(1) 非球形顆粒在流體中受力模型的建立與完善. 在目前主流的非球形顆粒兩相流的點顆粒模型中, 顆粒的阻力模型、力矩模型以及顆粒對流體反饋作用模型主要來自20 世紀Jeffery,Batchelor, Brenner 以及他們合作者的相關系列工作. 近年Dabade 等(2015)的流體慣性力矩模型是對非球形顆粒理論模型的重要補充, 但該模型的正確性與普適性還有待數(shù)值模擬以及實驗的驗證. 不過與球形顆粒相比, 非球形顆粒的模型依然還不夠完善, 比如Maxey-Riley 模型中提及的壓力梯度力、附加質量力以及Basset 歷史力等模型在非球形顆粒問題中往往被人為地忽略. 對于非球形顆粒而言, 這些力的準確表達式以及其對非球形顆粒運動的具體影響尚不明確.目前, 除了從理論層面推導, 未來可主要通過全分辨數(shù)值模擬的手段研究單顆粒在簡單流場中的運動來驗證與分析相關模型的正確性, 或進行相關模型的建立與改進工作. 非球形顆粒模型的建立與完善, 一方面能幫助研究者對非球形顆粒行為進行精準捕捉, 進而改進目前的顆粒模擬方法, 從而進行大規(guī)模的非球形顆粒兩相流的數(shù)值模擬, 另一方面能幫助研究者解釋非球形顆粒復雜行為背后的機理, 為相關工業(yè)設計提供理論依據(jù).

(2) 中高雷諾數(shù)非球形顆粒湍流兩相流的研究. 目前絕大部分非球形顆粒兩相流的數(shù)值模擬主要集中在較低雷諾數(shù)的情況, 有關非球形顆粒在中高等雷諾數(shù)的情況的研究屈指可數(shù), 以非球形顆粒槽道湍流為例, 目前, 摩擦雷諾數(shù)Reτ超過1000 的直接數(shù)值模擬的工作僅Jie 等(2019)與Ouchene 等(2018)兩項. 然而在實際的科學和工業(yè)問題中, 流場往往具有較高的雷諾數(shù), 在高雷諾數(shù)下, 流動中會存在明顯的大尺度結構. 大尺度結構對非球形顆粒的作用依然有待進一步研究. 與此同時, 對于高雷諾數(shù)的流場研究已經(jīng)有比較成熟的大渦模擬的方法取代直接數(shù)值模擬.但大渦模擬是對小尺度的流動進行?；? 而非球形顆粒的運動行為又往往對小尺度的流動結構比較敏感, Chen 等(2016)的工作展示了不同流體的亞格子模式就會對非球形顆粒的轉動與取向行為產生顯著影響. 因此, 未來對非球形顆粒相關的亞格子模型以及對現(xiàn)有的湍流亞格子模型的改進依然是一個比較重要且十分具有挑戰(zhàn)性的難題. 而最近的Pujara 等(2021)的工作表明無慣性橢球顆粒的運動行為具有形狀與尺度的無關性, 這可能展示了湍流的某種拉格朗日不變量隨湍流尺度無關. 這或許可以為非球形顆粒大渦模擬的研究工作提供新的研究思路.

(3) 非球形顆粒運動行為的機理分析. 非球形顆粒在流動中的復雜運動行為一方面由顆粒自身復雜的外形、顆粒慣性、顆粒運動學模型引起, 另一方面也受到顆粒周圍的復雜的流動環(huán)境的影響. 目前有關非球形顆粒的機理研究存在兩種基本方式, 第一種就是用當?shù)氐牧鲌鲂畔⑴c顆粒運動行為建立對應關系, 比如渦量與變形率張量與顆粒行為的聯(lián)系. 該方法的最大優(yōu)點是直觀, 但不可避免地忽略了顆粒拉格朗日運動的影響. 因此, 該類方法的最大缺點在于缺少拉格朗日信息, 且利用渦量與變形率解釋顆粒取向并不具備普適性. 盡管該類方法在各向同性湍流中取得非常大的成功, 但是我們還是能看到其在具有強剪切的壁面以及準二維的流動中是基本上失效的. 而第二種方式就是利用拉格朗日拉伸與壓縮方向來解釋顆粒的取向行為, 這是從拉格朗日的觀點出發(fā)進行解釋, 且該方法包括了必要的拉格朗日信息, 因此橢球顆粒的取向與流體拉格朗日拉伸或壓縮方向表現(xiàn)出極強的相關性. 然而, 該方法通常計算復雜, 需要得到大量的流體拉格朗日軌跡線, 且一般僅對于完全跟隨流體質點的非球形顆粒有比較好的效果. 非球形顆粒運動行為的機理研究的主要目的之一就是能夠希望通過周圍的流場信息推測出非球形顆粒的可能運動行為, 從而可以實現(xiàn)顆粒行為的控制或者預測. 然而, 在實際問題中流場的情況非常復雜, 使用第二種方式分析和預測顆粒的取向行為需要花費較大的計算成本. 因此, 第一種方式將是比較經(jīng)濟且有效的. 那么如何在第一種方式的基礎上或者其他方式, 提出更加有效的分析方法將可能是未來非球形顆粒兩相流研究一個重要的方向.

(4) 復雜流場中復雜顆粒的運動行為. 本綜述主要介紹了非球形顆粒在不可壓牛頓流體中的運動行為, 且流動的幾何邊界非常簡單. 未來可以從兩個方面來拓展非球形顆粒兩相流的研究范圍. 一個方面, 可以通過改變流場屬性、幾何等特征, 比如非牛頓流、密度分層流、熱對流或者可壓縮流動, 或研究粗糙壁面湍流、射流等復雜流動問題中的顆粒行為. 另一方面, 可以通過改變顆粒模型進而研究不同非球形顆粒模型對顆粒運動行為的影響, 比如游動、顆粒形狀反對稱、質心偏置等. 另外, 有關結合機器學習的智能顆粒的研究在近些年也得到了廣泛地關注, 相關綜述可見(邱敬然和趙立豪 2021).

致謝

國家自然科學基金(11 702 158, 11 911 530 141)和清華大學國強研究院(2019GQG1012)資助項目.