李 剛，孫東亮，齊亞強，陳 帥，宇 波

(1.北京化工大學(xué)機電工程學(xué)院，北京 100029；2.北京石油化工學(xué)院機械工程學(xué)院，北京 102617)

流動傳熱是一種非常復(fù)雜的物理現(xiàn)象，廣泛存在于石油、化工、冶金、能源等工業(yè)中，如化學(xué)、制藥工業(yè)中的蒸發(fā)器和反應(yīng)器；制冷工業(yè)中的冷凝器；核電站的蒸汽發(fā)生器、換熱器；海上平臺石油開采和精煉等都存在流動傳熱現(xiàn)象。目前，針對流動傳熱現(xiàn)象研究的方式主要采用試驗的方法，但是該方法受到投資大、周期長、優(yōu)化困難等一系列缺點的限制，難以滿足工程領(lǐng)域不斷增長的需求[1]。隨著科學(xué)技術(shù)的創(chuàng)新、計算機技術(shù)的進步以及現(xiàn)代流動傳熱計算方法的不斷完善，計算流體力學(xué)(CFD)成為研究流體力學(xué)各種物理現(xiàn)象的重要手段。流動傳熱的數(shù)值模擬已經(jīng)取代了一些昂貴的、危險的、難以實現(xiàn)的試驗，極大地降低了研制成本，縮短了研制周期。數(shù)值模擬方法在研究流動傳熱問題方面得到廣泛地應(yīng)用。然而，隨著所涉及的流動傳熱問題復(fù)雜程度不斷提高，方程數(shù)量變多，方程求解的難度變大，計算耗時大幅增加。如何快速求解這些復(fù)雜的流動傳熱問題成為學(xué)者研究的熱點。

近年來，圖像處理器(graphics processing unit，GPU)強大的浮點運算能力逐漸得到開發(fā)，與傳統(tǒng)數(shù)值計算所用的CPU相比，GPU有著計算能力強、精度高、并行能力強以及可移植等特點。為了快速求解這些復(fù)雜的流動傳熱問題，國內(nèi)外學(xué)者結(jié)合高效的GPU并行計算，完成了部分基于GPU并行計算的求解算法的研究。早在2008年，Brandvik等[2]利用GPU實現(xiàn)了三維Euler方程的快速求解，獲得了16倍的加速效果，該研究最初驗證了GPU應(yīng)用在CFD領(lǐng)域的有效性。Williams等[3]研究出了一種基于GPU共軛梯度求解器，專門求解大型的線性方程稀疏系統(tǒng)，進而減少CFD分析時間。Jin等[4]開發(fā)了一個基于GPU的三維非結(jié)構(gòu)化幾何多網(wǎng)格解算器，有效地解決了非結(jié)構(gòu)化三維網(wǎng)格的熱傳導(dǎo)問題。在160萬的網(wǎng)格規(guī)模下，獲得了24倍的加速倍數(shù)。彭峻等[5]通過改進WENO算法使其適合GPU加速，實現(xiàn)了高精度可壓縮流的GPU計算，并獲得了18倍的加速效果。鞠鵬飛等[6]為了實現(xiàn)隱式LUSGS時間推進，采用層內(nèi)GPU并行方法得到的加速倍數(shù)為8.89倍。

