湍流槽道流并行直接求解方法

習(xí)令楚, 謝嘉斌, 包 蕓

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

在流體力學(xué)的湍流研究中，壁湍流是一個(gè)重要的問(wèn)題，其在理論研究和工業(yè)應(yīng)用方面都具有關(guān)鍵意義，引起了人們極大的興趣. 在壁湍流的問(wèn)題當(dāng)中，槽道流由于具有簡(jiǎn)單的幾何模型和邊界條件，被視為一種非常有效的研究模型.

槽道湍流的實(shí)驗(yàn)研究與數(shù)值模擬研究成果非常豐富. 直接數(shù)值模擬被認(rèn)為是能夠得到最準(zhǔn)確數(shù)據(jù)的數(shù)值模擬方法. 1987年，Kim等[1]通過(guò)直接數(shù)值模擬的方式計(jì)算了Reτ≈180的槽道流，并將計(jì)算的結(jié)果與早前Eckemlmann等[2]的結(jié)果進(jìn)行了對(duì)比. 不過(guò)按照湍流研究的慣例，通常認(rèn)為槽道流的內(nèi)區(qū)為0.1個(gè)半槽道高度，而外區(qū)要大于50個(gè)壁面單位，Reτ必須足夠高才會(huì)在存在內(nèi)區(qū)與外區(qū)的交錯(cuò)段[3]. 同時(shí)Smits等[4]提出當(dāng)Reτ達(dá)到103時(shí)湍流研究的結(jié)果接近工業(yè)上所需要的雷諾數(shù)，研究具有實(shí)際應(yīng)用價(jià)值. Monty等[5]進(jìn)行了最大Re=2 000的槽道湍流實(shí)驗(yàn), 為Hoyas等[6]計(jì)算的Reτ=2000和Alamo等[7]計(jì)算的Reτ=1 000的結(jié)果提供實(shí)驗(yàn)驗(yàn)證數(shù)據(jù). Schultz等[8]進(jìn)行了1 000≤Reτ≤6 000的槽道流實(shí)驗(yàn)，主要的工作是測(cè)出不同雷諾數(shù)下的統(tǒng)計(jì)量，并且將新的數(shù)據(jù)與之前的DNS或者實(shí)驗(yàn)的數(shù)據(jù)進(jìn)行對(duì)比. 在數(shù)值模擬方面，Lozano等[9]計(jì)算了Reτ=4 200的槽道流，研究計(jì)算區(qū)域?qū)?shù)值模擬的影響，結(jié)論為計(jì)算區(qū)域(2πh×πh)在靠近壁面區(qū)域所得到的結(jié)果足夠準(zhǔn)確. Lee等[10]計(jì)算了Reτ=5 200的湍流槽道流，并且對(duì)之前的DNS與實(shí)驗(yàn)的結(jié)果進(jìn)行了驗(yàn)證與總結(jié). 上面的計(jì)算均采用譜方法進(jìn)行. Bernardini等[11]用二階有限差分法計(jì)算了Reτ=4 000的湍流槽道流，對(duì)他人的計(jì)算結(jié)果進(jìn)行了驗(yàn)證. 除了Bernadini等人采用二階有限差分法之外，Vela-Martin等[12]利用GPU對(duì)湍流槽道流開(kāi)展計(jì)算，采用的方法為二階有限差分法，計(jì)算的雷諾數(shù)達(dá)到Reτ=5300；Yamamoto等人[13]采用高階精度的有限差分法，計(jì)算了Reτ=8 000的槽道流，但是由于計(jì)算的網(wǎng)格分辨率不夠高，在速度脈動(dòng)的計(jì)算上面與目前所認(rèn)可的規(guī)律相比相差較大，目前對(duì)這組數(shù)據(jù)仍然存疑. 目前最大雷諾數(shù)的槽道流算例是Hoyas等人[14]計(jì)算的雷諾數(shù)達(dá)到Reτ=10 000的槽道流，這組數(shù)據(jù)計(jì)算區(qū)域較小(2πh×πh)，目的是為了得到槽道流湍流的單點(diǎn)統(tǒng)計(jì)規(guī)律.

