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

包蕓　葉豐　張義招

摘要： 對(duì)不可壓NS方程的數(shù)值計(jì)算，當(dāng)計(jì)算規(guī)模增大時(shí)，不論是采用湍流模型計(jì)算還是直接數(shù)值模擬（Direct Numerical Simulation，DNS），大規(guī)模的并行計(jì)算都難以實(shí)現(xiàn).該問題的關(guān)鍵是求解全場(chǎng)聯(lián)立的壓力泊松方程的并行計(jì)算技術(shù).利用并行近似解求解方案，創(chuàng)建高效大規(guī)模并行計(jì)算的不可壓NS方程的直接求解方法.三維窄方腔熱對(duì)流的DNS計(jì)算結(jié)果表明，該直接求解并行計(jì)算方法具有很好的并行效率，并且計(jì)算的三維湍流熱對(duì)流的特性是合理的.

關(guān)鍵詞： 泊松方程； 并行計(jì)算； 不可壓流動(dòng)； 湍流熱對(duì)流； 直接數(shù)值模擬

中圖分類號(hào)： O357.5文獻(xiàn)標(biāo)志碼： 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

基金項(xiàng)目： 國(guó)家自然科學(xué)基金（11372362）；中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金（14lgjc02）

作者簡(jiǎn)介： 包蕓（1960—），女，上海人，教授，博士，研究方向?yàn)橛?jì)算流體力學(xué)，（Email）stsby@mail.sysu.edu.cn0引言

在航空航天等高科技工程的推動(dòng)下，計(jì)算流體力學(xué)在可壓縮流動(dòng)的數(shù)值模擬計(jì)算技術(shù)領(lǐng)域進(jìn)步非凡.不可壓流動(dòng)的數(shù)值模擬技術(shù)也在不斷進(jìn)步.超級(jí)計(jì)算機(jī)硬件技術(shù)的快速發(fā)展為計(jì)算流體力學(xué)數(shù)值模擬的進(jìn)一步發(fā)展提供技術(shù)支持，高效并行計(jì)算技術(shù)的發(fā)展為進(jìn)一步擴(kuò)大不可壓NS方程的數(shù)值計(jì)算規(guī)模提供新的平臺(tái)，并使計(jì)算結(jié)果數(shù)據(jù)能更好地反映流體的流動(dòng)特性.

熱對(duì)流問題廣泛存在于自然界和工業(yè)界中，研究其對(duì)全球氣候變化、海洋環(huán)流、反應(yīng)堆設(shè)計(jì)、工業(yè)冷卻設(shè)計(jì)等問題的影響具有重要意義.[12]在Boussinesq假設(shè)下，湍流熱對(duì)流的描述方程為不可壓NS方程聯(lián)立溫度對(duì)流擴(kuò)散方程，因此熱對(duì)流問題也是典型的不可壓流動(dòng)問題.高瑞利數(shù)Ra的湍流RayleighBénard（RB）熱對(duì)流的直接數(shù)值模擬（Direct Numerical Simulation，DNS）面臨重大挑戰(zhàn).[3]隨著Ra的提高，熱對(duì)流進(jìn)入湍流狀態(tài)，DNS模擬的規(guī)模不斷增大導(dǎo)致計(jì)算所需要的成本巨大，數(shù)值計(jì)算難以實(shí)現(xiàn).目前，湍流熱對(duì)流的DNS只達(dá)到Ra=1012水平.[4]大規(guī)模可并行的高Ra湍流熱對(duì)流DNS及其海量數(shù)據(jù)結(jié)果分析已成為熱對(duì)流研究工作者們特別關(guān)注的問題.

在不可壓流動(dòng)NS方程的數(shù)值模擬計(jì)算中，不論采用何種計(jì)算模型或是DNS，其壓力泊松方程的求解都是大規(guī)模并行計(jì)算的難題：直接求解方法不易于并行.迭代求解壓力泊松方程通常采用局部算法從而較容易實(shí)現(xiàn)并行，但迭代計(jì)算過程占用大量的計(jì)算時(shí)間，所以很難實(shí)現(xiàn)高效率的大規(guī)模計(jì)算.這使得不可壓流動(dòng)的大規(guī)模數(shù)值模擬難以在有效時(shí)間內(nèi)滿足需求，因此妨礙大規(guī)模不可壓流動(dòng)NS方程的湍流模型計(jì)算或DNS的進(jìn)一步發(fā)展.一種新的泊松方程塊三對(duì)角近似求解方案[56]可解決在超級(jí)計(jì)算機(jī)上的高效并行直接求解問題.這使得建立不可壓流動(dòng)NS方程模擬的大規(guī)模高效并行計(jì)算方法成為可能，并可在超級(jí)計(jì)算機(jī)上實(shí)現(xiàn)三維高Ra湍流熱對(duì)流特性研究的DNS.

