沿流向微結(jié)構(gòu)溝槽流場直接數(shù)值模擬

李超群，李易,2,*，張晨曦，唐碩,2

1. 西北工業(yè)大學(xué) 航天學(xué)院，西安 710072 2. 陜西省空天飛行器設(shè)計重點實驗室，西安 710072

對于航空器和水下航行器等，減阻設(shè)計是一項重要內(nèi)容。尤其是對于飛機來說，減阻增升的氣動設(shè)計也是評價其先進性的指標之一[1]。減阻方法分為主動減阻和被動減阻兩種方式。相比于需要消耗能量的主動減阻方式，被動減阻不需要能量的消耗；其中，微型溝槽便是被動減阻的有效方式之一。NASA Langley研究中心的Walsh等[2-4]作為研究沿流向溝槽減阻原理的先驅(qū)者，先后對覆有不同形狀溝槽的平板阻力進行了大量實驗研究；研究發(fā)現(xiàn)，當沿流向溝槽無量綱高度h+<25并且溝槽無量綱寬度s+<30時，溝槽具有減阻效果，并當采用對稱V型且s+=h+=15的溝槽時，減阻效果最佳，可達8%。溝槽的存在雖然增加了浸潤面積，但是削弱了平板表面的湍流脈動，使阻力得到了降低。Sundaram等[5-6]對攻角為0°～12°的NACA0012翼型進行了實驗，發(fā)現(xiàn)溝槽減阻效果為10%～15%。Subashchandar等[7-9]對GAW(2)翼型和大攻角下的NACA 0012翼型進行了實驗，發(fā)現(xiàn)翼型表面附上溝槽壁面之后，其阻力明顯減小，減阻效果最大可達12%。Coustols和Schmitt[10]對表面覆有溝槽薄膜的空客A320 1/11縮比模型進行了實驗研究，發(fā)現(xiàn)當馬赫數(shù)為0.7時減阻效果最佳，黏性阻力降低了4.85%，相應(yīng)地，總阻力降低了1.6%。

胡海豹等[11]對不同來流速度下的溝槽表面湍動能分布律進行了實驗測量，指出溝槽結(jié)構(gòu)主要影響邊界層流場的黏性底層和過渡區(qū)。王晉軍等[12-13]利用激光多普勒測速計與氫氣泡流動顯示技術(shù)對溝槽表面湍流邊界層的帶條結(jié)構(gòu)進行了研究，結(jié)果發(fā)現(xiàn)溝槽平板的低速帶條具有更好的直線性，橫向流動得到抑制，壁面的湍流強度得到削弱，流動穩(wěn)定性增強。同時，崔光耀等[14]對溝槽上方的湍流邊界層進行了測量，發(fā)現(xiàn)減阻溝槽使得流向上的發(fā)卡渦數(shù)量明顯減少。

Zhang等[15]使用了隱式大渦模擬方法對不同攻角下表面覆有溝槽的翼型進行了數(shù)值模擬，結(jié)果表明雖然翼型表面的壓差阻力略微增大，但雷諾應(yīng)力大幅減小，二者的綜合作用使得阻力減小。周健等[16]結(jié)合已有實驗數(shù)據(jù)對雷諾平均Navier-Stokes(RANS)方法的k-ω模型進行修改，以得到新的?；瘮?shù)，使得新計算模型可在工程中對溝槽面減阻效果進行計算。

Choi等[17]采用二階中心格式對空間進行離散，同時結(jié)合Crank-Nicolson時間推進格式對沿流向溝槽湍流進行了直接數(shù)值模擬(Direct Numerical Simulation, DNS)，模擬結(jié)果對比了溝槽和平板的速度型；結(jié)果表明溝槽限制了流向渦的展向運動，使得渦引起的洗流影響范圍縮小，因而阻力得到了減小。同時Chu和Karniadakis[18]采用具有二階精度的譜方法對微結(jié)構(gòu)的層流和湍流流動進行了直接數(shù)值模擬，結(jié)果發(fā)現(xiàn)流向微溝槽僅針對湍流有較好的減阻效果。

