一种不可压缩二维流动的显式逐次超松弛并行算法

张晓慧， 柏君励， 顾解忡， 马 宁

(上海交通大学 海洋工程国家重点实验室； 高新船舶与深海开发装备协同创新中心；船舶海洋与建筑工程学院， 上海 200240)

计算流体力学方法已逐渐成为流动及流体传热问题研究的主流工具.研究表明[1-4]，与势流方法相比，采用计算流体力学方法能够更准确地预报船舶的耐波性.其中，主要方法是采用有限体积法(FVM)离散Navier-Stokes方程而形成代数方程组，再进行迭代求解.然而，采用计算流体力学方法求解一个算例需要耗费大量时间，尤其是迭代计算.因为每个时间层涉及多次外迭代，每次外迭代又涉及速度、压力或压力修正值的内迭代，所以提高迭代计算的效率已经成为计算流体力学的一个关键技术问题，而实行并行计算是提高计算效率的一条行之有效的途径.随着分组显式思想和区域分解方法[5]的提出，为并行计算创造了条件.目前，并行求解线性代数方程组的迭代方法主要包括Schwarz法、多重网格法、共轭梯度(CG)法以及经典的定常迭代法等.Schwarz法基于区域分解的思想，易于并行化但计算的收敛速度较慢.例如，文献[6]中通过乘性Schwarz方法以及预处理过程来加速迭代的收敛速度，但其很难程序化.多重网格方法[7]是通过将计算过程转换到粗网格中以加快迭代收敛速度.另外，Krylov子空间迭代法[8]包含了共轭梯度法及其扩展方法，是另一类较为重要的迭代方法.考虑到计算的效率和可操作性，本文选取经典的定常迭代法.虽然Jacobi迭代法容易实行并行操作，但其比Gauss-Seidel和逐次超松弛(SOR)算法的收敛速度慢.而相比于Gauss-Seidel算法，SOR算法的收敛速度更快.文献[9]中基于区域分解思想并结合SOR方法对Stokes方程进行并行计算.

最近，许秋燕等[10]结合区域分解和分组显式方法将计算域分割成多个子域，构造了子域并行迭代计算格式，以用于有限差分求解特定的泊松方程，最终求解的是常系数代数方程组.而对于流动以及流体传热问题，代数方程组的系数将随着外迭代和时间层的不断推进而不断变化.

本文采用船舶水动力学研究中常用的FVM求解不可压缩Navier-Stokes方程，从求解二维流体流动入手，在许秋燕等[10]研究的基础上，提出一种基于FVM的显式逐次超松弛并行算法(FV-pSOR).其中，采用4阶紧致格式离散流项和扩散项、4阶龙格-库塔法步进非定常项以及同位网格SIMPLE(Semi-Implicit Method for Pressure-Linked Equation)算法求解压力-速度耦合方程[11].基于区域分解，将计算域分割成4个子域而构造所需分组显式SOR迭代公式，并通过3种网格密度下典型二维驱动方腔流的计算来验证所提FV-pSOR算法的有效性.

1 数学模型

1.1 Navier-Stokes方程组及其数值离散

不可压缩流体的流动服从如下控制方程组：

(1)

(2)

本文并不是直接求解微分形式的控制方程组(式(1)、(2))，而是求解在同位网格中由FVM导出的积分形式的控制方程[11].同位网格是将速度、压力、黏度等物理量均定义在网格中心点处，如图1所示.其中，网格P的积分方程为

(3)

式中：Fc为对流项；FD为扩散项.

对流项可用前一时刻的速度进行线性化，即

(4)

扩散项为

(5)

图1 同位网格Fig.1 Collocated grid

式(4)、(5)中的变量包括速度和速度导数，采用4阶紧致格式插值[11].式(3)中时间层的步进采用4阶龙格-库塔法[11]，龙格-库塔法的每一步都采用 SIMPLE 算法求解速度-压力耦合方程.

经上述离散后，最终所得压力修正方程为

(6)

下面将具体说明式(6)的显式SOR并行迭代算法的构造过程.

1.2 FV-pSOR算法的基础迭代格式

FV-pSOR算法有4种基础SOR迭代格式，如图2所示.图中，横向代表第1个下标i，垂向代表第2个下标k，箭头表示迭代计算的行进方向.对于式(6)，4种基础SOR迭代格式的具体表达式如下：

基础迭代格式A1为

(7)

基础迭代格式A2为

(8)

基础迭代格式A3为

(9)

基础迭代格式A4为

(10)

图2 4种基础迭代格式Fig.2 Four iterative schemes

根据分组显式的思想，由上述4种基础SOR迭代格式构造的迭代格式A5见图3.在节点(i,k+1)处采用基础格式A3，节点(i,k)处采用基础格式A4，节点(i+1,k)处采用基础格式A2，节点(i+1,k+1)处采用基础格式A1.因为当前第n次迭代的值已知，所以只有4个节点处第n+1次迭代的值为

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

未知，可由以下4个基础格式构造，即

(11)

联立求解可得

