高雷諾數湍流風場大渦模擬的并行直接求解方法1)

包 蕓 習令楚

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

引言

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

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

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

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

1 大渦模擬控制方程和亞網格模式

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

其中，u為速度，P為壓力，Re為雷諾數，τij是亞網格應力.

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

其中Smagorinsky 的渦黏模式

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

大渦模擬的控制方程特性與NS方程DNS模擬在求解過程中沒有本質的區別，可以用相同的數值求解方法進行數值計算.

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

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

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

本文擬將湍流熱對流DNS 模擬的并行直接求解方法，擴展到湍流風場的LES 模擬中，希望能顯著有效的提高壁湍流風場LES 模擬的計算效率，由此可進一步擴大壁湍流風場的計算規模，提高計算壁湍流風場的雷諾數.

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

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

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

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

在二維湍流熱對流DNS 模擬的并行直接求解方法PDM-DNS[26]中，壓力泊松方程是二維的，因此僅需要在水平方向進行FFT 變換解耦泊松方程，將方程變為三對角方程，再利用PPD 并行算法，綜合建立高效并行計算方法.

在“天河二號”超級計算機上，大規模的并行計算可以進行MPI 和OpenMP 混編，MPI 并行計算需要交換分區邊界的數據.

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

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

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

三維壓力泊松方程為

其對應的上下壓力無梯度邊界條件為

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

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

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

M對應x方向的網格數Nx，N對應y方向的網格數Ny.以上變換可以通過使用FFTW 軟件包實現.將上式變換應用到壓力泊松方程中，可得到沿z方向聯立的三對角方程

對于給定的p和q，變換后的壓力泊松方程在需要xy和z方向上解耦，變為系列的三對角方程組.完成所有三對角方程組的求解后，通過對應的變換即可完成壓力泊松方程的求解.

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

2.2 并行計算規模和計算效率

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

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

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

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

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

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

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

3 三維風場的規模并行計算

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

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

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

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

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

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

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

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

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

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

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

4 結論

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

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

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

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