對于沒有強剪切的流動，二階格式精度足夠；但是，對于復(fù)雜的近壁湍流和壁面剪切流，有必要采用高階格式[19]。為處理流場中存在間斷的問題，Liu等[20]依據(jù)ENO(Essentially Non-Oscillatory)格式提出了WENO(Weighted Essentially Non-Oscillatory)格式，WENO格式在處理非連續(xù)或有間斷的問題時具有良好的穩(wěn)定性和魯棒性。Zhang和Jackson[19]驗證了5階WENO格式結(jié)合低耗散、低色散Runge-Kutta(Low Dissipation and Dispersion Runge-Kutta scheme, LDDRK)方法求解不可壓縮流動的方法；數(shù)值結(jié)果表明，針對不同的算例，計算結(jié)果與精確解、已有的數(shù)值或?qū)嶒灲Y(jié)果吻合較好，證明了該方法具有較高的精度和魯棒性。

為更精確地解析壁面湍流，本文采用Jiang和Shu[21]提出的構(gòu)建WENO格式的方法構(gòu)建7階WENO格式，同時應(yīng)用分數(shù)步時間推進格式與LDDRK結(jié)合的方法求解不可壓縮流動，并解析微結(jié)構(gòu)溝槽壁面湍流，以探究其減阻的可能機理。

本文首先給出所采用的控制方程以及數(shù)值格式；然后，對比Kim等[22]的結(jié)果驗證算法的正確性與精度；接著，給出槽道的計算模型以及邊界條件設(shè)置；最后，列出計算結(jié)果及分析，并做相應(yīng)的總結(jié)和展望。

1 數(shù)值方法

1.1 控制方程

通常，流體的運動可以用三維非定常的Navier-Stokes方程和能量方程來描述，但對于低速的無黏不可壓流動，控制方程可以簡化為動量方程和連續(xù)方程。由于采用的數(shù)值方法為有限差分法，因此采用的控制方程的形式為微分形式。為便于求解，使用槽道的半高度δ、槽道中線處速度Ul、密度ρl以及動力黏度μl將Navier-Stokes方程進行無量綱化。

這樣，微分形式的無量綱化的非定常不可壓縮流動的控制方程為

(1)

式中：ui為瞬時速度矢量的3個分量(x、y和z分量)，使用槽道中線處速度Ul進行無量綱化；t為時間變量；xj為笛卡爾坐標，使用槽道半高度δ進行無量綱化；p為壓強，使用ρlUl2進行無量綱化；Re為雷諾數(shù)(基于槽道高度的1/2)，Re=ρlUlδ/μl。

1.2 數(shù)值格式

對于沒有較強的剪切或間斷的湍流流動，二階中心格式的精度已經(jīng)足夠，但因主要針對溝槽表面的壁湍流和壁面剪切流進行模擬，有必要采用高精度格式對近壁面處的湍流進行模擬[19]。對于處理存在剪切流和光滑流動的流場，魯棒性和計算效率優(yōu)秀的WENO格式是一種較好的選擇。

WENO格式是在高精度的ENO格式上改進的一種格式。Harten和Osher[23]提出ENO格式的思路是構(gòu)造“模板(Stencil)”，并在諸多模板中挑選出最光滑的模板構(gòu)造多項式；但這種方法舍棄了那些光滑性不佳的模板，造成了計算的浪費。在此基礎(chǔ)上，Liu等[20]對ENO格式進行了改進，提出了WENO格式，WENO格式利用這些模板的加權(quán)組合可以在光滑區(qū)構(gòu)造出更高階差分表達式，這樣提高了求解的精度和計算資源的利用率。隨后，Jiang和Shu[21]對該WENO方法進行了進一步改進，構(gòu)造了新的光滑度量因子，這樣，格式的計算效率得到了進一步提高，這也就是目前廣泛使用的WENO格式。

對于直接數(shù)值模擬的空間離散來說，5階WENO格式的耗散并不足夠低[24]，為了精確模擬微結(jié)構(gòu)壁面的湍流流動，同時考慮計算消耗，使用7階WENO格式對對流項進行離散，構(gòu)造方法使用Jiang和Shu在文獻[21]中提出的方法。

