浸入運動邊界-格子Boltzmann方法4種固含率計算方法對比研究*

夏 明，鄧柳泓，黃剛海，徐遠臻

(1.湘潭大學(xué) 巖土力學(xué)與工程安全湖南省重點實驗室，湖南 湘潭 411105；2. 中南大學(xué) 土木工程學(xué)院，湖南 長沙 410075；3. 西南交通大學(xué) 交通隧道工程教育部重點實驗室，成都 610031)

0 引言

格子Boltzmann方法(LBM)[1-2]是一種基于流體微觀模型和介觀動理論的數(shù)值方法，它視流體為離散粒子集合體，通過離散粒子的碰撞和遷移來實現(xiàn)流體流動的.LBM具有粒子方法的背景，能方便地模擬流體和固體顆粒的相互作用，流固耦合[3-4]模擬比傳統(tǒng)方法有優(yōu)勢.

關(guān)于LBM流固耦合邊界的處理，過去30多年已形成幾種方法.Ladd[5-6]在1994年提出修正的反彈格式法(LBB)用于模擬顆粒懸浮問題，當(dāng)顆粒速度較大時，LBB會產(chǎn)生數(shù)值振蕩.為克服該數(shù)值振蕩，隨后Chun等[7]提出了插值反彈邊界方法(IBB)，但該方法在固體邊界處有時不能滿足無滑移邊界條件.Feng等[8]引入浸入邊界法[9](IBM)處理流固耦合邊界，其基本思想是采用獨立的兩種網(wǎng)格，即歐拉網(wǎng)格和拉格朗日網(wǎng)格，分別模擬流場和固體邊界；IBM用動量交換法計算流固作用力，這種顯式求解流體力的做法無法完全滿足無滑移邊界條件.之后一些學(xué)者提出了流體力隱式計算方法[10]，但需要復(fù)雜矩陣運算和較高計算內(nèi)存要求，且高雷諾數(shù)的流動模擬中易出現(xiàn)流體力計算振蕩[11].Noble等[12]提出的浸入運動邊界法(IMB)克服了其他耦合方式動量不連續(xù)的問題，在較低分辨率網(wǎng)格下能夠提供具有足夠代表性的網(wǎng)格非一致邊界.IMB引入格子固含率和附加碰撞項修正LBM方程，保留碰撞算子的局部性和遷移操作的簡單性，其計算流體力更加平順，是一種非常適合流固耦合的數(shù)值方法，因而被運用到諸多現(xiàn)象的模擬[13-14].

在用IMB法實現(xiàn)LBM流固耦合時，需要引入一個新的物理量，即固含率.所謂固含率，是指被固體覆蓋的格點單元比例.固含率是IMB法修正LBM方程的重要參數(shù).在固體內(nèi)部，固含率是1；在純流體格點，固含率是0；在邊界點，固含率為固體占據(jù)格點格子控制體的面積，因此固含率在0和1之間.在處理流固耦合格點時，需要精確地算出每一個被顆粒覆蓋格點的固含率.固含率的計算方法通常有4種，即蒙特卡洛法(MCM)、單元分解法(UDM)、近似多邊形法(APM)和閉合邊界法(CBM).這些方法在理論上各有優(yōu)劣，但它們的精度和耗時的衡量標準鮮有研究.本文針對常見固含率計算方法的精度和耗時做系統(tǒng)地測量，為優(yōu)選固含率計算方法提供借鑒.

1 浸入運動邊界-格子Boltzmann模型

二維問題常采用D2Q9模型[15]進行數(shù)值模擬，不考慮外力項的LBM方程為：

(1)

式中：fa(r，t)為格點的分布函數(shù)；r為格點坐標；ea為速度的離散方向；δt為時間步長；t為當(dāng)前時間；τ為無量綱松弛時間；feq為速度空間的局部平衡態(tài)分布函數(shù).

式(1)各變量的取值如下：

(2)

式中：c=δx/δt，其中δx和δt分別表示網(wǎng)格步長和時間步長，通常x和y方向的網(wǎng)格步長相同，為便于統(tǒng)一計算，一般取δx=δy=1.

(3)

IMB法處理流固耦合邊界時，被顆粒覆蓋的格點LBM方程變?yōu)椋?/p>

(4)

(5)

對比式(1)和式(4)可以看出，用IMB法修正受顆粒影響的格點LBM方程本質(zhì)上是在方程右邊加入固含率和附加碰撞項.

2 固含率計算方法

2.1 格點固含率

