沿流向微结构沟槽流场直接数值模拟

李超群，李易,2,*，张晨曦，唐硕,2

1. 西北工业大学 航天学院，西安 710072 2. 陕西省空天飞行器设计重点实验室，西安 710072

对于航空器和水下航行器等，减阻设计是一项重要内容。尤其是对于飞机来说，减阻增升的气动设计也是评价其先进性的指标之一[1]。减阻方法分为主动减阻和被动减阻两种方式。相比于需要消耗能量的主动减阻方式，被动减阻不需要能量的消耗；其中，微型沟槽便是被动减阻的有效方式之一。NASA Langley研究中心的Walsh等[2-4]作为研究沿流向沟槽减阻原理的先驱者，先后对覆有不同形状沟槽的平板阻力进行了大量实验研究；研究发现，当沿流向沟槽无量纲高度h+<25并且沟槽无量纲宽度s+<30时，沟槽具有减阻效果，并当采用对称V型且s+=h+=15的沟槽时，减阻效果最佳，可达8%。沟槽的存在虽然增加了浸润面积，但是削弱了平板表面的湍流脉动，使阻力得到了降低。Sundaram等[5-6]对攻角为0°～12°的NACA0012翼型进行了实验，发现沟槽减阻效果为10%～15%。Subashchandar等[7-9]对GAW(2)翼型和大攻角下的NACA 0012翼型进行了实验，发现翼型表面附上沟槽壁面之后，其阻力明显减小，减阻效果最大可达12%。Coustols和Schmitt[10]对表面覆有沟槽薄膜的空客A320 1/11缩比模型进行了实验研究，发现当马赫数为0.7时减阻效果最佳，黏性阻力降低了4.85%，相应地，总阻力降低了1.6%。

胡海豹等[11]对不同来流速度下的沟槽表面湍动能分布律进行了实验测量，指出沟槽结构主要影响边界层流场的黏性底层和过渡区。王晋军等[12-13]利用激光多普勒测速计与氢气泡流动显示技术对沟槽表面湍流边界层的带条结构进行了研究，结果发现沟槽平板的低速带条具有更好的直线性，横向流动得到抑制，壁面的湍流强度得到削弱，流动稳定性增强。同时，崔光耀等[14]对沟槽上方的湍流边界层进行了测量，发现减阻沟槽使得流向上的发卡涡数量明显减少。

Zhang等[15]使用了隐式大涡模拟方法对不同攻角下表面覆有沟槽的翼型进行了数值模拟，结果表明虽然翼型表面的压差阻力略微增大，但雷诺应力大幅减小，二者的综合作用使得阻力减小。周健等[16]结合已有实验数据对雷诺平均Navier-Stokes(RANS)方法的k-ω模型进行修改，以得到新的模化函数，使得新计算模型可在工程中对沟槽面减阻效果进行计算。

Choi等[17]采用二阶中心格式对空间进行离散，同时结合Crank-Nicolson时间推进格式对沿流向沟槽湍流进行了直接数值模拟(Direct Numerical Simulation, DNS)，模拟结果对比了沟槽和平板的速度型；结果表明沟槽限制了流向涡的展向运动，使得涡引起的洗流影响范围缩小，因而阻力得到了减小。同时Chu和Karniadakis[18]采用具有二阶精度的谱方法对微结构的层流和湍流流动进行了直接数值模拟，结果发现流向微沟槽仅针对湍流有较好的减阻效果。

对于没有强剪切的流动，二阶格式精度足够；但是，对于复杂的近壁湍流和壁面剪切流，有必要采用高阶格式[19]。为处理流场中存在间断的问题，Liu等[20]依据ENO(Essentially Non-Oscillatory)格式提出了WENO(Weighted Essentially Non-Oscillatory)格式，WENO格式在处理非连续或有间断的问题时具有良好的稳定性和鲁棒性。Zhang和Jackson[19]验证了5阶WENO格式结合低耗散、低色散Runge-Kutta(Low Dissipation and Dispersion Runge-Kutta scheme, LDDRK)方法求解不可压缩流动的方法；数值结果表明，针对不同的算例，计算结果与精确解、已有的数值或实验结果吻合较好，证明了该方法具有较高的精度和鲁棒性。