不妨令f代表速度矢量的任意分量ui，以正通量為例，構(gòu)建7階WENO格式，需要使用8個點并將其分成4個模板Sk(k=0,1,2,3)，每個模板中包含4個點：

(2)

(3)

表1 7階WENO格式中系數(shù)的值Table 1 Values of in the 7th WENO scheme

注:l為模板內(nèi)節(jié)點編號。

(4)

光滑度量因子ISk表達式為

(5)

將式(4)和式(5)代入式(3)中，再將式(3)代入式(2)中，即可得到關(guān)于正通量的7階WENO格式。

對于時間推進格式，采取基于分數(shù)步的時間推進格式[26](式(6)和式(7))，同時結(jié)合文獻[27]中給出的具有4階精度的最優(yōu)4-6 LDDRK格式(4和6分別表示第1和第2個時間步的放大因子)。這樣，非定常流動的時間精度大幅提高[19]。

(6)

同時：

(7)

對于黏性項采用6階中心差分格式即可，即

(8)

2 計算程序驗證

為驗證程序的正確性，使用了Kim等[22]的槽道流動結(jié)果進行驗證，計算條件與文獻[22]中的相同。槽道壁面流場的平均速度u+和均方根(Root-Mean-Square, RMS)速度脈動——urms、vrms和wrms對比如圖1所示，對比結(jié)果表明二者吻合較好，證明該程序在求解該類問題時具有較好的精度和可靠性。

圖1 本文程序與文獻[22]計算結(jié)果對比Fig.1 Comparison of calculation results of proposed design with results in Ref.[22]

同時，為體現(xiàn)高精度格式的優(yōu)勢，圖1還對比了高階格式和2階中心格式的速度型和均方根脈動速度曲線，可以發(fā)現(xiàn)：對于速度型，低階格式在過渡區(qū)和對數(shù)區(qū)精度略低于高階格式；對于脈動速度，低階格式所得的脈動結(jié)果要低于高階格式和文獻[22]中的結(jié)果。

3 微結(jié)構(gòu)溝槽建模

計算區(qū)域的槽道流動外形如圖2所示，該模型中，下壁面為無滑移的V型溝槽平面，上壁面為無滑移平板，其中，溝槽平面和平板壁面所在位置分別為y=-δ和y=δ，這樣，對比上下壁面的阻力即可得到表面溝槽的減阻效果[17-18,28]。

圖2 微尺度溝槽槽道計算區(qū)域Fig.2 Computational domain of channel mounted with riblets

溝槽近壁面網(wǎng)格如圖4所示，不同情況下網(wǎng)格的詳細信息如表2所示。其中在流向(x方向)上為均勻網(wǎng)格，Δx+=0.8；在垂直方向(y方向)上采用非均勻的雙曲正切分布，在壁面處Δy+=0.9；在展向(z方向)上采用均勻網(wǎng)格。

圖3 計算區(qū)域尺寸Fig.3 Size of computational domain

圖4 溝槽網(wǎng)格示意圖Fig.4 Schematic of grids near riblets

表2 不同情況下網(wǎng)格參數(shù)Table 2 Parameters of grids in different cases

注:N1、N2和N3為沿x、y和z方向的格點數(shù)；Δx+、Δy+和Δz+分別為沿x、y和z方向網(wǎng)格的無量綱長度

雖然并沒有進行網(wǎng)格無關(guān)性驗證，但3個方向上的網(wǎng)格密度已經(jīng)比文獻[17]的DNS模擬采用的網(wǎng)格高，因此可以認為該模擬采用的網(wǎng)格量是足夠的。

對于對稱的V型溝槽壁面，控制溝槽形狀的參數(shù)有兩個——溝槽的寬度s+以及溝槽斜壁與底面的夾角α。模擬的溝槽的寬度s+范圍為13～44，夾角α固定為60°，因此控制溝槽外形的參數(shù)只有寬度s+。流動沿x正方向，其關(guān)于槽道半高度的雷諾數(shù)Re=5 000。模擬時，基于最小網(wǎng)格尺度，將全局Courant-Friedrichs-Lewy(CFL)數(shù)控制為0.2，將無量綱時間步長控制為0.001；不同情況下，單個算例所計算的時間步至少為5×106。

