基于積分法的粗糙壁面湍流邊界層發(fā)展預(yù)測(cè)

張濤，朱曉軍，彭飛，閔少松

海軍工程大學(xué)艦船工程系，湖北武漢430033

0 引 言

粗糙表面上的流動(dòng)在生產(chǎn)、生活中隨處可見(jiàn)，包括地貌上風(fēng)的運(yùn)動(dòng)、河床上水的流動(dòng)，以及人造表面（例如：船舶、飛機(jī)、管道系統(tǒng)等）上的流動(dòng)。粗糙度對(duì)壁面流動(dòng)的影響研究一直是流體動(dòng)力學(xué)中的熱點(diǎn)問(wèn)題，其準(zhǔn)確影響至今尚未完全確定［1］。從上世紀(jì)的Nikuradse 粗糙砂粒管道試驗(yàn)到著名的Moody 圖，許多學(xué)者都曾嘗試給出粗糙度影響的高精度工程計(jì)算方法。Schultz 等［2］通過(guò)試驗(yàn)，研究了砂粒以及船體污損中的生物粘膜［3］、絲狀海藻［4］和藤壺［5］等粗糙形態(tài)對(duì)壁面流動(dòng)的影響，并以此為基礎(chǔ)評(píng)估了海洋污損生物對(duì)艦艇航速、軸功率和經(jīng)濟(jì)性的影響［6］。Jimenéz［7］對(duì)以往有關(guān)粗糙度對(duì)湍流影響的研究成果予以了評(píng)述，詳細(xì)分析了粗糙形態(tài)的分類和粗糙度影響的機(jī)理，認(rèn)為粗糙壁面對(duì)湍流的影響受2 個(gè)參數(shù)的控制：粗糙度雷諾數(shù)ks+和粗糙度函數(shù)ΔU+。近年來(lái)，國(guó)內(nèi)開(kāi)展了一些有關(guān)脊?fàn)畋砻鏈p阻方面的試驗(yàn)［8］和數(shù)值計(jì)算［9］研究，取得了一些成果。

目前，求解湍流邊界層的方法有微分法和積分法2 種。微分法是直接對(duì)N-S 方程進(jìn)行求解，即CFD 方法。在目前的CFD 軟件中，粗糙壁面流動(dòng)的計(jì)算是采用壁面函數(shù)或者低雷諾數(shù)湍流模型等近壁面處理方法來(lái)對(duì)粗糙形態(tài)進(jìn)行建模，這2種方法或依賴經(jīng)驗(yàn)關(guān)系式，或增加粗糙元附近的網(wǎng)格密度，但由于經(jīng)驗(yàn)關(guān)系式的通用性較差，且增加網(wǎng)格密度并不意味著能提高計(jì)算精度；因此，粗糙壁面的流動(dòng)計(jì)算被認(rèn)為是CFD 的致命弱點(diǎn)之一［10］。

積分法發(fā)展較早，它利用動(dòng)量或能量積分方程展開(kāi)求解，在喪失一些湍流結(jié)構(gòu)信息的同時(shí)也避開(kāi)了粗糙壁面邊界層復(fù)雜的湍流尺度計(jì)算，較為簡(jiǎn)便，而且對(duì)于壓力梯度不是很大的湍流邊界層，積分法的結(jié)果并不比最好的微分法差［11］。Dvorak［12］利用積分法對(duì)均勻砂粒粗糙的翼型邊界層進(jìn)行了計(jì)算，計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)吻合較好，但其所使用的關(guān)聯(lián)方程中的表面摩擦力經(jīng)驗(yàn)公式只適用于二維均勻分布的砂粒，而且沒(méi)有考慮外層速度剖面對(duì)壁面律的偏離。本文將利用積分法對(duì)零壓力梯度下的棱錐形粗糙壁面湍流邊界層開(kāi)展研究，以粗糙度函數(shù)來(lái)表征粗糙度對(duì)平均流動(dòng)的影響，并結(jié)合動(dòng)量積分方程和壁面律進(jìn)行計(jì)算分析。

1 求解方法

1.1 粗糙度影響的表示

粗糙度對(duì)壁面流動(dòng)最主要的影響體現(xiàn)在改變近壁面流動(dòng)的平均速度分布［7］上，因此，可以把平均速度分布的下降量作為粗糙度對(duì)壁面流動(dòng)影響的表征，稱之為粗糙度函數(shù)ΔU+。

光滑壁面湍流邊界層的平均速度分布可用Coles 公式表示為

式中：y 為邊界層內(nèi)以壁面為原點(diǎn)的法向坐標(biāo)；U為y 處順流方向的平均速度；uτ為摩擦速度；Π為尾流參數(shù)，其與順流方向的壓力梯度和壁面粗糙情況有關(guān)；W 為尾流函數(shù)，代表著外層速度剖面對(duì)對(duì)數(shù)律的偏離；δ 為邊界層厚度；B0為光滑壁面對(duì)數(shù)律的截距；κ 為卡門(mén)常數(shù)，通常取為0.41；υ為流體的運(yùn)動(dòng)粘性系數(shù)。