在流動傳熱問題的數(shù)值模擬中，最重要的部分是對Navier-Stokes(N-S)方程組的求解。非原始變量法和原始變量法是2種最常用的求解算法。在非原始變量法中最具代表性的方法為渦量-流函數(shù)法[7]。該類方法求解方程組時要引入新的變量，所求的結(jié)果并不是速度和壓力場，邊界處理較為復(fù)雜，應(yīng)用相對較少。原始變量方法的最大優(yōu)點是不僅邊界處理相對簡單，而且可以直接計算出所需要的壓力場和速度場，所以應(yīng)用較為廣泛，發(fā)展也相對較快。針對非穩(wěn)態(tài)流動傳熱的數(shù)值模擬，在求解N-S方程組原始變量的主要算法中投影算法應(yīng)用較為成熟，該方法被普遍應(yīng)用于有關(guān)CFD問題的求解。但是隨著非穩(wěn)態(tài)流動傳熱問題研究的深入，如針對復(fù)雜區(qū)域的三維問題，網(wǎng)格數(shù)量通常高達百萬或千萬，通常需要幾天甚至幾個月的時間才能獲得一個收斂結(jié)果，遠遠不能適應(yīng)工程上快速設(shè)計的要求。為了實現(xiàn)非穩(wěn)態(tài)流動傳熱問題的快速求解，筆者開發(fā)了基于投影法的GPU并行計算策略，并以二維頂蓋驅(qū)動流為例研究了GPU并行計算求解性能。

1 控制方程及離散

1.1 控制方程無量綱化

針對二維頂蓋驅(qū)動不可壓縮流，以方腔的長度和高度H、流體密度ρ、頂蓋拖動速度utop作為無量綱標尺，將原有N-S控制方程組無量綱化[8]：

質(zhì)量守恒方程：

(1)

動量守恒方程：

(2)

(3)

1.2 方程離散

對于投影算法，中心差分格式被用于離散方程中的導(dǎo)數(shù)[9-10]。針對壓力梯度項采用隱式格式，速度項采用顯示格式，Δτ為時間步長，ΔX、ΔY為空間步長。

對速度分量U進行離散：

對于內(nèi)點處理：

非穩(wěn)態(tài)項：

(4)

對流項：

(5)

(6)

擴散項：

(7)

(8)

壓力項：

(9)

(10)

(11)

離散化后的動量方程：

(12)

同理可得：

(13)

連續(xù)性方程離散：

(14)

(15)

其中：

2 投影法

為了快速準確地求解N-S方程組，投影法實現(xiàn)的方式是對該方程組的速度壓力進行解耦求解，且速度和壓力解耦后無需要求速度-壓力的插值函數(shù)滿足Ladyzhenskaya-Babuska-Brezzi(LBB)相容性條件[11]。該條件目的為實現(xiàn)速度-壓力混合有限元求解的穩(wěn)定，但是實現(xiàn)該條件的求解過程較為復(fù)雜。而投影法的主要思想為在原有方程中引入一個無需要滿足散度為零的中間變量，分開計算速度和壓力，從而無需滿足該條件。

(16)

其中：η滿足·η=0，因此投影法實施過程中通過在N-S方程中引入1個無散度適量，該變量可以表示為1個無旋(有勢)矢量場與無散度速度場之和：

(17)

首先，通過一階離散的方式將N-S方程中的時間導(dǎo)數(shù)項進行離散，得到1個半離散化方程：

(18)

(19)

(20)

(21)

式(21)的右端項由步驟①計算，而壓力pn+1項在給定的邊界條件下能夠求出。

(22)

總之，當利用投影法求解下一時刻的壓力場和速度場時，動量方程和連續(xù)性方程不需要重復(fù)迭代，只需要通過上述3步就能完成求解。

3 CUDA編程模型

CUDA是一個用于CPU+GPU異構(gòu)系統(tǒng)的通用并行計算架構(gòu)。該架構(gòu)的最大優(yōu)勢是降低了GPU通用計算的編程門檻，研究人員不需要掌握復(fù)雜的圖形知識就可以直接使用C/Fortran等語言為CUDA架構(gòu)編寫程序，因此該架構(gòu)在計算流體力學(xué)(CFD)并行計算方面得到廣泛地應(yīng)用[13-14]。筆者采用NVIDIA的CUDA編程模型實現(xiàn)基于投影法的GPU并行計算研究。

