胡海豹，宋保維，劉占一，黃橋高，潘 光

(西北工業(yè)大學(xué)航海學(xué)院，陜西 西安 710072)

0 引言

從仿生學(xué)角度探索降低航行體阻力的高效途徑，一直是減阻研究領(lǐng)域內(nèi)學(xué)者最為熱衷的研究課題之一。中國科學(xué)院童秉綱院士也曾在香山科學(xué)會議的報告中稱，魚類等水生動物和有翼昆蟲等飛行動物經(jīng)歷了億年的進(jìn)化過程，為了攫取食餌、逃避敵害、生殖繁衍和集群活動等生存需要，在漫長的適應(yīng)環(huán)境的自然選擇過程中，發(fā)展了各具特色的在水中游動和空中飛行的非凡能力，其整體功能漸趨優(yōu)化，為當(dāng)前的人造航行器和飛行器望塵莫及［1］。

人類在對海洋魚類的研究中驚訝地發(fā)現(xiàn)，鯊魚速度一般在50km/h左右，箭魚速度甚至可超過130km/h。海洋動物以其柔弱的身軀，卻能達(dá)到配備強(qiáng)大動力的艦船所難以企及的速度，在廣闊的大海里遨游穿梭……。學(xué)者對魚類外形的深入觀測發(fā)現(xiàn)，鯊魚、箭魚等魚類之所以具有良好阻力特性，一方面在于其擁有良好的流體動力外形，還有一個重要的原因就在于魚類表皮上遍布著脊?fàn)罱Y(jié)構(gòu)(圖1為不同種類鯊魚局部表皮的電子顯微鏡掃描圖［2］)。正是這些獨(dú)特的脊?fàn)罱Y(jié)構(gòu)的存在，有效降低了其在水中游動時所受的阻力。另外，在長期經(jīng)受大風(fēng)洗禮的廣袤沙漠里，也遍布著與海洋魚類表皮結(jié)構(gòu)類似的沙丘結(jié)構(gòu)(如圖2所示)。

受此啟發(fā)，國內(nèi)外學(xué)者開始關(guān)注脊?fàn)畋砻鏈p阻技術(shù)的研究。自20世紀(jì)80年代初，NASA蘭利研究中心Walsh、Wilkinson、Choi等先后報道了大量不同截面形狀縱向脊?fàn)畋砻娴臏p阻試驗結(jié)果，并得到了最大約8%的減阻量。目前國外在該技術(shù)領(lǐng)域已進(jìn)入工程測試階段，如空中客車將A320試驗機(jī)表面積的約70%貼上脊?fàn)畋砻姹∧ぃ_(dá)到了節(jié)油1% ～2%的效果;NASA蘭利中心在Learjet型飛機(jī)上開展的類似飛行試驗顯示，脊?fàn)畋砻娴臏p阻量約為6%左右［3］。國內(nèi)從20世紀(jì)90年代初才開始關(guān)注溝槽表面減阻技術(shù)的研究，如西工大宋保維等完成的縱向脊?fàn)畋砻嫠醇帮L(fēng)洞實驗獲得了超過6%的實驗室減阻量［4］;北航王晉軍等相繼報道的縱向脊?fàn)畋砻鏈y力及流場測試實驗表明，縱向脊?fàn)畋砻婢植孔枇p少高達(dá) 13% ～ 26%［5］;西交大宮武旗［1，6］等報道的縱向脊?fàn)畋砻骘L(fēng)洞試驗也得到了約7%的減阻量。與縱向脊?fàn)畋砻鏈p阻研究相比，橫向脊?fàn)畋砻鏈p阻屬減阻研究的新領(lǐng)域，在國內(nèi)外均無較全面的報道。

圖1 不同種類鯊魚局部表皮的放大照片F(xiàn)ig.1 Partial skin photograph of sharks

圖2 沙漠里的沙丘Fig.2 Photograph of sandbank in a desert

為促進(jìn)對脊?fàn)畋砻嫖⒂^流場的認(rèn)識，論文則著重探討脊?fàn)畋砻媪鲌鰯?shù)值計算方法研究，為脊?fàn)畋砻鏈p阻機(jī)理的研究提供參考。