減阻控制是湍流槽道流的研究熱點(diǎn)和重點(diǎn)之一，實(shí)驗(yàn)與數(shù)值模擬都進(jìn)行了大量的研究. 湍流減阻具有代表性的工作包括Quadrio等[15]采用的壁面展向振動(dòng)的方式、Luchini等[16]通過(guò)壁面溝槽限制了流向渦的展向流動(dòng)以及Maden等[17]通過(guò)等離子體激勵(lì)器控制的方式減小了槽道湍流的阻力等研究.

隨著超級(jí)計(jì)算機(jī)性能與并行技術(shù)的提高，高Re數(shù)的槽道流DNS計(jì)算成為可能. 本文結(jié)合所在課題組基于湍流熱對(duì)流所建立的并行直接求解方法(PDM-DNS)，將較高效的并行求解方法引入湍流槽道流計(jì)算當(dāng)中，實(shí)現(xiàn)較高效并行的湍流槽道流DNS模擬計(jì)算. 由于本文采用的計(jì)算方法是有限差分法，結(jié)合浸沒(méi)邊界法等計(jì)算技術(shù)，方便處理復(fù)雜的邊界條件，較容易將槽道流的計(jì)算研究擴(kuò)展到湍流流動(dòng)控制等計(jì)算研究領(lǐng)域.

1 槽道流的計(jì)算方法

由于計(jì)算規(guī)模巨大，高Re數(shù)湍流槽道流的DNS模擬往往需要使用高效并行計(jì)算技術(shù)并依托超級(jí)計(jì)算機(jī)進(jìn)行. 湍流槽道流的計(jì)算結(jié)果數(shù)據(jù)可以為壁湍流理論研究以及湍流流動(dòng)控制等研究提供基礎(chǔ)依據(jù).

1.1 控制方程與其求解

不可壓流動(dòng)的控制方程為

(1)

(2)

槽道流計(jì)算中，速度與壓力在流向與展向方向采用周期邊界條件, 上下壁面邊界條件為速度無(wú)滑移和零壓力梯度. 計(jì)算網(wǎng)格在流向與展向方向采用等距網(wǎng)格，垂直壁面方向采用壁面附近加密的非等距網(wǎng)格. 空間采用二階精度中心差分格式，時(shí)間采用二階精度的Runge-Kutta法進(jìn)行計(jì)算. 本文采用的計(jì)算步驟為投影法.

實(shí)際的槽道流是由壓力驅(qū)動(dòng)引起的流動(dòng). 本文采用定流量計(jì)算(CFR,Constant Flow Rate)，定流量指設(shè)定系統(tǒng)的流量恒定，系統(tǒng)的總流量不發(fā)生改變[18]. 控制流量的方法為在每個(gè)計(jì)算的單元內(nèi)引入外加力Fx與Fy. 設(shè)初始的流量為Q0,計(jì)算n+1時(shí)刻的流量為Qn+1，所以有

(3)

對(duì)于流向與展向的速度u和v，通過(guò)引入外加力F修正以保證流量恒定.

1.2 直接并行求解方法與并行效率

對(duì)于槽道流的DNS計(jì)算來(lái)說(shuō)，當(dāng)Re數(shù)增大就需要提高計(jì)算網(wǎng)格精度，所以采用高效的并行計(jì)算在高Re數(shù)湍流槽道流的DNS模擬中是不可缺少的.

課題組在之前為解決極高Ra數(shù)湍流熱對(duì)流的DNS數(shù)值模擬，建立了直接并行求解方法PDM-DNS(Parallel Direct Method of Direct Numerical Simulation)[19]，Direct Method指的是直接求解法，在求解的過(guò)程對(duì)于壓力項(xiàng)的求解沒(méi)有采用迭代的方式而是直接求解壓力泊松方程得到壓力. 將這一計(jì)算方法引入到槽道流的計(jì)算中來(lái)，初步實(shí)現(xiàn)具有較高并行效率的湍流槽道流DNS模擬并行計(jì)算. PDM-DNS方法的計(jì)算特點(diǎn)是，采用OpenMP和MPI混編并行，可在“天河二號(hào)”超級(jí)計(jì)算機(jī)上運(yùn)行. 其中關(guān)鍵技術(shù)是壓力泊松方程的高效并行求解，采用FFT解耦和PDD三對(duì)角方程并行求解技術(shù).