在CUDA架構(gòu)中，GPU主要處理并行計算任務(wù)。而主控制邏輯任務(wù)和串行計算部分則由CPU完成。CUDA的編程模型如圖1所示。圖中的內(nèi)核(kernel)函數(shù)表示在算法中并行的步驟。由圖1可以看出，并行性包含2個層次(或粒度)：一個是在線程網(wǎng)格(grid)中線程(thread)塊與線程塊之間的并行；另一個是線程塊中線程之間的并行。對于大多數(shù)CFD算法而言，網(wǎng)格的劃分模式將顯著影響GPU上程序的效率，最高效的方式是1個網(wǎng)格節(jié)點對應(yīng)1個線程。此外，內(nèi)存組織和通信重疊也對程序的效率有重要影響，因此良好的網(wǎng)格劃分模式也應(yīng)該滿足對線程束(warp)全局內(nèi)存的聯(lián)合訪問。但是，沒有一種通用的優(yōu)化方法適用于GPU上的所有并行算法，關(guān)于如何劃分網(wǎng)格、組織內(nèi)存等，需要針對不同的問題進行不同的分析和優(yōu)化。

圖1 CUDA編程模型Fig.1 CUDA programming model

4 投影法的GPU并行化

4.1 投影法的分析

GPU強大的浮點運算能力能夠使投影法快速而高效地實現(xiàn)，但這并不等于算法中所有的運算過程都適合在GPU上進行并行計算，因此該算法的關(guān)鍵問題在于找出合適的部分交由GPU來并行計算。GPU與CPU架構(gòu)不同，前者的設(shè)計目的是面向矩陣類型的數(shù)值計算，而后者的設(shè)計目的是面向指令的高效執(zhí)行，因此前者更擅長于高度并行的復(fù)雜計算，而后者卻擅長執(zhí)行程序以及復(fù)雜邏輯事物。通過研究投影法的具體步驟可以發(fā)現(xiàn)，在進行計算中間速度U、V以及迭代求解壓力P時，各個子塊之間相關(guān)性不大，可以同時進行運算，符合并行計算的前提；另外，算法的性能瓶頸主要在迭代計算過程中，尤其是當網(wǎng)格規(guī)模較大時，傳統(tǒng)的CPU串行計算無法滿足快速計算的要求。為此，采用CUDA架構(gòu)線程塊通過GPU的并行計算能力實現(xiàn)投影法的求解，從而提高算法的運行速度。

4.2 并行過程設(shè)計

適合并行計算的部分為中間速度U、V以及壓力P的求解。對于控制方程的離散方式是基于交錯網(wǎng)格系統(tǒng)[15]，即中間速度U、V以及壓力P分別對應(yīng)1套網(wǎng)格，具體如圖2所示。

圖2 交錯網(wǎng)格系統(tǒng)Fig.2 Staggered grid system

為了達到在GPU中并行的目的，通常是1個線程對應(yīng)1個節(jié)點，當計算網(wǎng)格數(shù)目較大時，1個線程處理若干個節(jié)點。CUDA中線程組成線程塊，所以自然把計算區(qū)域分解為若干個與線程塊對應(yīng)的區(qū)域，線程塊的大小與維度由用戶指定。在GPU上求解中間速度U、V時，首先需要將所需數(shù)據(jù)從CPU內(nèi)存?zhèn)鬏數(shù)紾PU的全局儲存器，然后根據(jù)網(wǎng)格大小數(shù)來確定所需最優(yōu)線程塊數(shù)以及每個線程塊最優(yōu)線程數(shù)，一般32個線程組成1個線程粒度(warp)。因此，盡可能地使每個線程塊的線程數(shù)是線程粒度的整數(shù)倍，以實現(xiàn)高效的并行計算。然后將離散后的方程直接放入對應(yīng)的核函數(shù)中進行并行求解，值得注意的是，當計算求解邊界點時，由于二維數(shù)組中的1行或1列數(shù)值是已知的，因此可以直接轉(zhuǎn)化定義為一維數(shù)組求解，以此來提高計算效率。在投影算法中，雅可比迭代法用于求解壓力泊松方程。需要注意的是，當求解1個時層壓力值后，一定要進行同步操作，以免造成不必要的錯誤。