为更精确地解析壁面湍流，本文采用Jiang和Shu[21]提出的构建WENO格式的方法构建7阶WENO格式，同时应用分数步时间推进格式与LDDRK结合的方法求解不可压缩流动，并解析微结构沟槽壁面湍流，以探究其减阻的可能机理。

本文首先给出所采用的控制方程以及数值格式；然后，对比Kim等[22]的结果验证算法的正确性与精度；接着，给出槽道的计算模型以及边界条件设置；最后，列出计算结果及分析，并做相应的总结和展望。

1 数值方法

1.1 控制方程

通常，流体的运动可以用三维非定常的Navier-Stokes方程和能量方程来描述，但对于低速的无黏不可压流动，控制方程可以简化为动量方程和连续方程。由于采用的数值方法为有限差分法，因此采用的控制方程的形式为微分形式。为便于求解，使用槽道的半高度δ、槽道中线处速度Ul、密度ρl以及动力黏度μl将Navier-Stokes方程进行无量纲化。

这样，微分形式的无量纲化的非定常不可压缩流动的控制方程为

(1)

式中：ui为瞬时速度矢量的3个分量(x、y和z分量)，使用槽道中线处速度Ul进行无量纲化；t为时间变量；xj为笛卡尔坐标，使用槽道半高度δ进行无量纲化；p为压强，使用ρlUl2进行无量纲化；Re为雷诺数(基于槽道高度的1/2)，Re=ρlUlδ/μl。

1.2 数值格式

对于没有较强的剪切或间断的湍流流动，二阶中心格式的精度已经足够，但因主要针对沟槽表面的壁湍流和壁面剪切流进行模拟，有必要采用高精度格式对近壁面处的湍流进行模拟[19]。对于处理存在剪切流和光滑流动的流场，鲁棒性和计算效率优秀的WENO格式是一种较好的选择。

WENO格式是在高精度的ENO格式上改进的一种格式。Harten和Osher[23]提出ENO格式的思路是构造“模板(Stencil)”，并在诸多模板中挑选出最光滑的模板构造多项式；但这种方法舍弃了那些光滑性不佳的模板，造成了计算的浪费。在此基础上，Liu等[20]对ENO格式进行了改进，提出了WENO格式，WENO格式利用这些模板的加权组合可以在光滑区构造出更高阶差分表达式，这样提高了求解的精度和计算资源的利用率。随后，Jiang和Shu[21]对该WENO方法进行了进一步改进，构造了新的光滑度量因子，这样，格式的计算效率得到了进一步提高，这也就是目前广泛使用的WENO格式。

对于直接数值模拟的空间离散来说，5阶WENO格式的耗散并不足够低[24]，为了精确模拟微结构壁面的湍流流动，同时考虑计算消耗，使用7阶WENO格式对对流项进行离散，构造方法使用Jiang和Shu在文献[21]中提出的方法。

不妨令f代表速度矢量的任意分量ui，以正通量为例，构建7阶WENO格式，需要使用8个点并将其分成4个模板Sk(k=0,1,2,3)，每个模板中包含4个点：

(2)

(3)

表1 7阶WENO格式中系数的值Table 1 Values of in the 7th WENO scheme

注:l为模板内节点编号。

(4)

光滑度量因子ISk表达式为

(5)

将式(4)和式(5)代入式(3)中，再将式(3)代入式(2)中，即可得到关于正通量的7阶WENO格式。

对于时间推进格式，采取基于分数步的时间推进格式[26](式(6)和式(7))，同时结合文献[27]中给出的具有4阶精度的最优4-6 LDDRK格式(4和6分别表示第1和第2个时间步的放大因子)。这样，非定常流动的时间精度大幅提高[19]。

(6)

同时：

(7)

对于黏性项采用6阶中心差分格式即可，即

(8)

2 计算程序验证

为验证程序的正确性，使用了Kim等[22]的槽道流动结果进行验证，计算条件与文献[22]中的相同。槽道壁面流场的平均速度u+和均方根(Root-Mean-Square, RMS)速度脉动——urms、vrms和wrms对比如图1所示，对比结果表明二者吻合较好，证明该程序在求解该类问题时具有较好的精度和可靠性。

