基于VOF方法的可調(diào)參數(shù)對氣泡聚并過程計算精度與成本的影響

周萍，蔣怡，廖義香，李家棟

(1. 中南大學 能源科學與工程學院，湖南 長沙，410083；2. 亥姆霍茨德累斯頓-羅森多夫研究中心，德國 德累斯頓，01328)

氣泡聚并現(xiàn)象廣泛存在于自然界中，并且在石油、化工、冶金、生物制藥以及食品加工等行業(yè)中得到廣泛應(yīng)用[1]。氣泡聚并過程通過改變氣液兩相間界面面積，影響相間的傳熱和傳質(zhì)，從而影響鼓泡塔、流化床等工業(yè)生產(chǎn)裝置的產(chǎn)能與效率[2-4]。因此，對氣泡聚并過程如氣泡形狀變化、聚并時間等進行詳細了解，可以更好地理解氣泡流動行為，對工業(yè)生產(chǎn)設(shè)備的設(shè)計和操作優(yōu)化至關(guān)重要。

氣泡聚并過程從機理上可分為氣泡碰撞、薄膜排干和液膜破裂3 個階段。在氣泡碰撞階段，2個氣泡由于速度不同而互相靠近，如在同軸上升情況下，后面的氣泡受前面氣泡尾渦的影響而加速，碰撞后氣泡之間形成一層液膜；在薄膜排干階段，2個氣泡之間的表面張力和毛細管力造成的虹吸使得液膜中液體流出，液膜變薄，氣泡更加接近；在液膜破裂階段，2個足夠接近的氣泡間液膜破裂，2個氣泡聚并成一個更大的氣泡[5-8]。

對氣泡聚并過程進行研究時，重點關(guān)注的是當2 個氣泡之間的液膜厚度為1.0×10-8～1.0×10-4m時液膜的排干速率[9]，這個階段的時間間隔在毫秒甚至微秒級。這種微尺度氣液兩相流現(xiàn)象的實驗研究對實驗檢測設(shè)備的采樣頻率與分辨率要求非常高，因此，受實驗設(shè)備性能的限制，對氣泡聚并過程的實驗研究較少[8,10-11]。另一方面，隨著計算流體力學(CFD)技術(shù)的快速發(fā)展，數(shù)值模擬方法在氣泡聚并研究中的應(yīng)用日益廣泛。

氣液兩相相界面的準確描述是實現(xiàn)氣泡聚并過程數(shù)值模擬的關(guān)鍵。相界面處理的方法分為兩大類：歐拉界面捕捉(Front-Capturing)和拉格朗日界 面 追 蹤(front-tracking)[12-13]。 前 者 包 括VOF[1,14-17]、Level-Set[18-20]以及Phase Field[21]等方法，而后者則包括Cell-Type、SPH(smoothedparticle-hydrodynamics)， PIC(particle-in-cell) 等 方法。近年來，還有學者將VOF 方法和Level-Set 方法這兩種具有代表性的界面捕捉方法進行耦合，得到了CLSVOF 模型[22]。值得一提的是，格子玻爾茲曼方法(LBM)是一種將宏觀流體力學和微觀分子動力學二者關(guān)聯(lián)起來的介觀理論，當其被應(yīng)用于氣泡聚并行為研究時，在處理相界面附近的混合流體方面表現(xiàn)出明顯優(yōu)勢[23]。由于LBM 方法受到計算量的限制，實際工程中的氣泡聚并行為研究仍以宏觀參數(shù)的求解為主。而VOF 方法因其固有的質(zhì)量守恒特性、處理表面拓撲結(jié)構(gòu)改變問題的合理性等特點，得到了廣大學者的青睞[24-25]。

采用VOF 方法對氣泡聚并的非穩(wěn)態(tài)問題進行數(shù)值模擬時，常常會遇到計算量及精度之間的平衡問題，即如何設(shè)置計算流程中可調(diào)參數(shù)(網(wǎng)格尺寸、時間步長以及動量方程的迭代次數(shù)等)，可以實現(xiàn)在精度不受太大影響的前提下減少計算量。此外，在現(xiàn)有的基于VOF 方法的氣泡聚并過程研究中，有關(guān)最大庫朗數(shù)、動量方程循環(huán)次數(shù)、控制方程循環(huán)次數(shù)等可調(diào)參數(shù)的協(xié)調(diào)設(shè)置的研究較少[26-29]。

