基于CFD技術(shù)的脊?fàn)畋砻嫱牧髁鲌?chǎng)特性研究

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

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

0 引言

源自于鯊魚皮的脊?fàn)畋砻鏈p阻技術(shù)已被人們成功應(yīng)用于石油輸運(yùn)、航空、航海、泳衣制造等廣大的領(lǐng)域。人們對(duì)脊?fàn)畋砻娴难芯恐饾u深入，在通過改變脊?fàn)罱Y(jié)構(gòu)形狀尺寸以追求減阻效果最大化的同時(shí)，對(duì)脊?fàn)畋砻鏈p阻機(jī)理的研究也從來沒有停止過。

美國(guó)航空航天局蘭利研究中心的Walsh是最早開始系統(tǒng)研究脊?fàn)畋砻娴膶W(xué)者之一，他通過大量的風(fēng)洞試驗(yàn)得出結(jié)論:在脊?fàn)罱Y(jié)構(gòu)無量綱尺寸小于30時(shí)，脊?fàn)畋砻娲嬖跍p阻效果;脊?fàn)罱Y(jié)構(gòu)附近形成的“二次渦”與流向渦的相互作用是脊?fàn)畋砻鏈p阻的重要原因［1－2］。日本的Yuji Suzuki利用三維粒子測(cè)速設(shè)備在水洞進(jìn)行了脊?fàn)畋砻鏈p阻試驗(yàn)，發(fā)現(xiàn)脊?fàn)罱Y(jié)構(gòu)能抑制順流向的能量向法向和展向轉(zhuǎn)化，從而實(shí)現(xiàn)減阻［3］，而El-Samni也通過對(duì)矩形管道進(jìn)行的數(shù)值模擬研究得到了類似的結(jié)論［4］。Haecheon Choi等利用直接數(shù)值模擬方法發(fā)現(xiàn):脊?fàn)罱Y(jié)構(gòu)無量綱尺寸小于20時(shí)，由于流向渦尺寸大約在30左右，無法深入脊?fàn)罱Y(jié)構(gòu)，只能停留于脊?fàn)罱Y(jié)構(gòu)頂部，因此脊?fàn)罱Y(jié)構(gòu)內(nèi)(除頂部區(qū)域)流動(dòng)較為安靜，壁面剪應(yīng)力相比光滑表面減小了［5］。本文對(duì)一系列表面布置不同尺寸順流向脊?fàn)罱Y(jié)構(gòu)的表面利用CFD方法進(jìn)行了湍流狀態(tài)下的數(shù)值模擬，并對(duì)脊?fàn)畋砻嫱牧髁鲌?chǎng)特性進(jìn)行了分析。

1 基于CFD技術(shù)的計(jì)算方法

本文選取的是一段三維的計(jì)算域，如圖1所示。

圖1 計(jì)算域Fig.1 Computation region

計(jì)算域底面沿x方向分為兩段，前段為平板，作為發(fā)展段，后段為測(cè)試段，用來對(duì)比研究脊?fàn)畋砻婧凸饣砻?。脊?fàn)罱Y(jié)構(gòu)沿流向布置，截面形狀均為高寬相等的V型。計(jì)算選取的脊?fàn)罱Y(jié)構(gòu)無量綱尺寸s+從6.7至43.7 不等。

s+的定義為:，其中s為脊?fàn)罱Y(jié)構(gòu)實(shí)際尺寸，u*為壁面摩擦速度，v為流體介質(zhì)的運(yùn)動(dòng)粘性系數(shù)。

圖2 脊?fàn)罱Y(jié)構(gòu)附近網(wǎng)格(垂直流向的截面圖)Fig.2 Grids near riblet structure(section perpendicular to flow direction)

脊?fàn)罱Y(jié)構(gòu)附近網(wǎng)格如圖2所示，脊?fàn)罱Y(jié)構(gòu)內(nèi)采用三棱柱網(wǎng)格，其他區(qū)域采用六面體網(wǎng)格。此次計(jì)算模擬的是平板表面脊?fàn)罱Y(jié)構(gòu)的流場(chǎng)特性，因此計(jì)算域上表面設(shè)置為對(duì)稱邊界條件。為避免側(cè)壁干擾，計(jì)算域兩個(gè)側(cè)面也都設(shè)置為對(duì)稱邊界。入口面設(shè)置為速度入口，為了保證測(cè)試段處于湍流狀態(tài)，筆者編制程序在入口面設(shè)定了一個(gè)已充分發(fā)展為湍流的速度分布。出口面為自由流。