图1 本文程序与文献[22]计算结果对比Fig.1 Comparison of calculation results of proposed design with results in Ref.[22]

同时，为体现高精度格式的优势，图1还对比了高阶格式和2阶中心格式的速度型和均方根脉动速度曲线，可以发现：对于速度型，低阶格式在过渡区和对数区精度略低于高阶格式；对于脉动速度，低阶格式所得的脉动结果要低于高阶格式和文献[22]中的结果。

3 微结构沟槽建模

计算区域的槽道流动外形如图2所示，该模型中，下壁面为无滑移的V型沟槽平面，上壁面为无滑移平板，其中，沟槽平面和平板壁面所在位置分别为y=-δ和y=δ，这样，对比上下壁面的阻力即可得到表面沟槽的减阻效果[17-18,28]。

图2 微尺度沟槽槽道计算区域Fig.2 Computational domain of channel mounted with riblets

沟槽近壁面网格如图4所示，不同情况下网格的详细信息如表2所示。其中在流向(x方向)上为均匀网格，Δx+=0.8；在垂直方向(y方向)上采用非均匀的双曲正切分布，在壁面处Δy+=0.9；在展向(z方向)上采用均匀网格。

图3 计算区域尺寸Fig.3 Size of computational domain

图4 沟槽网格示意图Fig.4 Schematic of grids near riblets

表2 不同情况下网格参数Table 2 Parameters of grids in different cases

注:N1、N2和N3为沿x、y和z方向的格点数；Δx+、Δy+和Δz+分别为沿x、y和z方向网格的无量纲长度

虽然并没有进行网格无关性验证，但3个方向上的网格密度已经比文献[17]的DNS模拟采用的网格高，因此可以认为该模拟采用的网格量是足够的。

对于对称的V型沟槽壁面，控制沟槽形状的参数有两个——沟槽的宽度s+以及沟槽斜壁与底面的夹角α。模拟的沟槽的宽度s+范围为13～44，夹角α固定为60°，因此控制沟槽外形的参数只有宽度s+。流动沿x正方向，其关于槽道半高度的雷诺数Re=5 000。模拟时，基于最小网格尺度，将全局Courant-Friedrichs-Lewy(CFL)数控制为0.2，将无量纲时间步长控制为0.001；不同情况下，单个算例所计算的时间步至少为5×106。

对于边界条件，流向(x方向)和展向(z方向)上设置为周期性边界，上下壁面设置为无滑移壁面。

所有计算均在计算集群上进行，计算采用4个节点，集群每个节点CPU型号为E5-2670 v2 @ 2.50 GHz，24核，128 GB内存，节点间采用InfiniBand(IB)技术互连。

4 计算结果及分析

给出了不同尺度下沟槽的减阻效果，并着重对减阻(s+=21.9)和增阻(s+=44.0)情况进行细致分析。针对这两种情况，当速度场的统计特性趋于稳定时开始进行数据采集，共采集了500个无量纲时间单元(tUl/δ)内的数据并对其进行了分析。为演示其收敛特性，图5绘制了s+=21.9的沟槽与平板壁面在500个无量纲时间内壁面剪切率随时间的变化关系，可以发现此时二者表面的流动已经处于统计稳定状态。

分析数据后给出了减阻效果、沟槽表面和无滑移平板表面的平均速度型、自相关系数、湍流强度、瞬时速度场以及流向涡分布规律。

图5 沟槽壁面和平板壁面剪切率的时间变化曲线Fig.5 Time history of wall-shear rates at both riblets and flat plate

4.1 减阻效果

对于沟槽的减阻效果，用式(9)衡量：

(9)

式中：η为阻力减小比例；Dr为沟槽壁面的阻力；Df为对应的平板壁面的阻力。若η<0，则沟槽壁面具有减阻效果；若η>0，则沟槽壁面具有增阻效果。

经过数值模拟后，得不同尺寸沟槽的减阻效果如图6所示，其中所对比的实验数据来自Bechert等[30]的实验结果。可以发现阻力减小的最大值为9%左右，此时s+=16.0；当s+=44.0时，阻力增大了6.5%左右。

图6 不同尺寸沟槽减阻效果对比Fig.6 Comparison of drag reducing ratios in different cases

4.2 平均速度型