為此，本文以同軸2個氣泡的聚并過程為研究對象，采用OpenFOAM 開源軟件、VOF 方法以及PIMPLE算法對氣泡的聚并過程進行數(shù)值模擬。在網(wǎng)格尺寸滿足界面厚度精確度與網(wǎng)格獨立性要求的前提下，探索計算流程中的可調(diào)參數(shù)(最大庫朗數(shù)Comax、相方程循環(huán)次數(shù)nα和控制循環(huán)迭代次數(shù)npimple)對計算量與精度間平衡性問題的影響規(guī)律，以期為氣泡聚并過程數(shù)值模擬時可調(diào)參數(shù)的設(shè)置提供指導。

1 數(shù)學模型與求解方法

VOF 方法基于歐拉法描述輸運方程，假設(shè)多種流體(或相)互不相溶，引入相體積分數(shù)來實現(xiàn)對每一個計算單元相界面的追蹤。

1.1 數(shù)學模型

不可壓縮牛頓流體的連續(xù)性方程和動量方程分別如式(1)、(2)所示。

式中：?為哈密頓算子；U為速度矢量，m/s；ρ為密度，kg/m3；P為壓力，Pa；μ為黏度，Pa·s；t為時間，s；g為重力加速度，m/s2；F為表面張力，N/m3。

式中：σ為表面張力系數(shù)，N/m；α為相體積分數(shù)；下標l表示液體；k為界面的曲率。

式中：n為相界面單位法向量。

式中：δ為極小的非零項，以確保分母不為零，可由計算得到；N為網(wǎng)格數(shù)量；為所有網(wǎng)格的體積之和，下標i表示第i個網(wǎng)格。

在VOF 方法中，相分數(shù)αg和αl分別表示氣相與液相在網(wǎng)格單元中的體積分數(shù)。建立并求解相分數(shù)輸運方程獲得氣液兩相體積分數(shù)，可以實現(xiàn)對氣液相界面的捕捉。當系統(tǒng)中只存在氣液兩相時，僅考慮其中一相的相體積分數(shù)方程即可，另一相體積分數(shù)通過關(guān)系式αg+αl=1 求得。以液相的相體積分數(shù)輸運方程為例，其計算式為

當網(wǎng)格單元內(nèi)同時存在氣液兩相時，混合流體的密度和黏度計算式為：

1.2 迭代算法及可調(diào)參數(shù)

