基于积分法的粗糙壁面湍流边界层发展预测

张涛，朱晓军，彭飞，闵少松

海军工程大学舰船工程系，湖北武汉430033

0 引 言

粗糙表面上的流动在生产、生活中随处可见，包括地貌上风的运动、河床上水的流动，以及人造表面（例如：船舶、飞机、管道系统等）上的流动。粗糙度对壁面流动的影响研究一直是流体动力学中的热点问题，其准确影响至今尚未完全确定［1］。从上世纪的Nikuradse 粗糙砂粒管道试验到著名的Moody 图，许多学者都曾尝试给出粗糙度影响的高精度工程计算方法。Schultz 等［2］通过试验，研究了砂粒以及船体污损中的生物粘膜［3］、丝状海藻［4］和藤壶［5］等粗糙形态对壁面流动的影响，并以此为基础评估了海洋污损生物对舰艇航速、轴功率和经济性的影响［6］。Jimenéz［7］对以往有关粗糙度对湍流影响的研究成果予以了评述，详细分析了粗糙形态的分类和粗糙度影响的机理，认为粗糙壁面对湍流的影响受2 个参数的控制：粗糙度雷诺数ks+和粗糙度函数ΔU+。近年来，国内开展了一些有关脊状表面减阻方面的试验［8］和数值计算［9］研究，取得了一些成果。

目前，求解湍流边界层的方法有微分法和积分法2 种。微分法是直接对N-S 方程进行求解，即CFD 方法。在目前的CFD 软件中，粗糙壁面流动的计算是采用壁面函数或者低雷诺数湍流模型等近壁面处理方法来对粗糙形态进行建模，这2种方法或依赖经验关系式，或增加粗糙元附近的网格密度，但由于经验关系式的通用性较差，且增加网格密度并不意味着能提高计算精度；因此，粗糙壁面的流动计算被认为是CFD 的致命弱点之一［10］。

积分法发展较早，它利用动量或能量积分方程展开求解，在丧失一些湍流结构信息的同时也避开了粗糙壁面边界层复杂的湍流尺度计算，较为简便，而且对于压力梯度不是很大的湍流边界层，积分法的结果并不比最好的微分法差［11］。Dvorak［12］利用积分法对均匀砂粒粗糙的翼型边界层进行了计算，计算结果与试验数据吻合较好，但其所使用的关联方程中的表面摩擦力经验公式只适用于二维均匀分布的砂粒，而且没有考虑外层速度剖面对壁面律的偏离。本文将利用积分法对零压力梯度下的棱锥形粗糙壁面湍流边界层开展研究，以粗糙度函数来表征粗糙度对平均流动的影响，并结合动量积分方程和壁面律进行计算分析。

1 求解方法

1.1 粗糙度影响的表示

粗糙度对壁面流动最主要的影响体现在改变近壁面流动的平均速度分布［7］上，因此，可以把平均速度分布的下降量作为粗糙度对壁面流动影响的表征，称之为粗糙度函数ΔU+。

光滑壁面湍流边界层的平均速度分布可用Coles 公式表示为

式中：y 为边界层内以壁面为原点的法向坐标；U为y 处顺流方向的平均速度；uτ为摩擦速度；Π为尾流参数，其与顺流方向的压力梯度和壁面粗糙情况有关；W 为尾流函数，代表着外层速度剖面对对数律的偏离；δ 为边界层厚度；B0为光滑壁面对数律的截距；κ 为卡门常数，通常取为0.41；υ为流体的运动粘性系数。

由粗糙度函数ΔU+的定义，粗糙壁面湍流边界层的速度分布可表示为

粗糙度函数ΔU+取决于粗糙的本质特征和雷诺数，关于粗糙度函数的形式，一般认为它与粗糙度参数成对数关系，可表示为

式中：h'为粗糙度参数；B 和C 为常系数，可根据具体的粗糙形式而定。

粗糙度参数的表示方法有多种，例如粗糙元的高度、间隔以及其形状特征等。很多文献将粗糙元的高度作为粗糙度唯一的表示参数，例如，Nikuradse 经典试验和Moody 图等，这种表示方法适用于粗糙元密集均匀分布的表面。然而，工程实际应用中的粗糙表面通常都不是均匀分布的，因此仅使用粗糙元高度这一个参数明显不足以对粗糙形态进行充分描述。一般认为，对于非均匀随机分布的粗糙表面，可以使用统计方法进行描述，统计量包括概率密度函数、均方根偏差、偏态和峰度等，例如，Townsin［13］在研究船体表面粗糙的阻力影响时，是使用粗糙元高度的谱矩作为粗糙度参数。本文是采用Newcastle University upon Tyne 的船舶性能组提出的粗糙度参数［14］：

式中：Rq 为采样长度上平均剖面的偏差均方根；Sa 为采样长度上轮廓平均斜率的绝对值。

1.2 粗糙边界层的特性参数关系式

当y=δ 时，U=Ue（其中Ue为自由流速），粗糙壁面的速度分布式（2）可转化为

将式（5）减去式（2），可得粗糙壁面的亏损律

