紅黑著色的相場(chǎng)兩相流并行投影算法

摘要：提出一種交錯(cuò)并行的有限體積投影算法求解基于相場(chǎng)法的兩相流控制方程組。該算法的實(shí)施主要依賴于壓力泊松方程的顯式推進(jìn)設(shè)計(jì)，從而突破投影算法求解不可壓Navier-Stokes方程的效率瓶頸。同時(shí)，提出了一種交錯(cuò)掃描策略來更新節(jié)點(diǎn)上的變量，以實(shí)現(xiàn)更緊湊的時(shí)空耦合。本算法與相場(chǎng)模型相結(jié)合，能夠高效、準(zhǔn)確地捕獲相界面的動(dòng)態(tài)拓?fù)渥兓y(cè)試算例結(jié)果表明：網(wǎng)格量為131 072，采用8線程CPU并行時(shí)，新提出并行算法的效率達(dá)到串行標(biāo)準(zhǔn)投影算法的80倍以上。

關(guān)鍵詞：兩相流；相場(chǎng)法；有限體積法；投影算法；并行計(jì)算

中圖分類號(hào)：U448.213 " " " " "文獻(xiàn)標(biāo)志碼：A " " "文章編號(hào)：1000-582X（2024）06-075-11

A novel red-black coloring parallel projection algorithm for two-phase flow using the phase field method

WANG Xiaoshuang， ZHANG Liangqi， XIAO Yao， ZENG Zhong

（College of Aerospace Engineering， Chongqing University， Chongqing 400044， P. R. China）

Abstract： In this study， an innovative parallel finite volume projection algorithm with a novel red-black coloring approach is proposed to solve the two-phase flow control equations using the phase field method. This strategy relies on the explicit advancement of the pressure Poisson equation， thus effectively overcoming the efficiency limitation inherent in the traditional projection algorithms for solving the incompressible Navier-Stokes equations. Moreover， we implement an interleaved scanning strategy for updating nodal variables， which significantly enhances the spatiotemporal coupling in a compact manner. The integration of this technique with the phase field method facilitates more accurate capture of interface dynamics and topology at a lower cost. Test results show that， utilizing a grid size of 131 072 and an 8-thread CPU parallelization， the proposed parallel algorithm is more than 80 times more efficient than the serial standard projection algorithm.

Keywords： two-phase flows; phase field; finite volume method; projection algorithm; parallel computing

兩相流在自然界和工業(yè)界中廣泛存在，例如雨滴[1]、噴霧霧化[2]和微流控芯片[3]等。對(duì)兩相流界面動(dòng)力學(xué)的深入研究，不僅可以幫助人們更好地理解自然現(xiàn)象，還可以為工程技術(shù)和生物醫(yī)學(xué)領(lǐng)域的發(fā)展提供重要的理論和實(shí)踐支撐。然而，兩相流精確且高效的數(shù)值計(jì)算仍然是計(jì)算流體動(dòng)力學(xué)的重要挑戰(zhàn)。當(dāng)采用基于壓力修正的投影算法[4]求解不可壓Navier-Stokes方程時(shí)，往往將其解耦為3個(gè)方程，依次計(jì)算得到中間速度 ，當(dāng)前壓力 和當(dāng)前速度

， （1）

， （2）

， （3）

式中： 是密度；μ是動(dòng)力黏度； 是速度矢量；P是壓力； 是質(zhì)量通量，即 ； 是外力項(xiàng)（例如重力、表面張力等）。

