基于LES模擬的微溝槽湍流減阻特性數(shù)值研究

朱 涵，黎義斌，李小斌，郭東升，張紅娜，李鳳臣

(1.蘭州理工大學能源與動力工程學院，甘肅 蘭州 730050；2.天津大學機械工程學院，中低溫熱能高效利用教育部重點實驗室，天津 300350；3.中山大學中法核工程與技術學院，廣東 珠海 519082)

湍流減阻對水上船舶、潛艇、飛機以及長輸油管道等運輸工具的節(jié)能減排具有重要意義.在流體中運動的物體會受到摩擦阻力和壓差阻力，消耗了絕大部分能量，并且在超音速的飛行器還會受到激波阻力，這三種阻力占比中，其中摩擦阻力約占總阻力的40%～50%[1].由于阻力與航速呈平方關系，與推進功率呈三次方關系[2]，減小摩擦阻力對提高航速、航程和節(jié)約能源具有重要的意義.

溝槽減阻技術，即在物體和流體接觸的表面上布置微小的縱向溝槽，以達到湍流減阻的目的.該技術實際也是效仿魚類的仿生減阻技術之一，縱向溝槽在物體表面的布置，將一定程度上改變與粘性阻力緊密相關的湍流擬序結構，抑制了渦結構的發(fā)展，在適當條件(如一定的雷諾數(shù)Re范圍、溝槽尺度)下具有明顯的減阻效果，該技術在管道輸送、航天、機械設備及體育運動中已有運用的實例.李育斌等[3]在1∶12的機翼模型表面具有湍流流動的區(qū)域順流向粘貼肋條薄膜后，進行試驗后發(fā)現(xiàn)可減小阻力5%～8%；南京航空航天大學的潘家正[4]提出了“微型滾動軸承”理論，在湍流邊界層底部按照一定間距橫向分布的溝槽能鎖住流動的小渦，并獲得了約10.2%的減阻效果；KSB公司[5]也在多級泵的葉片表面加工了一定形狀的溝槽后，發(fā)現(xiàn)綜合效率提高了1.5%.另外，在管道流動中，內壁面帶微肋條/溝槽結構的計算中，也獲得了流動阻力減小[6-7]及傳熱效果改善的結論[8-9].可見，在相應平面上布置溝槽形狀，可以起到減阻增輸?shù)男Ч?

在溝槽減阻的機理研究方面，大部分理論認為溝槽面上湍流邊界層的粘性底層增厚，近壁區(qū)的湍流強度、雷諾切應力降低.粘性底層增厚，使得局部表面摩擦阻力(粘性切應力)降低，實現(xiàn)減阻.另外，從湍流結構上來講，部分理論認為湍流擬序結構在溝槽面上發(fā)生了變化，大尺度渦結構減弱，從而降低了阻力.溝槽結構減阻的研究盡管有了較大發(fā)展，但仍有許多問題需要解決，如大Re數(shù)條件下的流動減阻機理、溝槽適用性等內容需要深入研究.

早期已有許多學者對縱向微溝槽在牛頓流體中的性能進行了數(shù)值模擬研究，但對縱向微溝槽在牛頓流體中的減阻機理依然不完善，且絕大多數(shù)研究針對于較大尺度的溝槽寬度及很小的尖峰寬度，如V型、薄肋片式矩形等縱向微溝槽[10].但上述溝槽結構易受到流體動力學破壞且減阻范圍限制較大，大尺度的溝槽不能適應高速流動的減阻需求.本文通過參考魚類等仿生結構，建立了微米尺度的縱向溝槽微結構，溝槽構型為半圓形.分別對光滑平板和微溝槽面的湍流流場進行了大渦模擬，觀測了溝槽壁面湍流邊界層的擬序結構和流動特征，并定量分析了溝槽的減阻效果及其減阻機理.

1 數(shù)值計算方法及結果驗證

為了分析溝槽表面的減阻效果，對光滑表面和溝槽表面分別進行數(shù)值模擬，對比兩者的流動特性.

1.1 模型建模及邊界條件

光滑平板幾何結構如圖1(a)所示，其流向長度為0.1 m，法向高度為0.04 m，展向寬度為0.04 m.對于帶溝槽表面，參考仿生結構，設計了溝槽為半圓形構型，其溝槽寬度s為60 μm，深度h為30 μm，間隔t為20 μm，單個溝槽形狀見圖1(b).考慮到溝槽寬度很小，且欲構建的溝槽數(shù)量較多，整體建模困難，故其整體結構通過單溝槽陣列獲得.此溝槽流向長度和法向高度均與光滑平板結構相同.

圖1 流體域幾何模型示意圖