由粗糙度函數(shù)ΔU+的定義，粗糙壁面湍流邊界層的速度分布可表示為

粗糙度函數(shù)ΔU+取決于粗糙的本質(zhì)特征和雷諾數(shù)，關(guān)于粗糙度函數(shù)的形式，一般認(rèn)為它與粗糙度參數(shù)成對(duì)數(shù)關(guān)系，可表示為

式中：h'為粗糙度參數(shù)；B 和C 為常系數(shù)，可根據(jù)具體的粗糙形式而定。

粗糙度參數(shù)的表示方法有多種，例如粗糙元的高度、間隔以及其形狀特征等。很多文獻(xiàn)將粗糙元的高度作為粗糙度唯一的表示參數(shù)，例如，Nikuradse 經(jīng)典試驗(yàn)和Moody 圖等，這種表示方法適用于粗糙元密集均勻分布的表面。然而，工程實(shí)際應(yīng)用中的粗糙表面通常都不是均勻分布的，因此僅使用粗糙元高度這一個(gè)參數(shù)明顯不足以對(duì)粗糙形態(tài)進(jìn)行充分描述。一般認(rèn)為，對(duì)于非均勻隨機(jī)分布的粗糙表面，可以使用統(tǒng)計(jì)方法進(jìn)行描述，統(tǒng)計(jì)量包括概率密度函數(shù)、均方根偏差、偏態(tài)和峰度等，例如，Townsin［13］在研究船體表面粗糙的阻力影響時(shí)，是使用粗糙元高度的譜矩作為粗糙度參數(shù)。本文是采用Newcastle University upon Tyne 的船舶性能組提出的粗糙度參數(shù)［14］：

式中：Rq 為采樣長(zhǎng)度上平均剖面的偏差均方根；Sa 為采樣長(zhǎng)度上輪廓平均斜率的絕對(duì)值。

1.2 粗糙邊界層的特性參數(shù)關(guān)系式

當(dāng)y=δ 時(shí)，U=Ue（其中Ue為自由流速），粗糙壁面的速度分布式（2）可轉(zhuǎn)化為

將式（5）減去式（2），可得粗糙壁面的虧損律

將式（7）代入邊界層位移厚度δ*和動(dòng)量厚度θ 的定義式中，可得邊界層特性參數(shù)關(guān)系式：

其中，

1.3 動(dòng)量積分方程

零壓力梯度下的二維動(dòng)量積分方程為

將式（8）代入式（10）中，并考慮到對(duì)于零壓力梯度下平板邊界層，可得

其中：

1.4 粗糙壁面方程

根據(jù)粗糙壁面邊界層的速度分布來(lái)推導(dǎo)粗糙壁面方程，由式（5），可得

將式（8）代入式（13）中并對(duì)x 求導(dǎo)，可得

其中：

聯(lián)立式（11）和式（14），可得到以dδ/dx 和dw/dx 為變量的線性方程組，利用Crammer 法則，可解得δ 和w 的二元常微分方程組，在已知粗糙度函數(shù)和流向初始點(diǎn)的速度分布的情況下，使用數(shù)值積分法便可求解出δ 和w 在邊界層順流方向上的變化情況。

2 計(jì)算模型

Schultz 等［15］對(duì)三維正四棱錐體的粗糙形式進(jìn)行了試驗(yàn)研究，其由3 種棱錐體高度和3 種棱錐體傾斜角共同組合出了9 種粗糙形態(tài)，如圖1 所示。每種粗糙形態(tài)分別在4 種流速（1，3，5，7 m/s）下進(jìn)行試驗(yàn)，在流向x=1.35 m 位置處進(jìn)行36 次速度分布測(cè)量。試驗(yàn)在美國(guó)海軍學(xué)院的高速水洞中進(jìn)行，采用激光多普勒測(cè)速儀（LDV）進(jìn)行測(cè)量。本文選用前6 種粗糙形態(tài)進(jìn)行計(jì)算分析。

圖1 粗糙形態(tài)示意圖（左上為主視圖，左下為俯視圖，右側(cè)為單個(gè)粗糙元的形狀）Fig.1 Demonstration of roughness texture（upper left：the front view;bottom left：the top view;right：the geometry of single roughness element）

根據(jù)粗糙度參數(shù)h'的定義式，圖1 所示粗糙形態(tài)的h'可表示為

式中，kt和α 分別為棱錐體的高度與傾角。