二維問題的固體顆粒常用圓盤表示，因而本文以圓盤顆粒展開固含率計算方法的論述.

如圖1所示，流場被劃分成格子，格子控制體為4個格子中心圍成的區(qū)域，格點固含率為格子控制體中固體區(qū)域所占百分比.如圖2所示，流場格點分為4種，正常格點(不在固體內(nèi)部且任何方向不與固體格點連接)、固體內(nèi)部格點(在固體內(nèi)部且任何方向不與流體格點連接)、固體邊界格點(在固體內(nèi)部且某一方向與流體格點連接)和流體邊界格點(不在固體內(nèi)部且某一方向與固體格點連接).正常格點按式(1)進行碰撞遷移，無須計算固含率；固體內(nèi)部格點按式(4)和式(5)進行碰撞遷移，但其固含率為1，無須額外計算；固體邊界格點和流體邊界格點也按式(4)和式(5)進行碰撞遷移，該兩類格點固含率在0和1之間，需要專門計算.由于格子控制體的面積是1，因此就圖2中的格子控制體(a)和(b)，它們的固含率實際上就等于陰影部分的面積.

圖1 格子控制體示意圖Fig.1 Schematic diagram of lattice control body

圖2 計算域中的四種格子控制體示意圖Fig.2 Schematic diagram of four lattice control body in computational domain

2.2 固含率4種計算方法

假設(shè)邊界格點坐標為(dx，dy)，下面闡述常用的4種邊界格點固含率計算方法.

2.2.1 蒙特卡洛法

該法利用統(tǒng)計思想計算格點固含率.在格點的格子控制體內(nèi)，用隨機數(shù)生成一系列隨機點，固體區(qū)域內(nèi)部的隨機點數(shù)與隨機點總數(shù)的比值即為格點固含率.MCM計算固含率步驟如下：

1)利用線性同余法分別生成x、y兩組隨機數(shù)組合成隨機點坐標.

先在區(qū)間[0，1]，生成均勻分布的隨機數(shù)，線性同余法的遞推公式為：

(6)

式中：M=2 147 483 648；p=314 159 269；q=453 806 245；X0取以系統(tǒng)時間為種子的隨機正整數(shù)；mod(M)表示對M取余；Rn為產(chǎn)生的隨機數(shù).

為了得到很好的偽隨機數(shù)，M取值是整數(shù)型的最大值；p和q必須足夠大，且不相等.

利用下式的抽樣公式，可生成區(qū)間[a，b]均勻分布的隨機數(shù).

t=a+(b-a)R，

(7)

式中：R為區(qū)間[0，1]生成的隨機數(shù)；t為區(qū)間[a，b]生成的隨機數(shù).

x坐標對應(yīng)的區(qū)間為[dx-0.5，dx+0.5]，y坐標對應(yīng)的區(qū)間為[dy-0.5，dy+0.5].

2)判斷生成的點與固體顆粒的位置關(guān)系.

假設(shè)生成的隨機點坐標為(dxm，dym)，圓盤的圓心坐標為(dCx，dCy)，半徑為r，則位置關(guān)系的判斷式為：

(8)

3)采用式(9)計算格點固含率.

(9)

式中：iNum為位于圓盤顆粒內(nèi)部的隨機點數(shù)量；iTotal為隨機點總數(shù).

2.2.2 單元分解法

UDM與MCM類似，只是把隨機點換成子單元進行統(tǒng)計，其具體步驟如下：

1)把格子控制體劃分成邊長為1/n的n2個正方形子單元，第i列j行的子單元中心坐標為：(dCx-0.5+(i+0.5)/n，dCy-0.5+(j+0.5)/n)；

2)判斷子單元中心與固體顆粒的位置關(guān)系，判斷公式同式(8)；

3)按式(9)計算固含率.

2.2.3 近似多邊形法

如圖3所示，邊界格點有5種可能情況，情況①～情況⑤的格子控制體依次有1個、2個、3個、4個和0個頂點位于圓盤顆粒內(nèi)部.APM將位于格子控制體內(nèi)部的圓弧線視為線段，即認為格子控制體被分割為2個多邊形.此時，圖3中5種情況的固含率依次為：

(10)

圖3 APM示意圖Fig.3 Schematic diagram of APM

式中：q1和q2的含義如圖3所示，具體為被圓弧切割的格子控制體棱邊位于圓盤外部的長度.

下面講述q1的計算方法，q2的計算與q1相同.

