一種不可壓縮二維流動的顯式逐次超松弛并行算法

張曉慧， 柏君勵， 顧解忡， 馬 寧

(上海交通大學(xué) 海洋工程國家重點實驗室； 高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心；船舶海洋與建筑工程學(xué)院， 上海 200240)

計算流體力學(xué)方法已逐漸成為流動及流體傳熱問題研究的主流工具.研究表明[1-4]，與勢流方法相比，采用計算流體力學(xué)方法能夠更準(zhǔn)確地預(yù)報船舶的耐波性.其中，主要方法是采用有限體積法(FVM)離散Navier-Stokes方程而形成代數(shù)方程組，再進(jìn)行迭代求解.然而，采用計算流體力學(xué)方法求解一個算例需要耗費大量時間，尤其是迭代計算.因為每個時間層涉及多次外迭代，每次外迭代又涉及速度、壓力或壓力修正值的內(nèi)迭代，所以提高迭代計算的效率已經(jīng)成為計算流體力學(xué)的一個關(guān)鍵技術(shù)問題，而實行并行計算是提高計算效率的一條行之有效的途徑.隨著分組顯式思想和區(qū)域分解方法[5]的提出，為并行計算創(chuàng)造了條件.目前，并行求解線性代數(shù)方程組的迭代方法主要包括Schwarz法、多重網(wǎng)格法、共軛梯度(CG)法以及經(jīng)典的定常迭代法等.Schwarz法基于區(qū)域分解的思想，易于并行化但計算的收斂速度較慢.例如，文獻(xiàn)[6]中通過乘性Schwarz方法以及預(yù)處理過程來加速迭代的收斂速度，但其很難程序化.多重網(wǎng)格方法[7]是通過將計算過程轉(zhuǎn)換到粗網(wǎng)格中以加快迭代收斂速度.另外，Krylov子空間迭代法[8]包含了共軛梯度法及其擴(kuò)展方法，是另一類較為重要的迭代方法.考慮到計算的效率和可操作性，本文選取經(jīng)典的定常迭代法.雖然Jacobi迭代法容易實行并行操作，但其比Gauss-Seidel和逐次超松弛(SOR)算法的收斂速度慢.而相比于Gauss-Seidel算法，SOR算法的收斂速度更快.文獻(xiàn)[9]中基于區(qū)域分解思想并結(jié)合SOR方法對Stokes方程進(jìn)行并行計算.

最近，許秋燕等[10]結(jié)合區(qū)域分解和分組顯式方法將計算域分割成多個子域，構(gòu)造了子域并行迭代計算格式，以用于有限差分求解特定的泊松方程，最終求解的是常系數(shù)代數(shù)方程組.而對于流動以及流體傳熱問題，代數(shù)方程組的系數(shù)將隨著外迭代和時間層的不斷推進(jìn)而不斷變化.

本文采用船舶水動力學(xué)研究中常用的FVM求解不可壓縮Navier-Stokes方程，從求解二維流體流動入手，在許秋燕等[10]研究的基礎(chǔ)上，提出一種基于FVM的顯式逐次超松弛并行算法(FV-pSOR).其中，采用4階緊致格式離散流項和擴(kuò)散項、4階龍格-庫塔法步進(jìn)非定常項以及同位網(wǎng)格SIMPLE(Semi-Implicit Method for Pressure-Linked Equation)算法求解壓力-速度耦合方程[11].基于區(qū)域分解，將計算域分割成4個子域而構(gòu)造所需分組顯式SOR迭代公式，并通過3種網(wǎng)格密度下典型二維驅(qū)動方腔流的計算來驗證所提FV-pSOR算法的有效性.

1 數(shù)學(xué)模型

1.1 Navier-Stokes方程組及其數(shù)值離散

不可壓縮流體的流動服從如下控制方程組：

(1)

(2)

本文并不是直接求解微分形式的控制方程組(式(1)、(2))，而是求解在同位網(wǎng)格中由FVM導(dǎo)出的積分形式的控制方程[11].同位網(wǎng)格是將速度、壓力、黏度等物理量均定義在網(wǎng)格中心點處，如圖1所示.其中，網(wǎng)格P的積分方程為

(3)

式中：Fc為對流項；FD為擴(kuò)散項.

對流項可用前一時刻的速度進(jìn)行線性化，即

(4)

擴(kuò)散項為

(5)

圖1 同位網(wǎng)格Fig.1 Collocated grid

