柔性仿羽毛结构抑制边界层转捩的初步探索

叶正寅, 洪 正, 武 洁

(西北工业大学, 西安 710072)

0 引 言

对于常规民用飞机而言，阻力增加1%就会在一次飞行中额外多消耗3 T的燃油，在相同起飞重量的情况下，会缩短120 km的航程[1]；对于相同起飞距离，升阻比提高1%，意味着可以增加14个乘客[2]。在构成民用飞机阻力的成分中，表面摩擦阻力占到了50%[3]。层流边界层的摩擦阻力远小于湍流边界层。研究表明[3]: 飞机机翼、平尾、垂尾和短舱分别贡献了总阻力的18%、4%、3%和3%。如果层流边界层在这些部件表面上能够达到20%、30%甚至40%的比例，则飞机阻力分别能够减少8%、12%甚至16%。因此，推迟边界层的转捩在飞行器减阻方面具有很大的实际价值。

在目前已知的边界层控制手段中，有吹/吸气控制[4]、表面沟槽控制[5-7]、振荡射流控制[8-9]、振荡洛伦兹力控制[10]、展向变形的行波法控制[11-13]、流向变形控制[14-15]、展向粗糙带控制[16]、边界层燃烧[17]等。当然，在航空领域的工程应用方面，目前最主要的途径仍然是通过设计层流机翼的方法推迟转捩的发生[18]。为了实现层流的有效控制，边界层的控制方法是当前的重点研究方向。但是，一些手段在向实际应用推广的过程中存在很大的难度和代价，如展向行波法、洛伦兹力控制技术[13]。

层流的被动控制方法具有成本低的特点，上述沟槽控制转捩的方法、粗糙表面的方法都属于被动控制措施。事实上，人们甚至在高速流动中，也仍然在探索有效的被动控制方法[19-20]。

是否存在更直接的边界层控制思路？自然界是人类寻求灵感的源泉，正如文献[21]中所讲，人类在发展飞行器的历史进程中，从鸟类飞行中得到过多种启示,鸟类的羽毛经过千万年的进化，仍然蕴藏着许多奥秘。Tucher等人[22]测量了哈里斯“鹰”翼羽修剪前后在风洞中自由滑翔时的最小阻力，发现在7.3 m/s到15 m/s的速度之间，未修剪的“鹰”阻力只有修剪后“鹰”阻力的70%～90%。Parfitt等[23]发现，许多企鹅物种喙周围突出的羽毛区域平均减少了31%的阻力。

从我们对鸟类羽毛的直接观测中不难发现，当气流顺羽毛方向时流动是顺利的，可当气流从侧面流过羽毛时，羽毛的侧缘就会卷起来阻止当地展向流动(由于羽毛侧缘卷起时基本上沿流向排列，在来流方向的迎风面积很小，而沿展向方向的迎风面积大)。因此，对于流经表面的气流而言，当局部气流沿展向流动时，在展向方向阻力很大，流向阻力小得多。而且离壁面高度越低，阻力会越大，当高度超过一定值(对应羽毛侧缘卷起来的高度)时，阻力就没有了。也就是说，羽毛对近物面气流具有各向异性的阻力特征。而羽毛的流固耦合特性就是羽毛产生这种各向异性阻力的“激发开关”。

通过对边界层的大量研究[24-29]，人们发现存在着一种共性特点，即：边界层从层流转化为湍流的转捩过程，必然出现Λ涡和发卡涡，而这两种涡的“腿”就是一种倾斜的局部流向涡，它会诱导气流的展向流动，如果羽毛规则的柔性边缘因为展向流动而卷起，就会对局部展向流动有阻滞作用并抑制Λ涡和发卡涡的发展，从而抑制转捩的发生。

如前面所提到，研究者们利用被动或主动的措施来减小边界层的阻力，但大多是针对湍流区的减阻。然而，由于湍流边界层的阻力远大于层流边界层，减阻的效果没有层流控制来得更直接有效。鸟类高效的飞行效率必然对应着其羽毛表面有着极低的飞行阻力，从而可以推断层流流动在羽毛表面是占主导的。那么，观察到的羽毛表面的各向异性阻力特征是否对边界层从层流到湍流的转捩有抑制效果呢？在国家数值风洞项目的资助下，作者对柔性羽毛进行了流固耦合特性分析，旨在探索羽毛的这种柔性结构特征在抑制边界层转捩中的效果。

1 羽毛模型