對于邊界條件，流向(x方向)和展向(z方向)上設(shè)置為周期性邊界，上下壁面設(shè)置為無滑移壁面。

所有計算均在計算集群上進行，計算采用4個節(jié)點，集群每個節(jié)點CPU型號為E5-2670 v2 @ 2.50 GHz，24核，128 GB內(nèi)存，節(jié)點間采用InfiniBand(IB)技術(shù)互連。

4 計算結(jié)果及分析

給出了不同尺度下溝槽的減阻效果，并著重對減阻(s+=21.9)和增阻(s+=44.0)情況進行細致分析。針對這兩種情況，當速度場的統(tǒng)計特性趨于穩(wěn)定時開始進行數(shù)據(jù)采集，共采集了500個無量綱時間單元(tUl/δ)內(nèi)的數(shù)據(jù)并對其進行了分析。為演示其收斂特性，圖5繪制了s+=21.9的溝槽與平板壁面在500個無量綱時間內(nèi)壁面剪切率隨時間的變化關(guān)系，可以發(fā)現(xiàn)此時二者表面的流動已經(jīng)處于統(tǒng)計穩(wěn)定狀態(tài)。

分析數(shù)據(jù)后給出了減阻效果、溝槽表面和無滑移平板表面的平均速度型、自相關(guān)系數(shù)、湍流強度、瞬時速度場以及流向渦分布規(guī)律。

圖5 溝槽壁面和平板壁面剪切率的時間變化曲線Fig.5 Time history of wall-shear rates at both riblets and flat plate

4.1 減阻效果

對于溝槽的減阻效果，用式(9)衡量：

(9)

式中：η為阻力減小比例；Dr為溝槽壁面的阻力；Df為對應(yīng)的平板壁面的阻力。若η<0，則溝槽壁面具有減阻效果；若η>0，則溝槽壁面具有增阻效果。

經(jīng)過數(shù)值模擬后，得不同尺寸溝槽的減阻效果如圖6所示，其中所對比的實驗數(shù)據(jù)來自Bechert等[30]的實驗結(jié)果?？梢园l(fā)現(xiàn)阻力減小的最大值為9%左右，此時s+=16.0；當s+=44.0時，阻力增大了6.5%左右。

圖6 不同尺寸溝槽減阻效果對比Fig.6 Comparison of drag reducing ratios in different cases

4.2 平均速度型

對于減阻(s+=21.9)和增阻(s+=44.0)情況，溝槽表面上方速度型如圖7所示，其速度使用中線處速度進行無量綱化處理?？梢园l(fā)現(xiàn)，在相同y坐標上，溝槽底部上方的速度值大于溝槽頂部上方的速度值，但是隨著y值的增大，二者之差越來越小，并且當y大于某個值時，二者速度型重合，這也意味著在此范圍之外，流場各處受到溝槽不同部分的影響程度相同。

為詳細顯示溝槽上方的平均速度型，圖8繪制了溝槽頂部、溝槽斜壁中點和溝槽底部上方的速度型情況。對比平板表面速度型可知：減阻情況下，溝槽表面速度型在對數(shù)區(qū)出現(xiàn)了明顯的上移現(xiàn)象；增阻情況下，溝槽表面的平均速度型則有明顯的下移。這也一定程度上印證了文獻[31]中的相關(guān)結(jié)論。

圖7 s+=21.9和s+=44.0情況下溝槽表面的平均速度型Fig.7 Mean velocity profiles over riblets in cases where s+=21.9 and s+=44.0

圖8 s+=21.9和s+=44.0情況下溝槽頂部、斜壁中點以及底部上方的平均速度型Fig.8 Mean velocity profiles over peaks, midpoints and valleys of riblets in cases where s+=21.9 and s+=44.0