1 脊?fàn)畋砻媪鲌鰯?shù)學(xué)模型

1.1 控制方程的確定

目前湍流數(shù)值模擬方法有三種:直接數(shù)值模擬(DNS)，大渦數(shù)值模擬(LES)，雷諾平均數(shù)值模擬(RANS)［7］。三種湍流數(shù)值模擬方法應(yīng)用于相同雷諾數(shù)條件下流場理論計算時，直接數(shù)值模擬的網(wǎng)格尺度最小，所以要求計算機(jī)的內(nèi)存最大，計算時間最長。雷諾平均數(shù)值模擬方法的網(wǎng)格尺度允許較大，因此要求計算機(jī)的內(nèi)存小，計算時間短;大渦數(shù)值模擬則介于兩者之間［8］。

脊?fàn)畋砻鏈p阻理論研究中，主要涉及脊?fàn)罱Y(jié)構(gòu)內(nèi)部及附近流場內(nèi)主要旋渦的發(fā)生和發(fā)展。三種湍流數(shù)值模擬算法中，直接數(shù)值模擬需要超高計算量、算法本身尚不成熟;大渦模擬與雷諾平均數(shù)值模擬同樣需要在近壁面采用壁面律來近似，但雷諾平均數(shù)值模擬方法計算量相對較小，且能夠提供脊?fàn)罱Y(jié)構(gòu)內(nèi)部主要旋渦分布等微觀流場信息［9］。經(jīng)綜合權(quán)衡，本文采用了雷諾平均的數(shù)值模擬方法。脊?fàn)畋砻媪鲌龅目刂品匠倘缦?

連續(xù)方程:

雷諾平均Navier-Stokes方程:

由于所模擬的時均流場屬定常流動，將上述公式化簡并在各個坐標(biāo)軸上展開:

1.2 湍流模型的選擇

目前常用的湍流模型很多，但每種湍流模型均有一定的使用范圍要求。經(jīng)過對標(biāo)準(zhǔn)k-ε模型、RNG k-ε模型等多種湍流模型的試用發(fā)現(xiàn)，相比其它湍流模型，RNG k-ε模型具有兩特點(diǎn):(a)通過修正湍流動粘度考慮了平均流動中的旋轉(zhuǎn)及旋流流動的影響;(b)在ε方程中增加了反映主流時均應(yīng)變率的項。因此，該模型可以更好地處理高應(yīng)變率及流線彎曲程度較大的流動［10］，也更適合于脊?fàn)畋砻娼趨^(qū)旋渦流動的仿真。該模型的方程形式如下:

其中，η0=4.377，β =0.012，Cμ=0.0845，αk= αε=1.39，C2ε=1.68。其它各項的具體計算公式如下:

2.3 邊界條件的設(shè)置

為解決網(wǎng)格數(shù)量與硬件條件、計算精度與收斂速度之間的矛盾，論文根據(jù)所選控制方程的類型，聯(lián)合使用了對稱邊界條件與周期性邊界條件等［11］。典型的計算域邊界定義如下:

(1)縱向脊?fàn)畋砻鏁r，沿z軸方向左、右兩邊為一對周期性入口和出口;

(2)橫向脊?fàn)畋砻鏁r，沿z軸方向左、右兩邊為速度入口和自由流出口;

(3)沿y軸方向上，上表面設(shè)置為光滑平板表面，下表面設(shè)置為脊?fàn)畋砻?

(4)沿x軸方向兩側(cè)為對稱面邊界條件。

2 脊?fàn)罱Y(jié)構(gòu)表面網(wǎng)格剖分方法

2.1 脊?fàn)畋砻嬗嬎阌虻慕?/h3>

2.1.1 縱向脊?fàn)畋砻嬗嬎阌虻慕?/p>

