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

關(guān)鍵詞：不可壓縮流動(dòng)；SIMPLE算法；BICGSTAB算法；并行計(jì)算；數(shù)值模擬

中圖分類號(hào)：V211.3" " " 文獻(xiàn)標(biāo)志碼：A" " " " " 文章編號(hào)：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

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

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

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

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

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

1" 數(shù)值求解

1.1" 不可壓縮流動(dòng)控制方程

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

1.2" BICGSTAB矩陣求解方法

通過對(duì)方程進(jìn)行離散，目的就是為了形成一個(gè)線性方程組，流場(chǎng)的求解最后都會(huì)歸結(jié)于求解線性方程組

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

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

2" 并行計(jì)算

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

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

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

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

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

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

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

3" 算例驗(yàn)證

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

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

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

3.2" 非定常圓柱繞流

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

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

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

4" 結(jié)束語

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

參考文獻(xiàn)：

[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)格的復(fù)雜流動(dòng)問題數(shù)值模擬[D].南京：南京航空航天大學(xué)，2023.

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