Wu等[31]指出，無論溝槽表面產(chǎn)生減阻效果還是增阻效果，其湍流-非湍流界面乃至整個流場均會受到溝槽的影響。以減阻情況為例，如圖8(a) 所示，溝槽不同位置處上方的速度型曲線在y+=50附近匯合，這表明在y+<50的區(qū)域之內(nèi)，溝槽的不同部分對流場中同一y+處流場的影響程度不盡相同；而在此區(qū)域之外，流場各處雖仍會受到溝槽的影響，但其受影響程度卻趨于一致。

4.3 自相關(guān)系數(shù)

討論的自相關(guān)性是脈動值的自相關(guān)性，自相關(guān)性代表了統(tǒng)計意義上的定量關(guān)系，并且如果速度變量關(guān)于位置是完全獨立的，那么他們的自相關(guān)系數(shù)為0[32]。自相關(guān)性是用統(tǒng)計方法表示速度在不同位置之間的相關(guān)關(guān)系，不同位置處速度分量的聯(lián)系程度可以用自相關(guān)系數(shù)表示，以x方向的速度分量u為例，其沿展向不同位置處的相關(guān)性系數(shù)定義為

(10)

空間自相關(guān)性函數(shù)的另一個性質(zhì)為

R(x,y,∞)=0

(11)

該性質(zhì)表明湍流場中相距很遠的兩點的隨機變量幾乎是統(tǒng)計獨立的，幾乎不存在相關(guān)性。因此，如果沿展向兩個位置處速度相關(guān)性系數(shù)趨近0，那么就可以說明展向距離足夠?qū)?可以使湍流流動得到充分發(fā)展。

通過對湍流流場內(nèi)的數(shù)據(jù)計算可得其相關(guān)性系數(shù)曲線，以減阻情況(s+=21.9)下流場速度沿展向的自相關(guān)系數(shù)(如圖9所示)為例?？梢园l(fā)現(xiàn)在溝槽近壁面處、槽道中線處和平板近壁面處，速度分量的相關(guān)性沿展向逐漸降低，最后降低到0附近。這表明流動的區(qū)域足夠大,可以使得沿展向湍流脈動不具有相關(guān)性，同時，這也說明湍流流場已經(jīng)充分發(fā)展。

圖9 速度沿展向的相關(guān)性系數(shù)Fig.9 Correlations of spanwise two-point velocity

4.4 湍流強度

溝槽頂部和溝槽底部上方的湍流強度示意圖如圖10所示，RMS速度脈動使用壁面剪應(yīng)力進行無量綱化處理，對比了減阻(s+=21.9)和增阻(s+=44.0)情況下的3個方向上的湍流強度。

在減阻情況下，溝槽上方沿x、y和z方向的湍流強度均有所降低，在展向上湍流強度最高均降低了約10%，對于流向湍流強度降低了約7%；與之相比，在增阻的情況下流向上溝槽頂部上方湍流強度明顯高于平板，垂直方向上溝槽表面湍流強度與平板壁面接近，展向上溝槽表面湍流強度略高于平板壁面?？梢哉J為溝槽產(chǎn)生減阻效果時，能較明顯降低3個方向上的湍流強度。

圖10 s+=21.9和s+=44.0情況下溝槽壁面湍流強度Fig.10 Turbulent intensity over riblets surface in cases where s+=21.9 and s+=44.0

4.5 瞬時速度場

給出減阻(s+=21.9)和增阻(s+=44.0)情況下某時刻流場速度u分量云圖和瞬時速度矢量場，如圖11和圖12所示。

從圖11(b)可以看出，當溝槽產(chǎn)生減阻效果時，大尺度流向渦位于溝槽之上且非?？拷鼫喜垌敳?，同時與溝槽頂部頻繁碰撞，可以形象地描述為流向渦“浮”于溝槽之上，并且由于溝槽的約束作用，流向渦的展向運動也受到了抑制，因此展向的湍流脈動也會得到削弱。

從圖12(b)可以看出，當溝槽產(chǎn)生增阻效果時，流向渦大部分位于溝槽內(nèi)部，并且與溝槽底部

圖11 s+=21.9情況下溝槽表面流動的瞬時流場Fig.11 Instantaneous flow field over riblets in case s+=21.9