在瞬態(tài)氣泡聚并數(shù)值模擬時，常用的控制方程迭代算法包括PISO(pressure-implicit splitoperator，壓力隱式分裂算子)算法[30-31]和PIMPLE(結(jié)合PISO 與SIMPLE 算法[32-33]，SIMPLE 算 法為壓力耦合方程組的半隱式方法算法[34-35]。PISO 算法采用速度預(yù)測—壓力校正的計算策略，即在求解動量方程(速度預(yù)測)之后，基于壓力泊松方程對速度場進行迭代校正(壓力校正)。由于PISO 算法只進行一次動量方程求解，為避免對流項線性化處理時(面)通量信息滯后帶來的影響，通常采用較小的時間步長來減少時間步進后速度場的變化幅度。PIMPLE算法采用類似的計算策略，但在每個時間步長內(nèi)對整套控制方程進行多次迭代計算以更新通量信息，因此，該算法允許使用更大的時間步長。

圖1 所示為基于OpenFOAM 的VOF 求解器interFoam 中PIMPLE 算法的計算流程。其中，對計算精度以及計算量有顯著影響的主要可調(diào)參數(shù)包括最大庫朗數(shù)Comax、控制方程循環(huán)次數(shù)npimple以及相方程循環(huán)次數(shù)nα。

圖1 PIMPLE算法流程Fig.1 Flow chart of PIMPLE algorithm

Comax建立了時間步長與網(wǎng)格內(nèi)流體速度及網(wǎng)格尺寸的約束關(guān)系[36]：

式中：Co為庫朗數(shù)；Δt為時間步長，s；ΔV為網(wǎng)格單元體積，m3；uf為網(wǎng)格面心速度，m/s；下標f為網(wǎng)格單元面編號；Sf為網(wǎng)格單元法向面積矢量，m2。在模擬氣泡聚并時，氣液界面處的網(wǎng)格分辨率往往要求較高(即網(wǎng)格尺寸較小)，通過Comax限制Δt，有助于提高計算精度，但過小的Δt會增大計算量。

對于PIMPLE循環(huán)，首先根據(jù)前一時間步獲得的速度更新面通量，之后對相方程進行求解：將Δt劃分為nα個子時間步長Δtsub，即Δtsub=Δt/nα；在Δtsub至nαΔtsub時間范圍內(nèi)對相體積分數(shù)進行求解(以Δtsub為計算時間步長)，并更新氣、液相體積分數(shù)；對Δtsub至nαΔtsub內(nèi)網(wǎng)格單元的面質(zhì)量通量ρUf·Sf進行積分，用于后續(xù)動量方程計算?？梢钥吹?，在不影響其他方程計算時間步長的情況下，相方程在Δt范圍內(nèi)嵌套了時間步長更小的瞬態(tài)計算；nα為該瞬態(tài)計算的迭代次數(shù)，nα越大，相方程的數(shù)值穩(wěn)定性及精度越高，但計算成本也會相應(yīng)增加。

PIMPLE 算法的動量預(yù)測采用PISO 算法的策略，通過一步預(yù)測、多次校正實現(xiàn)對速度場以及壓力場的計算。npimple控制了PIMPLE 算法中整套控制方程的循環(huán)求解次數(shù)，當npimple=1 時，PIMPLE 算法等同于PISO 算法。若控制方程的循環(huán)次數(shù)達到npimple時，則PIMPLE 算法結(jié)束，計算進入下一時間步長。與nα類似，npimple的取值也需要從平衡計算精度與計算量的角度考慮。

2 物理模型與定解條件

本研究選取同軸氣泡聚并過程開展數(shù)值實驗，定量評價上述可調(diào)參數(shù)，包括Comax、nα及npimple對數(shù)值模擬的影響，同時得到較優(yōu)的可調(diào)參數(shù)組合。本研究模擬的聚并過程為圓柱體容器中軸線上2個等徑氣泡的聚并過程，由于球形氣泡及圓柱體容器均具有軸對稱性，因此，可以假設(shè)氣液兩相的流動具有軸對稱性，簡化模型為二維計算域[37]。圖2所示為同軸等徑氣泡聚并的物理模型。為了減小壁面對模擬結(jié)果的影響(壁面效應(yīng))，物理模型的寬度設(shè)置為dc=10d0，同時高度取h=20d0，以保證足夠的氣泡上升空間。在初始階段，將直徑d0=0.56 mm 的2 個靜止氣泡置于計算區(qū)域底部，兩氣泡中心的距離為1.18d0，下部氣泡中心距離底面d0。

圖2 計算區(qū)域示意圖Fig. 2 Schematic diagram of computational domain

主要邊界條件設(shè)置見表1，壁面為無滑移條件，頂部為壓力出口。重力加速度g為9.81 m/s2，與y方向相反。

表1 主要邊界條件Table 1 Main boundary conditions

表2所示為數(shù)值模擬工況的主要物性參數(shù)，系統(tǒng)處于大氣壓(101 325 Pa)條件中。在計算的初始階段，設(shè)置時間步長為2×10-5s，同時采用自適應(yīng)時間步長，根據(jù)網(wǎng)格尺寸、流體速度以及庫朗數(shù)的相關(guān)設(shè)置，在計算過程中自動調(diào)整時間步長；此外，設(shè)置殘差為1×10-8來保證控制方程的收斂性。

表2 氣液兩相主要物性參數(shù)Table 2 Main physical properties of gas and liquid

3 網(wǎng)格獨立性分析與自適應(yīng)網(wǎng)格

3.1 網(wǎng)格獨立性分析

在網(wǎng)格獨立性分析的計算工況中均采用均勻布置的結(jié)構(gòu)化網(wǎng)格(uniform mesh, UM)。針對基準工況，設(shè)置4組不同尺寸的網(wǎng)格來驗證網(wǎng)格的無關(guān)性，網(wǎng)格的尺寸及數(shù)量如表3 所示，4 組方案的分辨率均達到10-5m量級。本節(jié)模擬工況中可調(diào)參數(shù)取值分別為Comax=2、nα=1及npimple=2。

表3 均勻網(wǎng)格算例的網(wǎng)格尺寸及數(shù)量Table 3 Mesh size and number of cases with uniform grid

將體積分數(shù)云圖中2個氣泡邊界連接成為一個封閉整體的這一時刻定義為氣泡聚并時刻。圖3所示為均勻網(wǎng)格算例的氣泡聚并時刻與網(wǎng)格數(shù)量的關(guān)系。當網(wǎng)格數(shù)量從156 800 個(UM3)增加到278 258 個(UM4)時，2 個工況的氣泡形狀基本相同，且均在0.011 s 發(fā)生聚并，氣泡聚并時刻趨于穩(wěn)定。

圖3 均勻網(wǎng)格算例的氣泡聚并時刻與網(wǎng)格數(shù)量的關(guān)系Fig. 3 Relationship between bubble coalescence time of cases with uniform grid number

圖4所示為不同網(wǎng)格數(shù)量下聚并時刻前上、下2 個氣泡的上升速度隨時間的變化曲線。從圖4 可以看出：當網(wǎng)格尺寸從3.0×10-5m(UM1)減小到2.0×10-5m(UM3)，上、下2 個氣泡上升速度的變化趨勢以及大小都趨于一致；使用網(wǎng)格方案UM2時，上方氣泡在0.006 s 后的上升速度與網(wǎng)格方案UM3、UM4 存在明顯的區(qū)別。因此，綜合考慮網(wǎng)格精度與計算成本，選擇網(wǎng)格方案UM3 進行后續(xù)的數(shù)值模擬研究。

圖4 上、下2個氣泡聚并時刻前的上升速度Fig. 4 Rising velocity of the upper and lower bubbles before coalescence

3.2 自適應(yīng)網(wǎng)格

液膜排干階段是氣泡聚并過程中的關(guān)鍵階段，其發(fā)生的尺度在微米量級，對氣泡周圍網(wǎng)格的細化尺度提出了更高的要求，若采用均勻網(wǎng)格，則勢必導致計算量巨大。為了更好地捕捉這一過程且盡可能減少計算量，本節(jié)在均勻布置的最優(yōu)網(wǎng)格方案(網(wǎng)格方案UM3)基礎(chǔ)上，使用自適應(yīng)網(wǎng)格(adaptive mesh refinement, AMR)探究最佳網(wǎng)格細化次數(shù)，確定計算精度與計算資源間的平衡點。

自適應(yīng)網(wǎng)格的基本思想是根據(jù)需要動態(tài)改變計算域的網(wǎng)格結(jié)構(gòu)，在物理量變化劇烈的計算區(qū)域，采用空間尺度較小的精細網(wǎng)格予以計算；在物理量變化緩慢的區(qū)域，采用空間尺度較大的粗網(wǎng)格來計算，從而提高計算效率。針對氣泡聚并過程，在氣液界面處使用精細網(wǎng)格，氣泡內(nèi)部及氣泡周圍的液相區(qū)域使用較粗的網(wǎng)格。設(shè)計4組自適應(yīng)網(wǎng)格方案，通過多次細化最小網(wǎng)格尺寸可達微米級，同時網(wǎng)格數(shù)量達到近200萬個，具體參數(shù)見表4。

表4 自適應(yīng)網(wǎng)格算例配置參數(shù)Table 4 Configuration parameters of cases with AMR

圖5 所示為以液體體積分數(shù)αl≤0.5 為判別條件截取出的兩氣泡整體的中心位置處速度隨時間的變化關(guān)系。從圖5 可以看到，在計算初期(0～3 ms)，中心點速度為0 m/s，這是因為初期兩氣泡間存在一定距離，中心位置處于兩氣泡間的液相中，速度為0 m/s；到計算末期(17～20 ms)，中心位置速度再次降為0 m/s，這是因為氣泡在上升過程形變?yōu)椤懊睜睢?，到后期形變幅度大，中心點處于兩氣泡下方的區(qū)域中，速度仍為0 m/s。當網(wǎng)格細化次數(shù)達到3 次以上(AMR3 和AMR4)，兩氣泡中心位置上升速度在各個時刻都十分接近，而細化次數(shù)為1 時中心位置處速度存在較大的差異，工況AMR2的氣泡中心上升速度隨時間的變化趨勢大體上與AMR3、AMR4的相同。