采用式（3）和式（16）對(duì)試驗(yàn)數(shù)據(jù)進(jìn)行擬合，求得B=0.32，C=-2.78，擬合結(jié)果如圖2 所示。考慮到Coles 速度分布中的尾流參數(shù)Π 與壁面粗糙情況有關(guān)，其在壁面粗糙情況下要比光滑情況下大得多，根據(jù)文獻(xiàn)［16］的研究可知，在一般粗糙壁面情況下，Π 取0.7 較為合適，因此，在下文的計(jì)算中Π=0.7。

圖2 粗糙度函數(shù)ΔU+的擬合情況Fig.2 Fitting of roughness function ΔU+

3 計(jì)算結(jié)果分析

采用4 階Runge-Kutta 方法對(duì)δ 和w 的二元常微分方程組進(jìn)行求解，積分起點(diǎn)δ0和w0（當(dāng)x=0.2 m 時(shí)）通過(guò)以下方式進(jìn)行確定：δ0通過(guò)基于冪次律的光滑平板經(jīng)驗(yàn)關(guān)系式進(jìn)行估算，見(jiàn)式（17），而w0則通過(guò)牛頓迭代法求解式（13）的零點(diǎn)來(lái)確定。

通過(guò)編程實(shí)現(xiàn)上述算法，計(jì)算得到4 種流速下6 種粗糙形態(tài)在x=1.35 m 處的局部摩擦系數(shù)Cf和邊界層厚度δ 值。將試驗(yàn)測(cè)量點(diǎn)處的計(jì)算結(jié)果與文獻(xiàn)［15］中的試驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比，并與光滑壁面下基于冪次律的邊界層厚度和局部摩擦系數(shù)公式（式（18））進(jìn)行對(duì)比，其結(jié)果如圖3 和圖4 所示。圖中3 條曲線分別表示6 種粗糙形態(tài)在4 種流速下文獻(xiàn)［15］中的試驗(yàn)值、本文的計(jì)算值以及光滑壁面在相應(yīng)流速下的理論值。

圖3 在試驗(yàn)測(cè)點(diǎn)處局部摩擦系數(shù)Cf 的計(jì)算值與試驗(yàn)值和光滑壁面理論值的對(duì)比Fig.3 Comparison of the calculated values，the experimental data and the theory values of smooth walls for local friction factor Cf at the measurement location

圖4 在試驗(yàn)測(cè)點(diǎn)處邊界層厚度δ 的計(jì)算值與試驗(yàn)值和光滑壁面理論值的對(duì)比Fig.4 Comparison of the calculated values，the experimental data and the theory values of smooth walls for boundary layer thickness δ at the measurement location

由圖3 可以看出：

1）在粗糙情況下，Cf的計(jì)算值與試驗(yàn)值吻合較好，平均絕對(duì)誤差為4.9%，而文獻(xiàn)［15］的試驗(yàn)中使用總應(yīng)力方法確定的摩擦速度uτ的誤差為±6%，可見(jiàn)本文采用的積分法準(zhǔn)確度較高。

2）相對(duì)光滑壁面，粗糙壁面的Cf要大得多，可見(jiàn)粗糙度對(duì)壁面摩擦阻力的影響很大。

3）在粗糙度的影響下，Cf不再隨自由流速的增加而減小，而是受流速和粗糙度的共同影響。

由圖4 可以發(fā)現(xiàn)：首先，在粗糙情況下，δ 的計(jì)算結(jié)果與試驗(yàn)值的變化趨勢(shì)吻合較好，兩者間的平均絕對(duì)誤差為7.1%，誤差相對(duì)于局部摩擦系數(shù)較大，這可能與δ 的積分起點(diǎn)采用光滑壁面的邊界層厚度有關(guān)；其次，可以發(fā)現(xiàn)在粗糙度的影響下，δ 不再隨自由流速的增加而減小，而是呈現(xiàn)出隨自由流速的增加而增加的趨勢(shì)，這說(shuō)明粗糙度的攝動(dòng)隨自由流速的增加而增強(qiáng)了。

將第2 種粗糙形態(tài)（kt=450 μm，α=45°）在自由流速為7 m/s 情況（即第8 次測(cè)量）下的局部摩擦系數(shù)Cf與邊界層厚度δ 沿流向的變化情況及光滑壁面進(jìn)行了對(duì)比，如圖5 和圖6 所示。

圖5 粗糙壁面與光滑壁面局部摩擦系數(shù)Cf 沿流向變化情況的對(duì)比Fig.5 Comparison of the local friction factor Cf between the rough walls and the smooth walls along the flow direction

圖6 粗糙壁面與光滑壁面邊界層厚度δ 沿流向變化情況的對(duì)比Fig.6 Comparison of the boundary layer thickness δ between the rough walls and the smooth walls along the flow direction