圖12 s+=44.0情況下溝槽表面流動的瞬時流場Fig.12 Instantaneous flow field over riblets in case s+=44.0

的壁面相互碰撞。

4.6 順流向渦分布規(guī)律

為更清晰地描述順流向渦的分布規(guī)律，圖13繪制了溝槽附近150個時刻的流向渦分布，流向渦的位置取流向渦渦量最大值所在的空間位置。

從圖13中溝槽上方流向渦分布情況可以發(fā)現(xiàn)，減阻(s+=21.9)情況下，流向渦主要位于溝

圖13 s+=21.9和s+=44.0情況下流向渦位置Fig.13 Positions of streamwise vortices in cases where s+=21.9 and s+=44.0

槽上方并且在溝槽頂部附近有較多的渦聚集，流向渦僅能與溝槽頂部的尖峰頻繁碰撞而無法直接與溝槽內(nèi)壁互相作用，一方面，這種行為誘導(dǎo)出了二次渦；另一方面，由于沒有大量的渦進入溝槽底部，溝槽底部的流動便會比較平穩(wěn)，亦即降低了壁面處的湍流脈動。相反，增阻(s+=44.0)情況下，一方面，流向渦會進入溝槽內(nèi)部，造成流向渦直接與溝槽內(nèi)壁頻繁碰撞，增大了壁面的湍流脈動；另一方面，溝槽本身尺寸的增加導(dǎo)致浸潤面積增加，這二者均造成了阻力增大。

為表明流向渦沿垂直方向的分布規(guī)律，定義Πp(y)為流向渦的數(shù)目[33]，圖14繪制了平板表面上方減阻(s+=21.9)和增阻(s+=44.0)情況下溝槽頂部上方的流向渦的分布規(guī)律。

通過圖14的曲線對比可以看出，溝槽的存在會影響其上方流向渦的數(shù)目：當溝槽具有減阻效果時，其上方的流向渦的數(shù)目減少；而當溝槽具有增阻效果時，其上方的流向渦的數(shù)目增多。結(jié)合圖13的結(jié)果，可以認為對于減阻(s+=21.9)情況，溝槽的存在一方面減少了壁面上方的流向渦的數(shù)目，另一方面使得流向渦較集中地分布在溝槽頂部附近，與溝槽頂部頻繁碰撞。

圖14 s+=21.9和s+=44.0情況下流向渦數(shù)目沿垂直分布Fig.14 Population trend of streamwise vortices in wall units in cases where s+=21.9 and s+=44.0

5 結(jié) 論

采用高精度空間離散格式和時間推進格式求解溝槽表面的流動，同時研究了不同尺度溝槽對壁面湍流的影響以及減阻效果，得到如下結(jié)論：

1) 7階WENO格式結(jié)合用LDDRK格式改進的分數(shù)步時間推進格式對于求解該類問題具有較高精度。

2) 減阻情況下，大尺度流向渦的尺度略大于溝槽的寬度，這使得大尺度渦只能與溝槽尖峰碰撞而無法進入溝槽之內(nèi)，因而溝槽內(nèi)部的流動較為平緩，溝槽表面的湍流脈動也會減弱；同時，溝槽也使得近壁面處順流向渦的數(shù)目減少。因此，溝槽的存在雖然增大了壁面與流體之間的浸潤面積，但是使其表面的湍流脈動得到了削弱，并使近壁面處流向渦的數(shù)目得到降低，最后其綜合作用使得阻力減小。

3) 模擬發(fā)現(xiàn)，在溝槽上方的流場中存在某一高度閾值，當在此高度之內(nèi)時，流場受溝槽不同部分的影響程度不同；而在此高度之外時，流場受溝槽不同部分的影響程度趨于一致。同時，對比增阻情況，減阻情況下的這一高度閾值較小。

經(jīng)過數(shù)值模擬后，雖然得到了一些結(jié)論，但本文工作仍具有局限性，溝槽減阻機理仍然有待進一步研究，在未來的工作中將結(jié)合相應(yīng)的實驗等繼續(xù)對溝槽表面的流動進行深入分析。