引言

在航空航天领域, 飞行器边界层转捩会增大壁面摩擦阻力, 提高热交换, 影响流动分离位置及发动机燃料混合等, 准确预测边界层转捩能够有效改进飞行器的气动性能, 合理设计热防护系统, 降低燃料消耗, 因此边界层转捩位置的准确预测已经成为一个非常重要的工程问题.在翼型设计领域, 使用全湍流或者全层流计算将会带来较大的气动力误差, 因此, 急须建立针对飞行器表面转捩位置的预测方法, 从而提高飞行器气动力的预测精度.在工程实际中, 边界层流动通常是三维的, 典型的如后掠机翼, 涡轮叶片, 带有攻角的圆锥、椭锥等.诱导三维边界层失稳的模态主要有附着线不稳定性、Taylor-Görtler不稳定性、流向Tollmien-Schlichting波不稳定性和横流不稳定性.其中, Tollmien-Schlichting波不稳定性和横流不稳定性是诱导三维边界层转捩的主要机制^[1].工程实践中常用的转捩预测方法有e^N法^[2]、经验关系法^[3]、转捩准则^[4]和转捩模式等, 其中应用最为广泛的是e^N法和转捩模式. e^N法基于线性稳定性理论^[5], 通过求解O-S方程来预测低湍流度下沿流向的自然转捩和分离流转捩, 但在求解过程中, 需要沿壁面法向的1阶和2阶速度导数, 这些信息在三维流场中, 因为求解路径的不确定而难以得到.基于Reynolds平均方程(RANS)的转捩模式方法成为工程领域转捩预测的主要方法.在转捩模式的提出与构造方面, 目前文献公开发表的方法多达十几种, 如Langtry等^[6]的γ-Re_θt模式, Walters等^[7]的k_L-k-ω模式, Warren等^[8]的k-ζ模式以及Wang等^[9]的k-ω-γ模式等.其中, γ-Re_θt模式集合了转捩经验关系式和低Reynolds数湍流模型的优势, 并被耦合到SST k-ω湍流模型中, 转捩的判别基于当地变量, 因而大大提高了通用CFD求解复杂外形的能力, 特别是对于三维复杂情况下的流动, 在工程中应用更加广泛.但是, 该模型在构造中针对的是T-S波不稳定性诱发的转捩, 对于横流诱导的转捩预测效果较差.近几年来, 人们大量通过当地变量来构造一些参数, 以解决横流预测问题. Seyfert等^[10]将C1准则耦合到转捩模式中; Grabe等^[11]引入形状因子概念希望将横流的求解当地化, 但该类方法仅适用于翼型; Medida等^[12]给出了横流转捩起始位置函数F_onset, 但由于涉及边界层外缘信息, 给出的流场转捩参数分布并不光滑.此外, Müller等^[13]、Choi等^[14]于2014年分别给出了将γ-Re_θt转捩模式拓展到三维横流转捩预测的方法, 均取得了一定的效果.但这些判据都不是完全基于当地变量信息, 构造的横流转捩预测方法均受几何外形约束, 应用范围存在局限.考虑到γ-Re_θt转捩模式在自然转捩、旁路转捩和分离诱导转捩方面的优良表现, 本文对其进行基本构造分析, 引入横流间歇因子概念, 将横流转捩预测添加到原有的γ-Re_θt转捩模式中, 对后掠机翼的三维边界层转捩进行分析, 验证改进的模式对于三维边界层转捩的预测精度和可靠性, 同时构建出计算固定前缘后掠角机翼横流不稳定转捩的判定方法.

1 横流转捩机理

