基于并行BICGSTAB的不可壓流場計算方法

摘" 要：不可壓縮流動廣泛存在于日常生活及工程應用之中，隨著數(shù)值模擬規(guī)模不斷擴大，CFD算法越來越重視并行計算。該文采用SIMPLE算法對不可壓流場進行求解，為求解算法中的大型稀疏矩陣，該文發(fā)展一套基于MPI框架的并行BICGSTAB矩陣求解方法。最后采用定常頂蓋驅(qū)動流算例及非定常圓柱繞流算例驗證并行算法的可行性及準確性。

關鍵詞：不可壓縮流動；SIMPLE算法；BICGSTAB算法；并行計算；數(shù)值模擬

中圖分類號：V211.3" " " 文獻標志碼：A" " " " " 文章編號：2095-2945（2024）25-0149-05

Abstract： Incompressible flows are prevalent in daily life and engineering applications. With the continuous expansion of numerical simulation scales， CFD algorithms increasingly emphasize parallel computing. This paper utilizes the SIMPLE algorithm to solve incompressible flow. We developed a parallel BICGSTAB matrix solver based on the MPI framework to address the large sparse matrices in the algorithm. Finally， the feasibility and accuracy of the parallel algorithm are validated using the steady lid-driven cavity flow and unsteady flow around cylinder test cases.

Keywords： incompressible flow; SIMPLE; BICGSTAB; parallel computation; numerical simulation

不可壓縮流動廣泛存在于日常生活及工程應用中。一般情況下，對于液體及低馬赫數(shù)（馬赫數(shù)Malt;0.3）的氣體認為是不可壓縮流體。

對于不可壓縮流場求解方法最常見的有人工壓縮法及壓力修正法。壓力修正法中最具有代表性的算法為SIMPLE[1]系列算法，包括SIMPLE、SIMPLER、SIMPLEC和PISO算法等。在SIMPLE算法中，需要分別求解動量方程及壓力修正方程。方程通過空間及離散后會形成線性方程組，對方程的求解最后都將歸結于求解一個大型稀疏矩陣。

目前，對于大型稀疏矩陣的求解有多種方法，包括經(jīng)典迭代法及Krylov子空間方法。其中，Krylov子空間方法是求解大型非對稱稀疏矩陣的最有效方法之一，其收斂速度快且不會因矩陣規(guī)模的增大而降低效率。而雙共軛梯度穩(wěn)定法[2]（Biconjugate Gradient Stabilized Method，BICGSTAB）則是Krylov子空間方法中一種經(jīng)典算法，目前被廣泛應用于求解大型稀疏矩陣。

隨著計算機的發(fā)展及數(shù)值模擬規(guī)模的不斷擴大，串行計算顯然無法滿足工程需要。因此，CFD算法越來越看重并行計算，可以說并行計算必不可少。

為了提高求解效率，本文在SIMPLE算法的基礎上發(fā)展了一套基于并行BICGSTAB的不可壓流場求解方法。對于線程之間的數(shù)據(jù)傳輸過程，本文通過調(diào)用MPI庫函數(shù)來完成。采用定常頂蓋驅(qū)動流算例及非定常圓柱繞流算例驗證算法的可行性。

1" 數(shù)值求解

1.1" 不可壓縮流動控制方程

空間離散采用有限體積法，對流項采用二階迎風格式，時間離散采用二階迎風歐拉形式。整體采用同位網(wǎng)格SIMPLE算法進行求解?？紤]到SIMPLE算法流程已有諸多文獻涉及，這里不再做過多贅述。

1.2" BICGSTAB矩陣求解方法

通過對方程進行離散，目的就是為了形成一個線性方程組，流場的求解最后都會歸結于求解線性方程組

大型非對稱稀疏矩陣的求解一直以來都是計算力學中的重點。本文采用雙共軛梯度穩(wěn)定法BICGSTAB進行大型稀疏矩陣的迭代求解。BICGSTAB的迭代過程主要涉及一系列矩陣與數(shù)組的計算，流程如下。

本文以式（4）中r0的二范數(shù)|r0|作為方程收斂判斷，如果方程收斂，那么|r0|一定是趨近于0的。因此在每個迭代步內(nèi)進行矩陣迭代求解之前都會計算|r0|并將其輸出，以此作為判斷方程收斂的參考依據(jù)。

2" 并行計算

并行計算就是將整體網(wǎng)格分為若干個部分，再分配給多個線程，每個線程負責一部分網(wǎng)格，分別進行單獨計算。并行計算大體可分為2部分：一是網(wǎng)格并行分區(qū)及并行邊界的處理，二是矩陣的并行計算。對于線程之間的數(shù)據(jù)傳輸過程，本文通過調(diào)用MPI庫函數(shù)來完成。

網(wǎng)格劃分形成的邊界稱為并行邊界。并行邊界本質(zhì)上屬于網(wǎng)格內(nèi)部面，只是由于并行分區(qū)而失去另一邊的單元。為了能讓并行邊界正常參與離散計算，通常會在并行邊界另一側(cè)設置虛擬網(wǎng)格單元，用以彌補因并行分區(qū)而缺失的單元，如圖2所示。每個并行邊界上的虛擬網(wǎng)格單元都能在其他線程上找到與之對應的單元，因此，在網(wǎng)格劃分之時需要通過網(wǎng)格單元編號建立起一一對應的關系。每一個虛擬網(wǎng)格單元的編號都唯一對應了一個物理網(wǎng)格單元編號，在需要用到虛擬單元上的數(shù)據(jù)時，就可以根據(jù)編號從其他線程找出對應的單元，并通過調(diào)用MPI庫函數(shù)完成數(shù)據(jù)的傳輸。