4.3 投影法GPU并行計算流程

基于CUDA核函數(shù)的異步執(zhí)行方式，實現(xiàn)了投影法的GPU并行計算。類似于Chen等[16]的實施方式，投影法的GPU并行計算流程如圖3所示。

圖3 投影法的GPU并行計算流程Fig.3 GPU parallel computing flow chart of projection method

整個工作分為兩部分：一部分由CPU完成；另一部分由GPU完成。在CPU上完成不適合并行運算的工作，主要包括：定義變量、初始化數(shù)據(jù)、網(wǎng)格生成、獲取網(wǎng)格坐標數(shù)據(jù)、誤差收斂判斷等。具體實現(xiàn)CUDA異構(gòu)運算的過程為：當變量聲明及存儲空間分配后，將數(shù)據(jù)初始化，接著將數(shù)據(jù)拷貝到GPU的全局內(nèi)存，依次使用3個內(nèi)核函數(shù)完成中間速度U和V的求解、壓力P的求解以及下一時層U和V的更新，最后將數(shù)據(jù)拷貝到CPU，進行收斂判斷。如果滿足收斂判斷條件，則結(jié)束循環(huán)，輸出結(jié)果。如果不滿足收斂條件，則繼續(xù)進行循環(huán)計算。

5 計算結(jié)果與分析

以二維方腔頂蓋驅(qū)動流為例進行實驗仿真。在數(shù)值試驗中，GPU為NVIDIA Quadro K1200，擁有512個流處理器(SP)，顯存帶寬最高可達80GBps。CPU為Intel(R)Xeon(R)E5-1660 v4，CPU處理單元的頻率是3.2 GHz。

當雷諾數(shù)Re=1 000、網(wǎng)格數(shù)為250×250時，基于投影法的GPU并行計算所得流場如圖4(a)所示，該計算結(jié)果與文獻[7](渦量-流函數(shù))和文獻[17](同位網(wǎng)格的SIMPLER算法，CPU并行)中給出的仿真結(jié)果[如圖4(b)所示]一致，驗證了該并行算法的正確性。以下對基于投影法的GPU并行計算方法求解性能進行分析研究。

5.1 雷諾數(shù)的影響

為了分析不同流動工況對GPU并行計算性能的影響規(guī)律，對比研究了不同雷諾數(shù)下CPU串行計算時間和GPU并行計算時間，這里雷諾數(shù)為無量綱準則數(shù)，決定了頂蓋驅(qū)動流的流動工況。在不同網(wǎng)格規(guī)模下，Re=100、Re=1 000、Re=5 000時CPU串行計算時間如表1所示。從表1中可以看出，隨著雷諾數(shù)的增加，CPU串行計算時間逐漸減少。相應(yīng)的GPU并行計算時間如表2所示。GPU并行計算與CPU串行計算具有相同的規(guī)則，隨著雷諾數(shù)的增加，GPU并行計算時間也逐漸減少。

圖4 計算所得頂蓋驅(qū)動流的流場圖Fig.4 Velocity field of lid-driven cavity flow

表1 不同雷諾數(shù)下CPU串行計算時間

Table 1 Serial computation time of CPU under different Reynolds numbers

網(wǎng)格數(shù)串行計算時間/sRe=100Re=1000Re=500050×500.50.30.23100×1001053250×2501697644500×50053231853869

表2 不同雷諾數(shù)下GPU并行計算時間

Table 2 Parallel computation time of GPU under different Reynolds numbers

網(wǎng)格數(shù)并行計算時間/sRe=100Re=1000Re=500050×502.71.61.0100×100742250×250462012500×5001846530

