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

海軍工程大學艦船工程系，湖北武漢430033

0 引 言

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

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

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

1 求解方法

1.1 粗糙度影響的表示

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

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

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

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

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

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

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

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

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

當y=δ 時，U=Ue（其中Ue為自由流速），粗糙壁面的速度分布式（2）可轉化為

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

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

其中，

1.3 動量積分方程

零壓力梯度下的二維動量積分方程為

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

其中：

1.4 粗糙壁面方程

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

將式（8）代入式（13）中并對x 求導，可得

其中：

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

2 計算模型

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

圖1 粗糙形態(tài)示意圖（左上為主視圖，左下為俯視圖，右側為單個粗糙元的形狀）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）對試驗數(shù)據(jù)進行擬合，求得B=0.32，C=-2.78，擬合結果如圖2 所示。考慮到Coles 速度分布中的尾流參數(shù)Π 與壁面粗糙情況有關，其在壁面粗糙情況下要比光滑情況下大得多，根據(jù)文獻［16］的研究可知，在一般粗糙壁面情況下，Π 取0.7 較為合適，因此，在下文的計算中Π=0.7。

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

3 計算結果分析

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

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

圖3 在試驗測點處局部摩擦系數(shù)Cf 的計算值與試驗值和光滑壁面理論值的對比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 在試驗測點處邊界層厚度δ 的計算值與試驗值和光滑壁面理論值的對比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的計算值與試驗值吻合較好，平均絕對誤差為4.9%，而文獻［15］的試驗中使用總應力方法確定的摩擦速度uτ的誤差為±6%，可見本文采用的積分法準確度較高。

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

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

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

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

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

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

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

4 結 語

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

［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］胡海豹，宋保維，劉占一，等.基于湍流邊界層時均速度分布的脊狀表面減阻規(guī)律研究［J］. 航空動力學報，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 型溝槽尖峰形狀對減阻效果及流場特性影響的數(shù)值分析［J］.水動力學研究與進展（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］陳懋章.粘性流體動力學基礎［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.