本文所用流體介質(zhì)為水，控制方程為雷諾平均N-S方程。計(jì)算選取的湍流模型為Reynolds應(yīng)力方程模型(Reynolds Stress equation Model，簡(jiǎn)稱 RSM)，意在著重分析脊?fàn)罱Y(jié)構(gòu)對(duì)近壁區(qū)雷諾應(yīng)力分布造成的影響。該湍流模型更加嚴(yán)格地考慮了流線型彎曲、旋渦、旋轉(zhuǎn)和張力的快速變化，對(duì)于復(fù)雜流動(dòng)有更高的精度預(yù)測(cè)潛力。

流場(chǎng)數(shù)值計(jì)算采用SIMPLE方法，SIMPLE算法是目前工程上應(yīng)用最廣泛的一種流場(chǎng)計(jì)算方法，主要用于求解不可壓流場(chǎng)。其核心是采用“猜測(cè)—修正”的過程，在交錯(cuò)網(wǎng)格的基礎(chǔ)上來計(jì)算流場(chǎng)，從而達(dá)到求解動(dòng)量方程的目的［6］。文獻(xiàn)［7］曾利用SIMPLE方法對(duì)三維湍流流場(chǎng)進(jìn)行數(shù)值分析，獲得了與試驗(yàn)相吻合的湍流特性［7］。

作者也曾利用該方法對(duì)具有不同間隔的脊?fàn)畋砻孢M(jìn)行了湍流減阻研究，發(fā)表于《系統(tǒng)仿真學(xué)報(bào)》2009年19 期［8］。

2 結(jié)果及分析

2.1 與實(shí)驗(yàn)值的比較

圖3給出了脊?fàn)畋砻鏈p阻效果隨脊?fàn)罱Y(jié)構(gòu)無量綱尺寸s+的變化規(guī)律，并與 Walsh試驗(yàn)所得結(jié)果［1］進(jìn)行了比較。Walsh所采用的脊?fàn)罱Y(jié)構(gòu)截面形狀同樣是寬高相等的V型。其試驗(yàn)在一個(gè)矩形截面的小型風(fēng)洞進(jìn)行，主流風(fēng)速為 7.6m/s，雷諾數(shù)為 9.26 ×104，測(cè)得最大減阻量為8%。本文的計(jì)算中，雷諾數(shù)為9.95×105，所得最大減阻量為8.63%，此時(shí)s+=18.2，減阻量隨s+變化規(guī)律與Walsh所測(cè)結(jié)果也基本一致，只是減阻量最大值對(duì)應(yīng)的s+略有差別。

圖3 脊?fàn)畋砻鏈p阻效果與無量綱尺寸的關(guān)系Fig.3 Relation between drag reduction effect and non-dimension size

2.2 脊?fàn)罱Y(jié)構(gòu)附近流場(chǎng)

圖4和圖5分別給出了減阻和增阻情況下y-z面脊?fàn)罱Y(jié)構(gòu)附近速度矢量圖。兩種情況均可以觀察到以下現(xiàn)象:在脊?fàn)罱Y(jié)構(gòu)頂部形成了一對(duì)反向旋轉(zhuǎn)的流向渦對(duì);這些渦將動(dòng)量從脊?fàn)罱Y(jié)構(gòu)底部沿脊?fàn)罱Y(jié)構(gòu)壁面輸運(yùn)至頂部;渦在z方向的尺度決定于脊?fàn)罱Y(jié)構(gòu)的寬度，約為脊?fàn)罱Y(jié)構(gòu)寬度的一半。減阻和增阻情況下y-z面脊?fàn)罱Y(jié)構(gòu)附近流向渦對(duì)中速度最大值分別為外流場(chǎng)速度的 1.4%和 2.1%。

圖4 y-z面脊?fàn)罱Y(jié)構(gòu)附近流線圖(s+=18.2)Fig.4 Stream line around riblet structure in y-z plane(s+=18.2)

圖6和圖7分別給出了減阻和增阻情況下脊?fàn)罱Y(jié)構(gòu)附近流向渦渦量云圖。將流向渦渦量利用公式無量綱化，則在減阻和增阻情況脊?fàn)罱Y(jié)構(gòu)附近流向渦中心處的該值分別為0.26和0.36，表明增阻情況(s+=39.4)下脊?fàn)罱Y(jié)構(gòu)附近流向渦強(qiáng)度比減阻情況(s+=18.2)下大。

圖6 脊?fàn)罱Y(jié)構(gòu)附近流向渦渦量云圖(s+=18.2)Fig.6 Contour of streamwise vorticity near riblet structure(s+=18.2)

圖7 脊?fàn)罱Y(jié)構(gòu)附近流向渦渦量云圖(s+=39.4)Fig.7 Contour of streamwise vorticity near riblet structure(s+=39.4)

2.3 雷諾正應(yīng)力分布曲線

