基于转捩模型的位移边界层研究

邓家钰，王成恩

(上海交通大学机械与动力工程学院，上海 201100)

1 引言

在航空航天领域中，边界层的层流到湍流的转捩过程的研究是计算流体力学(CFD)领域的一个关键且困难的问题。工程上对转捩流动的模拟大多是采用由实验总结出来的经验关系式，国外研究人员提出的一些转捩经典模型，如Abu-Ghanam&Ahaw(AGS)模型、Steelant & Dik模型、Suzen & Huang(S-H)模型、Menter & Langry(M-L)[1-5]模型等。其中M-L转捩模型对之前提出的转捩模型进行修正，在基于SST k-ω[6]湍流模型的基础上，融入基于流场当地变量的转捩模型。该模型结合了转捩经验关系式和低雷诺数湍流模型的优势，提高了模型的通用性和准确性。

转捩现象与边界层外的流动密切相关，对于平板边界层来说，当自由流湍流强度低于1%发生自然转捩，即边界层失去稳定导致层流瓦解；当湍流强度大于1%的，边界层对于自由流脉动的感受而导致bypass转捩[7]。Jacobs&Du-bin[8]以及Brandt[9]、Schaltter[10]关于自由流湍流影响下的bypass转捩的数值分析研究做了详细的工作。

关于平板边界层的研究，G.B.Schubauer[11]和A.M.Sacill[12]进行了平板边界层转捩实验。国内学者董平[13]对经典的T3系列[12]平板边界层实验进行实验和仿真对比，顾金生[14]提出了不依赖经验系数的湍流模型，郭玉波等人[15]进行了粘性流场的数值模拟实验。

在平板边界层转捩的研究中，何克敏[16]所研究的低湍流度风洞中湍流度对平板边界层转捩影响的实验、徐国亮[17]和符松[18]等人对平板的bypass转捩进行的实验研究。J.H.M.Fransson等人[19]研究的受自由流湍流影响的平板转捩实验、Schrader[20]对自由流为湍流的横掠平板转捩研究的数值仿真、陈奕[21]对平板边界层转捩的模拟仿真，引发了本文的研究点：对于来流为自由流湍流的平板的位移边界层厚度随着自由流湍流强度改变有着什么样的变化规律。

本文利用商用CFD软件FLUENT的SST tran-sition模型，模拟了经典的平板边界层转捩实验T3A和T3B[13]，以及Schrader[20]所做的平板边界层数值模拟实验。分析模拟结果得出以下结论：在自由流湍流强度大于1%的前提下，随着湍流强度的减小，位移边界层厚度沿着流向方向的变化曲线在转捩区域出现下凹现象，且随着湍流强度减少，现象越来越明显；在自由流湍流强度小于1%的提前下，在湍流强度接近1%，位移边界层厚度沿着流向方向的变化曲线在转捩区域出现下凹现象，但随着湍流强度减少，该现象逐渐消失。

2 理论方法

SST transition模型是SST k-ω传输方程和其它两个传输方程耦合的结果，其中一个运输方程用于控制间歇，另外一个方程为控制转捩起始条件，用动量厚度雷诺数表示。

2.1 流体一般控制方程

质量守恒方程，即连续性方程

(1)

动量守恒方程，即Navier-Stokes方程

(2)

能量守恒方程

(3)

(4)

其中，Φ是粘性导致的能量耗散函数。

2.2 Transition SST模型的传输方程

传输方程的间歇系数γ定义如下

(5)

传输源项定义为：

PJ1=2FlengthρS[γFlength]CJ3

(6)

(7)

其中，S为应变速率大小，Flength是一个控制转捩区域长度经验相关的因子。消减项/再层流化项被定义为

PJ2=(2CJ1)ρΩγFturb

(8)

Er3=Cr3

(9)

其中，Ω是涡量大小。转捩的开始由以下的系数控制

(10)

(11)

(12)

(13)

(14)

Fonset=max(Fonset2-Fonset3，0)

(15)

(16)

CJ1=0.03；CJ1=50；CJ3=0.5；σJ=1.0

2.3 转捩修正

对分离引起的转捩的修正如下

(17)

(18)

γeff=max(γ，γsep)

(19)

(20)

源项定义为

(21)

(22)

(23)

(24)

(25)

(26)

(27)

(28)

Reθt=f(Tu，λ)

(29)

(30)

(31)

第一个经验相关式子为当地湍流强度方程，Thawaite压力梯度因子λθ定义为

λθ=(θ2/v)dU/dS

(32)

2.4 耦合转捩模型和SST传输方程

给出SST传输方程后，为了计算更为准确，需要将SST传输方程和转捩模型耦合

(33)

Pk=γeffPkP

(34)

(35)

(36)

(37)

Ft=max(Florig，F3)

(38)

其中，Dk和Pk为SST模型的原始生成项和原始销毁项，F1orig为原始SST混合方程。

2.5 指定入口湍流等级

在入口给定的湍流强度根据粘性比不断衰减。粘性比越大，则湍流衰减率越小。湍动能的衰减由以下公式计算

