夏 明，鄧柳泓，黃剛海，徐遠(yuǎn)臻

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

0 引言

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

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

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

1 浸入運(yùn)動(dòng)邊界-格子Boltzmann模型

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

(1)

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

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

(2)

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

(3)

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

(4)

(5)

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

2 固含率計(jì)算方法

2.1 格點(diǎn)固含率

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

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

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

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

2.2 固含率4種計(jì)算方法

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

2.2.1 蒙特卡洛法

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

1)利用線性同余法分別生成x、y兩組隨機(jī)數(shù)組合成隨機(jī)點(diǎn)坐標(biāo).

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

(6)

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

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

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

t=a+(b-a)R，

(7)

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

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

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

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

(8)

3)采用式(9)計(jì)算格點(diǎn)固含率.

(9)

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

2.2.2 單元分解法

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

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

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

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

2.2.3 近似多邊形法

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

(10)

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

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

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

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

(11)

圖4 q1計(jì)算示意圖Fig.4 Schematic diagram of calculating q1

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

(12)

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

(13)

式中：φ=a2+b2-c2.

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

2.2.4 閉合邊界法

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

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

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

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

(14)

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

3 精度和效率測(cè)試

3.1 精度測(cè)試

3.1.1 隨機(jī)點(diǎn)數(shù)的確定

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

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

3.1.2 圓盤大小對(duì)精度的影響

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

(15)

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

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

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

3.2 精度測(cè)試

3.2.1 單個(gè)圓盤模型

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

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

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

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

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

圖10 流固耦合計(jì)算耗時(shí)Fig.10 The time consumption of fluid-solid coupling

4 結(jié)論

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

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

3)APM和CBM的計(jì)算效率比MCM高約2個(gè)數(shù)量級(jí)，比UDM高約1個(gè)數(shù)量級(jí)；CBM效率略小于APM；UDM比MCM計(jì)算效率高約1個(gè)數(shù)量級(jí).4種計(jì)算方法的效率排序?yàn)椋篈PM>CBM>UDM>MCM.

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