包 蕓 習(xí)令楚

(中山大學(xué)航空航天學(xué)院，廣州 510275)

引言

土地沙漠化是當(dāng)今人類面臨的一個重要生態(tài)環(huán)境問題，風(fēng)沙流是風(fēng)沙環(huán)境力學(xué)的核心研究內(nèi)容.風(fēng)沙流發(fā)生在高雷諾數(shù)大氣邊界層中，同時具有湍流和多相流兩個流體力學(xué)最具挑戰(zhàn)的復(fù)雜問題[1].在這類復(fù)雜流體力學(xué)問題中，風(fēng)場是沙塵運(yùn)動和輸運(yùn)的動力源和基礎(chǔ).蘭州大學(xué)在過去的二十年里，建立了世界上頂尖的野外風(fēng)沙特性觀察站，獲得了大量實際野外風(fēng)沙特性的第一手資料，并發(fā)性在大氣邊界層中存在超大型尺度湍流結(jié)構(gòu)[2].同時他們還展開了風(fēng)沙特性的風(fēng)洞實驗和數(shù)值模擬工作.在風(fēng)沙特性包括風(fēng)沙運(yùn)輸過程、電磁作用、沙塵暴等多方位的研究都取得大量的成果[3-9].風(fēng)沙流等研究成果很豐富[10-12].

采用高雷諾數(shù)的壁湍流模型研究風(fēng)場及其對沙塵特性的影響是當(dāng)今研究風(fēng)沙流特性的主要手段之一.計算機(jī)數(shù)值模擬是研究壁湍流及槽道流流動特性的重要方法[13-16].采用大渦模擬方法[17]，計算大氣邊界層的流動特性，已有很多研究工作和成果，并應(yīng)用到城市大氣污染擴(kuò)散和其他相關(guān)大氣風(fēng)場流動等研究[18-25].

對于自然環(huán)境大氣邊界層，數(shù)值模擬研究用于風(fēng)沙流的高雷諾數(shù)湍流風(fēng)場，需要足夠大的計算規(guī)模，才能有效的反映出實際野外環(huán)境中沙塵顆粒在氣流作用下的運(yùn)動規(guī)律.目前對高雷諾數(shù)湍流流動的數(shù)值模擬研究仍受到計算規(guī)模的限制.因此，充分利用我國的超級計算機(jī)資源，發(fā)展高效率并行計算技術(shù)，對高雷諾數(shù)湍流風(fēng)場以及進(jìn)一步的風(fēng)沙流數(shù)值模擬研究有重要的意義.

在對高瑞利數(shù)湍流熱對流進(jìn)行DNS 數(shù)值模擬的研究中，創(chuàng)建了高效的并行直接求解方法[26]，完成了大量的高和極高瑞利數(shù)的湍流熱對流DNS模擬，得到了豐富的計算數(shù)據(jù)，并在此基礎(chǔ)上取得了系列湍流熱對流物理特性的研究成果[27-30].在高雷諾數(shù)的壁湍流風(fēng)場的數(shù)值模擬研究中，大渦模擬(LES)是有效的計算工具.本文的工作是將用于湍流熱對流DNS 模擬的高效并行直接求解方法，擴(kuò)展到三維壁湍流風(fēng)場的大渦模擬計算中，建立新的用于三維壁湍流風(fēng)場大渦模擬的高效并行直接求解方法，為大渦模擬計算研究高雷諾數(shù)的壁湍流風(fēng)場提供有力的工具.

1 大渦模擬控制方程和亞網(wǎng)格模式

大氣邊界層中的空氣流動，在Bussinesqe 假定下可認(rèn)為是分層的不可壓縮流體.風(fēng)場大渦模擬控制方程是不可壓NS 方程.無量綱化的不可壓NS 方程為

其中，u為速度，P為壓力，Re為雷諾數(shù)，τij是亞網(wǎng)格應(yīng)力.

LES 方程和NS 方程的不同之處是LES 方程包含了亞網(wǎng)格通量.這些亞網(wǎng)格通量需要建立封閉模式，是大渦模擬的關(guān)鍵問題之一.關(guān)于復(fù)雜湍流中的亞網(wǎng)格通量模式仍在繼續(xù)研究和發(fā)展中，目前最常用的模式是渦黏和渦擴(kuò)散模型.Smagorinsky 渦黏模式是最早提出的亞網(wǎng)格模式，并被工業(yè)界廣泛應(yīng)用.渦黏模式的公式為

其中Smagorinsky 的渦黏模式

Cs是模式系數(shù)，?是網(wǎng)格尺度，,而是可解尺度速度變形率張量.各項同性湍流理論可確定模式系數(shù)Cs=0.10.但在實際應(yīng)用中發(fā)現(xiàn)，Smagorinsky 渦黏模式在底邊界壁面附近耗散過大，需要調(diào)整模式系數(shù).本文將采用Smagorinsky 渦黏模式作為亞網(wǎng)格模式進(jìn)行大渦模擬計算.