方程（1）和方程（3）都能顯式處理，但泊松方程（2）不含壓力瞬態(tài)項(xiàng)，只能隱式計(jì)算。當(dāng)采用有限體積法進(jìn)行離散求解時(shí)，隱式格式的每個(gè)時(shí)間步都需要重新構(gòu)造系數(shù)矩陣并求解，這往往是多相流數(shù)值計(jì)算效率瓶頸的主要原因。為了提高多相流數(shù)值模擬的計(jì)算效率，學(xué)者們開展了許多研究工作，主要可以分為下述3類方法：第1類是Dong等[5]提出的常系數(shù)矩陣時(shí)間步進(jìn)算法求解Navier-Stokes和Cahn-Hilliard方程，避免每個(gè)時(shí)間步的系數(shù)矩陣構(gòu)造，并結(jié)合預(yù)處理技術(shù)，從而實(shí)現(xiàn)大密度比兩相流的高效計(jì)算；Dodd等[6]也根據(jù)常系數(shù)矩陣的構(gòu)造思想和快速傅里葉變換，結(jié)合VOF模型加快了大密度比兩相流毛細(xì)波的計(jì)算；Shen等[7]通過引入一組輔助標(biāo)量，將系數(shù)矩陣轉(zhuǎn)變?yōu)槎ǔ＞仃?，?shí)現(xiàn)了高效且無條件能量穩(wěn)定的梯度流求解。第2類是加速線性方程組本身的收斂過程，文獻(xiàn)[8?10]采用多重網(wǎng)格法在粗網(wǎng)格和細(xì)網(wǎng)格之間反復(fù)轉(zhuǎn)換，有效地提高了迭代算法的收斂速度；Fuka[11]發(fā)布了基于快速傅里葉變換的開源Poisson方程求解器，可用于基于偽譜法和有限差分法離散線性方程組的高效求解。第3類是并行算法的設(shè)計(jì)，Dolean等[12]詳細(xì)介紹了區(qū)域分解算法，主要思想是將計(jì)算域分解成多個(gè)區(qū)域并行計(jì)算；Djanali等[13]結(jié)合區(qū)域分解的思想求解泊松方程，加速了投影算法的收斂速度；Adam等[14]通過為泊松方程引入一個(gè)瞬態(tài)項(xiàng)，在有限差分法框架下實(shí)現(xiàn)了不可壓Navier-Stokes方程的完全顯式并行計(jì)算。除了上述的3類方法外，自適應(yīng)網(wǎng)格[15?16]，GPU并行[17?18]和人工神經(jīng)網(wǎng)絡(luò)[19?20]等技術(shù)也常被用于改善多相流的計(jì)算效率。

總之，隱式或半隱式求解不可壓多相流時(shí)，線性方程組的存在既占用了大量內(nèi)存，也嚴(yán)重影響了計(jì)算效率，這個(gè)難題在三維模型數(shù)值計(jì)算中將更加突出。文中提出了并行投影的有限體積方法，實(shí)現(xiàn)不可壓Navier-Stokes方程的全解耦顯式計(jì)算。同時(shí)，結(jié)合紅黑著色的交替并行策略，提高了算法的守恒性。在氣泡上升等測(cè)試算例中，結(jié)果表明，在保證計(jì)算精度的同時(shí)，本文算法的效率遠(yuǎn)高于串行的標(biāo)準(zhǔn)投影算法。

1 界面捕捉方程

選擇相場(chǎng)法[21]作為界面捕捉算法，兩相界面的運(yùn)動(dòng)通過Cahn-Hilliard方程進(jìn)行描述。

， （4）

， （5）

式中： 為相分?jǐn)?shù)變量，初始化法向分布遵循反正切函數(shù)方程（6），范圍為[-1， 1]，其中，1表示一相，-1表示另一相，-1～1的分?jǐn)?shù)則表示兩相界面。M為遷移率，它定義了相界面區(qū)域的擴(kuò)散率，這里采用1個(gè)固定值，M=0.01。 是化學(xué)勢(shì)，由自由能泛函推導(dǎo)得來[22]。 表示相界面的厚度。

。 （6）

文中所有算例中相變量 "和化學(xué)勢(shì) "的邊界條件均為Neumann邊界

，

。 （7）

流場(chǎng)物性參數(shù)由相變量線性插值計(jì)算

，

。 （8）

相場(chǎng)法通過連續(xù)表面力模型[22-23]計(jì)算表面張力，相較于VOF模型避免了額外的界面重構(gòu)和曲率計(jì)算，σ為表面張力系數(shù)

。 （9）

2 求解算法

單一流體模型下不可壓縮兩相流的計(jì)算需要循環(huán)求解Navier-Stokes方程和界面捕捉方程。首先，進(jìn)行初始化，包括計(jì)算域、時(shí)間步設(shè)置、網(wǎng)格劃分以及相界面初始化等；然后，隨著時(shí)間步進(jìn)，循環(huán)求解Cahn-Hilliard方程和Navier-Stokes方程，求解流程如圖1所示。

