付少童，賈龍菲，王利民，*

(1.中國科學(xué)院 過程工程研究所 多相復(fù)雜系統(tǒng)國家重點實驗室，北京 100190；2.中國科學(xué)院大學(xué) 化學(xué)工程學(xué)院，北京 100049；3.沈陽化工大學(xué) 化學(xué)工程學(xué)院，沈陽 110142)

0 引言

氣固多相流系統(tǒng)廣泛存在于自然界和工業(yè)過程，自然界中如河流過障、泥石流、塵埃懸浮等，工業(yè)過程中如飛機(jī)獲得升力、橋梁振動、循環(huán)流化床反應(yīng)器等均涉及復(fù)雜的流固耦合運(yùn)動。流體與固體間相互作用為非線性的多物理現(xiàn)象[1]，體系的復(fù)雜性遠(yuǎn)超單相流問題，如何準(zhǔn)確解析移動的流固邊界是正確處理流固耦合的關(guān)鍵。

近年來，格子Boltzmann 方法[2]（lattice Boltzmann method,LBM）發(fā)展較為迅速，其屬于介于宏觀連續(xù)介質(zhì)模型和微觀分子動力學(xué)模型之間的介觀模型，物理背景清晰；相較于有限差分法、有限體積法、有限元法等常規(guī)的計算流體力學(xué)方法，LBM 具有求解簡單、容易并行等特點，受到國內(nèi)外學(xué)者的廣泛關(guān)注。目前基于LBM 研究者們構(gòu)建了多種描述移動邊界的方法，如邊界鏈法（link-bounce-back,LBB）、干顆粒耦合法（dry particle coupling method,DPC）和浸入邊界法（immersed boundary method,IBM）等。LBB 由Ladd等[3]提出，其將邊界視為階梯狀的球殼，并應(yīng)用于多顆粒沉降問題，但模擬邊界與物理邊界存在偏差。DPC 法[4]中內(nèi)部不再保留流體格點，即內(nèi)部不存在流體的慣性力，但需要頻繁生成和刪除格點，且無法保證局部質(zhì)量守恒。IBM 最早由Peskin[5]提出，固體邊界由一系列拉格朗日點描述，通過插值函數(shù)實現(xiàn)流體點與一系列固體邊界點的相互作用。

Nobel 等[6]提出了一種LBM 框架下處理運(yùn)動固體邊界的方法，即浸入運(yùn)動邊界法（immersed moving boundary,IMB），其引入格子控制體和格子固含率描述邊界，通過固含率實現(xiàn)流體到固體的平滑過渡，而權(quán)重系數(shù)作為固含率的函數(shù)用于分配流體碰撞算子和固體碰撞算子的作用比重。Cook 等[7]將IMB 與DEM 方法耦合，進(jìn)一步考慮固體間相互作用。Chen等[8]驗證該方法對運(yùn)動邊界的處理是完全有效的。Feng 等[9]在此基礎(chǔ)上引入大渦模擬（large eddy simulation,LES）對高雷諾數(shù)顆粒流體系統(tǒng)進(jìn)行模擬。Wang 等[10]將該方法與時間驅(qū)動硬球模型（timedriven hard-sphere,TDHS）相結(jié)合，成功復(fù)現(xiàn)了流化床中顆粒散式流態(tài)化和聚式流態(tài)化現(xiàn)象。此外，該方法還被應(yīng)用于處理變形凝膠輸運(yùn)[11]、黏性固體流動[12]以及水力壓裂[13]等問題。Xiong 等[14]基于該方法開展了129 024 個顆粒流動計算。Zhou 等[15]基于大規(guī)模氣固系統(tǒng)的直接模擬發(fā)現(xiàn)經(jīng)典的Wen &Yu阻力公式預(yù)測值偏高且方向存在偏差。Liu 等[16-17]對4.83 × 109個流體網(wǎng)格和115 200 個固體顆粒的周期懸浮進(jìn)行了直接數(shù)值模擬，發(fā)現(xiàn)了顆粒系統(tǒng)的統(tǒng)計量存在尺度依賴性，并提出了新的阻力關(guān)聯(lián)式修正。以上均反映了采用IMB 處理運(yùn)動固體邊界具備巨大潛力。