羽毛的结构和排列方式非常复杂，建立等价于真实羽毛的模型是非常困难甚至不可能的。因此，抓住所关注的主要特征对模型进行简化是十分必要的。如引言中提到的，在本文中，羽毛的柔性体现在存在展向流动时羽毛侧缘会卷起，从而阻碍该展向流动，而没有展向流动时羽毛则不会卷起即不起作用。此外，考虑到真实羽毛的特点，展向流动应当可以穿过羽毛，即羽毛不能完全抑制展向流动。对于羽毛表面的这种各向异性的阻力特征，本文采用体积力的形式对其进行建模。为简单起见，羽毛侧缘卷起时对物面法向的影响忽略不计。

1.1 展向阻力模型

对于卷起的羽毛侧缘给当地展向流动带来的阻力，有以下几点考虑：

(1)展向脉动速度越大，那么这种阻力也越大；

(2)越靠近表面，羽毛越密，则阻力也随之增大；

(3)这种阻力只能使得当地展向脉动速度减小，但不能改变其方向。

羽毛阻碍展向脉动流动产生的阻力，概念上与黏性阻力类似，则基于(1)、(2)两点且类比黏性力的公式，可将沿展向的阻力模化为：

(1)

为了确定σy的值，本文假设在离表面最近、高度为h0的控制点处，即展向脉动抑制作用最强的位置，展向脉动在一个计算时间步内恰好被衰减到0。于是，在该位置处利用动量定理近似有

(ρfy)h=h0VΔt=0-ρv′V

(2)

式中：V是离散控制体体积，Δt是计算的时间步长。将式(1)代入式(2)，可得:

(3)

更一般的，在任意高度h处，根据动量定理整理可得:

(4)

式中：下标tn、tn+1分别指代当前时刻和下一时刻。式(4)表明，使用式(3)来确定σy，则任意高度处的展向脉动速度相对衰减率只与高度成反比，与速度大小无关。

假设有一均匀流动的展向脉动速度固定为v′=1，网格和参数采用下文中的设置，但壁面处理为滑移壁面。这样，流动的变化仅受到模型公式(1)的影响。图1给出了1、10、100时间步后任一位置处v′沿壁面高度的分布情况。从图1中可以看到，模型对展向流动的衰减作用,越靠近壁面越强。即使在离壁面稍远的位置处，多次的衰减也将使展向流动趋于滞止。显然，模型只能削减速度的大小直到v′=0，然后模型不再起作用，速度方向也不会因此改变。

图1 展向脉动在模型作用下的衰减特征

1.2 流向阻力模型

羽毛卷起的侧缘除了主要对展向脉动流动起抑制作用外，也会对沿流向的流动产生或多或少的阻碍。因此，从完整性的角度考虑，也需要对羽毛卷起带来的流向阻力进行建模。

对于流向阻力的体积力模型，有以下几点考虑：

(1)流向阻力的公式与展向形式保持一致；

(2)展向脉动流动越强，则卷起的特征越显著，由此带来的流向阻力也越大。

与展向阻力模型保持形式一致，则流向阻力的模型可写为:

(5)

式中：fx为流向的体积力，σx是待确定的流向阻力系数，u是瞬时速度在流向的分量。与展向阻力系数σy相关联，σx的确定按如下公式:

(6)

2 数值方法和计算设置

2.1 控制方程

使用守恒形式的三维N-S方程来描述流动，形式如下:

(7)

式中：Q=(ρ,ρu,ρv,ρw,ρe)T是守恒变量，v、w是瞬时速度在展向和壁面法向方向上的分量，e是气体的总能，S=(0,ρfx,ρfy,0,uρfx+vρfy)T是方程的源项，羽毛的各向异性阻力模型正是在此处起作用。Fc、Gc、Hc是分别沿x、y、z方向的对流通量，Fv、Gv、Hv则是相应的黏性通量。

使用半离散的格心型有限差分方法来数值求解N-S方程。为求得控制点处的通量导数，先基于格心处已知的原始变量，通过五阶DCS(Dissipative Compact Scheme)插值得到界面上的原始变量的左迎风值和右迎风值。然后，对于对流通量，根据得到的左右迎风值，使用Roe格式计算界面处的对流通量。对于黏性通量，左右迎风值取算数平均作为界面值，然后由通量公式得到。得到界面上的对流通量和黏性通量后，采用六阶紧致中心格式来计算通量导数。最后，利用显式的四步二阶龙格-库塔法来计算下一时刻格心处的流场变量值。五阶DCS和格心型有限差分法的具体公式，可参阅文献[30]，本文中DCS中控制耗散的参数取为0.5。

