国家数值风洞风雷软件分布式粗糙壁面边界层流动模拟

Distributed rough surface boundary layer flow simulation using NNW-PHengLEI software

摘要

针对粗糙元均匀分布壁面，基于国家数值风洞风雷软件，在k-ω SST湍流模型和γ-Re_θt 转捩模型中分别引入粗糙壁面边界条件修正及粗糙壁面转捩动量厚度雷诺数，同时增加了部分粗糙表面设置，实现了局部/全局粗糙壁面边界层流动模拟能力。算例研究表明，对于全湍状态翼型，粗糙壁面导致失速攻角减小4°，1×10^-3倍弦长量级的等效粗糙元高度使最大升力系数下降约36%。同时，粗糙壁面会使边界层转捩提前，导致翼型的摩擦阻力增大。计算分析表明，Hellsten等提出的粗糙壁面比耗散率边界条件能够较准确地反映粗糙度造成的影响，气动力计算结果与实验值基本吻合，流场模拟结果符合粗糙壁面流动特征；同时，其能够保证与转捩模型及其粗糙壁面修正的稳定耦合计算，对于转捩位置及其随粗糙度和来流湍流度变化趋势的模拟结果符合实验结果，但转捩区长度还需通过更严谨的实验进行进一步验证和修正。

Abstract

Based on the NNW-PHengLEI software，for surfaces with uniformly distributed rough elements，rough surface boundary condition corrections and transition momentum thickness Reynolds numbers were introduced into the k-ω SST turbulence model and the γ-Re_θt transition model respectively.At the same time，some rough surface settings were added to achieve local/global rough surface boun-dary layer flow simulation capability.Test results show that when an airfoil is in a fully turbulent state，rough surfaces can reduce the stall angle of attack by 4°，and an equivalent roughness height of 1×10^-3 of chord length can cause a decrease in maximum lift coefficient by approximately 36%.Additionally，rough surfaces can advance the boundary layer transition，leading to an increase in skin friction drag.The computational analysis proves that the rough surface boundary condition proposed by Hellsten et al.can reflect the roughness effect more accurately than the dissipation rate boundary condition.The computational aerodynamic forces are consistent with the experimental data，and the flow field simulation results conform to the flow characteristics of rough surfaces. In addition，Hellsten’s model makes it possible to ensure stable coupling simulation with the transition model and its rough surface correlation function.The simulation results for the transition onset location and its trend with roughness level and freestream turbulent intensity are consistent with the experimental data. Still，the transition length needs further verification and correction through more rigorous experiments.

Graphical abstract

关键词

边界层流动 / 粗糙壁面 / 湍流模型 / 转捩模型 / 国家数值风洞风雷软件

Key words

boundary layer flow / rough surface / turbulence model / transition model / NNW PHengLEI software

在飞行器研制过程中，虽然目前的表面加工工艺已经能够保证较低的粗糙度水平^［1-2］，但由于机翼、叶片等部件的工作环境复杂多变，在使用和维护过程中其表面容易出现微小冰晶/沙尘/污垢堆积、漆皮剥落和腐蚀等情况，导致运行状态偏离设计状态，造成层流边界层转捩显著提前、湍流边界层厚度增大，影响飞行器的气动特性。与光滑表面相比，粗糙壁面不仅会导致机翼阻力增大，还增加了边界层下游分离的可能性，使层流等外形设计的收益降低^［3-4］，特别是对注重飞行效率的小型通用航空飞行器影响更大。日本本田公司研制的HondaJet高速公务机采用了大量层流设计技术，为避免喷漆等非设计因素的干扰，还专门研发了一套精细化喷漆工艺和配套专用车间^［5］。但由于粗糙度的影响未形成广泛的定量认识，使用中的壁面粗糙状态亦存在不确定性，目前基于计算流体力学（computational fluid dynamics，CFD）的气动设计工程应用大都是面向光滑壁面开展的^［6］。另一个重要限制因素是，现有的雷诺平均N-S（Reynolds-averaged navier-stocks，RANS）湍流和转捩模型基本是基于光滑壁面假设发展的^［7-9］，虽然提出较多修正模型，但由于不同粗糙分布的影响规律较为复杂，精细化、工程化的粗糙壁面流动实验并不易开展，导致粗糙壁面边界层流动模拟方法的可信度还没有得到广泛的验证和认可。

粗糙壁面按粗糙元分布情况的不同可分为离散式和分布式两类。从飞行器使用角度看，均匀分布式粗糙元更为常见，且粗糙壁面实验也常以黏贴砂纸的方式进行，因此本文以均匀分布式粗糙壁面作为研究对象。从流动模拟角度出发，对均匀分布的微小粗糙元直接进行几何建模和网格生成是不现实的，因此需要通过物理模型引入粗糙度的影响。根据湍流边界层分层结构，可以通过粗糙壁面函数来调整对数律^［10］。ANSYS FLUENT软件^［11］即采用这种方式，并针对网格相关性进行了修正处理。但这种壁面函数方法只适用于附着流动，同时对匹配点的选择非常敏感^［8］。另一种思路则是对RANS湍流/转捩模型进行修正。Schlichting等^［12］和Ligrani等^［13］在粗糙壁面流动方面进行了大量的实验研究，得到了粗糙度影响下的湍流边界层壁面律表达式，为粗糙壁面边界层研究及模型修正奠定了基础。