然而，IMB 中的權(quán)重系數(shù)作為實現(xiàn)流固耦合的重要參量，其形式依據(jù)經(jīng)驗確定[6]。文獻(xiàn)表明其在顆粒雷諾數(shù)小于5 的情況下可以保證精度[18]，但存在粘度依賴問題[19]。Wang 等[20]提出了一個廣義權(quán)重函數(shù)，其增大了各固含率下對固體權(quán)重的預(yù)測，結(jié)合雙松弛時間（two-relaxation-time,TRT）模型[21]用以修正粘度依賴；但在中等雷諾數(shù)下其和原始的權(quán)重函數(shù)一樣，存在固體受力預(yù)測值偏高的問題。

本文在原始權(quán)重函數(shù)的基礎(chǔ)上，提出了一個改進(jìn)的權(quán)重函數(shù)，通過引入零固含率處的權(quán)重因子多階導(dǎo)數(shù)為0 作為限制條件，改善中等雷諾數(shù)下固體受力的預(yù)測精度。通過靜止圓柱繞流、Taylor-Couette 流和振動圓柱繞流驗證該函數(shù)的有效性，表明改進(jìn)的權(quán)重函數(shù)可作為一種合理的浸入運(yùn)動邊界方案。

1 數(shù)值方法

1.1 格子Boltzmann 方法

標(biāo)準(zhǔn)的格子Boltzmann 方法將密度分解為多個方向離散的分布函數(shù)，并在格點網(wǎng)格內(nèi)將各分布函數(shù)碰撞遷移，通過對當(dāng)前格點上各分布函數(shù)的組合還原出密度、速度等宏觀量。其方程形式為：

其中，fi(x,t)為t時刻位于x格點的第i個分布函數(shù)；ci為格子速度，二維下一般采用9 個方向的分布函數(shù)，三維下一般使用19 或27 個方向的分布函數(shù)，如圖1 所示；?t為時間步長；ci?t為空間步長；?(f)為碰撞算子。單松弛下，方程（1）為：

圖1 D2Q9 模型Fig.1 D2Q9 model

通過對分布函數(shù)求零階矩和一階矩即可更新當(dāng)前格點下的流體速度u和流體密度ρ：

1.2 流固耦合實現(xiàn)

Nobel 等[6]提出的IMB 方法如圖2 所示，其通過引入格子控制體和格子固含率 εs統(tǒng)一描述流體、固體和流固邊界。

圖2 IMB 方法示意圖Fig.2 Schematic diagram of the IMB method

εs表示固體部分在格子控制體中的體積分率，εs=0 表示格點內(nèi)全部為流體，εs=1 表示格點內(nèi)全部為固體，固體邊界的固含率介于0 和1 之間，這樣便可實現(xiàn)對邊界的平滑描述。

進(jìn)一步地，基于非平衡態(tài)反彈的思想，在標(biāo)準(zhǔn)單松弛的LBM 中引入額外碰撞項描述固體作用：

式中，us為覆蓋在格點上的固體速度，u為格點的流體速度，?i代表與i反向；β(εs,τ)為權(quán)重函數(shù)，控制流體碰撞算子和額外碰撞項的計算比重，同樣滿足β ∈[0,1]。β為1 時完全由固體的額外碰撞項控制；β為0 時完全由流體控制，方程還原為標(biāo)準(zhǔn)的單松弛LBM 方程。β的形式較為固定，一種可直接令 β=εs，另一種應(yīng)用更廣泛的形式為[6]：

固體顆粒所受的力F和力矩T通過求和所有顆粒覆蓋固體格點的動量變化得出，分別為：

其中xj和xc分別為被顆粒覆蓋的固體格點位置和顆粒中心格子位置。

1.3 改進(jìn)的流固耦合權(quán)重函數(shù)

由上節(jié)可知，β在流固耦合中起到重要作用，而目前其形式通過經(jīng)驗確定。在實際應(yīng)用中發(fā)現(xiàn)原始的 β(εs,τ)在雷諾數(shù)提高后計算的固體受力值偏高，很有可能其對固體組分的權(quán)重做了過高預(yù)測，可對其做進(jìn)一步改進(jìn)，構(gòu)造時至少滿足以下原則：

在此基礎(chǔ)上，基于Noble 的形式引入以下假設(shè)：

即在固含率為0 處使 β關(guān)于固含率的b?1 階導(dǎo)均為0，相應(yīng)的改進(jìn)形式如下：

b=1 時方程（16）還原為原始的方程（10）形式。圖3 為 τ=0.6時不同b下權(quán)重函數(shù)隨固含率的變化。b＞1 時，方程在固含率為0 附近的過渡更加光滑，且整體上減小了原 β的預(yù)測值，b越高則修正效果越強(qiáng)。