為了衡量GPU并行計算性能，這里引入了GPU并行加速比的概念，定義為CPU串行計算時間與GPU并行計算時間的比值。在不同網(wǎng)格規(guī)模下，不同雷諾數(shù)時GPU并行加速比如圖5所示。從圖5中可以看出，雷諾數(shù)的變化也就是流動工況的變化，對GPU并行加速比幾乎沒有影響，這是由于在不同雷諾數(shù)下GPU并行計算時間與CPU串行計算時間具有相同的變化規(guī)律。

圖5 不同雷諾數(shù)下GPU并行加速比Fig.5 Speed-up ratio of GPU parallel computation under different Reynolds numbers

5.2 網(wǎng)格規(guī)模的影響

隨著所涉及的流動傳熱問題復(fù)雜程度不斷提高，網(wǎng)格數(shù)量不斷增加，對于三維復(fù)雜問題，網(wǎng)格規(guī)模通常高達百萬或千萬，計算耗時大幅增加，通常需要幾天甚至幾個月的時間才能獲得1個收斂結(jié)果，遠遠不能適應(yīng)工程上快速設(shè)計的要求，迫切需要開發(fā)適用于大規(guī)模網(wǎng)格的并行計算方法。因此，研究了網(wǎng)格規(guī)模對GPU并行計算性能的影響規(guī)律。

Re=1 000時不同網(wǎng)格大小的GPU并行加速比如表3所示。

表3 網(wǎng)格規(guī)模與GPU并行加速比的關(guān)系

Table 3 Relationship of speed-up ratio of GPU parallel computation and grid sizes

網(wǎng)格數(shù)CPU串行計算時間/sGPU并行計算時間/sGPU并行加速比50×500.31.60.2100×100541.2250×25076203.8500×50018536528.5

從表3中可以看出，當網(wǎng)格數(shù)量較小(50×50)時，GPU并行計算時間要大于CPU串行計算時間。這是由于網(wǎng)格規(guī)模較小時，GPU并行計算時間主要用于CPU與GPU之間的信息傳遞，從而降低其求解性能。隨著網(wǎng)格規(guī)模的增大，CPU與GPU之間信息傳遞耗時所占的比例逐漸減小，所以GPU并行求解性能會不斷提升。當網(wǎng)格規(guī)模從100×100增大到500×500時，GPU并行加速比從1.2大幅提高到28.5。由此可以看出，網(wǎng)格規(guī)模越大，GPU并行計算的性能優(yōu)勢越明顯，而且在計算機本身固有的計算能力下，GPU并行加速比會隨著網(wǎng)格數(shù)的增加持續(xù)增大。因此，GPU并行計算是一種非常適用于大規(guī)模網(wǎng)格的并行計算方法。

6 結(jié)論

為了實現(xiàn)非穩(wěn)態(tài)流動傳熱問題的快速求解，開發(fā)了基于投影法的GPU并行計算策略，并以二維頂蓋驅(qū)動流為例研究了GPU并行計算求解性能，得出以下主要結(jié)論：

(1)不同雷諾數(shù)下GPU并行計算時間與CPU串行計算時間具有相同的變化規(guī)律，因此雷諾數(shù)的變化對GPU并行加速比幾乎沒有影響，說明GPU并行計算適用于不同的流動工況，均具有相同的加速性能。

(2)隨著網(wǎng)格規(guī)模的增大，GPU并行求解性能不斷提升，當網(wǎng)格規(guī)模從100×100增大到500×500時，GPU并行加速比從1.2大幅提高到28.5。在計算機本身固有的計算能力下，GPU并行加速比會隨著網(wǎng)格數(shù)的增加持續(xù)增大，說明GPU并行計算非常適用于大規(guī)模網(wǎng)格的計算問題。

雖然僅以經(jīng)典的二維頂蓋驅(qū)動流為例展開GPU并行計算求解性能研究，但對于不同流動傳熱問題，由于具有相同的GPU并行計算流程和信息傳遞方式，因此所得結(jié)論具有普遍適用性。