針對縱向脊?fàn)畋砻娴男螤钐攸c(diǎn)，本文采用三維流場模擬計算方法，將計算域選定為脊?fàn)畋砻嫜亓飨虻囊欢稳S空間。由于脊?fàn)罱Y(jié)構(gòu)的尺寸微小，這就要求在脊?fàn)罱Y(jié)構(gòu)表面的網(wǎng)格尺寸應(yīng)足夠的小;而另一方面由于計算資源的限制，整個計算域的網(wǎng)格總數(shù)又不能太大。反復(fù)的試算結(jié)果表明，縱向脊?fàn)畋砻媪鲌稣承缘挠绊懺?0倍脊峰高度處已經(jīng)可以忽略。為盡量減少網(wǎng)格數(shù)量，本文取計算域高度等于10倍脊峰高度。因此，縱向脊?fàn)畋砻娴挠嬎阌蜃罱K尺寸選取為Lx=20h，Ly=10h，Lz=10s(如圖3所示)，其中h為脊峰高度，s為脊?fàn)罱Y(jié)構(gòu)寬度。

圖3 三維計算域示意圖Fig.3 Three-dimensional computation region sketch

2.1.2 橫向脊?fàn)畋砻嬗嬎阌虻慕?/p>

橫向脊?fàn)畋砻嬖谡瓜蚍较虻拿總€截面上具有完全相同的結(jié)構(gòu)，可近似為二維流場。鑒于脊?fàn)畋砻嬷饕谕牧鳡顟B(tài)下具有減阻效果，因此模擬計算時，必須保證脊?fàn)畋砻娴牧鲌鎏幱诔浞职l(fā)展的湍流狀態(tài)。為此，論文在所模擬的脊?fàn)畋砻嬷霸O(shè)置了一段流動發(fā)展段(對于平板表面流動，雷諾數(shù)超過5×105時才能達(dá)到湍流狀態(tài)。為此，論文算例中，針對不同流速范圍，在其計算域前端設(shè)置了不同長度的流動發(fā)展段)。同時為排除計算域出口對脊?fàn)畋砻媪鲌龅母蓴_，論文在所模擬脊?fàn)畋砻嬷笥衷O(shè)置了一段平板流動段。另外，本文將橫向脊?fàn)畋砻嬗嬎阌虻母叨纫踩?0倍脊峰高度。所以最終確定的計算域如圖4所示。

圖4 二維計算域示意圖Fig.4 Two-dimensional computation region sketch

2.2 脊?fàn)畋砻婢W(wǎng)格剖分策略

CFD技術(shù)中，網(wǎng)格質(zhì)量好壞將直接影響計算結(jié)果的精度。要真實模擬一個流動問題，必須結(jié)合具體問題建立滿足計算精度要求的網(wǎng)格。作者在開展脊?fàn)畋砻媪鲌鰯?shù)值模擬研究的過程中，針對脊?fàn)畋砻媪鲌鲇嬎阌虻奶攸c(diǎn)，總結(jié)了四條脊?fàn)畋砻婢W(wǎng)格剖分策略:

(1)脊?fàn)罱Y(jié)構(gòu)具有很好對稱性和周期性，可采用縮小計算域的方法來提高脊?fàn)罱Y(jié)構(gòu)表面上網(wǎng)格的分辨率。本文中，脊?fàn)畋砻娴挠嬎阌騼H選取相鄰的數(shù)個脊?fàn)罱Y(jié)構(gòu)。

(2)由于脊?fàn)畋砻嬗嬎阌虻膸缀涡螤畋容^復(fù)雜，且在脊?fàn)罱Y(jié)構(gòu)的脊峰處存在尖角，使得在脊?fàn)罱Y(jié)構(gòu)上難以剖分高質(zhì)量的結(jié)構(gòu)化網(wǎng)格。為此，可采用非結(jié)構(gòu)化網(wǎng)格剖分技術(shù)［12］。

(3)針對脊?fàn)罱Y(jié)構(gòu)底部流動速度較小而頂部較大的特點(diǎn)，采用網(wǎng)格尺寸函數(shù)和計算域分區(qū)兩種網(wǎng)格控制方法，可實現(xiàn)脊?fàn)畋砻婢W(wǎng)格分布呈底(脊谷)疏頂(脊峰)密的特點(diǎn)。