采用同位網(wǎng)格，如圖2所示，將所有變量都存儲(chǔ)在網(wǎng)格的體心。Cahn-Hilliard方程采用了易于并行化的二階Runge-Kutta TVD[24]算法顯式計(jì)算，本文的后續(xù)內(nèi)容主要關(guān)注不可壓Navier-Stokes方程的并行解耦，描述并行投影算法，再進(jìn)一步提出紅黑著色的并行投影算法。

2.1 并行投影算法

與標(biāo)準(zhǔn)的壓力修正投影算法方程（1）～（3）不同，并行投影算法[14]如下。

Function 1：預(yù)測(cè)步，由式（1）顯式計(jì)算中間速度 。

Function 2：校正步，考慮到當(dāng) "時(shí)，，則泊松方程可以修正為

。 （10）

將此方程用一階歐拉格式顯式推進(jìn)，循環(huán)求解K次（K=10～20），此處， ，h為網(wǎng)格寬度， ，因此，泊松方程得以顯式并行。

，

。 （11）

Function 3：最終步，更新得到 。

。 （12）

Adam等[14]基于有限差分法和中心差分格式實(shí)現(xiàn)了上述算法，有效提高了計(jì)算效率，但泊松方程中瞬態(tài)項(xiàng)的引入也帶來了額外的時(shí)間誤差，算法的守恒性和精度也有待提高。因此，為了提高并行投影算法的性能，采用守恒性更好的有限體積法，并結(jié)合5階WENO格式[25]離散對(duì)流項(xiàng)。此外，還提出了紅黑著色的交替并行策略以增強(qiáng)時(shí)空耦合。在本文的程序中，瞬態(tài)項(xiàng)采用一階歐拉格式，除了對(duì)流項(xiàng)以外的空間算子均采用中心差分格式。

2.2 紅黑著色的交錯(cuò)掃描策略

為了說明交錯(cuò)掃描策略，用F1（·）、F2（·）和F3（·）來表示2.1節(jié)所述的3個(gè)計(jì)算步驟，并將網(wǎng)格節(jié)點(diǎn)交替染成紅色和黑色，如圖2所示。

基于差分思想的數(shù)值方法（FVM， FDM）用差分代替微分來求解偏微分方程。以中心差分格式為例，每個(gè)節(jié)點(diǎn)的離散只與相鄰節(jié)點(diǎn)相關(guān)。因此，變量的更新可以交替執(zhí)行，第一步先由紅點(diǎn)的值計(jì)算黑點(diǎn)，然后再由上一步計(jì)算得到的黑點(diǎn)更新紅點(diǎn)的信息，如此反復(fù)執(zhí)行。結(jié)合紅黑著色的思想，對(duì)離散點(diǎn)進(jìn)行并行交替更新，分為以下6個(gè)步驟。

1）預(yù)測(cè)步（黑點(diǎn)）：黑色節(jié)點(diǎn)的 通過周圍紅色節(jié)點(diǎn)的值計(jì)算得到，分別用R和B作為下標(biāo)表示紅色節(jié)點(diǎn)和黑色節(jié)點(diǎn)。

。 （13）

2）校正步（紅點(diǎn)）：用上一步計(jì)算得到的黑點(diǎn)節(jié)點(diǎn)的 ，更新紅色節(jié)點(diǎn)的壓力。

。 （14）

3）最終步（黑點(diǎn)）：最后通過下式計(jì)算得到黑色節(jié)點(diǎn)的速度 。

。 （15）

至此完成了一半節(jié)點(diǎn)的計(jì)算，余下節(jié)點(diǎn)的計(jì)算采用類似的方式。

4）預(yù)測(cè)步（紅點(diǎn)）

。 （16）

5）校正步（黑點(diǎn)）

。 （17）

6）最終步（紅點(diǎn)）

。 （18）

