紅黑著色的相場兩相流并行投影算法

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

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

中圖分類號：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

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

， （1）

， （2）

， （3）

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

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

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

1 界面捕捉方程

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

， （4）

， （5）

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

。 （6）

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

，

。 （7）

流場物性參數由相變量線性插值計算

，

。 （8）

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

。 （9）

2 求解算法

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

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

2.1 并行投影算法

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

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

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

。 （10）

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

，

。 （11）

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

。 （12）

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

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

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

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

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

。 （13）

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

。 （14）

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

。 （15）

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

4）預測步（紅點）

。 （16）

5）校正步（黑點）

。 （17）

6）最終步（紅點）

。 （18）

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

3 算例測試

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

3.1 靜態液滴

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

。 （19）

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

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

。 （20）

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

3.2 氣泡上升

Hysing等[29]發布了一個二維氣泡上升的純數值標準算例，包含有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的物性參數如表1所示。

為了滿足CFL條件和適當的界面寬度，不同分辨率下的計算設置略有不同，網格寬度為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。對氣泡質心位置和上升速度隨時間變化的定量曲線進行比較，與參考文獻中的已發表數據吻合[31?32]，如圖8和圖9所示。由于不同的界面寬度下，Case 2中氣泡上升的速度略有差別，圖9（b）中，Xiao等[31]設置界面寬度為0.02 m時的氣泡上升速率曲線，和Aland等[32]設置界面寬度為0.01 m時的速率曲線，與本文所提出算法的測試結果一致。

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

3.3 潰壩模型

潰壩算例是海洋與海岸工程水動力學模擬中一個非常經典的基準算例，它與許多復雜的現象有關，比如，地表破裂、高壓射流和海洋造浪等。由于其密度比高，重力作用強，且流動會導致劇烈的氣液界面拓撲變化，是兩相流數值模擬中一個重要且具挑戰性的算例。Martin等[33]使用高速攝像機在實驗室進行了基準實驗，捕捉選定時間間隔的流體運動情況。許多研究人員已經通過數值方法[34?36]重復了這個案例，將其作為一個重要的基準算例。文中計算域為0.8 m×0.6 m，液柱長寬為H=W=0.2 m，使用400×300的網格，時間步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]的比較，結果表明，本文所提出的算法具有較好的精度，具備捕捉實際流動下大密度比兩相流界面變形的能力。圖11展示了自由液面的前端在t*=5.6時的形狀輪廓，數值測試與實驗的對比結果吻合。

4 結 "語

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

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

參考文獻

［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

網絡出版日期：2023-07-04

基金項目：國家自然科學基金資助項目（12102071，12172070）；中央高校基本科研業務費資助項目（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—），男，碩士研究生，主要從事不可壓多相流算法研究與程序開發方向研究，（E-mail）wangxs1995s@cqu.edu.cn。

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