2.2 计算设置

计算域如图2 所示。基本流为Blasius相似性解，在计算域入口处需引入扰动。初始扰动由3个T-S(Tollmien-Schlichting)波组成，具体形式如下

图2 计算域示意图

(8)

入口处扰动的参数选取参考文献[31]中的设置，如表1所示。

表1 入口扰动参数

计算域的三维尺寸(以入口处边界层的位移厚度无量纲化)和网格的参数如表2所示。展向和流向网格都是平均分布的，沿壁面法向，网格点按如下公式进行拉伸。

(9)

表2 计算域尺寸和网格参数

3 结果分析与讨论

首先，模拟了平板上的边界层转捩流动，一方面可以验证所采用的数值方法的可靠性，另一方面提供了一个没有加入羽毛模型的基准参考。然后，加入羽毛模型后再进行平板边界层转捩的模拟，对比基准结果就可以得到羽毛对边界层转捩的影响。

3.1 平板边界层转捩

尽管湍流是随机且复杂的，但仍是有迹可循的。无量纲的壁面高度定义为：

(10)

(11)

图3所示为无量纲的边界层速度分布。在入口x=0处，速度分布是典型的层流剖面，而在x=300处，速度分布已经演化为典型的湍流剖面。这说明使用的数值方法和计算设置可以反映层流到湍流转捩的过程。

图3 层流区和湍流区的速度剖面

壁面摩擦系数Cf的定义为：

(12)

对于层流，摩擦系数为

(13)

式中：δ是当地的边界层位移厚度。Ducros等人[33]提出了一种经验公式来描述湍流边界层摩擦系数的空间演化规律，具体形式为：

(14)

式中：l′=125δ0。图4中给出了平板上转捩流动，层流以及Ducros经验公式的壁面摩擦系数。从图中可以看到，大约从x=80开始，平板边界层的摩擦系数开始快速偏离层流对应的值，该点可视为层流向湍流转捩的起始点。摩擦系数在x=190达到峰值，之后流动进入完全湍流状态。湍流区的摩擦系数与经验公式的结果所在的区间是一致的。

图4 沿壁面的摩擦系数分布

平板边界层从层流到湍流的演化过程中，存在一些典型的大尺度相干结构。利用Q准则，图5呈现了转捩过程中不同阶段的涡结构，用距离壁面的高度进行了着色。为了更好地显示效果，展向通过复制延伸了一倍。可以看到转捩开始后先是形成Λ涡，然后是Λ涡演化成发卡涡，进而生成环状涡，并最终发展成环状涡链结构。对典型涡系结构的清晰捕捉也体现了数值模拟的可靠性。

(a)t=1.6T0

3.2 羽毛模型对平板转捩的影响

为了研究羽毛的各向异性阻力特征对转捩的影响，将第1节中建立的羽毛阻力模型加入N-S方程，其余设置均保持相同。流向阻力模型中的控制参数m考虑三种情况：(1)m=∞，对应羽毛带来的流向阻力忽略不计的情形；(2)m=1.5，对应流向阻力比展向阻力小的情形；(3)m=1，对应流向和展向阻力大小相等的情形。实际上，羽毛卷起的侧缘能够影响的高度是有限的。根据前面的平板边界层模拟结果，模型起作用的范围限制在z

图6给出了加入模型后计算得到的全流场的涡系结构图，上一小节中未加模型的结果也一并给出作为对比。可以看到在加入羽毛的阻力模型后，开始形成大尺度涡的位置较无模型情况明显靠后。无模型情况下，x=100附近就已经形成了发卡涡，而加入模型且m=∞时，x=300处发卡涡才出现，并且转捩到湍流过程中形成的涡系要稀疏得多。不同的m值得到的结果有显著的差异。m=1.5与m=∞相比，开始形成典型涡结构的位置相差不大，但后续过程中形成的涡系结构要丰富一些。m=1情况下，开始形成涡的位置较m=∞明显提前，x=200附近已经可以看到形成的发卡涡了。此外，之后形成的涡系结构相较之下也要丰富得多。

图6 沿平板的涡系结构(从上到下：无模型，m=∞,m=1.5,m=1)

图7 展向脉动速度均方根沿平板的分布(从上到下：无模型，m=∞,m=1.5,m=1)

图8 平均湍动能沿平板的分布(从上到下：无模型，m=∞,m=1.5,m=1)