(39)

对于在自由流中SST湍流模型的常数为

β=0.09，β*=0.0828

时间尺度由下面定义

(40)

x为进口流向距离，V为平均对流速度。涡流粘度被定义

(41)

湍动能衰减方程可由入口湍流强度和涡流粘性比重新定义为

(42)

2.6 位移边界层的定义

位移边界层厚度δ1的计算公式为

(43)

其中，u为流体速度，u∞为无穷远处自由流速度。

3 零压力梯度平板边界层仿真

使用商业CFD软件FLUENT的不可压缩非结构隐式RANS求解器对流动进行仿真，基本原理是数值求解流体力学基本方程Navier-Stokes(N-S)方程和SST k-ω湍流模型的控制方程、γ-Reθt转捩模型的控制方程。通过对经典的平板边界层转捩实验T3A和T3B[13]仿真，以及对Schrader[20]的数值仿真进一步深入探究，简化其几何模型，设置不同的入口湍流强度进行计算，分析零压力梯度平板边界层流动的位移边界层变化规律。

3.1 T3A和T3B仿真

通过对董平[13]T3A和T3B的实验，该实验的示意图如图1所示，进口自由流的湍流强度、进口自由流的速度以及粘性比如表1所示，将其看作二维流动。几何模型大小为1700mm×200mm，网格划分为1700×600。为了更为准确的预测转捩过程以及得到边界层内的流动速度分布，故对转捩发生的区域和壁面网格都进行了加密，距壁面的第一个网格点y+≤1。选择的物理模型为SST transition四方程转捩模型，差分格式均为二阶迎风格式，所有物理参考量均为入口参数。入口参数设置与实验进口条件相同，底面为No-Slip条件，顶部为Symmetry条件，出口设置为Outflow条件。

图1 平板转捩实验示意图

表1 实验的进口参数

3.2 仿真结果

对于T3A和T3B实验的计算结果，首先分析其壁面摩擦系数(Cf)能够很好的表征转捩发生的位置，通常认为Cf最小处为转捩的起始点，最大处为转捩的终止点)，与实验所测的值进行比较如图2和图3所示。对于自由流湍流强度为3%的T3A的计算结果，计算得到壁面摩擦系数与实验测量得到的值几乎相同；对于湍流强度为6%的T3B的计算结果，计算得到的壁面摩擦系数与实验值也基本一致。故可以分析得出使用SST transition模型模拟自由流横掠平板的边界层实验较为准确，能够比较准确捕捉转捩点的起始位置和终止位置。

图2 T3A 的摩擦系数

图3 T3B的摩擦系数

图4记录了T3A和T3B位移边界层厚度沿着流向方向的变化。仔细观察和分析图4，发现位移边界层沿着流向方向变化有一个奇怪的现象：对于T3A的仿真结果，位移边界层厚度沿着流向方向的变化趋势是先陡峭增大而又有略微趋近平缓增长再陡峭的增大；而对于T3B的仿真结果，位移边界层厚度沿着流向方向的变化趋势是一直增大。而且对于T3A来说，从它的壁面摩擦趋势图显示，x=0.5m附近为转捩的起始位置(壁面摩擦系数最小)，x=0.8m附近(壁面摩擦系数最大)为转捩的终止位置，恰巧位移边界层趋于平缓增长的趋势的位置大致在0.5m到0.8m这个区间。为了更加深入了解自由流横掠平板的位移边界层的变化规律，进一步对Schrader[20]所做的自由流湍流横掠平板的转捩数值模拟实验(该数值实验中，出现了位移边界层沿着流向方向出现下凹趋势)进行展开研究。

图4 T3A和T3B位移边界层厚度

3.3 Schrader数值仿真

将其实验简化为二维几何模型，模型长度为1200mm，宽度为80mm，选用网格1700×600。为了准确预测转捩过程和捕捉边界层内的流动速度分布，故对壁面网格进行局部加密，距壁面的第一个网格点y+≤1。选择物理模型为SST transition四方程转捩模型，差分格式均为二阶迎风格式，所有物理参考量均为入口参数。底面为No-Slip条件，顶部为Symmetry条件，出口设置为Outflow条件。入口条件设置不同的入口自由流湍流强度，湍流尺度为10mm。

3.3.1 仿真结果

自由流湍流强度大于1%的时候发生bypass转捩，自由流湍流强度小于1%的时候发生自然转捩。分别对自由流湍流强度大于1%和湍流强度小于1%计算和讨论，分析计算结果和进行比较。

对于自由流的湍流强度Tu=5.06%、Tu=3.73%、Tu=2.95%、Tu=2.53%、Tu=2.11%，为了能够较为清晰地捕捉转捩现象的发生和保证仿真的转捩发生在平板上，故选择自由流入口速度为V=5m/s。

壁面摩擦系数沿着流向方向的计算结果如图5所示，发现自由流的湍流强度越大，转捩位置越靠前。计算其位移边界层计算结果如图6所示，发现转捩位置越靠前，位移边界层有下凹趋势或者趋于平缓的位置越靠前，且该位置出现在转捩区域附近。