在邊界條件設置上，主要流動能量消耗在固體壁面上，將下壁面設定為無滑移壁面.因研究目標為充分發(fā)展狀態(tài)下的湍流流動結構，而湍流需要經(jīng)過層流區(qū)、過渡區(qū)才可發(fā)展成湍流，為節(jié)能計算資源，將流向方向(+z)上的前后面設置為周期性邊界，流體在流向上的運動通過設置流向上的質量流量來驅動.在展向方向(x)上，溝槽模型的網(wǎng)格數(shù)較多，亦將展向兩側面設置為周期性邊界.

1.2 計算域網(wǎng)格劃分及網(wǎng)格無關性驗證

光滑平板的網(wǎng)格劃分如圖2(a)所示，采用結構化網(wǎng)格劃分，流向均勻分布200個網(wǎng)格，展向上均勻分布100個網(wǎng)格，為了捕捉到近壁區(qū)的擬序條帶結構(如“上掃”與“下掠”)，在法向(+y)上采用邊界層加密，法向上對數(shù)分布100個網(wǎng)格，第一層網(wǎng)格為10 μm，增長比例為1.005，整體網(wǎng)格數(shù)為210萬.對于溝槽結構，如圖2(b)所示，也采用結構化網(wǎng)格，微溝槽模型流向均勻分布200個網(wǎng)格，法向上對數(shù)分布50個網(wǎng)格，在一個單元體結構內，展向分布8個網(wǎng)格.整個溝槽面陣列個數(shù)為50，總體網(wǎng)格數(shù)為410萬.

網(wǎng)格無關性驗證結果如圖3所示，以光滑平板為例，改變模型的網(wǎng)格劃分尺寸，分別生成170萬、190萬、200萬、210萬、220萬、230萬六種網(wǎng)格數(shù)量.通過監(jiān)測下壁面上的剪切應力，發(fā)現(xiàn)在網(wǎng)格數(shù)量210萬之后，剪切應力結果幾乎不再改變.210萬網(wǎng)格數(shù)與220萬、230萬數(shù)計算誤差分別等于0.094%、0.157%，同時210萬網(wǎng)格數(shù)與200萬、190萬、170萬網(wǎng)格計算計算誤差分別等于0.535%、1.794%、4.942%.可見當光滑平板網(wǎng)格數(shù)在210萬左右時，計算結果趨于穩(wěn)定，故以此作為計算標準；按照相同的分析方法對微溝槽結構模型網(wǎng)格進行無關性驗證，得到微溝槽網(wǎng)格數(shù)量為410萬.

圖3 光滑平板網(wǎng)格無關性驗證

1.3 湍流模型及數(shù)值計算方法

為了盡可能捕捉到流動中的細節(jié)結構并節(jié)能計算量，本文數(shù)值計算選取大渦模擬方法.大渦模擬理論有兩個基本假設，即(1)湍流的平均特性，由大尺度流動相干結構控制，而小尺度湍流結構的影響基本可以忽略；(2)高雷諾數(shù)下的小尺度湍流體現(xiàn)出各向同性的特點[11-13].故數(shù)值模擬重點放在比網(wǎng)格尺度大的大渦運動上，并通過引入模型來模擬小尺度的小渦運動，該模型稱為亞格子應力模型.這里選用WMLES(Algebraic Wall-Modeled LES Model)亞格子模型，在WMLES中，僅在對數(shù)層的內部激活模型的RANS部分，而邊界層的外部則被修改的LES公式覆蓋.由于邊界層的內部流動決定于LES模型的雷諾數(shù)，因此WMLES方法可以相同的網(wǎng)格分辨率應用于不斷增加的雷諾數(shù)，以進行流動模擬.WMLES模型是一種利用空間尺度過濾的結構函數(shù)模型.其代數(shù)表達式結合了混合長度模型、修正的Smagorinsky模型以及Piomelli的壁阻尼功能，其中渦黏度定義為

μt=min[(κdw)2，(CSmagΔ)2]·S·{1-exp[-(y+/25)3]}，

(1)

公式中：dw為計算點與壁面間距；S為應變率張量；參數(shù)κ=0.41、CSmag=0.2；濾波尺寸Δ根據(jù)特定的流場條件進行選取，有

Δ=min(max(Cw·dw；Cw·hmax；hwn)；hmax)，

(2)

公式中：hwn為沿壁面法向的網(wǎng)格步長；hmax為壁面網(wǎng)格的最大步長；Cw為通過實驗確定的常數(shù)，取值0.15.

在離散方法的選取上，梯度項采用Least Squares Cell Based進行離散，壓力項采用Body Force Weighted進行離散，動量項采用Bounded Second Order Differencing格式，瞬態(tài)項采用Bounded Second Order Implicit.速度和壓強之間的迭代關系，采用PISO相鄰校正算法進行求解.時間步長的選取對 PISO算法的精度非常重要，當在預測修正過程中，使用越小的時間步長，可取得越高的時間精度.在計算完成后，通過湍渦周轉時間τ=(ν/ε)1/2對時間步長進行后驗，其中，ν為運動粘度，ε為湍動能耗散率.步長選取為1×10-4s，小于一個湍渦周期，符合精度要求.

