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

張曉慧， 柏君勵(lì)， 顧解忡， 馬 寧

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

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

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

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

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

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

不可壓縮流體的流動(dòng)服從如下控制方程組：

(1)

(2)

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

(3)

式中：Fc為對(duì)流項(xiàng)；FD為擴(kuò)散項(xiàng).

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

(4)

擴(kuò)散項(xiàng)為

(5)

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

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

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

(6)

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

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

FV-pSOR算法有4種基礎(chǔ)SOR迭代格式，如圖2所示.圖中，橫向代表第1個(gè)下標(biāo)i，垂向代表第2個(gè)下標(biāo)k，箭頭表示迭代計(jì)算的行進(jìn)方向.對(duì)于式(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é)點(diǎn)(i,k+1)處采用基礎(chǔ)格式A3，節(jié)點(diǎn)(i,k)處采用基礎(chǔ)格式A4，節(jié)點(diǎn)(i+1,k)處采用基礎(chǔ)格式A2，節(jié)點(diǎn)(i+1,k+1)處采用基礎(chǔ)格式A1.因?yàn)楫?dāng)前第n次迭代的值已知，所以只有4個(gè)節(jié)點(diǎn)處第n+1次迭代的值為

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

未知，可由以下4個(gè)基礎(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é)點(diǎn)(i,k)處采用基礎(chǔ)格式A3，節(jié)點(diǎn)(i+1,k)處采用基礎(chǔ)格式A1，對(duì)應(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 迭代計(jì)算過程

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

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

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

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

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的迭代公式依次計(jì)算各區(qū)域中的點(diǎn)在下一個(gè)迭代步的數(shù)值.

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

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

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

3 算例

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

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

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

圖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ù)量下的云圖分辨率相同.為了便于對(duì)比，圖9和10中同時(shí)給出了SOR算法和文獻(xiàn)[12]中的速度分布結(jié)果.可以看出，本文提出的FV-pSOR算法的結(jié)果與SOR算法的結(jié)果完全相同(差值為10-4量級(jí))，而且與文獻(xiàn)[12]中的結(jié)果相符，表明本文提出的FV-pSOR算法與SOR算法具有同等計(jì)算精度.由圖9和10可以看出，本文提出的FV-pSOR算法與SOR算法的計(jì)算精度差異不大.而由表1可見，其計(jì)算效率是SOR算法計(jì)算效率的6～7倍，說明FV-pSOR算法能夠適應(yīng)數(shù)值求解Navier-Stokes方程的五對(duì)角線性方程組的并行迭代求解.

表1 網(wǎng)格數(shù)量分別為50×50、100×100及150×150時(shí)的計(jì)算耗時(shí)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算的計(jì)算效率明顯高于SOR算法，但在細(xì)網(wǎng)格條件下其計(jì)算效率的提高幅度比粗網(wǎng)格下的有所降低.

4 結(jié)語

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