(4)為比較相同條件下脊?fàn)畋砻婧推桨灞砻媪鲌鎏匦缘牟町?，同時又減小計算量，將與脊?fàn)畋砻?下表面)相對應(yīng)的平板面(上表面)設(shè)置為壁面，即可將脊?fàn)畋砻婧推桨灞砻媪鲌龅臄?shù)值模擬融合于同一算例中。

2.3 典型脊?fàn)罱Y(jié)構(gòu)的網(wǎng)格剖分

圖5～圖9為根據(jù)上述脊?fàn)畋砻婢W(wǎng)格剖分策略，生成的脊?fàn)罱Y(jié)構(gòu)表面網(wǎng)格圖。其中，圖5～圖7為橫向狀態(tài)下V型、圓弧型、U型脊?fàn)畋砻娑S計算域內(nèi)網(wǎng)格局部放大圖，圖8～圖9為縱向狀態(tài)下V型、U型脊?fàn)畋砻嫒S計算域內(nèi)網(wǎng)格局部放大圖。

圖5 V型脊?fàn)畋砻嬗嬎憔W(wǎng)格局部放大圖Fig.5 Partial enlargement grid figure of V-riblets

圖6 圓弧型脊?fàn)畋砻嬗嬎憔W(wǎng)格局部放大圖Fig.6 Partial enlargement grid figure of Semicircular riblets

圖7 U型脊?fàn)畋砻嬗嬎憔W(wǎng)格局部放大圖Fig.7 Partial enlargement grid figure of U-riblets

圖8 V型脊?fàn)畋砻嬗嬎憔W(wǎng)格Fig.8 Grid figure of V-riblets

圖9 U型脊?fàn)畋砻嬗嬎憔W(wǎng)格Fig.9 Grid figure of U-riblets

3 脊?fàn)畋砻媪鲌瞿M實例

3.1 橫向脊?fàn)畋砻媪鲌瞿M

圖10所示為水速5m/s時不同形狀脊?fàn)畋砻嫠俣确植际噶繄D，圖11為相同速度下不同形狀脊?fàn)畋砻娴撵o壓強(qiáng)分布圖，圖12則為不同速度下不同間距橫向U型脊?fàn)畋砻鏈p阻曲線。綜合分析圖10～圖12所示的橫向脊?fàn)畋砻媪鲌龇抡娼Y(jié)果可以發(fā)現(xiàn):

圖12 不同間距橫向脊?fàn)畋砻鏈p阻曲線Fig.12 Drag reduction profiles above spanwise riblet surfaces with different space

1)從圖10可以發(fā)現(xiàn)，脊?fàn)罱Y(jié)構(gòu)的存在致使在其脊谷中形成穩(wěn)定的大小、形狀和位置基本相同的低速旋渦。這些旋渦上部流動的方向與來流方向一致，而下部的流動方向與來流方向相反，且既沒有向外擴(kuò)散，也無明顯的相互影響［10］。因此，在橫向脊?fàn)罱Y(jié)構(gòu)的內(nèi)部其壁面剪應(yīng)力的方向與外部來流方向相反，即這些低速旋渦產(chǎn)生了粘性推力。

2)從圖11可見，在每個橫向脊?fàn)罱Y(jié)構(gòu)的迎流面處存在一片壓強(qiáng)增大區(qū)域(該區(qū)域內(nèi)靜壓力顯著比背流面大)，而在脊?fàn)罱Y(jié)構(gòu)間的臺階平面段靜壓力與外流場基本一致。由此可見，橫向脊?fàn)罱Y(jié)構(gòu)在產(chǎn)生粘性推力的同時，又伴隨產(chǎn)生了一部分粘性壓差阻力。

3)從圖12(a)、圖12(c)中可以看出，與光滑平板表面粘性阻力僅為粘性摩擦阻力不同，橫向脊?fàn)畋砻嬲承宰枇Πㄕ承詨翰钭枇驼承阅Σ磷枇刹糠?，其中粘性壓差阻力隨脊?fàn)罱Y(jié)構(gòu)間距的增大迅速降低，而粘性摩擦阻力則隨脊?fàn)罱Y(jié)構(gòu)間距的增大而增大(當(dāng)脊?fàn)罱Y(jié)構(gòu)間距a=0時，粘性摩擦阻力為負(fù)值，即摩擦阻力為推動力)。因此，隨脊?fàn)罱Y(jié)構(gòu)間距的增大，橫向脊?fàn)畋砻娴恼承宰枇χ袎翰钭枇Ρ戎亟档?，摩擦阻力比重則不斷提高。