诱导三维边界层转捩的失稳机制主要有第一、第二模态流向不稳定, Görtler离心不稳定, 附着线不稳定以及横流不稳定等, 其中横流不稳定是诱导三维边界转捩的重要因素^[14].以后掠翼边界层为例, 如图 1(a)所示^[14], 图中Λ为机翼后掠角; u和w分别为流向速度和横流速度, U_e为边界层外缘速度.由于机翼弦向压力梯度的存在, 边界层外缘无黏流线发生弯曲, 弯曲的流线产生离心力, 与径向的压力平衡.而边界层内, 由于压力保持为常数, 流向速度减小, 对应的离心力和压力不再平衡, 多出的压力使得沿垂直于边界层外缘无黏流线方向出现横向速度分量, 即横流.由于壁面和边界层外缘横流速度分量均为零, 因此边界层内横流速度剖面通常存在拐点, 流动易发生横流失稳.横流失稳属于无黏失稳, 相比于黏性失稳(T-S波), 其产生的扰动增长率要大很多, 且在很低Reynolds数条件下也易发生^[15].

Malik等^[17]的理论研究以及White等^[18]的实验研究表明, 横流扰动存在两种形式, 频率为零的定常驻涡模态(stationary crossflow mode)和相速度不为零的非定常行波模态(travelling crossflow mode). 图 1(b)^[15]给出了两类横流模态的示意图, 其中驻涡模态为边界层内一系列沿展向排列、同向旋转的定常涡, 方向近似平行来流, 在流场显示技术中表现为定常条纹结构; 而行波模态则表现为非定常的、以大的倾斜角沿流向下游的横向运动.低湍流度环境下, 横流不稳定主要表现为驻涡模态, 而在中高湍流度环境下横流不稳定则由行波模态主导.虽然稳定性理论预测表明, 最不稳定行波模态比最不稳定驻涡模态的增长率更大, 但由于实际飞行环境来流湍流度很低, 机翼表面通常出现驻涡模态, 可通过在机翼前缘布置一排微小粗糙单元进行控制.

2 数值模拟方法 2.1 原始γ-Re_θt转捩模式^[6]

为避免使用非当地变量, γ-Re_θt转捩模式引用应变率Reynolds数Re_ν代替动量厚度Reynolds数Re_θ, 从而隐含了边界层厚度的信息:

$ \begin{array}{l} \mathit{R}{\mathit{e}_\mathit{\nu }}{{\rm{(}}\mathit{x}{\rm{, }}\mathit{\tilde y}{\rm{)}}_{{\rm{max}}}}{\rm{ = 2}}{\rm{.193}}\mathit{R}{\mathit{e}_\mathit{\theta }}\left( \mathit{x} \right) \\ \;\;\;\;\;\mathit{R}{\mathit{e}_\mathit{\nu }}{\rm{ = }}\frac{{\mathit{\rho }{\mathit{y}^{\rm{2}}}}}{\mathit{\mu }}\frac{{\partial \mathit{u}}}{{\partial \mathit{y}}}{\rm{ = }}\frac{{\mathit{\rho }{\mathit{y}^{\rm{2}}}}}{\mathit{\mu }}\mathit{S} \end{array} $

式中, $ \mathit{\tilde y}$为Re_ν最大值的壁面法向位置, μ为分子黏性系数, S为应变率的模. y²S表征了边界层内扰动的增长, 而ν=μ/ρ反应了扰动的衰减, 因此Re_ν的使用具有一定的物理意义.当y²S随着边界层增厚逐渐增大而ν保持不变时, 则Re_ν达到某一临界值, 转捩发生.在使用应变率Reynolds数Re_ν的基础上, Menter和Langtry构造了间歇因子γ和转捩起始位置处动量厚度Reynolds数${\widetilde {\mathit{Re}}_{\mathit{\theta }{\rm{t}}}}$两个标准输运方程, 将Re_ν与实验得到的转捩拟合公式相耦合, 两个输运方程的具体形式为:

$ \begin{array}{l} \frac{{\partial \left( {\mathit{\rho \gamma }} \right)}}{{\partial \mathit{t}}} + \frac{{\partial \left( {\mathit{\rho }{\mathit{u}_\mathit{j}}\mathit{\gamma }} \right)}}{{\partial {\mathit{x}_\mathit{j}}}} = \frac{\partial }{{\partial {\mathit{x}_\mathit{j}}}}\left[{\left( {\mathit{\mu + }\frac{{{\mathit{\mu }_{\rm{t}}}}}{{{\mathit{\sigma }_{\rm{f}}}}}} \right)\frac{{\partial \mathit{\gamma }}}{{\partial {\mathit{x}_\mathit{j}}}}} \right] + {\mathit{P}_\mathit{\gamma }} - {\mathit{D}_\mathit{\gamma }}\\ \frac{{\partial \left( {\mathit{\rho }{{\widetilde {\mathit{Re}}}_{\mathit{\theta }{\rm{t}}}}} \right)}}{{\partial \mathit{t}}}{\rm{ + }}\frac{{\partial {\rm{(}}\mathit{\rho }{\mathit{u}_\mathit{j}}{{\widetilde {\mathit{Re}}}_{\mathit{\theta }{\rm{t}}}}{\rm{)}}}}{{\partial {\mathit{x}_\mathit{j}}}}{\rm{ = }}\frac{\partial }{{\partial {\mathit{x}_\mathit{j}}}}\left[{{\mathit{\sigma }_{\mathit{\theta }{\rm{t}}}}\left( {\mathit{\mu }{\rm{ + }}{\mathit{\mu }_{\rm{t}}}} \right)\frac{{\partial {{\widetilde {\mathit{Re}}}_{\mathit{\theta }{\rm{t}}}}}}{{\partial {\mathit{x}_\mathit{j}}}}} \right]{\rm{ + }}{\mathit{P}_{\mathit{\theta }{\rm{t}}}} \end{array} $

γ输运方程中转捩源项定义为:

$ \begin{array}{l} {\mathit{P}_\mathit{\gamma }}{\rm{ = }}{\mathit{F}_{{\rm{length}}}}{\mathit{c}_{\mathit{a}{\rm{1}}}}\mathit{\rho S}\sqrt {\mathit{\gamma }{\mathit{F}_{{\rm{onset}}}}} {\rm{(1 - }}{\mathit{c}_{\mathit{\varepsilon }{\rm{1}}}}\mathit{\gamma }{\rm{)}}\\ \;\;\;{\mathit{D}_\mathit{\gamma }}{\rm{ = }}{\mathit{c}_{\mathit{a}{\rm{2}}}}\mathit{\rho \Omega \gamma }{\mathit{F}_{{\rm{turb}}}}{\rm{(}}{\mathit{c}_{\mathit{\varepsilon }{\rm{2}}}}\mathit{\gamma - }{\rm{1)}} \end{array} $

式中, 函数F_length由经验拟合公式得到, 用于控制转捩区长度.函数F_onset与Re_ν相关, 用于触发转捩起始位置, 两者均为无量纲函数, 用于控制边界层内的间歇因子发展.

$ {\widetilde {\mathit{Re}}_{\mathit{\theta }{\rm{t}}}}$中的源项用于将输运标量$ {\widetilde {\mathit{Re}}_{\mathit{\theta }{\rm{t}}}}$与当地值Re_θt相对比, 定义为

$ {\mathit{P}_{\mathit{\theta }{\rm{t}}}}{\rm{ = }}{\mathit{c}_{\mathit{\theta }{\rm{t}}}}\frac{\mathit{\rho }}{\mathit{t}}{\rm{(}}\mathit{R}{\mathit{e}_{\mathit{\theta }{\rm{t}}}}{\rm{ - }}{\widetilde {\mathit{Re}}_{\mathit{\theta }{\rm{t}}}}{\rm{)(1}}{\rm{.0 - }}{\mathit{F}_{\mathit{\theta }{\rm{t}}}}{\rm{), }}\;\;\;\mathit{t}{\rm{ = }}\frac{{{\rm{500}}\mathit{\mu }}}{{\mathit{\rho }{\mathit{U}^{\rm{2}}}}} $

式中, Re_θt由经验拟合关系式得出.混合函数F_θt在边界层内为1, 在自由来流中为0, 因此, 该函数在边界层内将源项关闭, 而在自由流中允许输运标量$ {\widetilde {\mathit{Re}}_{\mathit{\theta }{\rm{t}}}}$通过扩散进入边界层中.这样, 经验公式和来流条件对于转捩的影响就与当地流场有机结合起来.