图9给出了这几种情况下的壁面摩擦系数。加入羽毛模型后，即使在转捩发生以后，壁面摩擦系数也没有增大很多，保持在与层流对应的摩擦系数相当的水平。特别地，在m=1的情况下，转捩之后的摩擦系数不但没有增大，反而下降得比层流时还要小。湍流区的摩擦系数比层流区的摩擦系数还要小，这与通常的认知是相悖的。实际上，在加入羽毛模型且不忽略流向阻力情况下，平板的阻力来源除了壁面的摩擦力外，还有羽毛阻碍流向流动产生的反作用力。因此，平板的总阻力系数应当写为:

图9 不同m值得到的摩擦系数分布

(15)

阻力系数的第I部分对应的正是壁面摩擦系数，第II部分对应的是羽毛阻碍流向流动带来的反作用力系数。

m=∞对应的是流向阻力小到可以忽略的情形，因此阻力仅来自壁面摩擦。对于m=1.5和m=1，图10给出了总阻力系数及其两部分系数的分布。转捩开始后生成的大尺度涡结构激发羽毛模型，可以看到，由此带来的反作用力阻力开始迅速增长并在湍流区维持较高水平。从总阻力系数来看，转捩之后阻力系数相比层流时显著增大，这样就与湍流区阻力比层流区阻力大的认知是一致的了。以阻力开始突增的点为转捩开始的位置，则m=1.5对应的转捩起始位置在x=240附近，而m=1对应x=160附近。此外，转捩之后与不加模型的平板阻力相比，m=1.5的总阻力仍然小一些，而m=1时反作用力带来的阻力占据了主导，总阻力不降反而增加了不止一倍。

图11 m=∞时的平均流向速度剖面

图12所示为m=1对应的不同位置的平均速度剖面。x=100、x=200处的速度分布仍与层流剖面相近，而在x=300、x=400处，速度剖面明显被修正。与m=∞不同的是，修正的速度分布并不是像湍流剖面那样逐渐变得饱满。特别地，修正的速度剖面在壁面处的梯度比层流剖面明显要小，这对应了图10中湍流区的壁面摩擦系数比层流还小的情形。壁面处缓慢的速度变化主要是由于m=1时模型带来的流向阻力削弱了平均流向流动，而这种削弱作用在壁面处最为强烈。另一方面，这也使得平均流动的速度分布与典型的层流速度分布和湍流速度分布有本质的差别。从图12中可以看到，层流速度分布或者湍流速度分布是不存在拐点的，而当m=1时，x=300、x=400处的速度分布显然是存在拐点的。拐点的出现使得流动的稳定性变弱，扰动更容易激发。这可能就是前面图6～图8m=1的结果中涡系结构更加丰富、脉动强度显著增强的原因所在。正因为湍流的增强，在大概z>1的高度范围内，x=300、x=400处平均速度分布的修正效果明显。

图12 m=1时的平均流向速度剖面

4 结 论

从观察发现，羽毛表面流向和展向流动的阻力特征并不相同。本文以体积力的形式对这种各向异性阻力进行了建模，并利用直接数值模拟的手段研究了其对平板边界层转捩的影响。研究得到如下主要结论：

1)准流向涡是边界层从层流转捩到湍流过程中典型的大尺度涡结构[34]，羽毛对展向流动的抑制作用可以削弱准流向涡，从而延缓了转捩的发生。即使边界层已经发展到了湍流状态，因为流向涡的削弱，导致湍流脉动强度显著降低，壁面摩擦系数降低到与层流相当的水平。

2)若考虑羽毛侧缘卷起对流向流动带来额外阻力，相较于无流向阻力时，模型带来的流向阻力使得转捩开始得更早且湍流脉动更强。但在湍流区的摩擦阻力却更小，甚至低于层流阻力。这是因为额外的流向阻力修正了来流方向的平均流动，使得壁面附近的速度变化得非常缓慢，速度梯度的减小了带来了壁面摩擦阻力的减小。此外，修正的速度型出现了拐点，导致流动的稳定性降低，从而促进了湍流的发展。

3)尽管模型的流向阻力带来了摩擦阻力的进一步降低，但总阻力却是增加的。流向阻力较大时，总的阻力甚至比典型湍流区阻力还大。因此，羽毛在展向和流向方向真实的力学特性关系对于减阻来说至关重要。

本文初步探究了羽毛对平板边界层转捩的影响，旨在为边界层层流控制提供一种新的思路，希望可以抛砖引玉，与广大研究者们共同深入探究。本文的模型及计算等方面仍存在很多不足之处，后续研究工作仍在进行当中。