1.4 模擬方法的可靠性驗證

為了驗證數(shù)值模擬方法的準確性，以光滑平板為模擬對象，計算3個不同來流速度(即不同Re數(shù))下的摩擦系數(shù)cf.

分別以流場法向高度H及邊界層厚度δ99為特征尺度定義的Re數(shù)為

(3)

公式中：ρ為水的密度；U為來流速度；μ為流體的動力粘度，這里流體的水介質.

壁面摩擦阻力為

F=τwΑ，

(4)

公式中：τw為壁面剪切力；A為壁面表面積.

摩擦系數(shù)cf定義為

(5)

根據(jù)邊界層理論，在光滑平板上，cf與Reδ99的理論關系式為一隱式關聯(lián)式：

(6)

在邊界層流動中，流向速度沿法向高度增加而增加，當速度達到平均來流的99%時，此時高度可認為是邊界層厚度δ99.以此計算出Reδ99，并反推出cf，獲得cf計算值與理論值的對比如圖4所示.可見計算值與理論值十分接近，說明LES的計算結果是準確的.在湍流流動結構方面，進一步給出來流速度7 m/s時的速度分布與速度脈動統(tǒng)計規(guī)律.該工況下，對應ReH=2.8×105，Reδ99=1.2×105，流向平均速度分布和速度脈動均方根分布，如圖5所示.從圖5可以看出，計算結果和壁湍流理論規(guī)律吻合很好，同時從湍流脈動強度隨法向高度的變化分布來看，分布合理.說明LES方法能獲得理想的計算結果，及本文方法適用于牛頓流體的壁湍流流動問題.

圖4 cf計算值與理論值對比

圖5 來流7 m/s時的速度分布與速度脈動統(tǒng)計規(guī)律

2 結果及分析

在驗證了計算結果的可靠性之后，下面給出溝槽減阻流動的特性及流場結構分析.

2.1 微溝槽減阻效果分析

本研究中，針對具有較寬尖峰寬度的半圓形縱向微溝槽進行分析，溝槽寬s=60 μm，深h=30 μm，間隔t=20 μm，t/s=0.33.定義無量綱溝槽寬度s+為

(7)

公式中：ν為運動粘度；cf為光滑平板的摩擦系數(shù).

通過對比溝槽結構與光滑平板所受到的摩擦阻力，減阻率η的計算如式(8)所示，若η>0，表示微溝槽表面實現(xiàn)了湍流減阻，反之則為增阻.

(8)

計算不同來流速度下的減阻效果η，結果如表1所示.可以看出，溝槽結構表面在來流10 m/s以下均出現(xiàn)了明顯的減阻效果，但隨著雷諾數(shù)的增大，減阻效果先增大后減小，并在ReH=2.8×105時，減阻效果達到最大，此時減阻效果為19.11%.來流20 m/s的工況，溝槽表面的阻力相比光滑表面增加，即此時溝槽出現(xiàn)了增阻效應.

表1 溝槽表面減阻效果

溝槽表面減阻率η隨s+變化的趨勢，如圖6所示.可以看出，減阻率先隨著s+的增大而增大，而在s+=14.64時達到最大值.在工況s+=34.97時出現(xiàn)負值，出現(xiàn)了增阻效果.從無量綱溝槽寬度s+定義出發(fā)，該參數(shù)實際和以s為特征尺度的Re具有同等的意義，說明在相同的溝槽尺寸下，減阻效果跟湍流流動狀態(tài)有著密切的聯(lián)系，這將在下面進行深入討論.

圖6 減阻率隨無量綱溝槽寬度變化趨勢

2.2 渦結構數(shù)值分析

首先從溝槽表面的速度統(tǒng)計進行分析，然后再進行渦結構分析.沿法向設置用于速度統(tǒng)計的監(jiān)測線，如圖7所示.該監(jiān)測線具體位于光滑平板壁面上方以及溝槽谷內外上方.通過上述監(jiān)測線，提取近壁面位置處的瞬時流向速度以及經(jīng)過時間統(tǒng)計后的平均速度分布.

圖7 法向監(jiān)測線布置示意圖

平均來流速度7 m/s時流向速度沿法向的分布，如圖8所示.此時溝槽表面s+=14.64.對溝槽工況來說，y+<0表示所在法向位置位于溝槽谷內，y+>0在溝槽之外，對于光滑平面來說，法向y+均大于0.

圖8 不同y+值的流向速度分布(來流速度7 m/s)