式(4)、(5)中的變量包括速度和速度導(dǎo)數(shù)，采用4階緊致格式插值[11].式(3)中時間層的步進(jìn)采用4階龍格-庫塔法[11]，龍格-庫塔法的每一步都采用 SIMPLE 算法求解速度-壓力耦合方程.

經(jīng)上述離散后，最終所得壓力修正方程為

(6)

下面將具體說明式(6)的顯式SOR并行迭代算法的構(gòu)造過程.

1.2 FV-pSOR算法的基礎(chǔ)迭代格式

FV-pSOR算法有4種基礎(chǔ)SOR迭代格式，如圖2所示.圖中，橫向代表第1個下標(biāo)i，垂向代表第2個下標(biāo)k，箭頭表示迭代計算的行進(jìn)方向.對于式(6)，4種基礎(chǔ)SOR迭代格式的具體表達(dá)式如下：

基礎(chǔ)迭代格式A1為

(7)

基礎(chǔ)迭代格式A2為

(8)

基礎(chǔ)迭代格式A3為

(9)

基礎(chǔ)迭代格式A4為

(10)

圖2 4種基礎(chǔ)迭代格式Fig.2 Four iterative schemes

根據(jù)分組顯式的思想，由上述4種基礎(chǔ)SOR迭代格式構(gòu)造的迭代格式A5見圖3.在節(jié)點(i,k+1)處采用基礎(chǔ)格式A3，節(jié)點(i,k)處采用基礎(chǔ)格式A4，節(jié)點(i+1,k)處采用基礎(chǔ)格式A2，節(jié)點(i+1,k+1)處采用基礎(chǔ)格式A1.因為當(dāng)前第n次迭代的值已知，所以只有4個節(jié)點處第n+1次迭代的值為

圖3 基本迭代格式A5示意圖Fig.3 Sketch of scheme A5

未知，可由以下4個基礎(chǔ)格式構(gòu)造，即

(11)

聯(lián)立求解可得

從而，由式(11)所得迭代格式A5的顯式表達(dá)式為

(12)

式中：

k=-1+b3b7+b1(b5-b3b5b7+b4b6b7)+

b6b8+b2(b4+b3b5b8-b4b6b8)

同樣地，由基礎(chǔ)SOR格式可以構(gòu)造圖4所示的另外4種基礎(chǔ)迭代格式.鑒于這4種基礎(chǔ)迭代格式的構(gòu)造方法類似，本文僅說明格式A6的構(gòu)造過程.在節(jié)點(i,k)處采用基礎(chǔ)格式A3，節(jié)點(i+1,k)處采用基礎(chǔ)格式A1，對應(yīng)的迭代方程組為

(13)

式中：

圖4 基本迭代格式A6、A7、A8及A9示意圖Fig.4 Sketch of scheme A6,A7,A8 and A9

求解式(13)，所得格式A6的顯式表達(dá)式為

(14)

同樣地，可以得到其他3種迭代格式.其中，格式A7為

(15)

式中：

格式A8為

(16)

式中：

格式A9為

(17)

式中：

至此，所需顯式SOR迭代格式全部構(gòu)造完成.

圖5 偶次迭代行進(jìn)過程示意圖Fig.5 Sketch of the process of even iteration

2 迭代計算過程

圖5和6示出了偶次迭代和奇次迭代的行進(jìn)順序.圖中，計算域被分割成4個子域.由圖5可見，偶次迭代順序是順次從中心向外行進(jìn)的.其迭代步驟如下：

圖6 奇次迭代行進(jìn)過程示意圖Fig.6 Sketch of the process of odd iteration

(1) 選定兩個坐標(biāo)方向的中點xp、xp+1和zq、zq+1，共有4個位置點，區(qū)域Λ由4個選定的點構(gòu)成.按照格式A5的迭代公式算出這4個點在下一個迭代步的數(shù)值.

(2) 根據(jù)箭頭的行進(jìn)方向，分別在區(qū)域Γ1、Γ2、Γ3和Γ4中按照格式A6、A7、A8和A9的迭代公式依次計算當(dāng)前區(qū)域中的點在下一個迭代步的數(shù)值.以Γ1為例，計算次序為

i=xp,xp+1
k=zq+2,zq+3,…,kmax

其中：kmax為k方向的最大網(wǎng)格數(shù).

(3) 在區(qū)域Ω1、Ω2、Ω3和Ω4中分別按照格式A1、A2、A3和A4的迭代公式依次計算各區(qū)域中的點在下一個迭代步的數(shù)值.