围绕工程常用的湍流模型，Aupoix等^［14］介绍了由法国宇航院（ONERA）和美国波音公司针对Spalart-Allmaras（S-A）湍流模型^［7］提出的粗糙壁面流动修正方法。针对k-ω湍流模型，Wilcox^［8］利用Schlichting等^［12］的结论提出了不同粗糙元高度下的ω边界条件。Hellsten等^［10］针对k-ω 切应力输运（shear stress transport，SST）湍流模型^［9］也提出了相应的关联函数表达式。Knopp等^［15］指出了Wilcox粗糙壁面边界条件在网格和收敛性方面的不足，并提出了同时修改湍动能和耗散率边界条件的改进方案。刘通等^［16］以边界层速度亏损量为指标，对该模型的经验参数进行了重新标定。但现有研究中对于不同粗糙壁面修正模型的性能对比研究较少，且算例多以平板为主，对评价不同粗糙壁面修正模型的可信度及工程应用中的选用缺乏指导意义。

对于粗糙壁面诱导的边界层转捩，许多学者基于RANS方法，尤其是间歇因子类建模方法进行了相关研究。Stripf等^［17］将离散粗糙单元模型与双层k-ε模型相结合构造了粗糙壁面湍流边界层模型，并引入粗糙壁面动量厚度雷诺数作为转捩触发判据，间歇因子作为转捩区模型，构造了当地化的转捩预测方法，后续还进一步进行了拓展^［18］。Lorenz等^［19］在此基础上构建了新的粗糙壁面诱导bypass转捩预测准则。目前，以Langtry等^［20］提出的γ-Re_θt 转捩模型为代表的基于当地化经验关系式的转捩建模方法应用广泛。如何将已有研究成果引入γ-Re_θt 转捩模型及与已有湍流模型粗糙壁面修正相耦合都是值得讨论的问题。

本文将Hellsten等^［10］和Knopp等^［15］粗糙壁面边界条件分别引入k-ω SST湍流模型，并尝试将Stripf等^［17］粗糙壁面转捩判据与γ-Re_θt 转捩模型^［20］相结合，基于国家数值风洞风雷软件，开发局部/全局粗糙壁面边界层流动模拟方法，并通过平板和翼型验证算例定量分析粗糙壁面对流动和气动特性的影响，对比分析不同修正方法的模拟能力及粗糙壁面条件下转捩模型与湍流模型耦合的适用性，为工程选用和实验模拟提供一定参考。

1 粗糙壁面边界层流动建模

对于微小粗糙元均匀分布的粗糙壁面，若不考虑粗糙元形状而只考虑平均（等效）高度k_s的影响，则可以采用湍流边界层内层中的黏性长度尺度ν/u_τ 对k_s进行无量纲化^［21］，其中：u_τ =（τ_w / ρ）^1/2为壁面摩擦速度；ν为运动黏性系数；τ_w为壁面切应力。得到的无量纲等效粗糙元高度k{}_{\mathrm{s}}^{+}=k_s·u_τ /ν与距壁面的无量纲高度y⁺的意义类似，可以用来描述壁面粗糙度水平。根据k{}_s^{+}对湍流边界层速度剖面的不同影响，Schlich-ting等^［12］将粗糙壁面大致分为水力学光滑壁面（hydraulically smooth surface，k{}_s^{+} ≤5）、过渡粗糙壁面（transitionally rough surface，5<k{}_s^{+} ≤70）和完全粗糙壁面（fully rough surface，k{}_s^{+} >70）三类。上述两个阈值也有不同的提法，例如Ligrani等^［13］采用的是2.25和90。实验结果表明，当粗糙度非常小时，粗糙元被淹没在黏性底层之内，不会对边界层流动造成大的影响，黏性底层仍然存在；随着k{}_s^{+}的增大，流动成为完全粗糙流，黏性底层将消失^［21］，壁面处的雷诺应力不再为0（对于RANS方法即涡黏性系数不为0）。

1.1 粗糙壁面湍流模型修正

k-ω SST湍流模型^［9］将k-ω和k-ε模型相融合，在边界层内外通过开关函数进行切换，并将Bradshaw假设引入涡黏性系数的定义中，得到了广泛的发展和应用。由于原始k-ω SST模型的ω物面边界条件{\omega }_{\mathrm{w}}与黏性有关，如式（1）所示。

{\omega }_{\mathrm{w}}=\frac{\mathrm{60}{\nu }_{\mathrm{w}}}{\beta {\left(\mathrm{\Delta }y\right)}^{\mathrm{2}}}（1）

式中：Δy为壁面法向首层网格高度；ν_w为壁面附近的黏性系数，因此可以采用边界层内层特征速度u_τ 和黏性长度尺度ν/u_τ 定义^［21］，根据量纲分析有ω_w~u_τ² /ν。若将壁面{\omega }_{\mathrm{w}}与粗糙元高度k{}_{\mathrm{s}}^{+}相关联，新的壁面{\omega }_{\mathrm{w}}可表示为

{\omega }_{\mathrm{w}}=\frac{u_{\tau }^{\mathrm{2}}}{v}{S}_{\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{g}\mathrm{h}}（2）

式中：S_rough为k{}_{\mathrm{s}}^{+}的关联函数，根据粗糙度对边界层的影响，两者应成反比关系：当S_rough→∞时（k{}_{\mathrm{s}}^{+}很小），ω_w→∞，此时式（2）与原始物面边界条件式（1）是一致的；而当S_rough→0时（k{}_{\mathrm{s}}^{+}很大），ω_w为有限值，以此驱动雷诺应力的生成。