圖5 自適應(yīng)網(wǎng)格算例兩氣泡中心位置速度Fig. 5 Velocity in the middle of two bubbles of cases with AMR

圖6所示為4種網(wǎng)格細化條件下，氣泡形狀隨時間變化的示意圖。從圖6可見，5 ms時4個工況的氣泡形態(tài)基本相同，細化次數(shù)多的工況氣液邊界相對來說更加銳利。在25 ms時，AMR1中氣泡已經(jīng)發(fā)生聚并，因而與其他工況下的氣泡形態(tài)區(qū)別較大。AMR3 和AMR4 的氣泡形狀無明顯區(qū)別，但與AMR2相比，前兩者的氣液邊界更加平滑。

圖6 自適應(yīng)網(wǎng)格細化次數(shù)對氣泡形態(tài)的影響Fig. 6 Effect of refinement times of the AMR on bubble shape

綜合考慮計算精度和計算成本，選擇3次網(wǎng)格細化操作，并在此基礎(chǔ)上開展可調(diào)參數(shù)對計算精度與計算成本影響的研究。

4 不同可調(diào)參數(shù)下聚并過程的數(shù)值模擬

在氣泡聚并過程中，氣泡間液膜厚度的變化是氣液相間相互作用的體現(xiàn)，直接反映了兩氣泡相對速度的變化。另外，液膜厚度是氣泡聚并模型中的關(guān)鍵參數(shù)[38]。因此，選擇液膜厚度作為評價可調(diào)參數(shù)影響的指標，其值取截出的2個氣泡整體的豎直中心線上液面間距離。計算效率采用數(shù)值模擬的CPU 時間進行衡量，即CPU 完成0.03 s(物理時間)的工況數(shù)值模擬計算及內(nèi)核系統(tǒng)調(diào)用所耗的時間總和。需要注意的是，本文4.1 節(jié)與4.2節(jié)的計算是在2個不同的計算平臺上完成的，計算平臺的性能參數(shù)分別如下：1) CPU 主頻3.4 GHz，內(nèi)存儲器容量為128 GB；2) CPU主頻2.5 GHz，內(nèi)存儲器容量192 GB。針對Comax、nα及npimple這3個可調(diào)參數(shù)的數(shù)值實驗方案如表5所示。