圖3 τ=0.6時不同b 下β 隨固含率的變化Fig.3 Variation of β with εs for different b at τ =0.6

2 數(shù)值驗證

2.1 靜止圓柱繞流

為了驗證權(quán)重函數(shù)的修正效果，本文采用上述方法，基于不同b的權(quán)重函數(shù)驗證不同雷諾數(shù)下均勻來流靜止圓柱繞流問題，并與文獻(xiàn)[23-26] 的結(jié)果對比。求解區(qū)域如圖4 所示。

圖4 圓柱繞流計算域及邊界條件Fig.4 Computational domain and boundary conditions for flow around a cylinder

計算區(qū)域為矩形，長和寬分別為800 和400，流體運(yùn)動黏度 υ=1.881×10?2，密度 ρ=1，圓柱直徑D=20。計算域左側(cè)為均勻來流入口邊界，速度u=U、v=0，采用Zou &He 邊界[22]；上下兩側(cè)均為周期邊界；右側(cè)為自由出流邊界，采用Neumann 邊界，即?u/?x=0、?v/?x=0。獲得圓柱的受力后，可通過下式分別計算阻力系數(shù)和升力系數(shù)：

式中Fx和Fy分別為流體對圓柱作用力在x和y方向的分量。

通過調(diào)整來流速度U實現(xiàn)不同雷諾數(shù)的計算，結(jié)果如表1 所示。各雷諾數(shù)下原始權(quán)重函數(shù)（b=1）所計算的阻力系數(shù)和升力系數(shù)均高于文獻(xiàn)[23-26]的結(jié)果，且雷諾數(shù)越大，偏離參考值越大，即公式（10）過高估計了固體的受力。隨著b的增加，不同雷諾數(shù)下所計算的升力系數(shù)和阻力系數(shù)均開始減小，b=3、Re=100、200 時阻力系數(shù)開始與文獻(xiàn)值接近。圖5為b=3、Re=100、200 時，達(dá)到動態(tài)穩(wěn)定后圓柱的阻力系數(shù)和升力系數(shù)隨時間的演化曲線，結(jié)果均呈現(xiàn)明顯的周期性。隨著b進(jìn)一步增大，各雷諾數(shù)下阻力和升力系數(shù)將進(jìn)一步降低。以上說明在較小雷諾數(shù)下公式（10）具備一定的預(yù)測性能，隨著雷諾數(shù)提高，公式（10）低估了流固邊界中流體組分的作用，需要相應(yīng)地降低固體權(quán)重，提升固體邊界的滲透性以增強(qiáng)流體的作用強(qiáng)度。

表1 不同雷諾數(shù)下圓柱的阻力系數(shù)和升力系數(shù)Table 1 Drag and lift coefficients of the cylinder at different Reynolds numbers

圖5 b=3 時，Re=100、200 對應(yīng)阻力系數(shù)和升力系數(shù)隨時間的演化Fig.5 Time evolution of the drag and lift coefficients at Re=100,200 for b=3

圖6 為b=3 時不同固含率下公式（16）隨黏度的變化曲線，并與公式（10）對比。改進(jìn)權(quán)重函數(shù)的β值均小于原始權(quán)重函數(shù)解，固含率越高修正效果越弱；固含率較低時黏度越大，改進(jìn)的權(quán)重函數(shù)修正效果越強(qiáng)。結(jié)合表1 分析表明，中等雷諾數(shù)下增加固體邊界的滲透性有利于獲得更加準(zhǔn)確的結(jié)果。

圖6 權(quán)重函數(shù)隨黏度的變化對比Fig.6 Comparison of the improved weighting function and the original weighting functionfor varying viscosity

2.2 Taylor-Couette 流

二維Taylor-Couette 流是流體力學(xué)中少數(shù)存在解析解的流動（僅限層流時），如圖7 所示。當(dāng)內(nèi)筒以角速度 ω1旋轉(zhuǎn)，外筒以角速度 ω2旋轉(zhuǎn)時，內(nèi)外筒間的速度分布為：

圖7 Taylor-Couette 流示意圖Fig.7 Schematic diagram of the Taylor-Couette flow

式中，R2為外筒半徑，γ為內(nèi)外徑之比：γ=R1/R2，r為該點與圓心的距離。該算例為曲線運(yùn)動邊界主導(dǎo)的流動，可用于檢驗運(yùn)動邊界的性能。