圖8和圖9分別給出了減阻和增阻情況下邊界層近壁區(qū)x向雷諾正應(yīng)力沿y方向分布曲線。由圖可見，兩種情況下脊?fàn)畋砻娼趨^(qū)x向雷諾正應(yīng)力較光滑表面均有明顯下降。Walsh、Choi等在試驗(yàn)中也觀察到了這一現(xiàn)象。但減阻情況下(s+=18.2)下降幅度更大，其峰值較光滑表面下降了16%，峰值位置與光滑表面基本一致。增阻情況(s+=39.4)下，相比光滑表面，峰值位置有所右移。兩種情況下x向雷諾應(yīng)力峰值大小、位置的差異反映出近壁區(qū)渦結(jié)構(gòu)大小、強(qiáng)度的不同。

圖8 x向雷諾正應(yīng)力分布(s+=18.2)Fig.8 Normal stress distribution in x direction(s+=18.2)

圖9 x向雷諾正應(yīng)力分布(s+=39.4)Fig.9 Normal stress distribution in x direction(s+=39.4)

圖10和圖11分別給出了減阻和增阻情況下邊界層近壁區(qū)y向雷諾正應(yīng)力沿y方向分布曲線。減阻情況(s+=18.2)下，測(cè)量位置在脊?fàn)罱Y(jié)構(gòu)頂部時(shí)，y向雷諾正應(yīng)力相比光滑表面在0＜y+＜70的區(qū)域明顯下降;測(cè)量位置在脊?fàn)罱Y(jié)構(gòu)底部時(shí)，在0＜y+＜10時(shí)，y向雷諾正應(yīng)力高于光滑表面，這是由脊?fàn)罱Y(jié)構(gòu)附近形成的流向渦引起的，而10＜y+＜70時(shí)，曲線與測(cè)量位置在頂部的曲線重合。而增阻情況(s+=39.4)下y向雷諾正應(yīng)力相比光滑表面明顯增大。

圖10 y向雷諾正應(yīng)力分布(s+=18.2)Fig.10 Normal stress distribution in y direction(s+=18.2)

圖11 y向雷諾正應(yīng)力分布(s+=39.4)Fig.11 Normal stress distribution in y direction(s+=39.4)

圖12和圖13分別給出了減阻和增阻情況下邊界層近壁區(qū)z向雷諾正應(yīng)力沿y方向分布曲線。兩種情況下z向雷諾正應(yīng)力相比光滑表面均明顯下降，峰值大小和位置沒有明顯差別。脊?fàn)畋砻鎧向雷諾應(yīng)力的減小表明脊?fàn)罱Y(jié)構(gòu)的存在確實(shí)對(duì)原本紊亂的近壁區(qū)流動(dòng)有了橫向限制，使流動(dòng)的橫向脈動(dòng)降低。但是減阻情況下向雷諾正應(yīng)力減小的區(qū)域?yàn)?＜y+＜70，增阻情況下這一區(qū)域?yàn)?＜y+＜50，原因在于增阻情況(s+=39.4)下，脊?fàn)罱Y(jié)構(gòu)內(nèi)流向渦規(guī)模、強(qiáng)度均比減阻情況(s+=18.2)下大。

圖12 z向雷諾正應(yīng)力分布(s+=18.2)Fig.12 Normal stress distribution in z direction(s+=18.2)

圖13 z向雷諾正應(yīng)力分布(s+=39.4)Fig.13 Normal stress distribution in z direction(s+=39.4)

2.4 雷諾剪應(yīng)力分布曲線

圖14和圖15分別給出了減阻和增阻情況下邊界層近壁區(qū)雷諾剪應(yīng)力沿y方向分布曲線?？梢悦黠@看出，減阻情況下，脊?fàn)畋砻胬字Z剪應(yīng)力較光滑表面在0＜y+＜70區(qū)域顯著下降;增阻情況下，測(cè)量位置在頂部時(shí)，脊?fàn)畋砻胬字Z剪應(yīng)力較光滑表面明顯變大，測(cè)量位置在底部時(shí)則在y+＞40的區(qū)域明顯變大。

2.5 壁面剪應(yīng)力分布云圖

圖14 雷諾剪應(yīng)力分布(s+=18.2)Fig.14 Shear stress distribution(s+=18.2)

圖15 雷諾剪應(yīng)力分布(s+=39.4)Fig.15 Shear stress distribution(s+=39.4)

圖16 脊?fàn)罱Y(jié)構(gòu)壁面剪應(yīng)力分布云圖(s+=18.2)Fig.16 Contour of wall shear stress distribution(s+=18.2)

圖17 脊?fàn)罱Y(jié)構(gòu)壁面剪應(yīng)力分布云圖(s+=39.4)Fig.17 Contour of wall shear stress distribution(s+=39.4)