1控制方程

無量綱化后的三維不可壓NS方程為Δ·V=0

Vt+（V·Δ）V=-Δp+1ReΔ2V （1）式中：V為速度矢量；p為壓力；Re為雷諾數(shù).計(jì)算邊界條件為無滑移邊界.

2并行直接求解數(shù)值計(jì)算方法

投影法是數(shù)值求解不可壓NS方程組的有效方法之一.[7]實(shí)際上，不論采用何種湍流模型或DNS，以及采用何種思路求解不可壓NS方程的速度壓力法，一般的動(dòng)量方程都可以采用顯式格式連續(xù)方程推導(dǎo)出壓力泊松方程并進(jìn)行迭代求解.大規(guī)模并行計(jì)算的主要困難為壓力泊松方程必須全場(chǎng)聯(lián)立求解.本文主要針對(duì)矩形計(jì)算域的特定情況，結(jié)合有效的泊松方程高效并行的直接求解方案，創(chuàng)建不可壓流動(dòng)DNS的可高效并行求解計(jì)算方法.

2.1網(wǎng)格布置和離散格式

計(jì)算區(qū)域取矩形，見圖1.網(wǎng)格數(shù)為nx×ny×nz.在設(shè)計(jì)大規(guī)模并行計(jì)算時(shí)，消息傳遞接口（Message Passing Interface，MPI）的計(jì)算區(qū)域沿xOy平面對(duì)z方向分割，即圖中x方向較粗的線.在該面上并行計(jì)算需要數(shù)據(jù)通信；區(qū)域內(nèi)部用OpenMP并行，無須數(shù)據(jù)通信.由于直接求解壓力泊松方程要用到二維快速傅里葉離散余弦變換[8]，因此在x和y方向必須采用等距網(wǎng)格且最好網(wǎng)格數(shù)是2k，在z方向上可根據(jù)計(jì)算的需要采用非等距網(wǎng)格.

圖 1計(jì)算網(wǎng)格及并行計(jì)算的分割

Fig.1Computational mesh and division for parallel computing

本文采用不可壓流動(dòng)計(jì)算時(shí)常用的交錯(cuò)網(wǎng)格，時(shí)間方向采用一階精度離散，空間采用二階精度離散格式.

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

在不可壓NS方程的數(shù)值求解過程中，采用投影法將計(jì)算分步驟進(jìn)行.動(dòng)量方程中的速度計(jì)算采用顯式格式，在求解中很容易實(shí)現(xiàn)并行.由連續(xù)方程推導(dǎo)出的壓力泊松方程在求解時(shí)需要全場(chǎng)聯(lián)立，是求解過程中計(jì)算工作量最大的部分.同時(shí)，聯(lián)立求解也給并行造成困難，是大規(guī)模高效并行計(jì)算的難點(diǎn).高效合理并可大規(guī)模并行計(jì)算的壓力泊松方程求解方案，是解決不可壓NS方程大規(guī)模并行計(jì)算的關(guān)鍵技術(shù).

建立三維泊松方程的直接求解方法.計(jì)算方法只用于矩形計(jì)算區(qū)域，x方向采用等距網(wǎng)格.具體做法為在xOy平面上使用二維快速傅里葉變換將空間3個(gè)方向上都要求聯(lián)立求解的壓力泊松方程解耦，使泊松方程變換為只在z方向上的三對(duì)角方程.將三對(duì)角方程變換為塊三角方程，設(shè)計(jì)高效且可并行計(jì)算求解方法，求解后再使用對(duì)應(yīng)的反變換得到原來壓力泊松方程的全場(chǎng)解.

使用二維離散余弦傅里葉變換將全場(chǎng)聯(lián)立的泊松方程在x和y方向上解耦.余弦傅里葉變換能使壓力泊松方程自動(dòng)滿足壓力梯度為0的邊界條件.變換通過FFTW軟件包[8]實(shí)現(xiàn).二維離散余弦傅里葉變換公式為