采用320 × 320 的矩形網(wǎng)格，在中心點放置同心圓筒，基于b=3 的改進(jìn)浸入運(yùn)動邊界對流場進(jìn)行驗證。流體運(yùn)動黏度 υ=0.2，密度 ρ=1。外筒半徑R2=150，ω2=0，內(nèi)筒邊界線速度為0.007 5。同時，定義L2 范數(shù)誤差：

式中uLB為采用該方法數(shù)值求解的流體格點速度。

圖8 為b=3 時，改進(jìn)浸入運(yùn)動邊界計算的流場與精確解對比，不同 γ下中心區(qū)域的流場均與精確解吻合較好。γ=0.8 時內(nèi)外筒邊界附近與精確解存在偏差，可能是該情況下解析流體的格點數(shù)量較少，邊界對周邊流場產(chǎn)生了擾動。

圖8 b=3 時，基于改進(jìn)權(quán)重函數(shù)得出的不同 γ下預(yù)測速度與精確解對比Fig.8 Comparison of the predicted velocities based on the improved weighting function (b=3) and the exact solutions for differentγ

表2 為不同 γ下b=1 和b=3 時，L2 范數(shù)的誤差對比。b=3 的改進(jìn)權(quán)重函數(shù)和b=1 的原始權(quán)重函數(shù)均不同程度地減小了誤差，γ=0.1 和0.2 時的效果最為明顯，誤差分別降低了61.5% 和76.7%，整體平均誤差降低了38.5%。以上表明改進(jìn)的浸入運(yùn)動邊界相較于傳統(tǒng)邊界可以提升流場的準(zhǔn)確性。

表2 不同γ 下改進(jìn)權(quán)重函數(shù)(b=3)與原始權(quán)重函數(shù)(b=1)計算Taylor-Couette 流與精確解誤差對比Table 2 Error comparison of the Taylor-Couette flow computed with the improved (b=3) and the original (b=1)weighting functions

進(jìn)一步地，基于 γ=0.5 的情況，令外徑物理長度為2.0，采用x=2.0/32、2.0/64、2.0/96、2.0/128、2.0/160的網(wǎng)格分辨率，以內(nèi)徑作為特征長度，對Re=10 的情況進(jìn)行精度驗證。由圖9 可知，算法整體保持了一階精度，b=1時的精度為1.041，b=3時的精度為1.039，b=5 時的精度為1.032。算法接近一階精度的原因是該算例基于固體邊界驅(qū)動流場，當(dāng)固含率εs=1.0 時，方程（8）和方程（9）恢復(fù)成具有一階精度的標(biāo)準(zhǔn)反彈格式。固體內(nèi)部格點雖不納入誤差計算，但依然基于方程（8）執(zhí)行碰撞遷移步驟會對流體區(qū)域產(chǎn)生一定影響。同時b=3 的誤差曲線低于b=1 和b=5的結(jié)果，表明改進(jìn)的浸入運(yùn)動邊界在改善流場準(zhǔn)確性的同時，不會降低邊界的收斂階。當(dāng)b繼續(xù)提高時，對比b=3 和b=5 的誤差線，可發(fā)現(xiàn)其誤差會有所提高，因此可將b=3 作為一個較優(yōu)的參數(shù)。

圖9 改進(jìn)浸入運(yùn)動邊界的收斂階Fig.9 Convergence order of the improved immersed moving boundary

2.3 振動圓柱繞流

本節(jié)基于b=3 時的權(quán)重函數(shù)，分析均勻來流下圓柱沿y方向周期性正弦振動的動力學(xué)現(xiàn)象。當(dāng)圓柱振動頻率與圓柱的自然渦脫落頻率相近時，圓柱升力系數(shù)波動頻率將與振動頻率一致，稱為“鎖頻”現(xiàn)象，其結(jié)果受模擬方法影響較為明顯。為了驗證圓柱振動的鎖頻區(qū)間，選用頻率比k=fe/f0作為無量綱參數(shù)，其中fe和f0分別為圓柱振動頻率和自然渦脫落頻率，由2.1 節(jié)Re=100 算例測得f0=7.9 × 10?5，模擬的計算域和邊界與圖4 一致。Koopmann[27]通過實驗得到圓柱鎖頻區(qū)間為k=0.75～1.25，且在k=1.0 兩側(cè)對稱分布。本文對Re=100、圓柱沿y方向振幅比為A/D=0.25 的情況進(jìn)行模擬，記錄k=0.5、0.9、1.0、1.1、1.5 下圓柱的升、阻力系數(shù)以及升力系數(shù)能量譜，并與文獻(xiàn)[27-29]結(jié)果對比。