对于减阻(s+=21.9)和增阻(s+=44.0)情况，沟槽表面上方速度型如图7所示，其速度使用中线处速度进行无量纲化处理。可以发现，在相同y坐标上，沟槽底部上方的速度值大于沟槽顶部上方的速度值，但是随着y值的增大，二者之差越来越小，并且当y大于某个值时，二者速度型重合，这也意味着在此范围之外，流场各处受到沟槽不同部分的影响程度相同。

为详细显示沟槽上方的平均速度型，图8绘制了沟槽顶部、沟槽斜壁中点和沟槽底部上方的速度型情况。对比平板表面速度型可知：减阻情况下，沟槽表面速度型在对数区出现了明显的上移现象；增阻情况下，沟槽表面的平均速度型则有明显的下移。这也一定程度上印证了文献[31]中的相关结论。

图7 s+=21.9和s+=44.0情况下沟槽表面的平均速度型Fig.7 Mean velocity profiles over riblets in cases where s+=21.9 and s+=44.0

图8 s+=21.9和s+=44.0情况下沟槽顶部、斜壁中点以及底部上方的平均速度型Fig.8 Mean velocity profiles over peaks, midpoints and valleys of riblets in cases where s+=21.9 and s+=44.0

Wu等[31]指出，无论沟槽表面产生减阻效果还是增阻效果，其湍流-非湍流界面乃至整个流场均会受到沟槽的影响。以减阻情况为例，如图8(a) 所示，沟槽不同位置处上方的速度型曲线在y+=50附近汇合，这表明在y+<50的区域之内，沟槽的不同部分对流场中同一y+处流场的影响程度不尽相同；而在此区域之外，流场各处虽仍会受到沟槽的影响，但其受影响程度却趋于一致。

4.3 自相关系数

讨论的自相关性是脉动值的自相关性，自相关性代表了统计意义上的定量关系，并且如果速度变量关于位置是完全独立的，那么他们的自相关系数为0[32]。自相关性是用统计方法表示速度在不同位置之间的相关关系，不同位置处速度分量的联系程度可以用自相关系数表示，以x方向的速度分量u为例，其沿展向不同位置处的相关性系数定义为

(10)

空间自相关性函数的另一个性质为

R(x,y,∞)=0

(11)

该性质表明湍流场中相距很远的两点的随机变量几乎是统计独立的，几乎不存在相关性。因此，如果沿展向两个位置处速度相关性系数趋近0，那么就可以说明展向距离足够宽,可以使湍流流动得到充分发展。

通过对湍流流场内的数据计算可得其相关性系数曲线，以减阻情况(s+=21.9)下流场速度沿展向的自相关系数(如图9所示)为例。可以发现在沟槽近壁面处、槽道中线处和平板近壁面处，速度分量的相关性沿展向逐渐降低，最后降低到0附近。这表明流动的区域足够大,可以使得沿展向湍流脉动不具有相关性，同时，这也说明湍流流场已经充分发展。

图9 速度沿展向的相关性系数Fig.9 Correlations of spanwise two-point velocity

4.4 湍流强度

沟槽顶部和沟槽底部上方的湍流强度示意图如图10所示，RMS速度脉动使用壁面剪应力进行无量纲化处理，对比了减阻(s+=21.9)和增阻(s+=44.0)情况下的3个方向上的湍流强度。

在减阻情况下，沟槽上方沿x、y和z方向的湍流强度均有所降低，在展向上湍流强度最高均降低了约10%，对于流向湍流强度降低了约7%；与之相比，在增阻的情况下流向上沟槽顶部上方湍流强度明显高于平板，垂直方向上沟槽表面湍流强度与平板壁面接近，展向上沟槽表面湍流强度略高于平板壁面。可以认为沟槽产生减阻效果时，能较明显降低3个方向上的湍流强度。

图10 s+=21.9和s+=44.0情况下沟槽壁面湍流强度Fig.10 Turbulent intensity over riblets surface in cases where s+=21.9 and s+=44.0

4.5 瞬时速度场

给出减阻(s+=21.9)和增阻(s+=44.0)情况下某时刻流场速度u分量云图和瞬时速度矢量场，如图11和图12所示。