从而，由式(11)所得迭代格式A5的显式表达式为

(12)

式中：

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

b6b8+b2(b4+b3b5b8-b4b6b8)

同样地，由基础SOR格式可以构造图4所示的另外4种基础迭代格式.鉴于这4种基础迭代格式的构造方法类似，本文仅说明格式A6的构造过程.在节点(i,k)处采用基础格式A3，节点(i+1,k)处采用基础格式A1，对应的迭代方程组为

(13)

式中：

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

求解式(13)，所得格式A6的显式表达式为

(14)

同样地，可以得到其他3种迭代格式.其中，格式A7为

(15)

式中：

格式A8为

(16)

式中：

格式A9为

(17)

式中：

至此，所需显式SOR迭代格式全部构造完成.

图5 偶次迭代行进过程示意图Fig.5 Sketch of the process of even iteration

2 迭代计算过程

图5和6示出了偶次迭代和奇次迭代的行进顺序.图中，计算域被分割成4个子域.由图5可见，偶次迭代顺序是顺次从中心向外行进的.其迭代步骤如下：

图6 奇次迭代行进过程示意图Fig.6 Sketch of the process of odd iteration

(1) 选定两个坐标方向的中点xp、xp+1和zq、zq+1，共有4个位置点，区域Λ由4个选定的点构成.按照格式A5的迭代公式算出这4个点在下一个迭代步的数值.

(2) 根据箭头的行进方向，分别在区域Γ1、Γ2、Γ3和Γ4中按照格式A6、A7、A8和A9的迭代公式依次计算当前区域中的点在下一个迭代步的数值.以Γ1为例，计算次序为

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

其中：kmax为k方向的最大网格数.

(3) 在区域Ω1、Ω2、Ω3和Ω4中分别按照格式A1、A2、A3和A4的迭代公式依次计算各区域中的点在下一个迭代步的数值.

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

如图6所示，奇次迭代的迭代顺序是从四角向中心行进.在区域Ω1、Ω2、Ω3和Ω4中分别按照格式A4、A3、A2和A1的迭代公式依次计算各区域中的点在下一个迭代步的数值.迭代完成后，所有点处的值均已得到更新.

至此，偶次迭代和奇次迭代的过程构建完成.计算时，奇次迭代和偶次迭代交替进行，直至精度满足控制要求，然后，进入下一个外迭代层或下一个时间层.

3 算例

以二维驱动方腔流为例来验证FV-pSOR算法的有效性.图7所示为二维方腔流在左右侧面和底面需满足的固壁无滑移条件，其顶面为给定的自左向右均匀剪切流动.水平速度设为u=1 m/s，垂向速度为v=0.

图7 二维方腔流边界条件Fig.7 Boundary conditions of 2D lid-driven cavity flow

取流动区域的边长为h=1 m，流动的雷诺数Re=uh/ν=1 000.计算时，将矩形计算域划分为均匀网格，网格数量N分别为50×50、100×100以及150×150.并行计算时，将矩形计算域分割成4个子域.本文分析方腔流发展到35 s左右时的计算结果.采用一台处理器为Intel酷睿i7-6700的计算机，调用DATE_AND_TIME函数所获两种算法在3种网格数量下的计算耗时见表1.可见：在网格数量N=50×50下，SOR算法的耗时超出FV-pSOR算法耗时12倍；在网格分别为100×100和150×150下，SOR算法的耗时分别超出FV-pSOR算法的耗时6和7倍.可见，对于相同的子域分割，采用细网格有利于提高计算效率.

图8～10为网格数量分别为50×50、100×100以及150×150下的水平速度分布云图、垂向速度v沿y=0.5 m 水平线的分布图、水平速度u沿x=0.5 m 垂直线的分布图.其中，3种网格数量下的云图分辨率相同.为了便于对比，图9和10中同时给出了SOR算法和文献[12]中的速度分布结果.可以看出，本文提出的FV-pSOR算法的结果与SOR算法的结果完全相同(差值为10-4量级)，而且与文献[12]中的结果相符，表明本文提出的FV-pSOR算法与SOR算法具有同等计算精度.由图9和10可以看出，本文提出的FV-pSOR算法与SOR算法的计算精度差异不大.而由表1可见，其计算效率是SOR算法计算效率的6～7倍，说明FV-pSOR算法能够适应数值求解Navier-Stokes方程的五对角线性方程组的并行迭代求解.

表1 网格数量分别为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算法，但在细网格条件下其计算效率的提高幅度比粗网格下的有所降低.

4 结语

本文在有限差分并行迭代算法的基础上提出了FV-pSOR算法，用于提高SOR方法求解Navier-Stokes方程的计算效率，并通过一个典型的二维驱动方腔流算例的计算，验证了FV-pSOR算法的有效性.结果表明，FV-pSOR算法的计算效率较高，且其计算精度与SOR算法差异不大.虽然FV-pSOR算法的计算格式是基于计算流体力学流动问题的数值求解提出的，但其同样适用于其他五对角线性方程组的并行迭代求解.