表1給出了在“天河二號(hào)”超級(jí)計(jì)算機(jī)上進(jìn)行槽道流DNS計(jì)算并行效率測(cè)試結(jié)果. 計(jì)算網(wǎng)格采用512×512×512，共1.34×109網(wǎng)格計(jì)算自由度. 采用PDM-DNS方法分別在超級(jí)計(jì)算機(jī)上應(yīng)用1～8個(gè)計(jì)算節(jié)點(diǎn)，即32～256計(jì)算核進(jìn)行計(jì)算，所有計(jì)算都能達(dá)到95%以上的并行效率，并且具有很好的并行效率的線性延展性. 因此，PDM-DNS方法可以給湍流槽道流DNS模擬提供較為高效的并行計(jì)算，使高Re數(shù)湍流槽道流的模擬得以實(shí)現(xiàn)并獲得大量有效的數(shù)據(jù).

表1 在nx×ny×nz=512×512×512的網(wǎng)格下并行相關(guān)參數(shù)

1.3 湍流的生成方法

槽道流以層流泊肅葉流動(dòng)作為初場(chǎng)，通過(guò)添加速度擾動(dòng)使其發(fā)展成為湍流槽道流. 擾動(dòng)大小幅值為初始速度大小的10%. 在外加擾動(dòng)的作用之下，流動(dòng)逐漸從層流向湍流進(jìn)行轉(zhuǎn)變.

圖1給出了槽道流壁面摩擦系數(shù)隨時(shí)間的變化過(guò)程. 在不斷添加人為的擾動(dòng)速度影響下，槽道壁面的摩擦系數(shù)隨著時(shí)間發(fā)生變化. 在無(wú)量綱時(shí)間T<20的時(shí)候，壁面摩擦系數(shù)幾乎是不變的；在20≤T≤35的階段，壁面的摩擦系數(shù)會(huì)出現(xiàn)一個(gè)突增的情況；在T>35之后，壁面摩擦系數(shù)發(fā)展到另一個(gè)較大的數(shù)值并趨于平穩(wěn)，此時(shí)流動(dòng)發(fā)展進(jìn)入湍流狀態(tài). 在流動(dòng)發(fā)展到湍流狀態(tài)之后撤去人為的速度擾動(dòng)，槽道流動(dòng)不會(huì)回到層流流動(dòng)而是穩(wěn)定在湍流流動(dòng)狀態(tài)，繼續(xù)計(jì)算最終得到充分發(fā)展的湍流槽道流，在此基礎(chǔ)上可得到湍流槽道流的研究數(shù)據(jù).

圖1 壁面摩擦系數(shù)隨著時(shí)間的變化

2 湍流槽道流DNS模擬結(jié)果與討論

設(shè)定半槽道的寬度為h，為了驗(yàn)證程序的準(zhǔn)確性，本文分別計(jì)算了Reτ≈180、550、1000、2 000時(shí)的4個(gè)算例，并且與相關(guān)的數(shù)值模擬進(jìn)行對(duì)比. 計(jì)算算例的主要參數(shù)具體如表2所示.

表2 槽道流算例的主要參數(shù)

對(duì)于網(wǎng)格精度，x方向控制在10個(gè)壁面單位以內(nèi)，y方向控制7個(gè)壁面單位以內(nèi). 根據(jù)Jimenez等[9]的研究成果，計(jì)算區(qū)域取2πh×πh中等范圍，可以保證槽道流內(nèi)區(qū)的計(jì)算具有足夠的準(zhǔn)確性. 本文所采用的對(duì)比數(shù)據(jù)是Lee[10]的DNS計(jì)算結(jié)果，他們的結(jié)果可靠性得到世界的認(rèn)可.

