0 引言

湍流模型是目前工程中湍流计算不可缺少的工具。工程上常用的湍流模型很多, BL模型是一种简单实用的零方程模型。它把湍流剖面分成内层与外层, 分别用混合长概念与尾流函数模化涡粘系数。这种模型没有普适性, 却可以较准确地模拟小曲率的湍流边界层。但对于超声速及高超声速湍流边界层, 这种模型能给出的精度却由于缺乏可靠的实验数据而无法确定。本文将利用直接数值模拟的结果来检验这种模型的精度。

实际的边界层流, 往往既有层流段, 又有湍流段。从层流到湍流的转捩位置的预测精度显然会影响到总的计算精度。本文虽只考虑湍流段的计算, 但转捩位置预测精度对湍流段计算精度的影响, 也是本文讨论的问题之一。

本文共选择了三个超声速湍流边界层的典型算例, 对比了BL模型与DNS的结果。分析了BL模型在计算壁面摩擦力、热传导系数、边界层厚度的增长及平均流剖面形状等方面的准确性, 以及转捩位置预测精度对计算结果的影响。另外, BL模型在提出时有对转捩位置预测的方法, 本文也将对这种预测方法做一次检验。

1 数值方法 1.1 BL模型的回顾

对N-S方程进行雷诺平均, 将出现二阶矩项, 使得方程不封闭。1872年, Boussinesq仿照分子粘性应力与速度变形率的关系, 引入了涡粘系数及湍流导热系数的概念, 即对于可压缩流的Favre平均, 有:

$ \begin{array}{*{20}{c}} {\tau _{ij}^{\rm T} = - \overline {\rho u_i^{''}u_j^{''}} = - \bar \rho \widetilde {u_i^{''}u_j^{''}} = {\mu _{\rm T}}(2{{\tilde S}_{ij}} - \frac{2}{3}\frac{{\partial {{\tilde u}_k}}}{{\partial {x_k}}}{\delta _{ij}})}\\ {q_j^{\rm T} = \overline {\rho u_j^{''}{h^{''}}} = \bar \rho \widetilde {u_j^{''}{h^{''}}} = - {k_{\rm T}}\frac{{\partial \tilde T}}{{\partial {x_j}}}\quad \quad \quad }\\ {{k_T} = \frac{{{c_p}{\mu _T}}}{{P{r_T}}}} \end{array} $

(1)

其中μ_T, k_T, Pr_T, S_ij分别为涡粘性系数, 湍流导热系数, 湍流Prandtl数与速度变形张量。

BL模型中, 涡粘系数μ_T的取值由下式来确定:

$ {\mu _{\rm T}} = \left\{ {\begin{array}{*{20}{l}} {{{({\mu _{\rm T}})}_{{\rm{inner}}}}\quad \quad y \le {y_{{\rm{crossover}}}}}\\ {{{({\mu _{\rm T}})}_{{\rm{outer}}}}\quad \quad y > {y_{{\rm{crossover}}}}} \end{array}} \right. $

(2)

其中, y为法向坐标, y_crossover为满足(μ_T) _inner= (μ_T) _outer的最小的法向距离。

内层用Prandtl-Van Driest建议的公式:

${({\mu _T})_{{\rm{inner}}}} = \rho l_{{\rm{mi}}x}^2|\omega |$

(3)

其中, l_mix=κy[1-exp (-y⁺/A⁺)], 表示混合长度; $|\omega | = \sqrt {{{(\frac{{\partial u}}{{\partial y}} - \frac{{\partial v}}{{\partial x}})}^2} + {{(\frac{{\partial v}}{{\partial z}} - \frac{{\partial w}}{{\partial y}})}^2} + {{(\frac{{\partial w}}{{\partial x}} - \frac{{\partial u}}{{\partial z}})}^2}} $, 表示涡量, 对我们的二维情况, 它可写为$|\omega | = |\frac{{\partial u}}{{\partial y}} - \frac{{\partial v}}{{\partial x}}|;{y^ + } = \frac{{{\rho _w}{u_\tau }y}}{{{\mu _w}}} = \frac{{y\sqrt {{\rho _w}{\tau _w}} }}{{{\mu _w}}}$, 为粘性长度。

外层公式为:

${({\mu _T})_{{\rm{outer}}}} = {\rm K}{C_{CP}}\rho {F_{{\rm{WA{\rm K}E}}}}{F_{{\rm{{\rm K}LEB}}}}(y)$

(4)

其中, K为Clauser常数, C_CP为附加常数。F_WAKE=min (y_MAX F_MAX, C_WK y_MAX u_DΙE²/F_MAX)为尾流函数, 其中, y_MAX与F_MAX由函数F (y) =y|ω|[1-exp (-y⁺/A⁺]决定, F_MAX为F (y)出现的最大值, y_MAX为与此最大值对应的y值, ${u_{{\rm{D{\rm I}E}}}} = {(\sqrt {{u^2} + {v^2} + {w^2}} )_{{\rm{{\rm M}AX}}}} - {(\sqrt {{u^2} + {v^2} + {w^2}} )_{{\rm{{\rm M}{\rm I}{\rm N}}}}}$。间歇因子${K_{{\rm{{\rm K}LEB}}}}(y) = {[1 + 5.5{(\frac{{{C_{{\rm{{\rm K}LEB}}}}~y}}{{{y_{{\rm{{\rm M}AX}}}}}})^6}]^{ - 1}}$。

各常数取值为:A⁺=26, C_CP=1.6, C_KLEB=0.3, C_WK=0.25, κ=0.4, K=0.0168, Pr_T=0.9。

1.2 计算方法

工程上计算湍流时一般要从层流开始, 此时用N-S方程计算, 涡粘系数设为0。而转捩点的判据为:

${({\mu _T})_{{\rm{{\rm M}AX}}}} \ge {C_{{\rm{{\rm M}U{\rm T}{\rm M}}}}}{\mu _\infty }$

(5)

其中, C_MUTM为经验常数, 文献[1]中建议取为14, μ_T由式(2)计算得到。随着计算向下游进行, 当某个剖面满足式(5)时, 认为该处转捩结束, 后面的剖面都用湍流模型来计算。

本文是与DNS数据做对比, 转捩位置已经知道, 故转捩判据不采用公式(5), 而从DNS数据中直接得到。

2 算例对比 2.1 超声速平板湍流边界层

对比数据来自文献[2], 为时间模式下的超声速平板湍流边界层的直接数值模拟结果。来流马赫数为4.5, 以初始流场的排移厚度δ^*为特征长度的雷诺数为17000, 来流温度为226.5K。流场中各物理量均以无穷远参数无量纲化, 长度单位以δ^*无量纲化。DNS计算网格为129×257×129, 计算域为31.4×40×20.9。壁面采用无滑移、绝热条件, 上边界用嵌边区, 流向与展向均为周期条件。

用湍流模型计算的条件与DNS相同。由于是时间模式, 原三维问题可化为平均场的一维问题进行求解, 但多出一个时间自变量, 相对于空间模式的流向坐标。考察DNS数据, 计算域的法向宽度取到20即可; 并且没有扰动传播的问题, 故上边界用固定边界条件即可。即流场的初始阶段为Blasius解, 随着时间的演化, 当t=t₀时(t₀为从DNS数据中得到的转捩时刻), 引入湍流模型进行求解。

首先作为网格的验证, 法向分别选取41、81及161个点, 比较计算了壁面摩擦系数、边界层厚度及形状因子的流向分布。发现除了法向41个点情况下对边界层厚度的计算, 其它结果基本吻合, 而引起这个差异的主要原因是数值积分误差。以下的计算中法向均取81个网格。

图 1 (a)为C_f曲线的对比结果。BL模型计算中转捩点t₀的选择是一个问题。图中BL1与BL2两条曲线分别是把t₀设为100与160处(DNS数据的转捩起止时刻)的模型计算结果。实际上, 这两条线几乎可以通过平移而使之重合。而且, 两条曲线在湍流区相差很小, 同时, 在湍流区域, 两条曲线均与DNS数据吻合较好。说明把t₀选在转捩区的哪个点对湍流区的C_f曲线的计算影响不大。图 1 (b)是边界层排移厚度(δ^*)与形状因子(H)随时间的演化曲线与DNS结果的对比, t₀选取为160, 对应于转捩结束时刻, 此时边界层排移厚度吻合最好。可以看出, BL模型结果中的δ^*和H与DNS结果都比较接近, 但是δ^*的增长率并不完全相同。从趋势来看, 随着时间的推进, 边界层排移厚度的差别会越来越大。选择不同的t₀不会改变曲线的增长率, 只会把曲线平移。以下的统计结果都是基于t₀取160的情况。

图 2 (a)给出的是t=854时刻壁面律的分布, 其中u_c⁺为经过Van Driest变换$({u_c} = \mathop \smallint _0^{\bar u} {(\bar \rho /{\bar \rho _w})^{1/2}}{\rm{d}}u)$并以摩擦速度${u_\tau } = \sqrt {{\tau _w}/{{\bar \rho }_w}} $无量纲化的流向速度。由于BL模型内层涡粘系数的建模依据是混合长度, 它必然会有很好的对数律区域, 这从图中可以看出。对比DNS, 发现在对数律区域DNS比模型结果高, 这是直接数值模拟时所取的网格不够密所致, 而这个缺点BL模型是没有的。另外, 对数律区域范围比DNS结果小。图 2 (b)与(c)分别为流向亏损速度$({\tilde U_e} - \tilde U)/{\tilde U_\tau }$与温度平均值的比较。我们可以看出, 不同时刻BL模型计算的剖面是相似的, 但是在边界层内大部分区域, BL模型的结果与DNS结果相差较大。

2.2 马赫数为2.5的超声速钝锥湍流边界层

对比数据来自文献[3], 为超声速零攻角钝锥湍流边界层的空间模式直接数值模拟结果。来流马赫数为2.5, 半锥角α为5°, 以锥头半径R为特征长度的雷诺数为5000, 来流温度为79K, 壁面取绝热条件。首先用NND格式求解包括激波的钝锥绕流问题, 再选取靠近壁面的一小块区域, 用高精度格式模拟湍流流场的空间演化。

流场中各物理量均以无穷远参数无量纲化, 长度单位以R无量纲化。计算域如图 3所示。DNS计算网格为1341×66×151, 计算域为484.5×23.7×12.6°。壁面采用无滑移绝热条件, 上边界给定, 展向为周期条件, 流向入口给定三组有限幅值T-S波, 出口用嵌边区边界条件。

用BL模型计算时实际上只需计算一个轴对称定常问题, 即求解一个x-o-y子午面的解。计算网格数为135×66, 计算域为484.5×23.7°。从DNS数据中可得到转捩位置x₀, 其后引入湍流模型求解。

图 4 (a)为C_f曲线的对比结果, 其中BL1与BL2曲线分别为以x₀=889与x₀=943处(分别为DNS数据的转捩起止位置)为转捩点的计算结果, 可以看出, 在充分发展湍流阶段, 两条曲线均与DNS吻合很好。图 4 (b)给出的是壁面摩擦力沿流向积分所得值, 纵坐标为${f_w} - {f_{w, {\rm{ inlet }}}} = \int_{{x_{{\rm{inlet }}}}}^x {{C_f}} \left( {2{\rm{ \mathsf{ π} }} \cdot x \cdot \sin \alpha } \right){\rm{d}}x$, 其中f_w为从锥头到当地的总摩擦力, x为流向坐标, 以圆锥母线与中轴线的交点为坐标原点, 下标inlet表示计算域入口处的值。可见, BL1曲线与DNS结果重合得非常好, 但是, 这种重合是x=900-1100之间的误差积分后被抵消所致, BL2曲线与DNS的差别也是对转捩过程的忽略引起的。

图 5 (a)与(b)分别为排移厚度及形状因子的流向分布的比较结果, 其中BL1与BL2曲线定义同图 4。可以看出, 排移厚度反而是BL2曲线与DNS吻合较好, 并且增长率也基本相同, 但是形状因子有一定差别。

可见, BL模型对湍流区壁面摩擦系数的模拟是与DNS结果吻合的, 转捩点选取在哪对它影响不大。但是壁面总摩擦力与边界层厚度都是与历史有关的整体量, 转捩点不同会直接影响其在湍流区域的值。值得说明的是, 在本算例中, 壁面总摩擦力是BL1 (转捩点选在转捩开始处)算得准, 边界层排移厚度是BL2 (转捩点选在转捩结束处)算得准, 说明这种差异并不能简单地通过改变转捩点的选取来消除。

以下计算均取x₀=943。图 6 (a)给出的是x=1214处壁面律的分布。可以看出, 在对数律以下, BL模型的结果与DNS结果很接近, 只是外区有所不同。图 6 (b)与(c)分别为流向亏损速度与温度平均值的比较。可以看出, BL模型计算的在不同位置的剖面是相似的, 但是在边界层内大部分区域, BL模型与DNS结果有一定差别。特别是, 二者给出的壁面温度有显著差别。图 6 (d)进一步给出了壁温沿流向的分布。再对比图 2, 说明BL模型因为以混合长度为出发点, 故能很好地模拟内层壁面律, 但却不能全面地反映湍流场。

2.3 马赫数为6的高超声速钝锥湍流边界层

以上数据都是绝热壁面的对比, 工程上常出现等温壁面的情况, 对湍流壁面传导热量的计算也是工程上常关注的一个指标。

本算例DNS数据为带冷却壁的高超声速零攻角钝锥湍流边界层空间模式直接数值模拟的结果, 来自文献[4]。来流马赫数为6, 半锥角α为5°, 以锥头半径R为特征长度的雷诺数为10000, 来流温度为79K, 壁面温度定为294K。计算方法同2.2, 计算网格为1461×151×101, 计算域为619×20×17.5°。

本文用BL模型计算时同样把原问题化为求解x-o-y子午面平均场的二维轴对称问题。计算网格数为147×151, 计算域为619×20。

图 7 (a)为C_f曲线与DNS的对比结果, 其中BL曲线为以x₀=844为转捩点的模型计算结果, 可以看出, 虽然从数值上看BL模型结果与DNS有所不同, 但差别并不大, 并且从趋势上看趋于一致。阻力积分值(图 7 (b))也很接近。

图 8为壁面热传导系数q_w的流向分布曲线, 可以看出, BL模型比DNS结果低, 并且两条曲线基本平行。从图 8 (b)可以看出, 模型预测出的热传导量积分值${Q_w} - {Q_{w, {\rm{ inlet }}}} = \int_{{x_{{\rm{inlet }}}}}^x {{q_w}} \left( {2{\rm{ \mathsf{ π} }} \cdot x \cdot \sin \alpha } \right){\rm{d}}x$(其中Q_w为从锥头到当地壁面传递的总热量)比DNS结果低。

可见, 同样是壁面处的局部统计量, BL模型计算出的平均温度梯度却不能像平均速度梯度那样准确。再结合前面已给出的对于绝热壁而言, BL模型和DNS计算出的壁面温度差别较大, 说明BL模型更多地关注流场中速度的求解, 而这样简单的模型是无法同时精确地预测动力学与热力学的物理量的。

2.4 对转捩位置预测的检验

对转捩位置的预测在工程上有重要的意义。但是, 目前尚无一种普适的预测转捩方法。工程上常借助一种半经验方法——e^N法预测转捩, 但它只适用于环境扰动很小的情况, 如在高空飞行的飞行器。BL模型在提出时同时也建议了一种转捩预测方法, 现对此方法做一次检验。

选取平板边界层的Blasius解为层流基本流, 用公式(5)预测转捩位置。来流马赫数取为4.5, 气体参数取30000m高空的空气参数。用文献[1]所建议的用BL模式转捩预测的方法, 得到以排移厚度δ^*为特征长度的转捩雷诺数(Re_δ) _tr=2160, 相当于该剖面到平板前缘距离3.63cm。而用线性稳定性理论得到的临界雷诺数(Re_δ) _c=2770 (第二模态T-S波)。

参考文献[1]中的算例1, 马赫数为2.85的情况下, 转捩雷诺数(Re_L) _tr≈30000, 此雷诺数是以当地到入口前缘的距离L为特征长度的, 转化为以排移厚度δ为特征长度, (Re_δ) _tr≈2630, 而它的临界雷诺数(Re_δ) _c≈11000 (由于马赫数较小, 此处为第一模态T-S波)。

因此, 用文献[1]所建议的方法预测转捩, 给出的是一种bypass转捩机制。而bypass转捩严重依赖于背景扰动的大小及频率。文献[1]中算例1的结果似乎和实验结果基本吻合, 而实验应是在风洞中做的。众所周知, 不同风洞的转捩实验结果可以相差悬殊, 其原因就是不同风洞的背景扰动差别很大。因此, 文献[1]中所建议的预测转捩方法无普适性。其结果与实验结果基本吻合的原因, 可能是其参数选择的结果。

3 结束语

BL模型是一种概念简单的湍流模型, 其计算量小, 并且对网格的依赖性不强, 广泛应用于工程上边界层湍流问题的计算。本文通过可压缩湍流边界层的三个典型算例, 对比了BL模型与DNS的计算结果, 得到以下结论:

(1) 对于湍流区的壁面摩擦系数, BL模型能计算出与DNS吻合得很好的结果, 并且对网格的依赖性不强; 但是, 转捩点的选取对其结果有一定影响, 并且这种影响不能通过改变转捩点的选取而消除。

(2) BL模型对热力学平均量的计算结果不像其对动力学平均量的计算结果那样精确, 它不能准确地预测定温壁面的热传导系数、总传导热量及绝热壁的壁面温度。其原因可能是湍流Prantdl数假定为常数与实际相差较大。

(3) BL模型计算的内层速度经过Van Driest变换后, 能很好地满足不可压缩湍流的壁面律, 这是由于该模型在计算内层涡粘系数时基于混合长概念; 但是, 对比平均流剖面, 它们虽然在湍流区同样具有相似性, 但是与DNS结果有一定的差别, 不能完全反映真实的湍流场。

(4) BL模型的转捩判定准则是针对某一风洞实验提出的, 它不适用于预测高空飞行器的转捩, 除非另外选择有关参数。