完成上述3步迭代后，所有點處的值均已得到更新.

如圖6所示，奇次迭代的迭代順序是從四角向中心行進(jìn).在區(qū)域Ω1、Ω2、Ω3和Ω4中分別按照格式A4、A3、A2和A1的迭代公式依次計算各區(qū)域中的點在下一個迭代步的數(shù)值.迭代完成后，所有點處的值均已得到更新.

至此，偶次迭代和奇次迭代的過程構(gòu)建完成.計算時，奇次迭代和偶次迭代交替進(jìn)行，直至精度滿足控制要求，然后，進(jìn)入下一個外迭代層或下一個時間層.

3 算例

以二維驅(qū)動方腔流為例來驗證FV-pSOR算法的有效性.圖7所示為二維方腔流在左右側(cè)面和底面需滿足的固壁無滑移條件，其頂面為給定的自左向右均勻剪切流動.水平速度設(shè)為u=1 m/s，垂向速度為v=0.

圖7 二維方腔流邊界條件Fig.7 Boundary conditions of 2D lid-driven cavity flow

取流動區(qū)域的邊長為h=1 m，流動的雷諾數(shù)Re=uh/ν=1 000.計算時，將矩形計算域劃分為均勻網(wǎng)格，網(wǎng)格數(shù)量N分別為50×50、100×100以及150×150.并行計算時，將矩形計算域分割成4個子域.本文分析方腔流發(fā)展到35 s左右時的計算結(jié)果.采用一臺處理器為Intel酷睿i7-6700的計算機(jī)，調(diào)用DATE_AND_TIME函數(shù)所獲兩種算法在3種網(wǎng)格數(shù)量下的計算耗時見表1.可見：在網(wǎng)格數(shù)量N=50×50下，SOR算法的耗時超出FV-pSOR算法耗時12倍；在網(wǎng)格分別為100×100和150×150下，SOR算法的耗時分別超出FV-pSOR算法的耗時6和7倍.可見，對于相同的子域分割，采用細(xì)網(wǎng)格有利于提高計算效率.

圖8～10為網(wǎng)格數(shù)量分別為50×50、100×100以及150×150下的水平速度分布云圖、垂向速度v沿y=0.5 m 水平線的分布圖、水平速度u沿x=0.5 m 垂直線的分布圖.其中，3種網(wǎng)格數(shù)量下的云圖分辨率相同.為了便于對比，圖9和10中同時給出了SOR算法和文獻(xiàn)[12]中的速度分布結(jié)果.可以看出，本文提出的FV-pSOR算法的結(jié)果與SOR算法的結(jié)果完全相同(差值為10-4量級)，而且與文獻(xiàn)[12]中的結(jié)果相符，表明本文提出的FV-pSOR算法與SOR算法具有同等計算精度.由圖9和10可以看出，本文提出的FV-pSOR算法與SOR算法的計算精度差異不大.而由表1可見，其計算效率是SOR算法計算效率的6～7倍，說明FV-pSOR算法能夠適應(yīng)數(shù)值求解Navier-Stokes方程的五對角線性方程組的并行迭代求解.

表1 網(wǎng)格數(shù)量分別為50×50、100×100及150×150時的計算耗時Tab.1 Time consumption under 50×50,100×100 and 150×150 grids

圖8 水平速度u的分布云圖Fig.8 Distribution of velocity component u

圖9 垂向速度v沿水平線y=0.5 m的分布情況Fig.9 Velocity v along horizontal direction y=0.5 m

圖10 水平速度u沿垂直線x=0.5 m的分布情況Fig.10 Velocity u along vertical direction x=0.5 m

綜上所述，本文提出的FV-pSOR算的計算效率明顯高于SOR算法，但在細(xì)網(wǎng)格條件下其計算效率的提高幅度比粗網(wǎng)格下的有所降低.

4 結(jié)語

本文在有限差分并行迭代算法的基礎(chǔ)上提出了FV-pSOR算法，用于提高SOR方法求解Navier-Stokes方程的計算效率，并通過一個典型的二維驅(qū)動方腔流算例的計算，驗證了FV-pSOR算法的有效性.結(jié)果表明，F(xiàn)V-pSOR算法的計算效率較高，且其計算精度與SOR算法差異不大.雖然FV-pSOR算法的計算格式是基于計算流體力學(xué)流動問題的數(shù)值求解提出的，但其同樣適用于其他五對角線性方程組的并行迭代求解.