不可压NS方程的高效并行直接求解

包芸　叶丰　张义招

摘要： 对不可压NS方程的数值计算，当计算规模增大时，不论是采用湍流模型计算还是直接数值模拟（Direct Numerical Simulation，DNS），大规模的并行计算都难以实现.该问题的关键是求解全场联立的压力泊松方程的并行计算技术.利用并行近似解求解方案，创建高效大规模并行计算的不可压NS方程的直接求解方法.三维窄方腔热对流的DNS计算结果表明，该直接求解并行计算方法具有很好的并行效率，并且计算的三维湍流热对流的特性是合理的.

关键词： 泊松方程； 并行计算； 不可压流动； 湍流热对流； 直接数值模拟

中图分类号： O357.5文献标志码： B

Efficient parallel direct solution on

incompressible NS equations

BAO Yun， YE Feng， ZHANG Yizhao

（Department of Mechanics， Sun Yetsen University， Guangzhou 510275， China）

Abstract： The largescale parallel computational techniques for incompressible NS equations solution is difficult to realize either for the turbulent models or the Direct Numerical Simulation （DNS） when the computational scale increases. It is key to implement the parallel computational technique for the pressure Poisson equation which needs the solutions for the whole flow field. An efficient parallel direct solution scheme for incompressible NS equations is developed using a largescale parallel approximate solution. The DNS calculation results of the heat convection in a narrow square cavity show that the parallel direct solution scheme has good parallel efficiency， and the characteristics of 3D turbulent heat convection obtained by the new scheme is reasonable.

Key words： Poisson equation； parallel calculation； incompressible flow； turbulent convection； direct numerical simulation