表5 可調(diào)參數(shù)數(shù)值實驗Table 5 Numerical experiment for adjustable parameters

4.1 可調(diào)參數(shù)對液膜厚度及CPU時間的影響

圖7 所示為最大庫朗數(shù)Comax對氣泡間液膜厚度以及CPU時間的影響。從圖7(a)可以看到，在計算剛開始時，液膜厚度為0.18d0，隨著計算過程進行，液膜厚度緩慢減??；隨后液膜厚度的減小速度逐漸增大，表明兩氣泡互相接近的速度也相應(yīng)變大；隨著液膜厚度進一步減小，其變化速度及氣泡接近速度逐漸放緩。對比4 個工況發(fā)現(xiàn)，Comax＞0.01 的3 組工況間液膜厚度隨時間的變化曲線區(qū)別較大，而當Comax從0.05 減小到0.01 時，2組工況的液膜厚度十分接近。另一方面，Comax與液膜變薄速度成反比，即Comax越小，相同時間內(nèi)液膜厚度越小(變薄速度大)。這主要是由于Comax較大時，時間步長也相應(yīng)較大，在數(shù)值擴散的影響下流體在網(wǎng)格單元上的移動出現(xiàn)滯后，從模擬結(jié)果上表現(xiàn)為氣泡聚并過程變慢，液膜變薄速度減小；當通過減小Comax來提高計算精度時，能有效改善流體移動的滯后現(xiàn)象。從圖7(b)可以發(fā)現(xiàn)，CPU 時間隨最大庫朗數(shù)增大而減小。綜上所述，Comax=0.05可以較好地兼顧計算精度及計算成本。

圖7 最大庫朗數(shù)對氣泡間液膜厚度及CPU時間的影響Fig. 7 Influence of the maximum Courant number on liquid film thickness and CPU time

圖8所示為不同相方程循環(huán)次數(shù)nα下氣泡間液膜厚度變化及所耗CPU 時間。觀察圖8(a)可以發(fā)現(xiàn)，當nα為1時，液膜厚度變化曲線明顯位于其他4個工況下的變化曲線的上方；當nα≥3時可以獲得穩(wěn)定的液膜厚度變化曲線。根據(jù)1.2 節(jié)的分析可知，增大nα相當于使用更小的時間步長求解相方程，有助于提高計算精度。nα越大，液膜變薄速度越大。結(jié)合圖8(b)可知，nα=3可以在保證計算精度的同時，節(jié)約計算成本。