紅黑著色并行算法的關(guān)鍵是節(jié)點(diǎn)信息的交替更新。完全顯式的并行投影算法解耦了節(jié)點(diǎn)之間的依賴關(guān)系，這使得計(jì)算易于并行化，文中使用OpenMP實(shí)現(xiàn)并行計(jì)算。

3 算例測(cè)試

通過3個(gè)經(jīng)典的兩相流算例測(cè)試了所提出算法的精度和效率，所有計(jì)算均為二維平面模型。

3.1 靜態(tài)液滴

初始時(shí)刻，一個(gè)直徑D = 0.4 m的靜態(tài)液滴放置于1 m×1 m二維平面方腔的中心。不考慮重力的作用，四壁均為無滑移邊界，兩相的密度和黏度比均為1，采用La數(shù)作為該算例的特征值。

。 （19）

通過改變密度來調(diào)整La數(shù)，而其他參數(shù)包括表面張力系數(shù)（1 N/m）和黏度（0.1 Pa?s）保持不變。設(shè)置網(wǎng)格量為64×64，時(shí)間步長dt=1×10-4。在實(shí)際的數(shù)值計(jì)算中，在相界面附近很難獲得精確的數(shù)值平衡，會(huì)產(chǎn)生所謂的“偽”或“寄生”速度，如圖3所示，速度矢量的分布符合Magnini等[26]和Mirjalili等[27]的預(yù)期。

為了使偽勢(shì)速度得到充分發(fā)展，計(jì)算了10 000步，用毛細(xì)數(shù)Ca來量化偽勢(shì)速度的強(qiáng)度。

。 （20）

圖4展示了網(wǎng)格細(xì)化后偽勢(shì)速度的演化過程（網(wǎng)格尺寸h分別為1/16、1/32、1/64、1/128），得到本算法的空間收斂階數(shù)在一階到二階之間。在本文所提出的算法框架下，相場(chǎng)模型比連續(xù)表面力模型[23]（一階）更快地收斂到尖銳界面極限，這與Mirjalili等[28]的測(cè)試是一致的。

3.2 氣泡上升

Hysing等[29]發(fā)布了一個(gè)二維氣泡上升的純數(shù)值標(biāo)準(zhǔn)算例，包含有2個(gè)測(cè)試用例，被廣泛用于兩相流算法的測(cè)試[30?31]。計(jì)算域和邊界條件如圖5所示。

該計(jì)算域是一個(gè)長度為1 m×2 m的全封閉腔體，氣泡初始位置為（x， y）=（0.5，0.5），單位m，初始直徑D=0.5 m。頂部和底部為無滑移固壁，左右為自由滑移固壁，重力指向定義域的底部。由于流體1（氣泡）的密度小于周圍流體2（液體），氣泡將在浮力作用下緩慢上升，圖6和圖7為不同時(shí)刻，氣泡的形狀。算例1和算例2的物性參數(shù)如表1所示。

為了滿足CFL條件和適當(dāng)?shù)慕缑鎸挾龋煌直媛氏碌挠?jì)算設(shè)置略有不同，網(wǎng)格寬度為h=1/64、h=1/128和h=1/256時(shí)，時(shí)間步長分別為1×10-4、5×10-5和2.5×10-5，界面寬度 "分別為0.02、0.01、0.005 m。對(duì)氣泡質(zhì)心位置和上升速度隨時(shí)間變化的定量曲線進(jìn)行比較，與參考文獻(xiàn)中的已發(fā)表數(shù)據(jù)吻合[31?32]，如圖8和圖9所示。由于不同的界面寬度下，Case 2中氣泡上升的速度略有差別，圖9（b）中，Xiao等[31]設(shè)置界面寬度為0.02 m時(shí)的氣泡上升速率曲線，和Aland等[32]設(shè)置界面寬度為0.01 m時(shí)的速率曲線，與本文所提出算法的測(cè)試結(jié)果一致。