先找到顆粒內(nèi)部的頂點和與之相關(guān)的外部點.如圖4所示，假設(shè)找到的內(nèi)部點和外部點分別為B(x2，y2)和A(x1，y1)，圓盤圓心為O(dCx，dCy)，圓盤半徑為r，點C為OA與圓弧的交點，點D為AB和圓盤邊界的交點，則q1為圖中AD的長度.由余弦定理得：

(11)

圖4 q1計算示意圖Fig.4 Schematic diagram of calculating q1

令OA=a，AB=b，OB=c，OD=r，則有

(12)

將式(12)代入式(11)，可求得兩根為：

(13)

式中：φ=a2+b2-c2.

因點D位于A和B之間，且AB=1，所以式(13)的兩個解只有滿足0≤q1≤1才是正解.

2.2.4 閉合邊界法

如圖5所示，采用CBM時，位于格子控制體內(nèi)的圓弧被視為真實圓弧，因而該法可求出準確的固含率.由此可以看出，CBM的固含率就是由APM算出的固含率加上圖5中弓形的面積，即S1+S2.

圖5 CBM示意圖Fig.5 Schematic diagram of CBM

下面講述圖5弓形面積S1的計算方法.

通過式(13)易算出圓弧與格子控制體的兩個交點坐標p1(x1，y1)和p2(x2，y2)，假設(shè)圓盤半徑為r，則三角形CiP1P2和扇形CiP1P2的面積計算式為：

(14)

弓形面積S1=S扇CiP1P2-SΔCiP1P2.

3 精度和效率測試

3.1 精度測試

3.1.1 隨機點數(shù)的確定

由固含率4種計算方法的計算步驟可看出，MCM和UDM的計算精度與隨機點數(shù)和子單元數(shù)有關(guān).因此，本文分別測試不同圓盤直徑下隨機點數(shù)和子單元數(shù)對計算精度的影響，達到優(yōu)選隨機點數(shù)和子單元數(shù)的目的.測試模型：流場格子數(shù)為2 000×2 000；測試圓盤圓心坐標為(1 000，1 000)，直徑分別為10、20、100、200、400、600、800、1 000和1 200；相應(yīng)測試的格子控制體中心坐標分別為(995，1 000)、(990，1 000)、(950，1 000)、(900，1 000)、(800，1 000)、(700，1 000)、(600，1 000)、(500，100)和(400，1 000).CBM算出的準確固含率分別為0.491 654、0.495 832、0.499 167、0.499 583、0.499 792、0.499 861、0.499 896、0.499 917和0.499 931.采用MCM和UDM分別計算上述格點固含率10次，所用隨機點數(shù)和子單元數(shù)依次為100、400、900、1 600、2 500、3 600、4 900、6 400、8 100和10 000，經(jīng)測試得：1)MCM的誤差不受直徑的變化影響，只與隨機點個數(shù)有關(guān)，10組圖的結(jié)果幾乎完全一致；如圖6所示，誤差隨隨機點數(shù)的增加而降低，但隨機點數(shù)為2 500時曲線會有轉(zhuǎn)折，當(dāng)隨機點數(shù)大于400時，誤差穩(wěn)定在1%及以下，最終誤差穩(wěn)定在0.4%～0.5%；2)UDM的誤差隨直徑的增大而降低，圖7(a)和圖7(b)給出了直徑分別為10和20時的計算誤差，其余模擬工況隨直徑的增大誤差不再變化，均在0.1%以下.UDM計算精度高，子單元數(shù)大于100可使誤差穩(wěn)定在1%以下.

圖6 MCM誤差結(jié)果圖Fig.6 The error of MCM

圖7 UDM誤差結(jié)果圖：(a)直徑為10；(b)直徑為20Fig.7 The error of UDM：(a) r=10；(b) r=20

為了探究格子控制體被圓弧切得的形狀對MCM和UDM計算精度的影響.測試了在圓盤直徑為20時，每一個邊界格點的固含率誤差.發(fā)現(xiàn)隨機點數(shù)大于1 000和子單元數(shù)大于100時，每個邊界格點的固含率計算誤差穩(wěn)定在2%左右.因此，下面測試選取的隨機點數(shù)為1 000，子單元數(shù)為100.

3.1.2 圓盤大小對精度的影響

為研究圓盤大小對4種計算方法計算固含率精度的影響，進行如下測試.測試工況與第3.1.1節(jié)一致.將各種模擬工況下各邊界點固含率的誤差值進行統(tǒng)計，采用平均值Er2表征，即：

(15)

式中：Xi(i=1，2，3，…，n)為第i個邊界點的固含率絕對誤差值(剔除絕對誤差為0的邊界點)；n為參與統(tǒng)計的邊界點數(shù).