由圖5 可知，Cf在壁面粗糙和光滑情況下沿流向的變化趨勢(shì)基本相同，只是兩者的截距不同。由圖6 可知，δ 在壁面粗糙和光滑情況下由于采用了相同的邊界層厚度經(jīng)驗(yàn)公式，因而兩者起始點(diǎn)的值相同，但在粗糙壁面，δ 沿流向的增加速度更快。

4 結(jié) 語(yǔ)

本文采用邊界層積分法對(duì)零壓力梯度下棱錐體形態(tài)的粗糙壁面湍流邊界層進(jìn)行計(jì)算，并與試驗(yàn)數(shù)據(jù)和光滑壁面的經(jīng)驗(yàn)公式進(jìn)行了對(duì)比，計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)吻合較好，說(shuō)明采用積分法求解粗糙壁面湍流邊界層精度較高?？梢?jiàn)，積分法可以較準(zhǔn)確、方便地預(yù)測(cè)出粗糙壁面湍流邊界層厚度和局部摩擦系數(shù)沿流向的發(fā)展，可為研究粗糙度對(duì)湍流邊界層的影響提供適于工程應(yīng)用的方法。在本文研究的基礎(chǔ)上考慮壓力梯度的影響因素，還可對(duì)粗糙螺旋槳葉切面湍流邊界層的發(fā)展展開(kāi)研究。

［1］TAYLOR J B，CARRANO A L，KANDLIKAR S G.Characterization of the effect of surface roughness and texture on fluid flow—past，present，and future［J］.In?ternational Journal of Thermal Sciences，2006，45（10）：962-968.

［2］SCHULTZ M P，F(xiàn)LACK K A. Turbulent boundary lay?ers over surfaces smoothed by sanding［J］. Journal of Fluids Engineering，2003，125：863-870.

［3］SCHULTZ M P，SWAIN G W. The influence of bio?films on skin friction drag［J］. Biofouling，2000，15（1/3）：129-139.

［4］SCHULTZ M P. Turbulent boundary layers on surfaces covered with filamentous algae［J］. Journal of Fluids Engineering，2000，122（2）：357-363.

［5］SADIQUE J，YANG X，MENEVEAU C，et al. Flow over barnacles characterization of barnacle geometry?and some initial flow characteristics［C］// 66th Annual Meeting of the APS Division of Fluid Dynamics，2013.

［6］SCHULTZ M P，BENDICK J A J A，HOLM E R，et al. Economic impact of biofouling on a naval surface ship［J］.Biofouling，2011，27（1）：87-98.

［7］JIMENéZ J. Turbulent flows over rough walls［J］. An?nual Review Fluid Mechanics，2004，36：173-196.

［8］胡海豹，宋保維，劉占一，等.基于湍流邊界層時(shí)均速度分布的脊?fàn)畋砻鏈p阻規(guī)律研究［J］. 航空動(dòng)力學(xué)報(bào)，2009，24（5）：1041-1047.HU Haibao，SONG Baowei，LIU Zhanyi，et al. Re?search about the characteristic of drag reduction over riblets surface based on average velocity profile in the turbulent boundary［J］. Journal of Aerospace Power，2009，24（5）：1041-1047.

［9］劉志華，董文才，夏飛. V 型溝槽尖峰形狀對(duì)減阻效果及流場(chǎng)特性影響的數(shù)值分析［J］.水動(dòng)力學(xué)研究與進(jìn)展（A 輯），2006，21（2）：223-231.LIU Zhihua，DONG Wencai，XIA Fei. The effects of the tip shape of V-groove on drag reduction and flow field characteristics by numerical analysis［J］. Chi?nese Journal of Hydrodynamics（Ser. A），2006，21（2）：223-231.

［10］PATEL V C. Perspective：flow at high Reynolds num?ber and over rough surfaces—Achilles heel of CFD［J］. Journal of Fluids Engineering，1998，120（3）：434-444.

［11］陳懋章.粘性流體動(dòng)力學(xué)基礎(chǔ)［M］.北京：高等教育出版社，2002.

［12］DVORAK F A. Calculation of turbulent boundary lay?ers on rough surfaces in pressure gradient［J］. AIAA Journal，1969，7（9）：1752-1759.

［13］TOWNSIN R L. The correlation of added drag with surface roughness parameters［J］. Fluid Mechanics and Its Applications，1991，6：181-191.

［14］SVENSEN T E，MEDHURST J S. A simplified meth?od for the accessment of propeller roughness penalties［J］.Marine Technology，1984，21：41-48.

［15］SCHULTZ M P，F(xiàn)LACK K A. Turbulent boundary layers on a systematically varied rough wall［J］. Phys?ics of Fluids，2009，21：1-9.

［16］CASTRO I P. Rough-wall boundary layers：mean flow universality［J］. Journal of Fluid Mechanics，2007，585（8）：469-485.