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

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

中圖分類號：U448.213 " " " " "文獻標志碼：A " " "文章編號：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]等。對兩相流界面動力學(xué)的深入研究，不僅可以幫助人們更好地理解自然現(xiàn)象，還可以為工程技術(shù)和生物醫(yī)學(xué)領(lǐng)域的發(fā)展提供重要的理論和實踐支撐。然而，兩相流精確且高效的數(shù)值計算仍然是計算流體動力學(xué)的重要挑戰(zhàn)。當(dāng)采用基于壓力修正的投影算法[4]求解不可壓Navier-Stokes方程時，往往將其解耦為3個方程，依次計算得到中間速度 ，當(dāng)前壓力 和當(dāng)前速度

， （1）

， （2）

， （3）

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

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

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

1 界面捕捉方程

選擇相場法[21]作為界面捕捉算法，兩相界面的運動通過Cahn-Hilliard方程進行描述。

， （4）

， （5）

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

。 （6）

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

，

。 （7）

流場物性參數(shù)由相變量線性插值計算

，

。 （8）

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

。 （9）

2 求解算法

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

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

2.1 并行投影算法

與標準的壓力修正投影算法方程（1）～（3）不同，并行投影算法[14]如下。

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

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

。 （10）

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

，

。 （11）

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

。 （12）

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

2.2 紅黑著色的交錯掃描策略

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

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

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

。 （13）

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

。 （14）

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

。 （15）

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

4）預(yù)測步（紅點）

。 （16）

5）校正步（黑點）

。 （17）

6）最終步（紅點）

。 （18）

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

3 算例測試

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

3.1 靜態(tài)液滴

初始時刻，一個直徑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，時間步長dt=1×10-4。在實際的數(shù)值計算中，在相界面附近很難獲得精確的數(shù)值平衡，會產(chǎn)生所謂的“偽”或“寄生”速度，如圖3所示，速度矢量的分布符合Magnini等[26]和Mirjalili等[27]的預(yù)期。

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

。 （20）

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

3.2 氣泡上升

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

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

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

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

3.3 潰壩模型

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

。 （21）

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

4 結(jié) "語

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

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

參考文獻

［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

基金項目：國家自然科學(xué)基金資助項目（12102071，12172070）；中央高?；究蒲袠I(yè)務(wù)費資助項目（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)師，主要從事計算流體力學(xué)特色數(shù)值方法研究及其在前沿流體力學(xué)問題中的應(yīng)用研究，（E-mail）zhangliangqi@cqu.edu.cn。