各種測試工況的Er2值如圖8(a)所示，從中可知：1)除APM外的3種方法，Er2值不受影響；原因在于，CBM誤差總是0，MCM與UDM均采用統(tǒng)計方法求解固含率，其求解結(jié)果不受格子控制體形狀影響；2)APM的誤差隨圓盤半徑增大而減??；這是因為APM的誤差為圖5中弓形的面積S1，當(dāng)圓盤顆粒增大時，S1逐漸減小.應(yīng)該指出，即便圓盤直徑僅為10倍格子長度，APM的計算誤差也僅為0.44%，表明APM容易達到很高的計算精度.

圖8 (a)固含率計算誤差；(b)固含率計算耗時Fig.8 (a)Solid ratio calculation error；(b)The time consumption of solid ratio calculation

3.2 精度測試

3.2.1 單個圓盤模型

記錄“3.1小節(jié)”下“3.2.1小節(jié)”各種測試工況的耗時(測試計算機為Intel(R) Xeon(R) CPU E5-2640 v3×2 @2.60 GHz、64G內(nèi)存、2T硬盤)，結(jié)果如圖8(b)所示.為體現(xiàn)不同計算規(guī)模下4種計算方法的效率，將圖8(b)的橫坐標變成邊界格點數(shù).從圖8(b)可以看出：1)APM和CBM效率比MCM高約2個數(shù)量級，比UDM高約1個數(shù)量級；原因在于，APM和CBM僅采用簡單的運算來求解固含率，因而計算效率高；而MCM和UDM需耗費時間生成隨機數(shù)和子單元，導(dǎo)致效率降低；2)APM計算效率較CBM高；這是因為CBM需多計算圖5中弓形的面積S1；3)UDM的效率比MCM高約1個數(shù)量級；原因在于UDM用少量的子單元可以達到大量隨機數(shù)的計算精度，子單元數(shù)小于隨機點數(shù)，因而計算效率高.

3.2.2 多顆粒沉降流固耦合模擬

為研究固含率計算方法在LBM流固耦合計算中的作用，本節(jié)將LBM與圓盤顆粒非連續(xù)變形分析方法(DDDA)進行耦合，分別用4種方法計算每種模擬工況的固含率，統(tǒng)計固含率計算環(huán)節(jié)的耗時.DDDA的基本原理參見參考文獻[16]和[17].如圖9所示，流固耦合模型模擬的是密閉方腔的多顆粒自由沉降，方腔尺寸為0.02 m×0.02 m，每個流體格子的真實物理長度為0.000 396 5 m，對應(yīng)格子區(qū)域為512×512；顆粒直徑為0.000 5 m，LBM松弛因子為0.991 5，顆粒和流體密度分別為1.01 g/cm3和1.0 g/cm3，流體計算時步為0.000 25 s，DDDA計算時步為0.000 5 s，1步固體計算對應(yīng)2步流體計算，分別用10、20、40、60、80和100個圓盤進行顆粒自由沉降模擬.圓盤顆粒生成可采用文獻[18]提出的高效算法，該算法可在0.9 s內(nèi)生成410萬顆粒.

圖9 流固耦合測試模型Fig.9 The test model of fluid-solid coupling

圖10給出了各種流固耦合模擬工況下4種方法計算固含率的耗時.從圖可以看出：1)APM和CBM的計算效率比MCM高約2個數(shù)量級，比UDM高約1個數(shù)量級；2)CBM效率略小于APM；3)UDM計算效率比MCM高約1個數(shù)量級.總體上看，4種計算方法的效率排序為：APM>CBM>UDM>MCM.

圖10 流固耦合計算耗時Fig.10 The time consumption of fluid-solid coupling

4 結(jié)論

1)MCM和UDM取1 000和100左右的隨機點和子單元時可保證固含率計算精度，此時誤差小于1%.

2)APM容易達到較高計算精度；當(dāng)圓盤直徑大于10倍格子長度時，APM的計算誤差小于0.44%.

3)APM和CBM的計算效率比MCM高約2個數(shù)量級，比UDM高約1個數(shù)量級；CBM效率略小于APM；UDM比MCM計算效率高約1個數(shù)量級.4種計算方法的效率排序為：APM>CBM>UDM>MCM.

4)用浸入運動邊界法進行DDDA-LBM流固耦合時，格點固含率計算應(yīng)采用APM，當(dāng)計算精度要求很高時應(yīng)采用CBM.