大渦模擬的控制方程特性與NS方程DNS模擬在求解過程中沒有本質(zhì)的區(qū)別，可以用相同的數(shù)值求解方法進(jìn)行數(shù)值計算.

高雷諾數(shù)湍流風(fēng)場的數(shù)值模擬研究受到計算規(guī)模的限制，一直是湍流數(shù)值模擬研究工作中的難題之一.因此，在超級計算機(jī)上研究和應(yīng)用高效的并行計算技術(shù)，對開展大規(guī)模高雷諾數(shù)湍流風(fēng)場的模擬以及為后續(xù)沙塵顆粒等運(yùn)動特性的聯(lián)立求解，有著現(xiàn)實的意義.

2 不可壓流動LES 模擬的并行直接求解方法

在過去的八年里，我們開展了極高瑞利數(shù)湍流熱對流的DNS 模擬研究.由于極高瑞利數(shù)湍流熱對流問題的特殊性，要求DNS 模擬的計算規(guī)模非常之大，使改進(jìn)和發(fā)展新的并行模擬求解方法的工作成為必然.2016 年終于在規(guī)模并行計算技術(shù)上有了突破性進(jìn)展，創(chuàng)建了不可壓流動湍流熱對流的并行直接求解方法(parallel direct method of DNS,PDMDNS)，并在“天河二號”超級計算機(jī)上實現(xiàn)了系列高瑞利數(shù)湍流熱對流的DNS模擬計算，并得到了系列的高瑞利數(shù)湍流熱對流特性的研究成果.

本文擬將湍流熱對流DNS 模擬的并行直接求解方法，擴(kuò)展到湍流風(fēng)場的LES 模擬中，希望能顯著有效的提高壁湍流風(fēng)場LES 模擬的計算效率，由此可進(jìn)一步擴(kuò)大壁湍流風(fēng)場的計算規(guī)模，提高計算壁湍流風(fēng)場的雷諾數(shù).

2.1 不可壓流動LES 的并行直接求解方法

由于該計算方法只限用于水平采用等距網(wǎng)格的矩形計算區(qū)域，適用于風(fēng)場的LES 模擬計算，如圖1 所示的風(fēng)場示意圖.

圖1 風(fēng)沙流中的風(fēng)場計算域示意圖Fig.1 The wind field calculation domain schematic

采用LES 模擬不可壓NS方程，本文采用投影法分步驟計算.動量方程中預(yù)估速度的計算采用容易實現(xiàn)并行的顯示格式.聯(lián)立壓力驅(qū)動速度項和連續(xù)方程，得到壓力泊松方程，需要全流場聯(lián)立求解，是數(shù)值求解過程中計算工作量最大的部分，同時全場聯(lián)立求解也給大規(guī)模并行造成困難.大規(guī)模并行計算的壓力泊松方程高效求解技術(shù)，是解決高雷諾數(shù)湍流風(fēng)場大規(guī)模并行計算的關(guān)鍵.

在二維湍流熱對流DNS 模擬的并行直接求解方法PDM-DNS[26]中，壓力泊松方程是二維的，因此僅需要在水平方向進(jìn)行FFT 變換解耦泊松方程，將方程變?yōu)槿龑欠匠?，再利用PPD 并行算法，綜合建立高效并行計算方法.

在“天河二號”超級計算機(jī)上，大規(guī)模的并行計算可以進(jìn)行MPI 和OpenMP 混編，MPI 并行計算需要交換分區(qū)邊界的數(shù)據(jù).

類似二維湍流熱對流DNS 的直接求解方法，在三維湍流風(fēng)場的LES 模擬中將矩形的風(fēng)場計算區(qū)域沿水平面對z進(jìn)行分割.相鄰的MPI 分割區(qū)域在并行計算中需要數(shù)據(jù)通訊，區(qū)域內(nèi)部可用OpenMP 并行，無需數(shù)據(jù)通訊.在此基礎(chǔ)上進(jìn)行壓力泊松方程的水平面FFT 變化時將不用并行計算，而是在縱向三對角方程求解時使用PDD 技術(shù)進(jìn)行并行計算.因而大大的提高了泊松方程并行求解的計算效率.

由此可見，二維流動的并行直接求解方法擴(kuò)展到三維計算，泊松方程的求解過程僅是在FFT 解耦的過程上有所不同.二維流動計算中用一維FFT，三維流動計算中用平面二維FFT 變換，其它求解過程基本一致.