将式（7）代入边界层位移厚度δ*和动量厚度θ 的定义式中，可得边界层特性参数关系式：

其中，

1.3 动量积分方程

零压力梯度下的二维动量积分方程为

将式（8）代入式（10）中，并考虑到对于零压力梯度下平板边界层，可得

其中：

1.4 粗糙壁面方程

根据粗糙壁面边界层的速度分布来推导粗糙壁面方程，由式（5），可得

将式（8）代入式（13）中并对x 求导，可得

其中：

联立式（11）和式（14），可得到以dδ/dx 和dw/dx 为变量的线性方程组，利用Crammer 法则，可解得δ 和w 的二元常微分方程组，在已知粗糙度函数和流向初始点的速度分布的情况下，使用数值积分法便可求解出δ 和w 在边界层顺流方向上的变化情况。

2 计算模型

Schultz 等［15］对三维正四棱锥体的粗糙形式进行了试验研究，其由3 种棱锥体高度和3 种棱锥体倾斜角共同组合出了9 种粗糙形态，如图1 所示。每种粗糙形态分别在4 种流速（1，3，5，7 m/s）下进行试验，在流向x=1.35 m 位置处进行36 次速度分布测量。试验在美国海军学院的高速水洞中进行，采用激光多普勒测速仪（LDV）进行测量。本文选用前6 种粗糙形态进行计算分析。

图1 粗糙形态示意图（左上为主视图，左下为俯视图，右侧为单个粗糙元的形状）Fig.1 Demonstration of roughness texture（upper left：the front view;bottom left：the top view;right：the geometry of single roughness element）

根据粗糙度参数h'的定义式，图1 所示粗糙形态的h'可表示为

式中，kt和α 分别为棱锥体的高度与倾角。

采用式（3）和式（16）对试验数据进行拟合，求得B=0.32，C=-2.78，拟合结果如图2 所示。考虑到Coles 速度分布中的尾流参数Π 与壁面粗糙情况有关，其在壁面粗糙情况下要比光滑情况下大得多，根据文献［16］的研究可知，在一般粗糙壁面情况下，Π 取0.7 较为合适，因此，在下文的计算中Π=0.7。

图2 粗糙度函数ΔU+的拟合情况Fig.2 Fitting of roughness function ΔU+

3 计算结果分析

采用4 阶Runge-Kutta 方法对δ 和w 的二元常微分方程组进行求解，积分起点δ0和w0（当x=0.2 m 时）通过以下方式进行确定：δ0通过基于幂次律的光滑平板经验关系式进行估算，见式（17），而w0则通过牛顿迭代法求解式（13）的零点来确定。

通过编程实现上述算法，计算得到4 种流速下6 种粗糙形态在x=1.35 m 处的局部摩擦系数Cf和边界层厚度δ 值。将试验测量点处的计算结果与文献［15］中的试验数据进行对比，并与光滑壁面下基于幂次律的边界层厚度和局部摩擦系数公式（式（18））进行对比，其结果如图3 和图4 所示。图中3 条曲线分别表示6 种粗糙形态在4 种流速下文献［15］中的试验值、本文的计算值以及光滑壁面在相应流速下的理论值。

图3 在试验测点处局部摩擦系数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 可以发现：首先，在粗糙情况下，δ 的计算结果与试验值的变化趋势吻合较好，两者间的平均绝对误差为7.1%，误差相对于局部摩擦系数较大，这可能与δ 的积分起点采用光滑壁面的边界层厚度有关；其次，可以发现在粗糙度的影响下，δ 不再随自由流速的增加而减小，而是呈现出随自由流速的增加而增加的趋势，这说明粗糙度的摄动随自由流速的增加而增强了。

将第2 种粗糙形态（kt=450 μm，α=45°）在自由流速为7 m/s 情况（即第8 次测量）下的局部摩擦系数Cf与边界层厚度δ 沿流向的变化情况及光滑壁面进行了对比，如图5 和图6 所示。

图5 粗糙壁面与光滑壁面局部摩擦系数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 可知，δ 在壁面粗糙和光滑情况下由于采用了相同的边界层厚度经验公式，因而两者起始点的值相同，但在粗糙壁面，δ 沿流向的增加速度更快。

4 结 语

本文采用边界层积分法对零压力梯度下棱锥体形态的粗糙壁面湍流边界层进行计算，并与试验数据和光滑壁面的经验公式进行了对比，计算结果与试验数据吻合较好，说明采用积分法求解粗糙壁面湍流边界层精度较高。可见，积分法可以较准确、方便地预测出粗糙壁面湍流边界层厚度和局部摩擦系数沿流向的发展，可为研究粗糙度对湍流边界层的影响提供适于工程应用的方法。在本文研究的基础上考虑压力梯度的影响因素，还可对粗糙螺旋桨叶切面湍流边界层的发展展开研究。

［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，FLACK 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］胡海豹，宋保维，刘占一，等.基于湍流边界层时均速度分布的脊状表面减阻规律研究［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 型沟槽尖峰形状对减阻效果及流场特性影响的数值分析［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，FLACK 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.