1.1.1 比耗散率边界条件修正模型

对比数值模拟^［8］和实验^［8，12］得到的粗糙壁面对数律式，可以得到S_rough与k{}_{\mathrm{s}}^{+}的关联关系。Wilcox^［8］针对k-ω湍流模型提出的S_rough关联函数为

S_{\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{g}\mathrm{h}}=\left\{\begin{array}{l}{\left(\mathrm{200}/k_{\mathrm{s}}^{+}\right)}^{\mathrm{2}},\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{k}_{\mathrm{s}}^{+}\le \mathrm{5}\\ \mathrm{100}/k_{\mathrm{s}}^{+}+\\ \left[{\left(\mathrm{200}/k_{\mathrm{s}}^{+}\right)}^{\mathrm{2}}-\mathrm{100}/k_{\mathrm{s}}^{+}\right]e^{\left(\mathrm{5}-k_{\mathrm{s}}^{+}\right)},\hspace{1em}{k}_{\mathrm{s}}^{+}>\mathrm{5}\end{array}\right.（3）

Hellsten等^［10］提出的关联关系为

S_{\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{g}\mathrm{h}}=\left\{\begin{array}{l}{\left(\mathrm{50}/k_{\mathrm{s}}^{+}\right)}^{\mathrm{2}},\hspace{1em}{k}_{\mathrm{s}}^{+}\le \mathrm{25}\\ \mathrm{100}/k_{\mathrm{s}}^{+},\hspace{1em}{k}_{\mathrm{s}}^{+}>\mathrm{25}\end{array}\right.（4）

同时，Hellsten等^［10］注意到了上述边界条件在完全粗糙壁面条件下使用，会在黏性底层区域启动涡黏性限制器（Bradshaw假设）这一不合理性，通过与Menter的讨论及参数优化向涡黏性限制器引入了额外的开关函数以确保其在黏性底层区域不会启动，获得了正确的粗糙壁面对数律^［21］。其中，式（5）为改进后的涡黏性系数ν_T限制器，式（6）为Hellsten等^［10］构造的开关函数。

{\nu }_{\mathrm{T}}=\frac{a_{\mathrm{1}}k}{\mathrm{m}\mathrm{a}\mathrm{x}(a_{\mathrm{1}}\omega ,\overline{S}{F}_{\mathrm{2}}{F}_{\mathrm{3}})}（5）

F_{\mathrm{3}}=\mathrm{1}-\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}{(\frac{\mathrm{150}\nu }{\omega y^{\mathrm{2}}})}^{\mathrm{4}}（6）

式中：a₁=0.31，a₁为Bradshaw假设常数；S为应变不变量；F₂为原始k-ω SST湍流模型^［9］的开关函数；y为距离壁面的高度。而湍动能k的物面边界条件仍保持原始k-ω SST湍流模型的设定，即k_w=0。

1.1.2 湍动能和比耗散率边界条件修正模型

Knopp等^［15］指出Wilcox的模型对于过渡粗糙表面的壁面摩擦系数预测并不十分准确，且Wilcox针对粗糙表面仍采用k_w=0的设定并不合理，导致了涡黏性系数与近壁面区域的理论值有偏差。为了实现对过渡粗糙表面的良好预测，Knopp等^［15］给出了另一个粗糙度湍流模型扩展，对k和ω的物面边界条件都进行了修正，如式（7）、（8）所示。

k=\frac{u_{\tau }^{\mathrm{2}}}{\sqrt[]{{\beta }^{\mathrm{*}}}}\mathrm{m}\mathrm{i}\mathrm{n}(\mathrm{1},\frac{k_{\mathrm{s}}^{+}}{\mathrm{90}})（7）

\omega =\mathrm{m}\mathrm{i}\mathrm{n}(\frac{u_{\tau }}{\sqrt[]{{\beta }^{\mathrm{*}}}\kappa {\tilde{k}}_{\mathrm{0}}},\frac{\mathrm{60}{\nu }_{\mathrm{w}}}{\beta {\left(\mathrm{\Delta }y\right)}^{\mathrm{2}}})（8）

其中

{\tilde{k}}_{\mathrm{0}}=\mathrm{0.03}{k}_{\mathrm{s}}{\phi }_{\omega }（9）

\begin{array}{l}{\phi }_{\omega }=\mathrm{m}\mathrm{i}\mathrm{n}\left[\mathrm{1},{\left(\frac{k_{\mathrm{s}}^{+}}{\mathrm{30}}\right)}^{\frac{\mathrm{2}}{\mathrm{3}}}\right]\mathrm{m}\mathrm{i}\mathrm{n}\left[\mathrm{1},{\left(\frac{k_{\mathrm{s}}^{+}}{\mathrm{45}}\right)}^{\frac{\mathrm{1}}{\mathrm{4}}}\right]•\\ \mathrm{m}\mathrm{i}\mathrm{n}\left[\mathrm{1},{\left(\frac{k_{\mathrm{s}}^{+}}{\mathrm{60}}\right)}^{\frac{\mathrm{1}}{\mathrm{4}}}\right]\end{array}（10）

式中：模型常数β^*=0.09；卡门常数κ=0.41。

1.1.3 程序实现

国家数值风洞风雷软件^［22］是由中国空气动力研究与发展中心面向全国发布的混合CFD计算开源平台，采用面向对象的大型软件设计理念，利用C++语言编写程序，具有良好的通用性和可扩展性。本文对风雷软件的修改主要集中在湍流结构网格求解器TurbSolverStruct中，即所有修改暂时只针对结构网格。参照风雷软件对公式进行无量纲化，例如修改k-ω SST湍流模型的比耗散率ω的边界条件，无量纲形式为

{\omega }_{\mathrm{w}}=\frac{u}{y{R}_e}{s}_{\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{g}\mathrm{h}}（11）

参考k-ω SST模型混合函数F₂的构造方式在程序中加入新的开关函数F₃，来实现粗糙壁面湍流模拟。

1.1.4 算例确认与参数校核

采用Mills等^［23］提出的粗糙平板摩擦系数半经验公式验证上述两种粗糙壁面边界条件植入程序的正确性，网格示意及边界条件设置如图1所示。平板长度为2 m，计算域高度为0.22 m，入口距平板前缘0.22 m，入口速度为100 m/s，来流单位雷诺数为1.5×10⁶。设置平板具有全长度粗糙壁面，无量纲粗糙元高度k{}_s^{+}分别为62、149、298、447。

使用3套网格对k{}_s^{+} =149的粗糙平板进行网格无关性验证，网格参数如表1所示。图2为不同网络计算的壁面摩擦系数对比。从图2可以看出，各条曲线基本重合。故下文粗糙壁面模拟均使用表1中的中等网格规模参数进行相应的网格构造。

由于CFD模型参数大多具有求解器相关性，而原始Hellsten模型对半经验公式曲线拟合效果不好，本文依据拟合程度将式（4）中的100调整为170，最终得到了图3的计算结果。从图3可以看出，不同粗糙度下两种粗糙壁面边界条件得到的摩擦系数分布与半经验公式趋势吻合较好，Knop模型在前半段与半经验公式更为贴合，而调整后的Hellsten模型在流动完全发展后更接近半经验公式，初步验证了两种粗糙壁面边界条件在风雷软件中实现的正确性，下文还将详细对比两种模型针对湍流边界层及耦合粗糙壁面转捩判据的使用效果。下文所有Hellsten粗糙壁面模型^［10］计算结果都是基于调整后的形式得到的。

1.2 粗糙壁面边界层转捩模拟

当边界层动量厚度雷诺数Re_θ 随流动增长到一个临界值Re_θt 时，边界层转捩就会发生，这是转捩模型构造的基本出发点。为了量化众多的转捩影响因素，转捩建模中需要引入经验关系式。原始γ-Re_θt 转捩模型主要考虑的是来流湍流度Tu和压力梯度λ_θ 的影响，因此转捩动量厚度雷诺数的经验关系式可写为Re_θt =f（Tu，λ_θ ），对于粗糙壁面需要引入粗糙度作为新的影响因素。

1.2.1 转捩模型粗糙壁面经验关系式修正

本文引入了Stripf等^［17］建立的粗糙壁面转捩动量厚度雷诺数，对于有量纲等效粗糙元高度k_s超过1%边界层位移厚度δ^*的情况，采用通过实验数据拟合的新关系式，即

\begin{array}{l}R{e}_{\theta t,\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{g}\mathrm{h}}\\ =\left\{\begin{array}{l}R{e}_{\theta t},k_{\mathrm{s}}/{\delta }^{\mathrm{*}}\le \mathrm{0.01}\\ {\left[\begin{array}{l}\frac{\mathrm{1}}{R{e}_{\theta t}}+\\ \mathrm{0.0061}{f}_{\mathrm{\Delta }}{(\frac{k_s}{{\delta }^{\mathrm{*}}}-\mathrm{0.01})}^{f_{{\sigma }_{\mathrm{t}\mathrm{u}}}}\end{array}\right]}^{-\mathrm{1}},k_{\mathrm{s}}/{\delta }^{\mathrm{*}}>\mathrm{0.01}\end{array}\right.\end{array}（12）

式中：Re_θt =f（Tu，λ_θ ）为原始γ-Re_θt 转捩模型的边界层转捩动量厚度雷诺数经验关系式。

有效湍流强度函数为

f_{{\sigma }_{\mathrm{t}\mathrm{u}}}=\mathrm{m}\mathrm{a}\mathrm{x}\left[\mathrm{0.9,1.61}-\mathrm{1.15}\mathrm{e}\mathrm{x}\mathrm{p}(-{\sigma }_{\mathrm{t}\mathrm{u}})\right]（13）

有效湍流度为

{\sigma }_{\mathrm{t}\mathrm{u}}=(R{e}_{\theta t}{/\mathrm{500})}^{-\frac{\mathrm{4}}{\mathrm{3}}}（14）

（15）

式（15）描述了不同粗糙单元间距（拓扑结构）的影响，式中：Δ _R =k_s /h_m，h_m 为粗糙元平均高程。由于在实际中遇到的粗糙度类型和分布可能会有显著变化，最大或等效粗糙元高度可能是唯一已知的参数，因此令f_Δ=1，由此可以得到粗糙元高度对边界层转捩的最大影响。

式（12）中还存在边界层位移厚度δ^*这一非当地变量，为使并行化CFD程序能够进行当地化求解，借鉴李虹杨等^［6］的思路，采用层流边界层速度剖面的经验公式^［24］进行求解。

\frac{\theta }{\delta }=\frac{n}{(\mathrm{1}+n)(\mathrm{2}+n)}（16）

\frac{{\delta }^{\mathrm{*}}}{\delta }=\frac{\mathrm{1}}{(\mathrm{1}+n)}（17）

式中：θ为边界层动量厚度；δ为边界层实际厚度；n为常数，对于Blasius层流边界层取为7。计算时首先采用Re_θt 的定义式和经验关系式迭代求解当地化边界层转捩动量厚度θ_t，再通过式（16）、（17）得到相应的边界层位移厚度δ^*。

1.2.2 程序实现

考虑粗糙度影响的转捩动量厚度雷诺数经验关系式为分段函数，接近光滑表面时仍使用原始经验关系式Re_θt，这也为设置局部粗糙度提供了便利。在原始γ-Re_θt 转捩模型程序基础上，将转捩判据改为分段函数Re_θt，rough，设置粗糙壁面起始与结束位置，通过网格节点坐标选择需要调用的转捩动量厚度雷诺数经验关系式及是否启动湍流模型粗糙壁面边界条件，以实现全局/局部粗糙壁面边界层流动模拟。最后，为了模拟粗糙壁面诱导转捩，需要同时考虑表面粗糙度对边界层转捩和湍流边界层发展的影响，这就涉及到粗糙壁面湍流模型边界条件和引入粗糙度影响后的转捩模型耦合使用的合理性。理论上看，两者对边界层流动模拟的作用分别体现在转捩后湍流的发展和转捩起始位置两个不同方面，即随着引入粗糙度影响的转捩模型发挥作用（触发转捩），湍流模型及其粗糙壁面边界条件的影响逐渐被引入转捩和湍流边界层；由于转捩过程中湍流模型未完全启动，因此粗糙壁面边界条件的作用也比较有限，主要影响充分发展的湍流边界层的流动特性。因此从模型原理出发，两种修正的作用阶段相对独立，耦合使用具有合理性，但仍需通过算例进一步验证说明。

2 算例应用与结果分析

本文计算均基于结构网格，采用风雷软件LU-SGS时间推进方法求解RANS方程和k-ω SST湍流模型、γ-Re_θt 转捩模型方程，空间离散采用有限体积法二阶Roe格式。有关控制方程、数值算法的具体形式及计算参数设置可参见风雷软件手册^［25］。本文所有算例均给定温度为室温15 ℃（即288 K）的来流条件。考察斯坦顿数的算例设为等温壁面（算例中将给出壁面温度），其他算例设为绝热壁面。

2.1 零压力梯度平板

本算例所用网格和边界条件与1.1.4小节相同，来流条件根据不同工况进行设置。

2.1.1 MSU粗糙平板湍流边界层

Hosni等^［26-27］研究了半球形粗糙沙粒对平板湍流边界层流动的影响，来流速度分别为12 m/s和58 m/s（对应来流单位雷诺数分别为0.8×10⁶和3.9×10⁶），壁面温度为300 K，等效粗糙度k_s分别为1.095 mm和0.293 mm，分别记为MSU（Mississippi State University）1~4平板。图4给出了采用两种不同湍流模型修正的计算结果与实验值的对比。从图4可以看出，粗糙度导致表面摩擦系数增大。与光滑平板相比，随着来流雷诺数（来流速度）的增大，边界层厚度减小，同一粗糙度水平（图4a和4b、图4c和4d）带来的摩擦系数增量有所降低；而对于同一来流状态（图4a和4c、图4b和4d），高粗糙度水平对应的壁面摩擦系数明显高于低粗糙度水平，k_s=1.095 mm时摩擦系数约为光滑平板的2倍。Hellsten模型得到的摩擦系数在各个条件下均大于Knopp模型，且整体来看与实验结果更为吻合，但在低速状态下（图4a和4c）湍流边界层发展初始阶段的摩擦系数模拟结果与实验数据有一定差距。

2.1.2 双尺度平板边界层转捩

Wang等^［28］、Abbott等^［29］开展了双粗糙度壁面边界层转捩实验，研究了上下游不同粗糙度水平对于边界层转捩的影响，如图5所示。其中平板前缘砂纸带长5 cm，下游砂纸带覆盖平板剩余表面，长度延伸至2 00 cm。表2给出了实验中使用的3种不同规格的砂纸参数。目数实际是筛网的开孔单位，指单位面积筛网上能够完整开孔的数量，即能开10个完整的孔为10目，能开100个完整的孔为100目。目数越大，孔的直径越小。对于砂纸上粗糙颗粒也可做类似定义。在取样长度内，R_a 为表面纵坐标绝对值的算术平均值，R_z 为表面纵坐标峰谷值的平均值^［30］。可见R_z 更接近于实际粗糙度水平的含义，故选择R_z 作为粗糙度输入参数。由于本文引入的粗糙壁面转捩动量厚度雷诺数经验关系式并未考虑双尺度粗糙度情形，因此暂不具备精确的变粗糙度模拟能力，接下来的算例仅针对实验中的光滑平板和局部粗糙平板进行计算。

平板边界层转捩计算条件如表3所示，其中U_inf为来流速度，Re_unit为来流单位雷诺数，Tu_inf为来流湍流度。计算中设置壁面温度为300 K，湍流涡黏性比为10。实验中给出的湍流度一般为物面前缘湍流度，而数值模拟设置的湍流度则为远场值。由于湍流模型在自由流中只有耗散机制没有生成机制，会造成来流湍流度衰减；加之计算远场较大，前缘湍流度将远低于实验水平，影响模拟精度。因此有必要开启自由流湍流度衰减抑制^［21］。

结果中每种工况都使用前缘砂纸目数/下游砂纸目数/来流湍流度×100进行标识。例如0/0/3.3表示前缘和下游（全部表面）均为光滑表面，来流湍流度为3.3%；60/0/0.5表示前缘为60目砂纸，下游为光滑表面，来流湍流度为0.5%；100/100/0.5表示前缘和下游（全部表面）均为100目砂纸，来流湍流度为0.5%。再依据表2和表3可分别查出不同砂纸目数和来流湍流度对应的粗糙度水平和来流条件。实验测量数据为斯坦顿数（Stanton number），其定义为

St=\frac{h}{\rho u{C}_p}（18）

式中：h为流体对流传热系数；ρ为密度；u为速度；C_p 为比热。斯坦顿数表征流动传递到流体的热量与流体热容量的比。由于层流和湍流的对流传热能力差别明显，斯坦顿数可以用来指示边界层转捩。不同湍流度下光滑平板边界层转捩计算结果与实验值的对比如图6a所示。从图6a可以看出，湍流度为0.2%时未发生转捩，随着湍流度的增大转捩位置前移。可见，采用原始k-ω SST湍流模型耦合γ-Re_θt 转捩模型得到的转捩位置和趋势与实验值基本一致，但层流和转捩区的斯坦顿数分布低于实验值。这可能是风雷层流边界层计算精度问题或斯坦顿数处理误差导致的。

来流单位雷诺数为3.34×10⁶、湍流度为0.5%状态下不同粗糙度平板（含光滑）的边界层转捩计算结果与实验值的对比如图6b所示，Hellsten模型和Knopp模型标识的结果均为与粗糙壁面转捩模型耦合模拟的结果，其中40和60目砂纸为前缘局部粗糙，100目砂纸为全表面粗糙。相比于光滑平板，粗糙壁面使得转捩位置明显前移，并随着前缘粗糙度增大（砂纸目数减小）前移量不断变大。两种耦合模拟在各粗糙度水平下均预测到了边界层转捩，但从转捩位置上看，耦合Hellsten模型的转捩模型模拟结果与实验吻合较好，但仍然存在层流和转捩区斯坦顿数计算值偏低的问题。相比之下，耦合Knopp模型得到的转捩位置雷诺数均延后了约2×10⁵，效果不佳。此外，由于Knopp模型对湍动能k的边界条件也进行了修正，导致与转捩模型耦合时计算稳定性变差。上述各状态均是在层流计算收敛结果的基础上进行续算得到的，Knopp模型耦合粗糙壁面转捩判据无法直接完成稳定计算。

该实验研究还表明在前缘粗糙带下游添加较小的第二粗糙带，在某些组合下会使边界层转捩延后，这可能与双尺度粗糙元高度变化对不稳定涡的抑制作用有关^［30］。

2.2 NACA65₂A215翼型

NACA65₂A215翼型实验由Ljungstroem^［4］开展，来流马赫数为0.182，翼型弦长1m，基于弦长的雷诺数为2.6×10⁶。实验中粗糙带覆盖了翼型整个上表面及下表面前缘0.15c（c为翼型弦长）的区域。自研计算网格（近壁区）如图7所示。周向布点513个并向前后缘加密，远场大小取为弦长的100倍，壁面法向首层网格高度为1×10^-5 m，保证在本算例计算状态下y⁺ <1，网格单元总数约为6.8万。选取k_s/c=1.54×10^-4和k_s/c=12.3×10^-4两种粗糙度状态进行全湍流状态模拟。

由于原始实验中光滑翼型的升力测量值偏低，因此改用Re=3.0×10⁶下NACA65₂A215翼型的实验结果^［29］，两种翼型低速条件下最大升力系数相差3%左右。两种粗糙壁面湍流模型修正的升力曲线计算结果与实验值的对比如图8所示。相比于光滑翼型，不同粗糙度水平均使翼型失速攻角减小4°左右，而较高的粗糙度水平会导致较大的升力损失（约36%）。对于过渡粗糙壁面（图8a），采用Hellsten模型的模拟结果更接近实验数据，Knopp模型虽然得到了与实验一致的失速攻角，但升力系数本身基本与光滑翼型计算结果一致，并未反映出粗糙度应有的影响；对于完全粗糙壁面（图8b），两种湍流模型修正的模拟结果接近且与实验吻合较好。图9进一步给出了8°攻角下不同粗糙翼型表面压力分布形态（粗糙壁面由Hellsten模型计算），可见过渡粗糙壁面翼型的压力分布形态与光滑翼型差别不大，但达到完全粗糙后，上下表面压力分布形态聚拢，导致升力系数显著减小。

图10为8°攻角下翼型附近涡黏性系数分布对比（粗糙壁面由Hellsten模型计算）。从图10可以看出，上表面前缘和下表面后缘附近近壁区涡黏性系数（即雷诺应力）随表面粗糙度变大而显著增大，直接导致摩擦阻力增大。同时，较大的雷诺应力向下游输运将使得边界层明显变厚，使流动在相同压力梯度下更容易发生分离。

2.3 NACA 0012翼型

采用风雷软件测试算例库中的结构化网格^［31］（图11），周向布点353个并向前后缘加密，远场大小取为弦长的10倍，壁面法向首层网格高度为1×10^-6m，保证在本算例计算状态下y⁺<1，网格单元总数约为5.8万。

2.3.1 前缘粗糙湍流边界层

Samadani等^［32］对NACA0012分布式粗糙前缘进行了实验研究，来流马赫数为0.2，翼型弦长为1 m，基于弦长的雷诺数为6×10⁶。实验在上下表面前缘0.08c的区域设置粗糙带，粗糙度k_s /c=1×10^-3。图12给出了全湍状态下的升力系数。相比于光滑翼型，前缘粗糙翼型失速攻角减小约4°，最大升力系数损失约37%。两种粗糙壁面湍流模型修正方法得到的升力曲线线性段结果基本相同，但在失速攻角附近，Hellsten模型与试验吻合更好。

2.3.2 局部粗糙边界层转捩

NACA0012翼型局部粗糙边界层转捩实验由Kerho等^［33］在伊利诺伊大学的亚音速风洞中完成，来流马赫数为0.1，弦长为0.534 m，基于弦长的雷诺数分别为0.75×10⁶和2.25×10⁶。翼型前缘不同位置处黏贴了粗糙带，粗糙颗粒高度为350 μm（k_s /c≈6.55×10^-4），粗糙带长度l、粗糙带起始位置s以及对应的粗糙带前后缘位置（沿弦长位置和沿表面位置）如表4所示。由于2.1.2节已经验证了Knopp模型与粗糙壁面转捩模型耦合效果不佳，且计算稳定性差，本节仅使用Hellsten模型耦合粗糙壁面转捩模型进行计算。计算中设置来流湍流度为0.1%，湍流涡黏性比为5，开启自由流湍流度衰减抑制^［21］。

通过图13给出的翼型上表面摩擦系数来判断转捩起始与结束位置，并通过箭头标注在实验结果堆积图中进行对比（图14）。从图14可以看出，对于低雷诺数状态（图14a），实验结果显示，表3中前3个粗糙带布置的边界层转捩位置和转捩区长度粗糙带长度和粗糙带位置变化较小。Hellsten模型耦合粗糙壁面转捩模型得到的转捩位置及其变化趋势与实验结果基本一致，但由于本文修正并未考虑粗糙度对转捩过程的影响，模拟得到的转捩区长度没有明显变化。对于算例4状态，计算转捩位置虽也较算例3略向上游移动，但与实验结果相差较大。与算例3相比，实验结果显示，粗糙带位置后移约0.03 c（s=24 mm）反而使转捩大幅度前移，这不但与边界层转捩源于上游不稳定性历史累计效应的认知相悖，也不符合实验本身^［33］给出的“转捩区域对粗糙带位置不敏感”的结论。实验结果有待进一步确认。

相比于低雷诺数状态，对于较高雷诺数状态（图14b），原始转捩模型捕捉到了光滑翼型转捩位置前移及转捩区长度变小。对于粗糙翼型，实验结果显示转捩发生在翼型前缘0.05c之前，且转捩区域长度均超过了40%c，转捩结束位置随粗糙带变长逐渐前移。计算结果反映出了转捩区域整体前移的趋势，但转捩位置和转捩区长度与实验结果差别仍然较大。

3 结论

本文基于国家数值风洞风雷软件，通过引入粗糙壁面湍流模型边界条件和转捩动量厚度雷诺数，实现了局部/全局粗糙壁面边界层流动模拟，主要结论如下：

（1）粗糙壁面会使边界层转捩提前发生，同时导致外形的摩擦阻力显著增大，近壁区雷诺应力的提前生成也使得湍流边界层变厚，加剧了较大攻角下的后缘分离。对于全湍状态翼型，粗糙度均造成了失速攻角4°左右的损失，k_s/c=1×10^-3量级的粗糙度将使最大升力系数下降约36%。

（2）整体来看，Hellsten等提出的针对k-ω SST湍流模型粗糙壁面修正能够较准确地反映粗糙度的影响（特别是对于过渡粗糙状态），全湍状态下气动力计算结果与实验值基本吻合，流场模拟结果符合粗糙壁面流动特征，同时仅改变比耗散率边界条件的方式也能够保证与转捩模型的稳定耦合计算，适合粗糙壁面边界层流动模拟工程应用。

（3）针对粗糙壁面的湍流模型和转捩模型修正的作用阶段和对模拟结果的影响相对独立，但需考虑耦合使用的计算稳定性问题。本文引入的Stripf等粗糙壁面转捩厚度雷诺数对于转捩位置及其变化趋势的模拟结果相对可信，但修正未涉及转捩区长度，同时基于风雷软件的实现可能存在层流区计算精度及斯坦顿数处理误差问题，还需后续研究解决。

粗糙壁面边界层流动模拟主要是为了处理过渡粗糙壁面和完全粗糙壁面，而在加工和实际使用中，应尽可能保证粗糙度水平达到水力学光滑。值得注意的是，双粗糙度平板流动实验表明^{［28，30］}，粗糙度变化可以对边界层转捩过程产生影响，达到流动控制的效果，相关机理和基于RANS的模拟方法值得深入探讨。

参考文献

[1] 丁明亮.航空钛合金加工方法及表面完整性控制研究［J］.中国机械，2023（26）：52-55.

[2] 张学军，陈冰清.航空典型金属材料增材制造组织、缺陷、表面、构型研究进展［J］.航空材料学报，2024，44（1）： 1-14.

[3] Roberts S K， Yaras M I.Boundary-layer transition over rough surfaces with elevated free stream turbulence［C］//ASME Turbo Expo：Power for Land，Sea，and Air.Vienna：ASME，2004：105-118.

[4] Ljungstroem B L G.Wind tunnel investigation of simulated hoar frost on a dimensional wing section with and without high lift devices［R］.Stockholm：Aeronautical Research Institute of Sweden，1972.

[5] Fujino M.Design and development of the hondaJet［J］.Journal of Aircraft，2005，42（3）：755-764.

[6] 李虹杨，郑赟，刘大响.粗糙壁面诱导的流动转捩数值模拟方法［J］.航空动力学报，2016，31（9）： 2251-2257.

[7] Spalart P， Allmaras S.A one-equation turbulence model for aerodynamic flows［C］//30th Aerospace Sciences Meeting and Exhibit.Reno：AIAA，1992：439.

[8] Wilcox D C.Turbulence Modeling of CFD［M］.3rd edition.California：DCW Indus-tries Inc.，2006.

[9] Menter F R.Two-equation eddy-viscosity turbulence models for engineering applications［J］.AIAA Journal，1994，32（8）：1598-1605.

[10] Hellsten A， Laine S.Extension of k-ω shear-stress transport turbulence model for rough-wall flows［J］.AIAA Journal，1998，36（9）：1728-1729.

[11] ANSYS，Inc.Ansys Fluent Theory Guide，Release 2024 R1［R］.Pennsylvania，USA：ANSYS，Inc.，2024.

[12] Schlichting H， Gersten K.Boundary-layer theory ［M］.9th edition.Berlin：Springer-Verlag，2017.

[13] Ligrani P M， Moffat R J.Structure of transitionally rough and fully rough turbulent boundary layers［J］.Journal of Fluid Mechanics，1986，162： 69-98.

[14] Aupoix B， Spalart P R.Extensions of the Spalart-allmaras turbulence model to account for wall roughness ［J］.International Journal of Heat and Fluid Flow，2003，24（4）： 454–462.

[15] Knopp T， Eisfeld B， Calvo J B.A new extension for k-ω turbulence models to account for wall roughness［J］.International Journal of Heat and Fluid Flow，2009，30（1）： 54-65.

[16] 刘通，蔡晋生，屈崑.壁面粗糙度湍流扩展模型及流动数值模拟［J］.航空动力学报，2018，33（8）： 1981-1989.

[17] Stripf M， Schulz A， Bauer H J.Modeling of rough wall boundary layer transition and heat transfer on turbine airfoils［J］.Journal of Turbomachinery，2008，130（2）：1139-1151.

[18] Stripf M， Schulz A， Bauer H J，et al.Extended models for transitional rough wall boundary layers with heat transfer part I： model formulations［J］.Journal of Turbomachinery，2009，131（3）：031017-031028.

[19] Lorenz M， Schulz A， Bauer H J.Predicting rough wall heat transfer and skin friction in transitional boundary layers：a new correlation for bypass transition onset［J］.Journal of Turbomachinery，2013，135（4）： 819-832.

[20] Langtry R B， Menter F R.Correlation-based transition modeling for unstructured parallelized computational fluid dynamics codes［J］.AIAA Journal，2009，47（12）： 2894-2906.

[21] 杜一鸣.涡黏性湍流模型修正与三维边界层转捩预测方法研究［D］.西安：西北工业大学，2021.

[22] 赵钟，何磊，何先耀.风雷（PHengLEI）通用CFD软件设计［J］.计算机工程与科学，2020，42（2）： 210-219.

[23] Mills A F， Hang X.On the skin friction coefficient for a fully rough hat plate［J］.Journal of Fluids Engineering，1983，105（3）： 364-365.

[24] 陈懋章.粘性流体力学基础［M］.北京： 高等教育出版社，2002.

[25] 孟丽媛，徐刚，万云博，等.风雷软件应用与开发指南［M］.绵阳：中国空气动力研究与发展中心，2023.

[26] Hosni M H， Coleman H W， Taylor R P.Measurements and calculations of rough-wall heat transfer in the turbulent boundary layer［J］.International Journal of Heat and Mass Transfer，1991，34（4/5）： 1067-1082.

[27] Hosni M H， Coleman H W， Garner J W，et al.Roughness element shape effects on heat transfer and skin friction in rough-wall turbulent boundary layers［J］.International Journal of Heat and Mass Transfer，1993，36（1）： 147-153.

[28] Wang T， Matthew C R.Effect of elevated free-stream turbulence on transitional flow heat transfer over dual scaled rough surfaces［J］.Journal of Heat Transfer，2005，127 （4）： 393-403.

[29] Abbott I， von Doenhoff A.Theory of Wing Sections［M］.2nd edition.New York： Dover Publications Inc，1958.

[30] Pinson M W， Wang T.Effect of two scale roughness on boundary layer transition over a heated flat plate part Ⅱ：boundary layer structure［J］.Journal of Turbomachinery，2000，122（2）： 308-316.

[31] A系列结构网格三维NACA0012翼型算例.红山开源平台风雷测试算例库 ［EB/OL］.（2022-09-05）［2023-03-20］

[32] ThreeD_NACA0012_SA_Struct_ 4 CPU/grid/%E7%BD%91%E6%A0 %BC%E 5%9C%B0%E5% 9D%80.txt.

[33] Samadani S， Morency F.Heat transfer correlations for smooth and rough airfoils［J］.Fluids，2023，8（2）： 66.

[34] Kerho M F， Bragg M B.Airfoil boundary-layer development and transition with large leading-edge roughness.AIAA Journal，1997，35（1）： 75-84.

基金资助

国家自然科学基金(12202284)

江西省飞行器设计与气动仿真重点实验室开放基金(EI202280263)

沈阳航空航天大学引进人才科研启动基金(22YB11)