紅色實線為層流情況時拋物線速度分布曲線，如圖8所示.該拋物線與溝槽谷內的瞬時速度分布(黑色實線)及時間平均速度分布(黑色虛線)均比較吻合，說明此時溝槽谷內大部分為層流流動，渦的侵入現(xiàn)象很少.在該工況下，微溝槽通過其縱向結構，使得在溝槽上方出現(xiàn)了“微型滾動軸承現(xiàn)象”，將流體與壁面的高速流動轉化為流體與流體間的高速流動，亦即造成了表面的部分滑移，從而減小了流動阻力.

來流速度7 m/s時，渦結構及其壁面其應力分布規(guī)律，如圖9所示.圖9中給出了光滑表面與溝槽表面流場渦結構特征.其中，使用Q準則作為渦結構識別的判據(jù)，兩者閾值相同.整體渦結構可見，溝槽結構出現(xiàn)的渦大多為升起的法向渦，而在光滑平板中有較多流向渦與壁面頻繁的沖刷，以致光滑平板產(chǎn)生了較高的剪切應力.

圖9 光滑壁面與溝槽結構表面流場渦結構(Q準則，Q=6×104)

針對不同來流速度的工況，渦結構的對比分析.同樣使用Q準則作為渦結構判據(jù)，如圖10所示.Q的閾值為Q/Qmax=0.007.可以看出，對應來流5 m/s和7 m/s，s+等于10.74和14.64，同樣渦結構大部分為升起的法向渦，且與溝槽壁面直接接觸的渦結構數(shù)量較少.相對比來流20 m/s的工況，s+=34.97，溝槽谷內有較多渦侵入，破壞了溝槽內固有的粘性底層，速度梯度增加明顯，故反而將出現(xiàn)增阻效果.上述結果也印證了s+<30，一般為減阻區(qū)間，而s+>30時，流動阻力相比光滑表面增加.

圖10 溝槽結構渦結構分布對比(閾值Q/Qmax=0.007)

2.3 低速條紋結構特性

在壁湍流中由于近壁面區(qū)不斷升起和下沉的渦結構，如圖11所示.會在不同法向高度的平面內呈現(xiàn)不同量級的低速區(qū)和高速區(qū)相間的速度分布，稱為低速條紋分布.圖中也給出了來流速度5 m/s時，截取y+=5光滑平板與溝槽結構的流向速度分布圖，可見其分布著低速條紋.低速條紋的出現(xiàn)，伴隨著渦結構的生長和消散，其密集程度也代表了漩渦產(chǎn)生的概率.下面進一步闡述在光滑表面及溝槽表面上的低速條紋結構分布.

針對兩種壁面，如圖12所示.在法向高度40 μm處建立了一條沿展向的監(jiān)測線，在光滑平板中，該線位于y+=10處，而在溝槽結構中，該線則位于y+=6.7處.

以來流速度7 m/s為例，在光滑平板與溝槽結構里，如圖13所示.不同流向位置處的流向速度w對展向距離的自相關分析結果Corr(w，w).通?？烧J定自相關函數(shù)Corr(w，w)中原點和第一個正峰值之間的距離可作為近壁低速條紋展向空間尺度的量度.顯然，該距離越小，說明低速條紋越密集，渦結構的演化越頻繁.

圖13 流向速度對展向距離的自相關函數(shù)

由圖13可以看出，光滑平板不同流向位置處的條紋間距在Δx+=98～176之間，而溝槽結構的條紋間距在Δx+=527～555之間，即光滑平板的條紋間距要遠小于溝槽結構中的條紋間距，說明在光滑平板的壁湍流結構中有著較多的湍流流向條帶分布，壁面和渦結構的相互作用十分頻繁，故易產(chǎn)生較大的摩擦阻力.

3 結 論

本文針對高雷諾數(shù)條件下微米級半圓形微溝槽結構的湍流減阻特性進行了數(shù)值模擬和結果分析，結果發(fā)現(xiàn)，在雷諾數(shù)ReH=1.99×105～7.96×105范圍內，半圓形溝槽結構有較高的減阻率，但在來流20 m/s條件下，相比光滑結構出現(xiàn)了增阻現(xiàn)象.在合適的雷諾數(shù)下，微溝槽結構在縱向上會形成“微型滾動軸承現(xiàn)象”，降低了壁面與流體之間的摩擦阻力，而在高雷諾數(shù)下，形成更細小的渦旋結構時，渦旋會侵入到溝槽谷內，故會出現(xiàn)增阻效果.溝槽結構的減阻效果還與溝槽尖端結構有關，尖端結構會阻礙展向渦的運動，使得溝槽結構在展向上分布較少的湍流條帶，提高了減阻效果.因此，微溝槽結構減阻效果是由其結構的橫縱向減阻效果共同作用的.