一種改進的浸入運動邊界算法

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

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

0 引言

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

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

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

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

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

1 數值方法

1.1 格子Boltzmann 方法

標準的格子Boltzmann 方法將密度分解為多個方向離散的分布函數，并在格點網格內將各分布函數碰撞遷移，通過對當前格點上各分布函數的組合還原出密度、速度等宏觀量。其方程形式為：

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

圖1 D2Q9 模型Fig.1 D2Q9 model

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

1.2 流固耦合實現

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

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

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

進一步地，基于非平衡態反彈的思想，在標準單松弛的LBM 中引入額外碰撞項描述固體作用：

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

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

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

1.3 改進的流固耦合權重函數

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

在此基礎上，基于Noble 的形式引入以下假設：

即在固含率為0 處使 β關于固含率的b?1 階導均為0，相應的改進形式如下：

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

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

2 數值驗證

2.1 靜止圓柱繞流

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

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

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

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

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

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

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

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

圖6 權重函數隨黏度的變化對比Fig.6 Comparison of the improved weighting function and the original weighting functionfor varying viscosity

2.2 Taylor-Couette 流

二維Taylor-Couette 流是流體力學中少數存在解析解的流動（僅限層流時），如圖7 所示。當內筒以角速度 ω1旋轉，外筒以角速度 ω2旋轉時，內外筒間的速度分布為：

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

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

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

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

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

圖8 b=3 時，基于改進權重函數得出的不同 γ下預測速度與精確解對比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 范數的誤差對比。b=3 的改進權重函數和b=1 的原始權重函數均不同程度地減小了誤差，γ=0.1 和0.2 時的效果最為明顯，誤差分別降低了61.5% 和76.7%，整體平均誤差降低了38.5%。以上表明改進的浸入運動邊界相較于傳統邊界可以提升流場的準確性。

表2 不同γ 下改進權重函數(b=3)與原始權重函數(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

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

圖9 改進浸入運動邊界的收斂階Fig.9 Convergence order of the improved immersed moving boundary

2.3 振動圓柱繞流

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

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

圖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 為穩定后，不同振動頻率下圓柱升力和阻力系數的隨時間演化曲線及升力系數的能量譜結果。由圖11（a）可知，當k=0.5 時，CL曲線中高幅值波和低幅值波交替出現，即拍頻現象，能量譜呈現雙峰形態，此時升力同時由fe和f0控制，主控頻率為f0，圓柱處于鎖頻區間之外；k=0.9 和k=1.1 時，升力系數隨時間的演化曲線不再由f0控制，而是鎖定在fe附近，此時處于鎖頻區間內。由圖11（d）可知，k=1.5 時，升力系數曲線再次出現拍頻現象，主控頻率為fe，處于鎖頻區間外。

圖11 A/D=0.25 時，不同振動頻率下圓柱阻力系數和升力系數隨時間的演化及升力系數的能量譜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 時，阻力系數均值CDmean隨頻率比k的變化規律，所得結果與Placzek等[29]的計算結果一致：CDmean隨著k的增加先增大后減小，在鎖頻區間內k=1.1 時取得極大值。由圖12可知，采用b=3 的改進權重函數計算的CDmean明顯與Placzek 等[29]的結果更為接近，表明本文提出的權重函數在中等雷諾數下可以對固體的動力學信息做出修正，可將該改進邊界算法應用于處理復雜的流固耦合問題。

圖12 阻力系數均值CDmean 隨頻率比k 的變化規律（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 結論

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