將式（4）代入壓力泊松方程，并令展開式兩邊對(duì)應(yīng)系數(shù)相等，可以得到一組只沿z方向聯(lián)立的三對(duì)角方程，使壓力泊松方程在x和y方向上解耦，求解過程變簡(jiǎn)單.

在以往的二維熱對(duì)流DNS中，利用追趕法求解三對(duì)角方程的泊松方程直接解法[9]，與采用迭代求解方法在計(jì)算機(jī)上單線程計(jì)算相比有效得多，但三對(duì)角方程的追趕法很難進(jìn)行大規(guī)模并行計(jì)算.

數(shù)學(xué)和計(jì)算機(jī)的研究者們嘗試建立塊三對(duì)角方程的大規(guī)模高效并行求解方案[56]，將變換得來的三對(duì)角方程寫成塊三對(duì)角的形式為A= （5）式中：A=Ai，Bi，Ci為塊三對(duì)角矩陣，Ai，Ci1≤i≤nz和Bi0≤i≤nz為M×N階矩陣；和為M×N維列向量，=xT1，…，xTnT，=T1，…，TnT.為已知方程的右端項(xiàng)，通過式（4）求得.通過變換分解塊三對(duì)角方程，得到在MPI分塊中縮減的塊三對(duì)角方程，其中絕大部分主對(duì)角絕對(duì)占優(yōu)，可采用只需與相鄰分塊通信的近似求解方法，減少并行計(jì)算中的數(shù)據(jù)通信量.剩下的少部分主對(duì)角不能絕對(duì)占優(yōu)的三對(duì)角方程，求解過程則需全局通信.在z方向分塊內(nèi)的計(jì)算網(wǎng)格數(shù)目越大，主對(duì)角占優(yōu)的三對(duì)角方程的數(shù)量越大.塊三對(duì)角方程的高效并行近似求解方法，是實(shí)現(xiàn)高效并行計(jì)算不可壓流動(dòng)NS方程中壓力泊松方程直接求解方法的關(guān)鍵步驟.

在計(jì)算得到塊三對(duì)角方程的解后，通過對(duì)應(yīng)的二維離散余弦傅里葉的反變換公式

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

可得到全場(chǎng)的壓力pn+1，完成泊松方程直接求解.

利用以上高效并行三維壓力泊松方程直接求解方法，聯(lián)合其他方便并行的動(dòng)量方程等計(jì)算，創(chuàng)建三維不可壓NS方程高效并行直接求解計(jì)算方法，使得在一些特定情況下，大規(guī)模高效并行的不可壓流動(dòng)NS方程湍流模型計(jì)算或者DNS成為可能.

3三維湍流熱對(duì)流大規(guī)模并行計(jì)算效率與計(jì)算結(jié)果3.1三維湍流熱對(duì)流方程

RB熱對(duì)流是研究流體對(duì)流傳熱的典型物理模型.在封閉的盒子內(nèi)，下底板加熱而上底板冷卻后形成對(duì)流傳熱的研究系統(tǒng).在Boussinesq假設(shè)下，無量綱化后的三維熱對(duì)流的描述方程為Δ·V=0

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

θt+（V·Δ）θ=1RaPrΔ2θ （7）式中：Ra為瑞利數(shù)；Pr為普朗特?cái)?shù)；θ為相對(duì)溫度，下底板為加熱，θ=0.5，上底板為冷卻，θ=-0.5.

通過熱對(duì)流方程組可以看出，整個(gè)計(jì)算過程實(shí)際上就是數(shù)值求解不可壓NS方程組聯(lián)立溫度的對(duì)流擴(kuò)散方程.本文對(duì)三維湍流熱對(duì)流進(jìn)行DNS.

3.2并行計(jì)算效率

采用本文建立的三維不可壓流動(dòng)的直接求解并行計(jì)算方法，在超級(jí)計(jì)算機(jī)“天河二號(hào)”上進(jìn)行加速比測(cè)試.每個(gè)計(jì)算節(jié)點(diǎn)包含24個(gè)計(jì)算物理核心.測(cè)試算例的計(jì)算網(wǎng)格數(shù)和物理計(jì)算核心數(shù)見表1.表 1測(cè)試算例

Tab.1Test casesnx×ny×nz節(jié)點(diǎn)數(shù)/個(gè)核心數(shù)/個(gè)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