从图11(b)可以看出，当沟槽产生减阻效果时，大尺度流向涡位于沟槽之上且非常靠近沟槽顶部，同时与沟槽顶部频繁碰撞，可以形象地描述为流向涡“浮”于沟槽之上，并且由于沟槽的约束作用，流向涡的展向运动也受到了抑制，因此展向的湍流脉动也会得到削弱。

从图12(b)可以看出，当沟槽产生增阻效果时，流向涡大部分位于沟槽内部，并且与沟槽底部

图11 s+=21.9情况下沟槽表面流动的瞬时流场Fig.11 Instantaneous flow field over riblets in case s+=21.9

图12 s+=44.0情况下沟槽表面流动的瞬时流场Fig.12 Instantaneous flow field over riblets in case s+=44.0

的壁面相互碰撞。

4.6 顺流向涡分布规律

为更清晰地描述顺流向涡的分布规律，图13绘制了沟槽附近150个时刻的流向涡分布，流向涡的位置取流向涡涡量最大值所在的空间位置。

从图13中沟槽上方流向涡分布情况可以发现，减阻(s+=21.9)情况下，流向涡主要位于沟

图13 s+=21.9和s+=44.0情况下流向涡位置Fig.13 Positions of streamwise vortices in cases where s+=21.9 and s+=44.0

槽上方并且在沟槽顶部附近有较多的涡聚集，流向涡仅能与沟槽顶部的尖峰频繁碰撞而无法直接与沟槽内壁互相作用，一方面，这种行为诱导出了二次涡；另一方面，由于没有大量的涡进入沟槽底部，沟槽底部的流动便会比较平稳，亦即降低了壁面处的湍流脉动。相反，增阻(s+=44.0)情况下，一方面，流向涡会进入沟槽内部，造成流向涡直接与沟槽内壁频繁碰撞，增大了壁面的湍流脉动；另一方面，沟槽本身尺寸的增加导致浸润面积增加，这二者均造成了阻力增大。

为表明流向涡沿垂直方向的分布规律，定义Πp(y)为流向涡的数目[33]，图14绘制了平板表面上方减阻(s+=21.9)和增阻(s+=44.0)情况下沟槽顶部上方的流向涡的分布规律。

通过图14的曲线对比可以看出，沟槽的存在会影响其上方流向涡的数目：当沟槽具有减阻效果时，其上方的流向涡的数目减少；而当沟槽具有增阻效果时，其上方的流向涡的数目增多。结合图13的结果，可以认为对于减阻(s+=21.9)情况，沟槽的存在一方面减少了壁面上方的流向涡的数目，另一方面使得流向涡较集中地分布在沟槽顶部附近，与沟槽顶部频繁碰撞。

图14 s+=21.9和s+=44.0情况下流向涡数目沿垂直分布Fig.14 Population trend of streamwise vortices in wall units in cases where s+=21.9 and s+=44.0

5 结 论

采用高精度空间离散格式和时间推进格式求解沟槽表面的流动，同时研究了不同尺度沟槽对壁面湍流的影响以及减阻效果，得到如下结论：

1) 7阶WENO格式结合用LDDRK格式改进的分数步时间推进格式对于求解该类问题具有较高精度。

2) 减阻情况下，大尺度流向涡的尺度略大于沟槽的宽度，这使得大尺度涡只能与沟槽尖峰碰撞而无法进入沟槽之内，因而沟槽内部的流动较为平缓，沟槽表面的湍流脉动也会减弱；同时，沟槽也使得近壁面处顺流向涡的数目减少。因此，沟槽的存在虽然增大了壁面与流体之间的浸润面积，但是使其表面的湍流脉动得到了削弱，并使近壁面处流向涡的数目得到降低，最后其综合作用使得阻力减小。

3) 模拟发现，在沟槽上方的流场中存在某一高度阈值，当在此高度之内时，流场受沟槽不同部分的影响程度不同；而在此高度之外时，流场受沟槽不同部分的影响程度趋于一致。同时，对比增阻情况，减阻情况下的这一高度阈值较小。

经过数值模拟后，虽然得到了一些结论，但本文工作仍具有局限性，沟槽减阻机理仍然有待进一步研究，在未来的工作中将结合相应的实验等继续对沟槽表面的流动进行深入分析。