對三維湍流風(fēng)場的LES 模擬，同樣并行計算的困難在于求解壓力泊松方程.三維泊松方程的直接求解方法的具體做法為，在三維空間的xoy水平平面上使用二維平面FFT 變換，將空間3 個方向上都要求聯(lián)立求解的壓力泊松方程解耦，使泊松方程變換成為只在z方向上的三對角方程.求解三對角方程后通過FFT 反變換就可得到全場泊松方程解.

三維壓力泊松方程為

其對應(yīng)的上下壓力無梯度邊界條件為

沿流向及展向方向的壓力邊界條件為周期邊界條件.

在x和y方向使用等距網(wǎng)格，z方向可使用變距網(wǎng)格.在點(i,j,k)對應(yīng)的3 個方向網(wǎng)格長度分別為?x,?y和?zk，壓力泊松方程的二階精度中心差分離散格式可寫成如下形式

根據(jù)以上離散泊松方程的特點，使用xy平面上的二維離散余弦傅里葉變換，可以對泊松方程的離散方程進(jìn)行解耦.為自動滿足壓力的邊界條件，選用二維離散余弦傅里葉變換，表達(dá)式為

M對應(yīng)x方向的網(wǎng)格數(shù)Nx，N對應(yīng)y方向的網(wǎng)格數(shù)Ny.以上變換可以通過使用FFTW 軟件包實現(xiàn).將上式變換應(yīng)用到壓力泊松方程中，可得到沿z方向聯(lián)立的三對角方程

對于給定的p和q，變換后的壓力泊松方程在需要xy和z方向上解耦，變?yōu)橄盗械娜龑欠匠探M.完成所有三對角方程組的求解后，通過對應(yīng)的變換即可完成壓力泊松方程的求解.

利用以上高效的壓力泊松方程并行直接求解方法，聯(lián)合其它易并行的動量方程等計算，本文將已創(chuàng)建的二維湍流熱對流并行直接求解方法，擴(kuò)展到三維不可壓湍流風(fēng)場的LES 模擬中，建立LES 模擬的并行直接求解方法，為大規(guī)模高效并行計算高雷諾數(shù)湍流流動，提供全新的數(shù)值計算技術(shù)和計算方法.

2.2 并行計算規(guī)模和計算效率

實際計算研究中，壁湍流風(fēng)場的大渦模擬已有許多成熟的研究成果.但在實際自然災(zāi)害問題的風(fēng)場研究中，空間尺度都很大，從而導(dǎo)致風(fēng)場具有很高的雷諾數(shù)，需要計算規(guī)模巨大.如何實現(xiàn)大規(guī)模高雷諾數(shù)風(fēng)場數(shù)值計算仍是一個瓶頸問題.本論文的研究特點主要在采用新的并行計算技術(shù)，建立了一個可進(jìn)行大規(guī)模高效并行的計算方法，可以顯著的提高風(fēng)場大渦模擬的計算效率，為實現(xiàn)高Re風(fēng)場以及由風(fēng)帶來的物質(zhì)輸運(yùn)的計算研究提供有價值的計算工具.

采用半槽道壁湍流大渦模擬計算，對新方法PDM-LES 的并行效率進(jìn)行驗證.

計算域為三維矩形.作為驗證計算，計算網(wǎng)格在(x,y,z)方向均取256×256×256.以x作為流動方向，y方向為與x方向垂直的另一水平方向，z方向為垂直地面方向.邊界條件為x方向上下游及y方向為周期邊界，z方向下底面為無滑移邊界，上邊界為對稱邊界條件.在計算過程中，采用定流量計算方法，以保證流動持續(xù)進(jìn)行，即在每個時間迭代計算中對由于粘性摩擦損失的流量進(jìn)行修正.

計算在“天河二號”超級計算機(jī)上進(jìn)行.每個計算機(jī)節(jié)點可以有32 線程無需交換數(shù)據(jù)的并行計算，采用OpenMP 方式進(jìn)行并行計算.節(jié)點間在并行計算是需要交換數(shù)據(jù)，采用MPI方式進(jìn)行并行計算.

表1 桿不同計算節(jié)點的計算效率Table 1 The computational efficiency of different compute nodes

表1 給出了采用不同節(jié)點數(shù)的計算結(jié)果，分別采用1，2 和4 個節(jié)點進(jìn)行了計算效率的對比計算.以1 個節(jié)點的計算作為基礎(chǔ)計算工作量，可以看到，采用1 個節(jié)點每天可迭代71.6 萬個時間步長，當(dāng)用4 個節(jié)點共128 核的并行計算，加速比為3.59，計算效率達(dá)到90%.

由此可見，本文建立的針對風(fēng)場的壁湍流大渦模擬的并行直接求解方法(PDM-LES)，具有較高的并行計算效率，可用于規(guī)模風(fēng)場的數(shù)值計算研究.

3 三維風(fēng)場的規(guī)模并行計算