基于3種不同的網(wǎng)格分辨率來比較本文所提出的并行算法與串行標(biāo)準(zhǔn)投影算法（使用預(yù)條件共軛梯度法求解線性方程組）的效率。計(jì)算采用的同樣的硬件條件（AMD Ryzen 7 4800H，Radeon Graphics 2.90 GHz）和設(shè)置。表2中列出了采用8線程OpenMP并行的C++程序在3種網(wǎng)格分辨率情況下的CPU時(shí)間。因?yàn)橥耆?guī)避了線性方程組的求解，本文所提出算法的計(jì)算效率遠(yuǎn)高于基于標(biāo)準(zhǔn)投影算法的串行程序，計(jì)算規(guī)模越大時(shí)加速效果也越顯著。

3.3 潰壩模型

潰壩算例是海洋與海岸工程水動(dòng)力學(xué)模擬中一個(gè)非常經(jīng)典的基準(zhǔn)算例，它與許多復(fù)雜的現(xiàn)象有關(guān)，比如，地表破裂、高壓射流和海洋造浪等。由于其密度比高，重力作用強(qiáng)，且流動(dòng)會(huì)導(dǎo)致劇烈的氣液界面拓?fù)渥兓?，是兩相流?shù)值模擬中一個(gè)重要且具挑戰(zhàn)性的算例。Martin等[33]使用高速攝像機(jī)在實(shí)驗(yàn)室進(jìn)行了基準(zhǔn)實(shí)驗(yàn)，捕捉選定時(shí)間間隔的流體運(yùn)動(dòng)情況。許多研究人員已經(jīng)通過數(shù)值方法[34?36]重復(fù)了這個(gè)案例，將其作為一個(gè)重要的基準(zhǔn)算例。文中計(jì)算域?yàn)?.8 m×0.6 m，液柱長寬為H=W=0.2 m，使用400×300的網(wǎng)格，時(shí)間步dt=1×10-7，模擬時(shí)間為0.9 s，無量綱時(shí)間

。 （21）

ρL=997，ρG=1.185，μL=8.9×10-3，μG=1.789 4×10-5，重力加速度g=9.81 m/s2，忽略表面張力。圖10是本文算法與已有實(shí)驗(yàn)[37]的比較，結(jié)果表明，本文所提出的算法具有較好的精度，具備捕捉實(shí)際流動(dòng)下大密度比兩相流界面變形的能力。圖11展示了自由液面的前端在t*=5.6時(shí)的形狀輪廓，數(shù)值測(cè)試與實(shí)驗(yàn)的對(duì)比結(jié)果吻合。

4 結(jié) "語

提出了求解Navier-Stokes和Cahn-Hilliard方程的紅黑著色并行有限體積投影算法，通過引入人工的壓力瞬態(tài)項(xiàng)，對(duì)泊松方程進(jìn)行不完全迭代，實(shí)現(xiàn)了不可壓Navier-Stokes方程的全解耦顯式計(jì)算。此外，結(jié)合紅黑著色的思想，交替更新離散點(diǎn)，實(shí)現(xiàn)了更緊湊的時(shí)空耦合。本算法完全避免了線性方程組的求解，并且可以對(duì)每個(gè)網(wǎng)格點(diǎn)分別進(jìn)行計(jì)算，只需簡單地使用幾行OpenMP代碼，就可以直接實(shí)現(xiàn)節(jié)點(diǎn)循環(huán)的并行化。

采用靜態(tài)液滴、氣泡上升和潰壩模型這3種經(jīng)典的兩相流基準(zhǔn)算例，對(duì)本文所提出的算法進(jìn)行了數(shù)值驗(yàn)證。結(jié)果表明，本文算法具備精確捕捉大密度比兩相流復(fù)雜界面變化的能力。更重要的是，在相同的設(shè)置和硬件環(huán)境下，本文算法的效率遠(yuǎn)高于串行投影算法，在十萬網(wǎng)格量時(shí)，加速比約80倍，且計(jì)算規(guī)模越大，加速效果越明顯。本文算法為工業(yè)流體中不可壓縮多相流的高性能計(jì)算提供了新的思路。

參考文獻(xiàn)

［1］ "Barros A P， Prat O P， Testik F Y. Size distribution of raindrops[J]. Nature Physics， 2010， 6（4）： 232.

［2］ "O’Sullivan J J， Norwood E A， O’Mahony J A， et al. Atomisation technologies used in spray drying in the dairy industry： a review[J]. Journal of Food Engineering， 2019， 243： 57-69.

