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

包蕓　葉豐　張義招

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

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

中圖分類號： 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

基金項目： 國家自然科學(xué)基金（11372362）；中央高校基本科研業(yè)務(wù)費(fèi)專項資金（14lgjc02）

作者簡介： 包蕓（1960—），女，上海人，教授，博士，研究方向為計算流體力學(xué)，（Email）stsby@mail.sysu.edu.cn0引言

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

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

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

1控制方程

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

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

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

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

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

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

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

Fig.1Computational mesh and division for parallel computing

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.2并行計算效率

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

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