2.1 平均速度

在壁湍流的理論當(dāng)中，一般認(rèn)為在流動(dòng)的內(nèi)區(qū)到外區(qū)的平均流向速度是符合對(duì)數(shù)律的. 對(duì)數(shù)律的公式可以寫(xiě)為

(4)

其中，κ指卡門(mén)系數(shù). 在不同組的實(shí)驗(yàn)當(dāng)中所測(cè)出來(lái)的κ值也有些許差異，Osterlund[20]測(cè)出κ=0.38，Monty[5]測(cè)出κ=0.37，Nagib[21]測(cè)出κ=0.39. 由此可以認(rèn)為κ在0.38附近，本文選取κ=0.38.

圖2給出了本文4個(gè)算例的平均速度剖面線，并給出相應(yīng)的Lee[10]的DNS計(jì)算結(jié)果. 從圖2可以看出,CH1的結(jié)果在對(duì)數(shù)段相比其它的結(jié)果略高，原因是在Reτ≈180下，當(dāng)z+≈50時(shí)，z/δ≈0.28，當(dāng)z+≈30時(shí)，z/δ≈0.17，所以在槽道流當(dāng)中并不存在交錯(cuò)段(z+>50,z/δ<0.1)，但是存在對(duì)數(shù)段(z+>30,z/δ<0.3). 在Reτ≈550 時(shí)，在z+=50下，z/δ≈0.09，所以在槽道流當(dāng)中既存在交錯(cuò)段(z+>50,z/δ<0.1)，也存在對(duì)數(shù)段(z+>30,z/δ<0.3). 同樣現(xiàn)象在Lee的數(shù)據(jù)當(dāng)中也有所體現(xiàn). 在CH2～CH4的結(jié)果中都存在基本一致的對(duì)數(shù)段，且與Lee的結(jié)果相比較為接近，表明計(jì)算結(jié)果基本合理.

圖2 不同雷諾數(shù)下的平均速度曲線

2.2 雷諾應(yīng)力張量

在比較+時(shí)同樣將所有算例的雷諾應(yīng)力與Lee的結(jié)果進(jìn)行了對(duì)比. 可以看出，雷諾正應(yīng)力的最大值基本都出現(xiàn)在z+= 15的位置，而且隨著Re數(shù)的增大，雷諾應(yīng)力也隨之增大. 本文的結(jié)果與Lee在2015年的DNS結(jié)果相近，可以證明本文結(jié)果的準(zhǔn)確性.

Lozano等人[9]在文章中指出，流向方向的雷諾應(yīng)力最大值與雷諾應(yīng)力之間滿足以下關(guān)系：

(5)

由于CH1所計(jì)算的Reτ較小，規(guī)律與CH2～CH4的結(jié)果相比有所區(qū)別. 從圖3(b)上來(lái)看，本文CH2～CH4的結(jié)果同樣也滿足方程(5).

(a) 與Δx+的關(guān)系

(b) 最大值隨著Reτ的變化

綜上所述，平均速度與雷諾正應(yīng)力結(jié)果的對(duì)比情況表明，本文湍流槽道流的DNS計(jì)算結(jié)果是可信的.

3 局部特殊壁面邊界條件

湍流流動(dòng)控制是目前最熱門(mén)的研究課題之一，壁面邊界是湍流流動(dòng)控制的研究重點(diǎn). 相較于常用的譜方法，速度壓力法壁面邊界條件易于處理是針對(duì)流動(dòng)控制數(shù)值研究的一個(gè)優(yōu)點(diǎn)，而計(jì)算精度問(wèn)題則可以通過(guò)盡可能地提高計(jì)算規(guī)模來(lái)解決. 借助類(lèi)似浸沒(méi)邊界法等計(jì)算技術(shù)，粗糙壁面、滑移壁面以及各種吹吸氣等壁面處技術(shù)都很容易實(shí)現(xiàn).