图5 不同湍流强度下的壁面摩擦系数

对于不同的自由流湍流强度，如果入口的自由流速度太小，则在有限长度的平板上可能不发生转捩；如果入口的自由流速度太大，则在平板 前缘部分可能很早就发生转捩，且不容易观测到转捩。对于自由流湍流强度Tu=2.11%、Tu=1.69%、Tu=1.26%设置来流入口速度V=10m/s，计算得到的壁面摩擦系数沿流向方向的变化如图5所示，位移边界层厚度沿流向方向的变化如图6所示。

图6 不同湍流强度下的位移边界层厚度

观察和分析自由流的湍流强度为Tu=2.11%、Tu=1.69%、Tu=1.26%的计算结果，发现自由流的湍流强度越大，转捩位置越靠前。随着自由流的湍流强度减少，相同入口速度下的位移边界层沿着流向方向的下凹现象越来越明显。

图7 不同湍流强度下的壁面摩擦系数

图8 不同湍流强度下的位移边界层厚度

综合自由流的湍流强度大于1%的计算结果，相对比较位移边界层厚度沿着流向方向的变化趋势发现，随着入口自由流的湍流强度减少，位移边界层沿着流向方向的变化趋势从沿着流向方向一直陡峭地增大过渡到先陡峭地增大而后平缓地增大再到相对陡峭地增大，最后过渡到出现先陡峭地增大而后减少再增大的趋势(位移边界层厚度沿着流向方向的变化曲线出现下凹的现象)。

而对于自由流的湍流强度低于1%的设置，需要进行略微调整。因为在大量数值计算后发现，对于给定低于1%的自由流的湍流强度，横掠平板的位移边界层厚度对自由流速度非常敏感，难以捕捉位移边界层趋于平缓增长或出现下降的变化趋势。故经过多组速度的仿真计算，得到位移边界层变化趋势如图9所示。

图9 低湍流强度下的位移边界层厚度

通过对比数据，发现在湍流强度小于1%的时候，位移边界层变化规律与自由流湍流强度大于1%的位移边界层变化规律是不同的，平板边界层位移厚度沿流向的变化曲线随着湍流强度变大而容易出现下凹或者趋于平缓的现象。

3.4 位移边界层现象分析

从图7和图8中任意选择一个湍流强度对应的壁面摩擦系数和位移边界层厚度的曲线考虑，壁面摩擦系数在层流区域是随流动发展减小的，在湍流区域也是随着流动发展减小的，但是在转捩区域却是增大的(这里取湍流强度为1.26%作为案例说明，大致在X=1.0m处开始转捩，X=1.2m处完成转捩)。位移边界层厚度δ1产生原因是因为壁面无滑移所导致的速度梯度，进而导致流体流过该区域的质量流量损失(图10中b的阴影区域)。观察不同X位置处对应的速度轮廓曲线，如下图11所示：从速度轮廓曲线中知道质量损失分为两块(质量损失为水平速度u在Y方向的积分，在图上表示为面积如图10(b)的阴影部分)，大约在Y=0.002以下和Y=0.002以上的两部分(以水平线Y=0.002分界)。随流动发展，下半部分的面积逐渐减少，上半部分的面积逐渐增大。在过渡区域，一度出现减少大于增大的情况(X=1.0m～1.2m)，上半部分增大的面积小于下半部分增大的面积，此时相当于总质量流量是减少的，故位移边界层厚度δ1减少。在湍流区，上半部面积显著增大，总的阴影面积也增大，此时位移边界层厚度δ1增大。这就解释了为什么边界层厚度δ1在转捩的过渡区域出现先减少后增大的现象。

图10 边界层内速度分布

图11 基于速度轮廓的原因分析

4 结论

通过对零压力梯度的平板边界层实验T3A、T3B[13]，以及Schrader[20]的仿真，模拟了对于不同自由流湍流度的自由流横掠平板层流到湍流的转捩过程，分析其计算结果和对比与之相关的实验得出以下的结论：

1)在来流自由流的湍流强度大于1%的时候，随着湍流强度的减少，平板的位移边界层厚度沿着流向方向先是增大而后稍微有下降而后又增大。并且该现象随着湍流度的降低，越来越明显。

2)在来流自由流的湍流强度小于1%的时候，随着湍流强度的减少，平板的位移边界层变化的趋势很难捕捉到，基本上是沿着流向方向一直是增大的趋势。在湍流强度接近1%的时候，位移边界层厚度沿着流向方向的变化会出现下降趋势。

3)湍流强度1%是一个分界点，湍流度大于1%，会发生bypass转捩，即边界层对自由流脉动的感应而导致的bypass转捩；而湍流度小于1%发生的是自然转捩，即边界层失稳而导致的层流瓦解。由于转捩的方式不同，而导致的位移边界层随湍流度的变化，而导致趋势变化不同。

4)位移边界厚度沿着流向方向出现下降或者趋于平缓的变化与转捩相关。