收稿日期： 2015[KG*9〗09[KG*9〗21修回日期： 2015[KG*9〗12[KG*9〗08

基金项目： 国家自然科学基金（11372362）；中央高校基本科研业务费专项资金（14lgjc02）

作者简介： 包芸（1960—），女，上海人，教授，博士，研究方向为计算流体力学，（Email）stsby@mail.sysu.edu.cn0引言

在航空航天等高科技工程的推动下，计算流体力学在可压缩流动的数值模拟计算技术领域进步非凡.不可压流动的数值模拟技术也在不断进步.超级计算机硬件技术的快速发展为计算流体力学数值模拟的进一步发展提供技术支持，高效并行计算技术的发展为进一步扩大不可压NS方程的数值计算规模提供新的平台，并使计算结果数据能更好地反映流体的流动特性.

热对流问题广泛存在于自然界和工业界中，研究其对全球气候变化、海洋环流、反应堆设计、工业冷却设计等问题的影响具有重要意义.[12]在Boussinesq假设下，湍流热对流的描述方程为不可压NS方程联立温度对流扩散方程，因此热对流问题也是典型的不可压流动问题.高瑞利数Ra的湍流RayleighBénard（RB）热对流的直接数值模拟（Direct Numerical Simulation，DNS）面临重大挑战.[3]随着Ra的提高，热对流进入湍流状态，DNS模拟的规模不断增大导致计算所需要的成本巨大，数值计算难以实现.目前，湍流热对流的DNS只达到Ra=1012水平.[4]大规模可并行的高Ra湍流热对流DNS及其海量数据结果分析已成为热对流研究工作者们特别关注的问题.

在不可压流动NS方程的数值模拟计算中，不论采用何种计算模型或是DNS，其压力泊松方程的求解都是大规模并行计算的难题：直接求解方法不易于并行.迭代求解压力泊松方程通常采用局部算法从而较容易实现并行，但迭代计算过程占用大量的计算时间，所以很难实现高效率的大规模计算.这使得不可压流动的大规模数值模拟难以在有效时间内满足需求，因此妨碍大规模不可压流动NS方程的湍流模型计算或DNS的进一步发展.一种新的泊松方程块三对角近似求解方案[56]可解决在超级计算机上的高效并行直接求解问题.这使得建立不可压流动NS方程模拟的大规模高效并行计算方法成为可能，并可在超级计算机上实现三维高Ra湍流热对流特性研究的DNS.

1控制方程

无量纲化后的三维不可压NS方程为Δ·V=0

Vt+（V·Δ）V=-Δp+1ReΔ2V （1）式中：V为速度矢量；p为压力；Re为雷诺数.计算边界条件为无滑移边界.

2并行直接求解数值计算方法

投影法是数值求解不可压NS方程组的有效方法之一.[7]实际上，不论采用何种湍流模型或DNS，以及采用何种思路求解不可压NS方程的速度压力法，一般的动量方程都可以采用显式格式连续方程推导出压力泊松方程并进行迭代求解.大规模并行计算的主要困难为压力泊松方程必须全场联立求解.本文主要针对矩形计算域的特定情况，结合有效的泊松方程高效并行的直接求解方案，创建不可压流动DNS的可高效并行求解计算方法.

2.1网格布置和离散格式

计算区域取矩形，见图1.网格数为nx×ny×nz.在设计大规模并行计算时，消息传递接口（Message Passing Interface，MPI）的计算区域沿xOy平面对z方向分割，即图中x方向较粗的线.在该面上并行计算需要数据通信；区域内部用OpenMP并行，无须数据通信.由于直接求解压力泊松方程要用到二维快速傅里叶离散余弦变换[8]，因此在x和y方向必须采用等距网格且最好网格数是2k，在z方向上可根据计算的需要采用非等距网格.

图 1计算网格及并行计算的分割

Fig.1Computational mesh and division for parallel computing

本文采用不可压流动计算时常用的交错网格，时间方向采用一阶精度离散，空间采用二阶精度离散格式.

2.2不可压NS方程的并行求解方案

在不可压NS方程的数值求解过程中，采用投影法将计算分步骤进行.动量方程中的速度计算采用显式格式，在求解中很容易实现并行.由连续方程推导出的压力泊松方程在求解时需要全场联立，是求解过程中计算工作量最大的部分.同时，联立求解也给并行造成困难，是大规模高效并行计算的难点.高效合理并可大规模并行计算的压力泊松方程求解方案，是解决不可压NS方程大规模并行计算的关键技术.

建立三维泊松方程的直接求解方法.计算方法只用于矩形计算区域，x方向采用等距网格.具体做法为在xOy平面上使用二维快速傅里叶变换将空间3个方向上都要求联立求解的压力泊松方程解耦，使泊松方程变换为只在z方向上的三对角方程.将三对角方程变换为块三角方程，设计高效且可并行计算求解方法，求解后再使用对应的反变换得到原来压力泊松方程的全场解.

使用二维离散余弦傅里叶变换将全场联立的泊松方程在x和y方向上解耦.余弦傅里叶变换能使压力泊松方程自动满足压力梯度为0的边界条件.变换通过FFTW软件包[8]实现.二维离散余弦傅里叶变换公式为

将式（4）代入压力泊松方程，并令展开式两边对应系数相等，可以得到一组只沿z方向联立的三对角方程，使压力泊松方程在x和y方向上解耦，求解过程变简单.

在以往的二维热对流DNS中，利用追赶法求解三对角方程的泊松方程直接解法[9]，与采用迭代求解方法在计算机上单线程计算相比有效得多，但三对角方程的追赶法很难进行大规模并行计算.

数学和计算机的研究者们尝试建立块三对角方程的大规模高效并行求解方案[56]，将变换得来的三对角方程写成块三对角的形式为A= （5）式中：A=Ai，Bi，Ci为块三对角矩阵，Ai，Ci1≤i≤nz和Bi0≤i≤nz为M×N阶矩阵；和为M×N维列向量，=xT1，…，xTnT，=T1，…，TnT.为已知方程的右端项，通过式（4）求得.通过变换分解块三对角方程，得到在MPI分块中缩减的块三对角方程，其中绝大部分主对角绝对占优，可采用只需与相邻分块通信的近似求解方法，减少并行计算中的数据通信量.剩下的少部分主对角不能绝对占优的三对角方程，求解过程则需全局通信.在z方向分块内的计算网格数目越大，主对角占优的三对角方程的数量越大.块三对角方程的高效并行近似求解方法，是实现高效并行计算不可压流动NS方程中压力泊松方程直接求解方法的关键步骤.

在计算得到块三对角方程的解后，通过对应的二维离散余弦傅里叶的反变换公式

pi，j，k=4MNMu=0Nv=0αuβvu，v，kcosπuiMcosπvjN （6）

可得到全场的压力pn+1，完成泊松方程直接求解.

利用以上高效并行三维压力泊松方程直接求解方法，联合其他方便并行的动量方程等计算，创建三维不可压NS方程高效并行直接求解计算方法，使得在一些特定情况下，大规模高效并行的不可压流动NS方程湍流模型计算或者DNS成为可能.

3三维湍流热对流大规模并行计算效率与计算结果3.1三维湍流热对流方程

RB热对流是研究流体对流传热的典型物理模型.在封闭的盒子内，下底板加热而上底板冷却后形成对流传热的研究系统.在Boussinesq假设下，无量纲化后的三维热对流的描述方程为Δ·V=0

Vt+（V·Δ）V=-Δp+1Ra/PrΔ2V+θ

θt+（V·Δ）θ=1RaPrΔ2θ （7）式中：Ra为瑞利数；Pr为普朗特数；θ为相对温度，下底板为加热，θ=0.5，上底板为冷却，θ=-0.5.

通过热对流方程组可以看出，整个计算过程实际上就是数值求解不可压NS方程组联立温度的对流扩散方程.本文对三维湍流热对流进行DNS.

3.2并行计算效率

采用本文建立的三维不可压流动的直接求解并行计算方法，在超级计算机“天河二号”上进行加速比测试.每个计算节点包含24个计算物理核心.测试算例的计算网格数和物理计算核心数见表1.表 1测试算例

Tab.1Test casesnx×ny×nz节点数/个核心数/个1 024×128×1 53644×241 024×128×1 53688×241 024×128×1 5361616×241 024×128×1 5363232×241 024×128×1 5366464×24