［3］ "Atencia J， Beebe D J. Controlled microfluidic interfaces[J]. Nature， 2005， 437（7059）： 648-655.

［4］ "Guermond J L， Minev P， Shen J. An overview of projection methods for incompressible flows[J]. Computer Methods in Applied Mechanics and Engineering， 2006， 195（44/45/46/47）： 6011-6045.

［5］ "Dong S， Shen J. A time-stepping scheme involving constant coefficient matrices for phase-field simulations of two-phase incompressible flows with large density ratios[J]. Journal of Computational Physics， 2012， 231（17）： 5788-5804.

［6］ "Dodd M S， Ferrante A. A fast pressure-correction method for incompressible two-fluid flows[J]. Journal of Computational Physics， 2014， 273： 416-434.

［7］ "Shen J， Xu J， Yang J. The scalar auxiliary variable （SAV） approach for gradient flows[J]. Journal of Computational Physics， 2018， 353： 407-416.

［8］ "Lee C S， Hamon F P， Castelletto N， et al. An aggregation-based nonlinear multigrid solver for two-phase flow and transport in porous media[J]. Computers amp; Mathematics With Applications， 2022， 113： 282-299.

［9］ "Liu T. A nonlinear multigrid method for inverse problem in the multiphase porous media flow[J]. Applied Mathematics and Computation， 2018， 320： 271-281.

［10］ "Weymouth G D. Data-driven Multi-Grid solver for accelerated pressure projection[J]. Computers amp; Fluids， 2022， 246： 105620.

［11］ "Fuka V. PoisFFT： a free parallel fast Poisson solver[J]. Applied Mathematics and Computation， 2015， 267： 356-364.

［12］ "Dolean V， Jolivet P， Nataf F. An introduction to domain decomposition methods[M]. Philadelphia， PA： Society for Industrial and Applied Mathematics， 2015.

［13］ "Djanali V， Armfield S W， Kirkpatrick M P， et al. Parallel speed-up of preconditioned fractional step navier-stokes solvers[J]. Applied Mechanics and Materials， 2014， 493： 215-220.

［14］ "Adam N， Franke F， Aland S. A simple parallel solution method for the navier-stokes cahn-hilliard equations[J]. Mathematics， 2020， 8（8）： 1224.

［15］ "Kuhn M B， Deskos G， Sprague M A. A mass-momentum consistent coupling for mesh-adaptive two-phase flow simulations[J]. Computers amp; Fluids， 2023， 252： 105770.

［16］ "Ding L， Shu C， Ding H， et al. Stencil adaptive diffuse interface method for simulation of two-dimensional incompressible multiphase flows[J]. Computers amp; Fluids， 2010， 39（6）： 936-944.

［17］ "Sakane S， Takaki T， Rojas R， et al. Multi-GPUs parallel computation of dendrite growth in forced convection using the phase-field-lattice Boltzmann model[J]. Journal of Crystal Growth， 2017， 474： 154-159.

［18］ "Ha S， Park J， You D. A GPU-accelerated semi-implicit fractional-step method for numerical solutions of incompressible Navier-Stokes equations[J]. Journal of Computational Physics， 2018， 352： 246-264.

［19］ "Qiu R D， Huang R F， Xiao Y， et al. Physics-informed neural networks for phase-field method in two-phase flow[J]. Physics of Fluids， 2022， 34（5）： 052109.

［20］ "Raissi M， Perdikaris P， Karniadakis G E. Physics-informed neural networks： a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[J]. Journal of Computational Physics， 2019， 378： 686-707.

［21］ "Ding H， Spelt P D M， Shu C. Diffuse interface model for incompressible two-phase flows with large density ratios[J]. Journal of Computational Physics， 2007， 226（2）： 2078-2095.

［22］ "Liu H H， Valocchi A J， Zhang Y H， et al. Lattice Boltzmann phase-field modeling of thermocapillary flows in a confined microchannel[J]. Journal of Computational Physics， 2014， 256： 334-356.

［23］ "Kim J. A continuous surface tension force formulation for diffuse-interface models[J]. Journal of Computational Physics， 2005， 204（2）： 784-804.