圖10 展示了動態(tài)穩(wěn)定后圓柱運(yùn)動到最下端的渦量場，圖10（b、c）中的尾渦分布較為均勻和規(guī)則。由圖10（a、d）可以看出：當(dāng)圓柱振動頻率偏離自然渦脫落頻率較遠(yuǎn)時，圓柱后方尾渦將不再對稱，振動頻率越高，脫落的渦尺寸越大。

圖10 圓柱振動到最下端時尾部渦量圖（Re=100、A/D=0.25）Fig.10 Vorticity distributions as the cylinder oscillates to the bottom at Re=100,A/D=0.25

圖11 為穩(wěn)定后，不同振動頻率下圓柱升力和阻力系數(shù)的隨時間演化曲線及升力系數(shù)的能量譜結(jié)果。由圖11（a）可知，當(dāng)k=0.5 時，CL曲線中高幅值波和低幅值波交替出現(xiàn)，即拍頻現(xiàn)象，能量譜呈現(xiàn)雙峰形態(tài)，此時升力同時由fe和f0控制，主控頻率為f0，圓柱處于鎖頻區(qū)間之外；k=0.9 和k=1.1 時，升力系數(shù)隨時間的演化曲線不再由f0控制，而是鎖定在fe附近，此時處于鎖頻區(qū)間內(nèi)。由圖11（d）可知，k=1.5 時，升力系數(shù)曲線再次出現(xiàn)拍頻現(xiàn)象，主控頻率為fe，處于鎖頻區(qū)間外。

圖11 A/D=0.25 時，不同振動頻率下圓柱阻力系數(shù)和升力系數(shù)隨時間的演化及升力系數(shù)的能量譜Fig.11 Time evolution of the drag and lift coefficients under different oscillation frequencies,and the power spectra density of the lift coefficient at A/D=0.25

圖12 為Re=100、A/D=0.25 時，阻力系數(shù)均值CDmean隨頻率比k的變化規(guī)律，所得結(jié)果與Placzek等[29]的計算結(jié)果一致：CDmean隨著k的增加先增大后減小，在鎖頻區(qū)間內(nèi)k=1.1 時取得極大值。由圖12可知，采用b=3 的改進(jìn)權(quán)重函數(shù)計算的CDmean明顯與Placzek 等[29]的結(jié)果更為接近，表明本文提出的權(quán)重函數(shù)在中等雷諾數(shù)下可以對固體的動力學(xué)信息做出修正，可將該改進(jìn)邊界算法應(yīng)用于處理復(fù)雜的流固耦合問題。

圖12 阻力系數(shù)均值CDmean 隨頻率比k 的變化規(guī)律（Re=100、A/D=0.25）Fig.12 Variation of the averaged drag coefficient CDmean with the frequency ratio k at Re=100,A/D=0.25

3 結(jié)論

本文提出了一種LBM 框架下的改進(jìn)浸入運(yùn)動邊界，通過假定零固含率處的權(quán)重因子多階導(dǎo)數(shù)為0 來減小固體權(quán)重的過高預(yù)測，且當(dāng)可調(diào)參數(shù)b=1 時權(quán)重函數(shù)還原為原始形式。由靜止圓柱繞流和Taylor-Couette 流測試，通過一定程度地降低固體組分的權(quán)重，不僅獲得更為合理的固體動力學(xué)參數(shù)，還提升了流場預(yù)測的準(zhǔn)確性，b=3 時的改進(jìn)效果最佳。將修正后的權(quán)重函數(shù)應(yīng)用于振動圓柱繞流問題，成功復(fù)現(xiàn)了圓柱振動時的鎖頻現(xiàn)象，且得到的阻力系數(shù)更準(zhǔn)確，表明該函數(shù)下的浸入運(yùn)動邊界具備普適性，對后續(xù)流固耦合問題預(yù)測的改進(jìn)具有較高價值。后續(xù)將開展修正參數(shù)的?；ぷ?，探索邊界處局部速度與修正參數(shù)的關(guān)聯(lián)，實現(xiàn)權(quán)重函數(shù)的自動調(diào)整以拓寬IMB的應(yīng)用場景。