4)從圖12(b)、圖12(d)中可見，橫向脊?fàn)畋砻婢哂辛己玫臏p阻效果，最大減阻量超過10%，且在不同流速下橫向脊?fàn)畋砻娴南鄬p阻量均隨脊?fàn)罱Y(jié)構(gòu)間距的增大呈現(xiàn)先增大后又減小的趨勢。以s=h=0.1mm的V型橫向脊?fàn)畋砻鏋槔?其中s為脊?fàn)罱Y(jié)構(gòu)寬度，h為脊?fàn)罱Y(jié)構(gòu)高度)，當(dāng)脊?fàn)罱Y(jié)構(gòu)間距等于脊?fàn)罱Y(jié)構(gòu)寬度(a=s)時，減阻效果相對較好。

3.2 縱向脊?fàn)畋砻媪鲌瞿M

圖13所示為4m/s水速下流向中間橫截面處V型脊?fàn)畋砻媪飨蛩俣确植紙D，圖14為相同工況下V型脊?fàn)畋砻嫘郎u分布圖，圖15為相同工況下V型脊?fàn)畋砻婕魬?yīng)力分布圖，圖16則為V型脊?fàn)畋砻鏈p阻量隨脊?fàn)罱Y(jié)構(gòu)無量綱寬度的變化曲線。綜合分析圖13～圖16所示的縱向脊?fàn)畋砻媪鲌龇抡娼Y(jié)果可以發(fā)現(xiàn):

1)從圖13(a)可以發(fā)現(xiàn)，在每一個脊?fàn)罱Y(jié)構(gòu)底部沿y向的速度變化較頂角處要緩慢的多，說明在其的底部存在著低速流帶;從圖14(b)可見，在縱向脊?fàn)畋砻娼趨^(qū)流場仍能夠分辨出粘性底層、過渡層及對數(shù)律層等分層結(jié)構(gòu)［4］，但縱向脊?fàn)畋砻嬲承缘讓釉龊瘢瑢?shù)律區(qū)上移，這與其它壁面湍流減阻技術(shù)研究中所得到的邊界層流場結(jié)構(gòu)變化規(guī)律相一致［3，11］。

圖13 流向中截面上V型脊?fàn)畋砻媪飨蛩俣确植紙DFig.13 Velocity distribution of the middle cross-section above V-riblets

圖14 流向中截面上V型脊?fàn)畋砻嫘郎u分布圖Fig.14 Vortex distribution of the middle cross-section above V-riblets

圖15 V型脊?fàn)畋砻婕魬?yīng)力分布圖Fig.15 Wall shear stress distribution above V-riblets

圖16 V型脊?fàn)畋砻鏈p阻量隨無量綱寬度的變化曲線Fig.16 Correlation between the none dimension riblet width and the drag reduction quantity of V-riblets

2)從圖14可以看出，縱向脊?fàn)罱Y(jié)構(gòu)的存在使得在每個縱向脊?fàn)罱Y(jié)構(gòu)頂點(diǎn)(頂部)兩側(cè)存在一對穩(wěn)定的反向“二次渦”。該“二次渦”對可能正是脊?fàn)畋砻娼趨^(qū)湍流猝發(fā)過程中起主導(dǎo)作用的順流向“反向旋轉(zhuǎn)渦對”與縱向脊?fàn)罱Y(jié)構(gòu)相互作用的產(chǎn)物，其產(chǎn)生必然會削弱順流向“反向旋轉(zhuǎn)渦對”的強(qiáng)度，進(jìn)而降低近壁區(qū)湍流猝發(fā)的強(qiáng)度［6，10］。

