不可压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

不同计算核数时直接求解并行计算方法的加速比见图2，其中加速比以计算节点24核心数为基数.计算中z方向上网格数1 536是24核心数的64倍.可知当使用32节点即32×24=768核心数时，具有约81.7%的计算效率，当使用64节点即64×24=1 536核心数时，具有约67%的计算效率.加速比随计算机核数的进一步增加仍有较好的可增长性.由于采用快速傅里叶变换解耦压力泊松方程的需要， 并行区域只在z方向上划分，由此本文创建的直接求解方法在大规模并行计算时z方向网格数与计算核心数之间有一定关联.总体上讲，本文的不可压流动直接求解方法在大规模并行计算上已得到较好的计算效率，可以进行较大规模的三维湍流热对流的DNS.新的直接求解计算方法与本课题组原有的迭代计算方法相比，节省1/2以上的总计算时间，计算技术的进步为系列三维热对流的DNS及物理特性的研究提供有力工具.

3.3三维湍流热对流计算结果

选取窄方腔，宽高比Γ=1/4，Pr=4.3，对不同Ra的三维热对流进行计算.低Ra计算采用512×64×768网格，高Ra计算采用1 024×128×1 200网格.首先计算流动的初场，根据传热特性计算时平均场统计的需要，当热对流中出现较稳定的大尺度环流流动后，继续迭代计算至少100万时间步.根据不同的Ra及流动特性，计算时间的迭代至少达到300万时间步.

RB热对流主要探讨由浮力产生的对流运动对流体传热特性的影响.Ra=1010的平均场温度等值面分布见图3，下底板加热为高温，上底板冷却为低温.温度作为热对流中的重要物理量，其平均场的分布反映传热过程中对流流动带动温度运动的情况.图3中方腔的中间部分没有温度的分布，表明中间没有流动也没有传热作用.图3显示由大尺度环流带动的温度分布形态，RB热对流中对流传热的热量主要沿着侧壁和棱边向上传输.

RB热对流研究的核心问题之一是对流传热效率，表征传热的物理参数是努塞尔数Nu，表示流体对流传热与热传导传热的比值.Nu与Ra之间存在一定的标度关系[3]，因此需要一系列Ra数的不同数值模拟结果进行研究.

图4中同时也给出计算得到的二维热对流Nu数的变化结果[10]，以及KACZOROWSKI等[11]的三维计算结果对比.由此可见，本文的三维湍流热对流的DNS模拟结果是合理的.图4中可以看到，二维和三维热对流的计算Nu随Ra的变化都存在很好的标度率.计算结果与理论预测和大量试验结果得到的结论一致.[3]Γ=1/4三维方腔的Nu总体偏大（向上平移）.在试验和计算中也发现同样的Nu平移现象，并且不同的Γ变化会引起不同的Nu向上平移量.[12]在传热Nu对Ra的标度率中，三维和二维流动计算结果产生差异的物理原因，还有待更深入的研究.

4结束语

在三维不可压流动特性的研究过程中，尤其是到湍流阶段，超大规模的数值模拟计算十分必要.依靠超级计算机技术，探索高效并行的计算方法和计算技术，进行大规模的高自由度三维不可压NS方程的数值计算，对更深入研究不可压流动的物理和流动特性具有很重要的意义.大规模并行数值模拟计算三维不可压NS方程，最难的计算方法和计算技术问题是压力泊松方程的高效并行求解.利用块三对角方程的大规模高效并行近似解求解方案，创建大规模高效并行计算三维不可压NS方程的直接求解方法.

DNS是研究高Ra湍流热对流的重要手段之一.热对流问题是典型的不可压流动问题.利用本文创建的高效并行不可压流动的直接求解方法，对高Ra三维湍流热对流进行DNS.通过并行效率的验证计算，证明本文建立的直接求解方法具有较好的并行效率.一系列不同Ra的三维窄方腔热对流计算得到的传热Nu具有合理的标度率.

本文创建的大规模高效并行计算的直接求解方法，为高Ra的湍流热对流大规模高效并行计算和数值模拟研究提供有价值的计算参考.参考文献：

[1]BERNARD C， GEMUNU G， FRANOIS H， et al. Scaling of hard thermal turbulence in RayleighBénard convection[J]. Journal of Fluid Mechanics， 1989， 204： 130. DOI： 10.1017/S0022112089001643.

[2]HE X， FUNFSCHILLING D， NOBACH H， et al. Transition to the ultimate state of turbulent rayleighbénard convection[J]. Physical Review Letters， 2012， 108： 024502. DOI： 10.1103/PhysRevLett.108.024502.

[3]AHLERS G， GROSSMANN S， LOHSE D. Heat transfer and large scale dynamics in turbulent RayleighBénard convection[J]. Reviews of Modern Physics， 2009， 81： 503537. DOI： 10.1103/RevModPhys.81.503.

[4]STEVENS R J A M， VERZICCO R， LOHSE D. Radial boundary layer structure and Nusselt number in RayleighBenard convection[J]. Journal of Fluid Mechanics， 2010， 643： 495507. DOI： 10.1017/S0022112009992461.

[5]ZHANG H， ZHANG W， SUN X H. a parallel algorithm for block tridiagonal systems[C]//Proc 2008 9th International Conference on Parallel and Distributed Computing， Applications and Technologies， 2008： 6265. DOI： 10.1109/PDCAT.2008.21.

[6]张武， 张衡. 初边值问题的块三对角可扩展并行算法[J]. 上海大学学报（自然科学版）， 2007， 13（5）： 497503.

ZHANG W， ZHANG H. Scalable parallel algorithm of block tridiagonal systems for initial boundary value problem[J]. Journal of Shanghai University（Natural Science）， 2007， 13（5）： 497503.

[7]ALEXANDRE J C. A numerical method for solving incompressible viscous flow problems[J]. Journal of Computation Physics， 1997， 2（1）： 1226. DOI： 10.1016/00219991（67）90037X.

[8]MATTEO F， STEVEN G J. The design and implementation of FFTW3[J]. Proceeding of IEEE， 2005， 93（2）： 216231. DOI： 10.1109/JPROC.2004.840301.

[9]徐炜， 包芸. FFT高效求解二维RayleighBénard热对流[J]. 力学学报， 2013， 45（4）： 16. DOI： 10.6052/0459187912334.

XU W， BAO Y. An efficient solution for 2D RayleighBénard convection using FFT[J]. Chinese Journal of Theoretical and Applied Mechanics， 2013， 45（4）： 16. DOI： 10.6052/0459187912334.

[10]宁浩. 维RayleighBénard热对流数值计算及大尺度环流特性[D]. 广州： 中山大学， 2013.

[11]KACZOROWSKI M， CHONG K L， XIA K Q. Turbulent flow in the bulk of RayleighBénard convection：aspectratio dependence of the smallscale properties[J]. Journal of Fluid Mechanics， 2014（747）： 73102. DOI： 10.1017/jfm.2014.154.

[12]HUANG S D， KACZOROWSKI M， NI R， et al. Confinementinduced heattransport enhancement in turbulent thermal convection[J]. Physical Review Letters， 2013（111）： 104501. DOI： 10.1103/PhysRevLett.111. 104501. （编辑于杰）