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

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

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

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

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

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

1 槽道流的計算方法

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

1.1 控制方程與其求解

不可壓流動的控制方程為

(1)

(2)

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

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

(3)

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

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

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

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

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

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

1.3 湍流的生成方法

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

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

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

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

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

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

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

2.1 平均速度

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

(4)

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

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

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

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

在比較+時同樣將所有算例的雷諾應(yīng)力與Lee的結(jié)果進(jìn)行了對比. 可以看出，雷諾正應(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所計算的Reτ較小，規(guī)律與CH2～CH4的結(jié)果相比有所區(qū)別. 從圖3(b)上來看，本文CH2～CH4的結(jié)果同樣也滿足方程(5).

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

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

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

3 局部特殊壁面邊界條件

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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