圖16和圖17分別給出了減阻和增阻情況下脊?fàn)罱Y(jié)構(gòu)壁面剪應(yīng)力分布云圖。由圖可見，兩種情況下脊?fàn)罱Y(jié)構(gòu)底部區(qū)域(深藍(lán)色)剪應(yīng)力都較小，剪應(yīng)力較大的部分都集中在頂部區(qū)域(紅色)，這表明脊?fàn)罱Y(jié)構(gòu)底部流動(dòng)較為安靜，而頂部流動(dòng)則十分劇烈，存在很強(qiáng)的動(dòng)量交換，這從圖4和圖5所示的脊?fàn)罱Y(jié)構(gòu)附近流線圖也可以看出。

3 脊?fàn)畋砻嫱牧髁鲌?chǎng)特性分析

眾多學(xué)者的理論及試驗(yàn)研究表明，湍流流動(dòng)中近壁面的流向渦對(duì)于壁面剪應(yīng)力有著決定性的影響。這些流向渦會(huì)導(dǎo)致強(qiáng)烈的下掃運(yùn)動(dòng)，使流動(dòng)在壁面法向和展向方向產(chǎn)生動(dòng)量交換，削弱了沿流動(dòng)方向的能量。

本文的數(shù)值模擬結(jié)果反映出，在脊?fàn)罱Y(jié)構(gòu)附近產(chǎn)生了相對(duì)規(guī)則的流向渦，一個(gè)脊?fàn)罱Y(jié)構(gòu)內(nèi)存在一對(duì)尺寸相當(dāng)、反向旋轉(zhuǎn)的流向渦對(duì)，單個(gè)流向渦尺寸大約是脊?fàn)罱Y(jié)構(gòu)尺寸的一半。筆者認(rèn)為這些流向渦是脊?fàn)畋砻鏈p阻或者增阻的重要原因，并從邊界層近壁面各方向雷諾應(yīng)力分布以及壁面剪應(yīng)力分布等方面對(duì)脊?fàn)畋砻媪鲌?chǎng)特性進(jìn)行了分析。

圖10和圖12反映出脊?fàn)罱Y(jié)構(gòu)的存在限制了流向渦的尺寸以及流動(dòng)的橫向脈動(dòng)，減小了壁面法向及展向的動(dòng)量交換，流動(dòng)方向能量轉(zhuǎn)化為法向及展向能量的比例減少，從而產(chǎn)生了減阻效果。從本文得出的減阻效果隨脊?fàn)罱Y(jié)構(gòu)無量綱尺寸變化圖中可以看出，脊?fàn)罱Y(jié)構(gòu)的無量綱尺寸超過30之后，脊?fàn)畋砻姹憩F(xiàn)出增阻。通過對(duì)增阻情況下脊?fàn)畋砻娓鱾€(gè)參數(shù)的分析發(fā)現(xiàn)，脊?fàn)罱Y(jié)構(gòu)無量綱尺寸超過30時(shí)，流向渦的尺寸及渦量明顯增加，加上沾濕面積的增大，是可能出現(xiàn)減阻效果減弱甚至增阻現(xiàn)象的。

［1］WALSH M J.Riblets as a viscous drag reduction technique［J］.AIAA Journal，1983，21(4)∶485-486.

［2］WALSH M J.Turbulent boundary layer drag reduction using riblets［R］.Orlando:AIAA 20th Aerospace Sciences Meeting，1982.

［3］SUZUKI Y，KASAGI N.Turbulent drag reduction mechanism above a riblet surface［J］.AIAA Journal，1994，32(1)∶1781-1790.

［4］El-SSMNI O A，CHUN H H，YOON H S.Drag reduction of turbulent flow over thin rectangular riblets［J］.International Journal of Engineering Science.2007，45∶436-454.

［5］CHOI H，MOIN P，KIM J.Direct numerical simulation of turbulent flow over riblets［J］.Journal of Fluid Mechanics，1993，255∶503-539.

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

［7］劉喜寧.基于SIMPLE方法下的三維湍流流場(chǎng)數(shù)值分析［D］.［碩士學(xué)位論文］.西北工業(yè)大學(xué)，2005.

［8］劉占一，胡海豹，宋保維.不同間隔脊?fàn)畋砻鏈p阻數(shù)值仿真研究［J］.系統(tǒng)仿真學(xué)報(bào)，2009，21(19)∶6025-6028.

［9］潘家正.湍流減阻新概念的實(shí)驗(yàn)探索［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，1996，14(3)∶304-310.

［10］李新華，董守平，趙志勇.減阻溝槽表面近壁湍流相干運(yùn)動(dòng)的實(shí)驗(yàn)研究［J］.山東大學(xué)學(xué)報(bào)，2006，41(2)∶106-110，143.