為了探討壁面邊界條件變化對(duì)槽道湍流的影響，本文在CH1算例的基礎(chǔ)上進(jìn)行初步嘗試. 在下壁面加入三排微小毛刺，采用浸沒(méi)邊界法實(shí)現(xiàn)，計(jì)算壁面局部變化對(duì)槽道湍流流動(dòng)的影響. 本文采用離散力法計(jì)算浸沒(méi)邊界法中的力源項(xiàng). 在N-S方程中加入虛擬力源項(xiàng)，在控制方程中加入體積力來(lái)替代邊界條件. 在計(jì)算過(guò)程中求解完投影法的速度壓力之后再加入虛擬力源項(xiàng)得到最終的速度.

如圖4所示，在周期邊界條件的槽道流中的1/4位置處，加入三排長(zhǎng)寬約為0.044h(8個(gè)壁面單位)、高度約為0.022h(4個(gè)壁面單位)交錯(cuò)排列的微小毛刺. 毛刺在整個(gè)展向具有微小間隔并均勻分布，在流向方向上對(duì)流動(dòng)的影響只有局部作用.

圖4 槽道流加入粗糙毛刺示意圖

圖5給出的是在毛刺下游16個(gè)到320個(gè)壁面單位處的展向平均速度剖面，分別為毛刺高度4個(gè)壁面單位的4、10、20、40和80倍距離，同時(shí)給出了無(wú)毛刺時(shí)全場(chǎng)平均的速度分布. 在毛刺的作用下，鄰近位置的平均速度剖面線變化很大，平均速度曲線有所升高. 當(dāng)遠(yuǎn)離毛刺位置時(shí)，毛刺對(duì)湍流速度剖面的影響逐漸減小，到80倍毛刺高度距離位置處，湍流速度剖面已恢復(fù)到與無(wú)毛刺流動(dòng)基本一致. 可見(jiàn)壁面條件的改變對(duì)湍流流動(dòng)影響顯著.

圖5 加入毛刺后不同位置的平均速度曲線

由此可見(jiàn)，本文的速度壓力法對(duì)湍流流動(dòng)控制的研究，無(wú)論是加入展向的振動(dòng)、壁面吹吸氣還是改變槽道流的壁面形狀都是易于實(shí)現(xiàn)的. 加上本文方法具有高效并行特點(diǎn)，可以開(kāi)展規(guī)模巨大的高Re數(shù)湍流槽道流DNS模擬計(jì)算，為湍流流動(dòng)控制研究提供有效的計(jì)算工具和有價(jià)值的數(shù)據(jù).

4 結(jié) 論

建立了槽道流DNS計(jì)算的PDM-DNS方法，在“天河二號(hào)”超級(jí)計(jì)算機(jī)上使用256個(gè)CPU核進(jìn)行計(jì)算，其規(guī)模超過(guò)109自由度，計(jì)算并行效率超過(guò)95%，并得到良好的線性加速比. 由此為開(kāi)展高Re數(shù)湍流槽道流DNS模擬提供了有力的計(jì)算工具.

4個(gè)不同Re數(shù)的槽道流計(jì)算結(jié)果顯示，速度剖面都具有粘性底層和對(duì)數(shù)段分布，并與湍流速度對(duì)數(shù)段的理論預(yù)測(cè)值分布一致，雷諾正應(yīng)力分布合理，所有計(jì)算結(jié)果與他人的DNS結(jié)果吻合. 本文提出的計(jì)算方法所得到的湍流槽道流計(jì)算結(jié)果和數(shù)據(jù)是合理可信的.

相較于譜方法，本文采用的速度壓力差分求解方法方便于復(fù)雜邊界條件的處理. 在下壁面加入局部毛刺的算例結(jié)果展示，結(jié)合浸沒(méi)邊界法等計(jì)算技術(shù)，本文方法可以擴(kuò)展到復(fù)雜壁面條件的湍流流動(dòng)控制等研究中.