采用建立的PDM-LES 方法，對三維壁湍流風(fēng)場進(jìn)行大渦模擬計算.取來流Re=10000 進(jìn)行試算，檢驗壁湍流計算結(jié)果的合理性.在以高度為1 的無量綱計算域中，計算域選取流動流向、展向和垂向長度分別為8×2×1，網(wǎng)格分別為Nx?Ny?Nz=1024 ?192 ?160.水平采用等距網(wǎng)格，垂向采用變距網(wǎng)格，計算網(wǎng)格規(guī)模約為3.800×107.

圖2 給出了來流Re=10 000 時三維瞬時速度分布.三維瞬時速度分布圖中可以看到，湍流流動的速度分布存在明顯的脈動特性.計算結(jié)果得出，在現(xiàn)有計算條件下Reτ≈4000.

圖2 三維瞬時速度云圖Fig.2 3D instantaneous speed cloud map

圖3 為不同高度xy平面的瞬時速度分布.圖3(a)為近壁面平行于xy平面上的瞬時速度分布(z+=12)，可以看到壁湍流近壁面區(qū)域分布著眾多的沿流動方向分布的低速條紋和高速條紋.圖3(b)為離壁面距離略高的xy平面上的瞬時速度分布(z+=200)，仍然可以看到明顯的條帶狀速度分布.隨著高度的增加，條帶狀速度分布的尺度增大.這些條紋反映出壁面附近的流動存在沿流向方向的條帶狀擬序結(jié)構(gòu).瞬時速度條帶狀擬序結(jié)構(gòu)是湍流流動的最主要流動特征，并且起到控制湍流流動特性的重要作用.

圖3 不同高度xy 平面的瞬時速度分布Fig.3 Instant velocity distribution of xy planes at different altitudes

計算湍流平均速度特性.選用多個不同時刻的瞬時速度場取平均值，再對水平平面的速度進(jìn)行平均計算，得到壁湍流流動在時間和空間上的平均速度場湍流特性縱向分布.

圖4 中給出了平均速度場U+和雷諾應(yīng)力的縱向分布特性.圖4(a)為速度U+的縱向分布，縱向距離z+取對數(shù).可以看到，第一個數(shù)值點的距離小于1 個z+，表明近壁面網(wǎng)格距離基本滿足計算的要求，本算例的計算網(wǎng)格選取是合理的.U+的縱向分布在z+≈20 處分為兩段，近壁面區(qū)域U+的分布曲線是彎曲的，離開壁面區(qū)域U+隨z+的分布在單對數(shù)坐標(biāo)中呈現(xiàn)直線關(guān)系，表明速度分布存在對數(shù)律.圖4(b) 是雷諾應(yīng)力的縱向分布情況.雷諾應(yīng)力的分布特性符合湍流邊界層的特性規(guī)律，但雷諾應(yīng)力的最大值偏小.

圖4 壁湍流平均速度場特性Fig.4 Average velocity field characteristics of wall turbulence

圖4 壁湍流平均速度場特性(續(xù))Fig.4 Average velocity field characteristics of wall turbulence(continued)

總體上說，本文計算的壁湍流風(fēng)場的湍流特性基本合理.

下一步工作是，在“天河二號”超級計算機(jī)上開展高雷諾數(shù)的壁湍流大渦模擬計算，增加計算網(wǎng)格規(guī)模，同時也增加有效的并行計算線程數(shù)，使得能在可接受的時間內(nèi)完成高雷諾數(shù)的壁湍流大渦模擬計算，為風(fēng)沙流的模擬計算等研究服務(wù).

4 結(jié)論

利用大渦模擬數(shù)值研究三維高雷諾數(shù)風(fēng)場特性，超大規(guī)模的數(shù)值模擬計算是十分必要的.依靠超級計算機(jī)技術(shù)的發(fā)展，建立高效并行的計算方法，對更深入研究計算實際野外風(fēng)沙流的物理和流動特性具有很重要的意義.

(1)本文將高效的并行直接求解方法擴(kuò)展到三維壁湍流風(fēng)場的大渦模擬計算中，建立了三維壁湍流風(fēng)場大渦模擬的高效并行直接求解方法PDM-LES.

(2)對并行直接求解方法的并行計算效率進(jìn)行的驗證計算結(jié)果表明，在“天河二號”超級計算機(jī)上進(jìn)行的1 至4 個節(jié)點，最大用到128 核的并行計算效率達(dá)到90%.

(3)計算Re=10 000 的結(jié)果顯示，湍流風(fēng)場瞬時速度分布在近壁區(qū)存在條帶狀的擬序結(jié)構(gòu)，并遠(yuǎn)離壁面條帶結(jié)構(gòu)尺度增加.平均場的速度分布符合湍流速度分布特性，即近壁面在z+≈20 之內(nèi)U+為線性分布，在外區(qū)U+滿足對數(shù)律分布.計算得到的風(fēng)場湍流特性基本合理.