［24］ "Gottlieb S， Shu C W. Total variation diminishing Runge-Kutta schemes[J]. Mathematics of Computation， 1998， 67（221）： 73-85.

［25］ "Friedrich O. Weighted essentially non-oscillatory schemes for the interpolation of mean values on unstructured grids[J]. Journal of Computational Physics， 1998， 144（1）： 194-212.

［26］ "Magnini M， Pulvirenti B， Thome J R. Characterization of the velocity fields generated by flow initialization in the CFD simulation of multiphase flows[J]. Applied Mathematical Modelling， 2016， 40（15/16）： 6811-6830.

［27］ "Mirjalili S， Ivey C B， Mani A L. Comparison between the diffuse interface and volume of fluid methods for simulating two-phase flows[J]. International Journal of Multiphase Flow， 2019， 116： 221-238.

［28］ "Mirjalili S， Khanwale M A， Mani A L. Assessment of an energy-based surface tension model for simulation of two-phase flows using second-order phase field methods[J]. Journal of Computational Physics， 2023， 474： 111795.

［29］ "Hysing S， Turek S， Kuzmin D， et al. Quantitative benchmark computations of two-dimensional bubble dynamics[J]. International Journal for Numerical Methods in Fluids， 2009， 60（11）： 1259-1288.

［30］ "Klostermann J， Schaake K， Schwarze R. Numerical simulation of a single rising bubble by VOF with surface compression[J]. International Journal for Numerical Methods in Fluids， 2013， 71（8）： 960-982.

［31］ "Xiao Y， Zeng Z， Zhang L Q， et al. A spectral element-based phase field method for incompressible two-phase flows[J]. Physics of Fluids， 2022， 34（2）： 022114.

［32］ "Aland S， Voigt A. Benchmark computations of diffuse interface models for two-dimensional bubble dynamics[J]. International Journal for Numerical Methods in Fluids， 2012， 69（3）： 747-761.

［33］ "Martin J C， Moyce W J， Penney W G， et al. Part IV. An experimental study of the collapse of liquid columns on a rigid horizontal plane[J]. Philosophical Transactions of the Royal Society of London Series A， Mathematical and Physical Sciences， 1952， 244（882）： 312-324.

［34］ "Yu C H， Sheu T W H. Development of a coupled level set and immersed boundary method for predicting dam break flows[J]. Computer Physics Communications， 2017， 221： 1-18.

［35］ "Meng W K， Liao L， Chen M， et al. An enhanced CLSVOF method with an algebraic second-reconstruction step for simulating incompressible two-phase flows[J]. International Journal of Multiphase Flow， 2022， 154： 104151.

［36］ "Li Y L， Wan L， Wang Y H， et al. Numerical investigation of interface capturing method by the Rayleigh-Taylor instability， dambreak and solitary wave problems[J]. Ocean Engineering， 2019， 194： 106583.

［37］ "Kamra M M， Mohd N， Liu C， et al. Numerical and experimental investigation of three-dimensionality in the dam-break flow against a vertical wall[J]. Journal of Hydrodynamics， 2018， 30（4）： 682-693.

（編輯 "鄭潔）

doi： 10.11835/j.issn.1000-582X.2023.259

收稿日期：2023-03-02

網(wǎng)絡(luò)出版日期：2023-07-04

基金項(xiàng)目：國家自然科學(xué)基金資助項(xiàng)目（12102071，12172070）；中央高?；究蒲袠I(yè)務(wù)費(fèi)資助項(xiàng)目（2021CDJQY-055）。

Foundation：Supported by National Natural Science Foundation of China （12102071， 12172070）， and the Fundamental Research Funds for the Central Universities （2021CDJQY-055）.

作者簡介：王小雙（1995—），男，碩士研究生，主要從事不可壓多相流算法研究與程序開發(fā)方向研究，（E-mail）wangxs1995s@cqu.edu.cn。

通信作者：張良奇，男，研究員，博士生導(dǎo)師，主要從事計(jì)算流體力學(xué)特色數(shù)值方法研究及其在前沿流體力學(xué)問題中的應(yīng)用研究，（E-mail）zhangliangqi@cqu.edu.cn。