圖8 相方程循環(huán)次數(shù)對氣泡間液膜厚度及CPU時間的影響Fig. 8 Influence of number of cycles of the phase equation on liquid film thickness and CPU time

圖9 所示為不同控制方程循環(huán)次數(shù)npimple對氣泡間液膜厚度變化及CPU 時間的影響。從圖9(a)可以看到，計算精度隨著npimple增大而提高，當npimple≥8 時，可以獲得穩(wěn)定數(shù)值解。從計算耗時來看(圖9(b))，CPU 時間與npimple成正比，因為增大npimple意味著PIMPLE 算法的外部循環(huán)次數(shù)增多，因此，消耗的CPU 時間增加。從平衡計算量與計算精度的角度出發(fā)，取npimple=8。

圖9 控制方程循環(huán)次數(shù)對氣泡間液膜厚度及CPU時間的影響Fig. 9 Influence of number of cycles of the governing equations on liquid film thickness and CPU time

4.2 較優(yōu)的可調(diào)參數(shù)組合

從4.1 節(jié)可知，3 種可調(diào)參數(shù)分別取Comax=0.05、nα=3 以及npimple=8 時，能夠在各自數(shù)值實驗組內(nèi)有效平衡計算精度與成本。本節(jié)選取計算精度最高時的3 個可調(diào)參數(shù)的最值(各參數(shù)取值范圍參考4.1節(jié))進行組合，即(Comax，nα，npimple)=(0.01，5，10)?；谕S氣泡聚并數(shù)值模擬，在求解精度方面將其與(Comax，nα，npimple)=(0.05，3，8)進行比較，數(shù)值實驗的相關(guān)設(shè)置見表6，其他設(shè)置與4.1節(jié)中的設(shè)置相同。

表6 可調(diào)參數(shù)組合數(shù)值實驗Table 6 Numerical experiment for optimal adjustable parameters

圖10、圖11 所示分別為Toptimal和Tmaxmin這2 個工況下液膜厚度隨時間的變化曲線及氣泡形狀的比較。以工況Tmaxmin為基準，液膜厚度比在1附近浮動，液膜厚度的平均相對誤差為3.21%。觀察氣泡形狀，發(fā)現(xiàn)2 組工況幾乎沒有區(qū)別。綜上所述，認為(Comax，nα，npimple)=(0.05，3，8)是針對同軸氣泡VOF數(shù)值模擬的較優(yōu)的可調(diào)參數(shù)組合。

圖10 Toptimal和Tmaxmin這2個工況中液膜厚度隨時間的變化規(guī)律Fig. 10 Variation law of liquid film thickness with time in cases of Toptimal and Tmaxmin

圖11 Toptimal和Tmaxmin這2個工況的氣泡形態(tài)比較Fig. 11 Comparison of bubble shape between two cases of Toptimal and Tmaxmin

5 結(jié)論

1)Comax與液膜變薄速度成反比，即Comax越小，數(shù)值模擬時Δt相應(yīng)減小，流體在網(wǎng)格單元上以更大速度進行輸運，導致液膜變薄速度增大，反之亦然；此外，Comax越小，計算精度越高，相應(yīng)的計算成本也會提高；基于數(shù)值實驗獲得其最優(yōu)取值為Comax=0.05。

2) 液膜減薄速度隨nα增大而增大，這主要是由于增加nα相當于使用更小的時間步長求解相方程，由此導致流體輸運速度增大，加速液膜厚度減?。挥嬎憔入Snα增大而提高，計算成本亦然；基于數(shù)值實驗獲得其最優(yōu)取值為nα=3。

3)npimple與液膜減薄速度成正比；計算精度隨著npimple增大而提高，計算成本也隨著控制方程的迭代次數(shù)增大而增加；基于數(shù)值實驗獲得其最優(yōu)取值為npimple=8。

4) 通過與計算精度最高的參數(shù)組合比較，獲得基于VOF 方法模擬氣泡聚并時的較優(yōu)參數(shù)組合為(Comax，nα，npimple)=(0.05，3，8)。