并行邊界本質(zhì)上屬于內(nèi)部面，其離散方法也與內(nèi)部單元面相同，唯一不同點在于另一側(cè)是虛擬單元，對應的網(wǎng)格編號也是虛擬單元的編號。虛擬單元不參與離散，只是用于為并行邊界內(nèi)側(cè)單元提供相應流場參數(shù)。因此，稀疏矩陣的規(guī)模不會發(fā)生變化。如圖3所示，每個線程內(nèi)都存有部分的系數(shù)矩陣，矩陣虛線右側(cè)部分則代表了虛擬網(wǎng)格單元。

BICGSTAB求解過程主要包括3種類型的數(shù)組運算：向量的數(shù)乘、向量的點乘、矩陣與向量相乘，即

第一類為向量的數(shù)乘，向量分段存儲于各個線程之中，只需要將各自的局部向量數(shù)組作點乘即可，無須進行線程間的數(shù)據(jù)傳輸。

第二類為向量的點乘，向量分段存儲于各個線程之中，需要先在每個線程內(nèi)作局部向量點乘，再將各個線程內(nèi)的點乘結果相加就是最終的結果，可以通過調(diào)用MPI_Allreduce（）函數(shù)進行歸約求和。

第三類為矩陣與向量相乘，系數(shù)矩陣與向量皆分段存儲于各個線程之中。根據(jù)有限體積法離散的系數(shù)矩陣每一行對應于一個網(wǎng)格單元，網(wǎng)格單元只與鄰單元有關，反映到系數(shù)矩陣中就是每一行中除了自身單元及鄰單元編號外的所有系數(shù)都為0。因此，在計算Ax時作單元循環(huán)，每一行只需要對該單元編號及鄰單元編號的系數(shù)進行對應乘積即可。在遇到并行邊界處單元時需要根據(jù)對應虛擬單元編號找到相應虛擬單元參與計算。需要注意的是，第三類計算在矩陣每一次迭代循環(huán)中都需要更新對應的虛擬網(wǎng)格上的向量參數(shù)。

3" 算例驗證

本文將以頂蓋驅(qū)動流及圓柱繞流這2個經(jīng)典算例驗證并行算法的可行性及精確性。

3.1" 定常頂蓋驅(qū)動流

頂蓋驅(qū)動流算例物理模型示意圖如圖4（a）所示，一個長寬皆為1的密閉方腔，四周都是壁面，頂部為移動壁面，其余皆為靜止壁面。圖4（b）為計算網(wǎng)格，采用混合網(wǎng)格，第一層網(wǎng)格高度為1×10-4，網(wǎng)格總量為17 700。

3.2" 非定常圓柱繞流

圖12展示了雷諾數(shù)Re=200下圓柱表面升阻力系數(shù)的變化曲線，從圖12中可以看出，流場從初始狀態(tài)開始發(fā)展，大約在t=60 s后，流場開始穩(wěn)定，升阻力系數(shù)呈現(xiàn)穩(wěn)定的周期性振蕩。

通過計算得出平均阻力系數(shù)為1.33，升力系數(shù)則是圍繞0一直振蕩。通常情況下，用Strouhal數(shù)描述升力系數(shù)的振蕩頻率，定義如下

本文在表1中將計算結果與文獻[5]及實驗[6]結果進行對比，并進行了誤差計算。從表1中可以看出，本文計算結果與文獻[5]及實驗結果十分接近，誤差處于可以接受的范圍內(nèi)，能夠滿足工程需求。

4" 結束語

本文采用SIMPLE算法對不可壓流場進行求解，并且在MPI并行框架下開發(fā)了一套并行BICGSTAB矩陣求解方法。最后以定常頂蓋驅(qū)動流及非定常圓柱繞流算例對算法進行了驗證，計算結果與文獻[5]及實驗結果進行了對比。從對比結果可以看出本文提出的并行算法是可行的，且計算結果具有較高的精度。

參考文獻：

[1] BARTON I E. Comparison of SIMPLE and PISO-type algorithms for transient flows[J]. International Journal for Numerical Methods in Fluids，1998，26（4）：459-483.

[2] VORST H A V D. Bi-CGSTAB： A fast and smoothly converging variant of bi-cg for the solution of nonsymmetric linear systems[J]. SIAM Journal on Scientific and Statistical Computing，1992，13（2）：0913035.

[3] BLAZEK J. Computational Fluid Dynamics： Principles and Applications[M]. Computational Fluid Dynamics Principl- es amp; Applications，2001.

[4] GHIA U， GHIA K N， SHIN C T. High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method[J]. Journal of Computational Physics，1982，48（3）：387-411.

[5] 干雨新.基于笛卡爾網(wǎng)格的復雜流動問題數(shù)值模擬[D].南京：南京航空航天大學，2023.

[6] WILLE R. Karman vortex streets[Z]. Advances in Applied Mechanics，1960.