通过以上公式求得有效间歇因子γ_eff, 应用到原始的SST湍流模型方程中，有效间歇因子γ_eff通过影响湍动能k输运方程的生成项P_k和耗散项D_k, 进而模拟转捩效应的增长, 预测边界层转捩.湍动能k输运方程的具体形式如下:

$ \begin{array}{l} \frac{{\partial \left( {\mathit{\rho k}} \right)}}{{\partial \mathit{t}}}{\rm{ + }}\frac{{\partial {\rm{(}}\mathit{\rho }{\mathit{u}_\mathit{j}}\mathit{k}{\rm{)}}}}{{\partial {\mathit{x}_\mathit{j}}}}{\rm{ = }}\frac{\partial }{{\partial {\mathit{x}_\mathit{j}}}}\left[{{\rm{(}}\mathit{\mu }{\rm{ + }}{\mathit{\sigma }_\mathit{k}}{\mathit{\mu }_{\rm{t}}}{\rm{)}}\frac{{\partial \mathit{k}}}{{\partial {\mathit{x}_\mathit{j}}}}} \right]{\rm{ + }}\mathit{\tilde P - }{{\mathit{\tilde D}}_\mathit{k}}\\ \;\;\;\;\;\;{{\mathit{\tilde P}}_\mathit{k}}{\rm{ = }}{\mathit{\gamma }_{{\rm{eff}}}}{\mathit{P}_\mathit{k}}{\rm{, }}{{\mathit{\tilde D}}_\mathit{k}}{\rm{ = min}}\left[{{\rm{max}}\left( {{\mathit{\gamma }_{{\rm{eff}}}}{\rm{, 0}}{\rm{.1}}} \right)} \right]{\mathit{D}_\mathit{k}}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\mathit{\gamma }_{{\rm{eff}}}}{\rm{ = max}}\left( {\mathit{\gamma }{\rm{, }}{\mathit{\gamma }_{{\rm{sep}}}}} \right) \end{array} $

有效间歇因子的使用主要是为了提高分离诱导边界层转捩预测的准确性.当流动出现边界层分离时, 模型允许当地间歇因子大于1(γ_eff=γ_sep>1), 从而引起湍动能迅速增长, 而湍动能输运方程中大的生成项使得边界层再附能被准确预测.

2.2 横流间歇因子构造^[19]

为了能够对三维边界层转捩进行准确模拟, 本文在原始的γ-Re_θt转捩模式上增加了横流间歇因子项, 用于预测横流不稳定性对转捩的诱导效果.该横流间歇因子具体构造如下^[19]:

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\mathit{\gamma }_{{\rm{eff}}}}{\rm{ = max(}}\mathit{\gamma }{\rm{, }}{\mathit{\gamma }_{{\rm{sep}}}}{\rm{, }}{\mathit{\gamma }_{{\rm{cf}}}}{\rm{)}}\\ {\mathit{\gamma }_{{\rm{cf}}}}{\rm{ = 1 - exp( - }}{\mathit{a}_{\rm{1}}}{\rm{\cdot(max(}}\mathit{R}{\mathit{e}_{{\rm{cf}}}}{\rm{/}}\mathit{R}{\mathit{e}_{{\rm{cf\_crit}}}}{\rm{ - 1}}{\rm{.0, 0}}{\rm{.0)}}{{\rm{)}}^{\rm{2}}}{\rm{\cdot}}{\mathit{F}_{{\rm{cf}}}}{\rm{)}} \end{array} $

式中, a₁为控制转捩区长度的参数, 此处, 取a₁=3.当地横流Reynolds数Re_cf和横流控制函数F_cf构造如下:

$ \begin{array}{l} \;\;\;\;\;\mathit{R}{\mathit{e}_{{\rm{cf}}}}{\rm{ = }}{\mathit{\rho }_{\rm{e}}}\mathit{\delta }{\mathit{U}_{{\rm{cf\_max}}}}{\rm{/}}{\mathit{\mu }_{\rm{e}}}\\ {\mathit{F}_{{\rm{cf}}}}{\rm{ = }}\left\{ \begin{array}{l} {\rm{1, }}\;\;{\mathit{U}_{{\rm{cf\_max}}}} \ge {\mathit{U}_{{\rm{cf\_crit}}}}\\ {\rm{0, }}\;\;{\mathit{U}_{{\rm{cf\_max}}}}{\rm{ < }}{\mathit{U}_{{\rm{cf\_crit}}}} \end{array} \right. \end{array} $

式中, δ, ρ_e, μ_e分别为边界层厚度、边界层外缘密度和边界层外缘分子黏性. U_{cf_max}为壁面法向的最大横流速度, 其大小可达到自由来流速度的4%~6%. Re_{cf_crit}为横流速度Reynolds数, 它决定了横流失稳诱导转捩的起始位置.将横流间歇因子应用于原始的γ-Re_θt间歇因子中, 就得到了改进的可以预测横流转捩的γ-Re_θt的转捩模式.

3 后掠翼型算例验证

为验证加入横流间歇因子的γ-Re_θt转捩模式的正确性以及考察模式对横流不稳定性诱导的边界层转捩的预测精度.选取NLF(2)-0415后掠翼型作为算例, 并将数值模拟结果与实验值进行比较.

3.1 NLF(2)-0415后掠翼型

该实验在亚利桑那州立大学的低湍流风洞完成, 该风洞来流湍流度(Tu_∞)可以达到0.8%, 风洞实验段为1.4 m×1.4 m×4.9 m. NLF(2)-0415翼型的无限翼展后掠机翼模型^[20]如图 2所示, 模型的后掠角为45°, 与来流呈-4°攻角.具体实验参数如表 1所示:基于弦长的自由来流Reynolds数Re=1.9×10⁶~6.5×10⁶, 自由来流湍流度Tu_∞=0.1%, 实验中采用萘(naphthalene)升华方法测定模型上表面的转捩位置. x_{tr, exp}为实验测定的相对转捩位置(转捩位置/机翼弦长).

计算网格如图 2所示, 为保证物面网格y⁺＜1, 第1层网格高度为2×10^-6 m, 法向网格增长率为1.1, 网格量为3×10⁶, 远场为弦长的20倍, 展长为弦长的2.5倍, 在展向两端设置周期边界条件, 以减少展向的长度与网格点数.

3.2 计算结果及分析

图 3对比了后掠机翼NLF(2)-0415的上翼面压强系数分布, 横坐标x/c表示弦长位置, 纵坐标-C_{p, 3}表示后掠翼型的压强分布. C_{p, 3}=C_p·cos²Λ^[20].在-4°攻角情况下, 流向方向弦长x/c=0.2~0.71的位置处, 出现了一大段较长的线性逆压梯度区间.在此区间内, 由稳定性理论可知, T-S波不稳定性被有效抑制, 而横流不稳定性被充分放大.由于NLF(2)-0415翼型头部钝度较小, 且上翼面没有壁面凹曲, 因此, 并未发生附着线不稳定和Görtler不稳定.所以, 横流不稳定性成为诱导上翼面流动失稳的最主要因素.

将上翼面的压强系数数值模拟结果和理论预测、实验值进行对比可发现, 与实验值相比, 数值结果与理论预测更加接近.同时, 不同Reynolds数下, 上翼面的压强系数基本不变, 该逆压条件下均可发展成有效的横流.

图 4为原始的γ-Re_θt转捩模式与加入横流项改进的γ-Re_θt转捩模式数值计算得到的转捩点位置随Reynolds数变化曲线及其与实验值的对比情况.从图中可以看出, 原始的γ-Re_θt转捩模式在不同的Reynolds数条件下, 预测的转捩位置变化不大, 均在弦长的80%左右, Reynolds数的变化并不会影响转捩点的位置.再看添加横流修正的转捩模式预测得到的转捩点位置, 与实验值吻合较好.分析其原因, 随着Reynolds数增大, 主导转捩的因素由T-S波不稳定性转化为横流不稳定性, 原始的γ-Re_θt转捩模式只能预测流向T-S波不稳定引起的转捩, 无法预测横流不稳定性诱发的转捩, 因而, 预测结果与实验值相差很大.同时, 可以看出改进的γ-Re_θt转捩模式能够有效地模拟横流效应, 预测出的转捩位置与实验值^[16]吻合较好.

为了进一步探究横流不稳定性对转捩位置的影响, 图 5给出了壁面流线和空间流线分布.其中, 灰色线代表壁面流线, 粉色线代表空间流线.从图 5可以看到, -4°攻角条件下, 上翼面表面流线沿流向方向呈现明显的“S”型偏转, 壁面附近, 流线弯曲较大, 往垂直壁面方向发展, 流线曲率不断降低, 由于机翼弦向压力梯度的存在, 边界层外缘无黏流线发生弯曲, 弯曲的流线产生离心力, 与径向压力平衡.而边界层内, 由于压力保持为常数, 但是速度减小, 导致离心力和压力不再平衡, 多出的压力使得沿垂直于边界层外缘无黏流线方向出现横向速度分量, 横流诱导边界层失稳, 转捩发生.

为了进一步分析NLF(2)-0415后掠机翼边界层内横流速度的分布情况, 截取展向中间位置(z=0.5 m), 流向x/c=0.2~0.7(线性逆压梯度部分)范围内8个位置点处的无量纲横流速度U_cf/U_e进行分析, 横轴表示无量纲横流速度, 纵轴d表示距离壁面的垂直距离.由图 6可以看到, 从x/c=0.2处到流向方向x/c=0.7位置处, 边界层逐渐变厚, U_cf/U_e先增大后减小.在x/c=0.65(图中黑色虚线图标)站位处, 横流速度达到最大值, 最大无量纲横流速度值超过边界层外缘速度的7%.根据Malik等^[17]的研究, 当横流速度达到6%~8%的强度时, 足以诱导边界层发生转捩.

图 7显示了在来流Reynolds数Re=3.73×10⁶下, 原始γ-Re_θt模式和加入横流后改进的γ-Re_θt模式后掠翼型上表面摩擦力系数分布云图的俯视图.壁面摩擦阻力主要由黏性剪切效应产生, 当流动由层流转捩为湍流后, 边界层变厚, 同时边界层内外的物质交换和能量交换更加显著, 壁面摩擦力也就越大.从图 7可以看到, 改进的转捩模式得到的转捩区域明显比原始的转捩模式预测得到的位置更加靠前, 同时, 整体的摩擦力系数也更大.这是因为, 对于该翼型, 在45°后掠角下, 横流不稳定性已经成为诱导转捩发生的主要原因, 而原始γ-Re_θt模式未考虑横流效应对转捩的影响, 所以, 添加横流间歇因子的γ-Re_θt转捩模式能够更加合理地预测诱导转捩的失稳因素, 也就能够得到更加合理的转捩位置和转捩区域分布.同时, 在横流效应较强的情况下须考虑横流对转捩位置的影响.

4 结论

边界层信息的准确获取对边界层转捩工程预测方法的发展意义重大.在所获得的边界层外缘及相关信息基础上, 本文引入了结合ρ_e, δ, U_{cf_max}和μ_e等参数构造的横流间歇因子γ_cf项, 将其加入到原始γ-Re_θt转捩模式的有效间歇因子输运方法中, 得到了能够预测横流效应的γ-Re_θt转捩模式.通过对NLF(2)-0415后掠翼算例的检验和验证, 改进的γ-Re_θt转捩模式能够模拟横流不稳定性, 对三维边界层转捩的预测精度大大提高.三维边界层流动中的横向压力梯度使得边界层内出现横流速度分量, 横流速度剖面通常存在拐点, 流动易发生横流失稳.对于NLF(2)-0415后掠机翼, 其上翼面压力在x/c=0.2~0.7范围内呈线性下降, 在逆压梯度的作用下, 上翼面出现大面积的横流作用区域.不同的Reynolds数下, 未考虑横流不稳定性的γ-Re_θt转捩模式预测得到的边界层转捩区域与实验结果差异较大, 而考虑横流不稳定性的γ-Re_θt转捩模式与实验结果更加吻合, 表明本文引入的方法对三维边界层横流转捩有较高的预测精度.