3)從圖15可見，縱向脊?fàn)畋砻姹诿婕魬?yīng)力沿流向基本保持恒定，且不同形狀脊?fàn)畋砻姹诿婕魬?yīng)力均從脊?fàn)罱Y(jié)構(gòu)的底部到頂部逐步增大，即脊?fàn)畋砻娴母呒魬?yīng)力區(qū)主要出現(xiàn)在其頂部處，而在脊?fàn)罱Y(jié)構(gòu)內(nèi)部剪應(yīng)力明顯較低。

4)從圖16可見，縱向脊?fàn)畋砻娴臏p阻量隨脊?fàn)罱Y(jié)構(gòu)無量綱寬度s+(s+=s·u/μ，其中μ為運(yùn)動粘性系數(shù)，u為當(dāng)?shù)啬Σ了俣?的變大而先增大后減小，且在s+=10～25的范圍內(nèi)均具有減阻效果，最大減阻量約6.5% ，這也與相關(guān)報道基本一致［3，11］。

4 結(jié)論

上述橫向和縱向脊?fàn)畋砻媪鲌瞿M結(jié)果說明，論文所總結(jié)的數(shù)值計算方法不僅能夠得到脊?fàn)畋砻娴臏p阻特性，還可以獲得脊?fàn)罱Y(jié)構(gòu)內(nèi)部及附近流場內(nèi)流速、旋渦等細(xì)微結(jié)構(gòu)的分布規(guī)律，且仿真結(jié)果與已有相關(guān)文獻(xiàn)報道基本相一致。因此，本文所總結(jié)的脊?fàn)畋砻媪鲌鰯?shù)值計算方法具有很強(qiáng)的適應(yīng)性，基于該方法的脊?fàn)畋砻媪鲌鰯?shù)值計算研究將可為脊?fàn)畋砻嫖⒂^流場結(jié)構(gòu)及減阻機(jī)理探索提供理論參考。

［1］童秉綱，陸夕云.關(guān)于飛行和游動的生物力學(xué)研究［J］.力學(xué)進(jìn)展，2004，34(1):1-8.

［2］BALL P.Engineering:shark skin and other solutions［J］.Nature，1999，400:507-508.

［3］王晉軍.溝槽面湍流減阻研究綜述［J］.北京航空航天大學(xué)學(xué)報，1998，24(1):31-34.

［4］王柯，宋保維，潘光.條紋溝槽表面水下航行器減阻實驗研究［J］. 力學(xué)與實踐，2005，27(2):18-21.

［5］王晉軍，李亞臣.溝槽面三角翼減阻特性實驗研究［J］.空氣動力學(xué)學(xué)報，2001，19(3):283-287.

［6］宮武旗，李新宏，黃淑娟.溝槽壁面減阻機(jī)理實驗研究［J］. 工程熱物理學(xué)報，2002，23(5):579-582.

［7］張兆順，崔桂香，許春曉.湍流理論與模擬［M］.北京:清華大學(xué)出版社，2005，9.

［8］KALITZIN G，MEDIC G，IACCARINO G，et al.Near-wall behavior of RANS turbulence models and implications for wall functions［J］.Journalof Computational Physics，2005，204(1):265-291.

［9］胡海豹，宋保維，潘光，等.鯊魚溝槽表皮減阻機(jī)理的仿真研究［J］. 系統(tǒng)仿真學(xué)報，2007，19(21):4901-4903.

［10］KOK J C.Resolving the dependence on freestream values for the k/ε turbulence model［J］.AIAA J.2000，38:1292-1295.

［11］王福軍.計算流體動力學(xué)分析［M］.北京:清華大學(xué)出版社，2004.

［12］劉超群.多重網(wǎng)格法及其在計算流體力學(xué)中的應(yīng)用［M］.北京:清華大學(xué)出版社，1995.

［13］胡海豹，宋保維，毛昭勇，等.隨行波表面減阻降噪機(jī)理探索［J］.火力與指揮控制，2007，36(10):1928-1932.

［14］王晉軍，蘭世隆，陳光.溝槽面湍流邊界層結(jié)構(gòu)實驗研究［J］.力學(xué)學(xué)報，2000，32(5):621-626.