数值模拟是研究涡分离流动重要方法之一。影响数值模拟结果的精度的因素很多，其中，湍流模型是影响计算结果的重要因素之一^[1]。按照不同湍流模型对N-S方程的求解方法差异，数值模拟可以分为直接数值模拟(direct numerical simulation，DNS)、大涡模拟(large eddy simulation，LES)^[2]和雷诺平均法模拟(reynolds averaged navier-stokes，RANS)^[3]三种方法。DNS与LES法具有较高的计算精度，但计算量巨大，计算时间过长，当前难以应用于大雷诺数下复杂流场的计算^[4]。对于分离流动的数值模拟，主要还是基于RANS方法开展。由于分离流动流动复杂，常规的k-ε与k-ω湍流模型，由于存在各自的缺陷，难以获得精确的计算结果^[5]。结合k-ε与k-ω模型优势的SST模型在当前应用最广，然而运用原始SST模型计算分离流动，其精度仍然无法令人满意。

已有的研究表明，运用SST湍流模型对亚临界雷诺数下的圆柱绕流以及圆柱体分离流动进行数值模拟，当雷诺数大于10⁴左右时，分离流动的横向位移以及圆柱绕流的升力系数的数值模拟结果均小于实验值。运用传统SST模型无法对分离流动现象进行准确的数值模拟，必须加以修正。

当前，很多学者对原始湍流模型进行了研究和改进^[6-11]。任芸针对离心泵的特点，利用RNG k-ε模型和Realizable k-ε模型对SST模型进行了改进，以考虑大旋转和曲率影响^[12]。甘文彪对原始SST模型进行了分离修正与可压性修正，提高了模型预测分离流动的能力^[13]。Zhang等对SST模型进行改进，用以计算海洋表面温度^[14]。这些改进针对各自的问题取得了较好的效果，但无法直接应用于分离流动的模拟。

本文基于SST模型的上述缺陷，根据分离流动的特性，在传统SST模型的基础上增加了能量传递项，提升了该模型对分离流动的模拟精度，并以开源软件OpenFOAM为平台，运用改进的SST湍流模型，对典型雷诺数下的圆柱绕流算例进行计算，验证其性能。

1 SST湍流模型 1.1 原始SST模型

剪应力输运(shear stress transport，SST)模型是由Menter于1994年提出的，其基本思路是通过变形，将k-ε湍流模型与Wilcox两方程k-ω湍流模型通过混合函数结合起来^[15]。在近壁面使用k-ω模型，在远离壁面充分发展的流场中采用k-ε湍流模型，充分结合了k-ω模型对分离流模拟精度高以及k-ε模型对湍流初始参数敏感度低的优势，是当前应用最广的湍流模型之一，其无量纲的形式如下

涡粘性

$ {\mathit{\boldsymbol{v}}_T} = \frac{{{a_1}k}}{{\max \left( {{a_1}\omega ,S{F_2}} \right)}} $

(1)

湍动能

$ \frac{{\partial k}}{{\partial t}} + {u_j}\frac{{\partial k}}{{\partial {x_j}}} = {P_k} - {\beta ^ * }k\omega + \frac{\partial }{{\partial {x_j}}}\left[ {\left( {\nu + {\sigma _k}{\nu _T}} \right)\frac{{\partial k}}{{\partial {x_j}}}} \right] $

(2)

比耗散率

$ \begin{array}{l} \frac{{\partial \omega }}{{\partial t}} + {U_j}\frac{{\partial \omega }}{{\partial {x_j}}} = a{S^2} - \beta {\omega ^2} + \frac{\partial }{{\partial {x_j}}}\left[ {\left( {\nu + {\sigma _k}{\nu _T}} \right)\frac{{\partial \omega }}{{\partial {x_j}}}} \right] + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;2\left( {1 - {F_1}} \right){\sigma _{\omega 2}}\frac{1}{\omega }\frac{{\partial k}}{{\partial {x_i}}}\frac{{\partial \omega }}{{\partial {x_i}}} \end{array} $

(3)

其中, 第二混合函数

$ {F_2} = \tanh \left[ {{{\left[ {\max \left( {\frac{{2\sqrt k }}{{{\beta ^ * }\omega y}},\frac{{500v}}{{{y^2}\omega }}} \right)} \right]}^2}} \right] $

(4)

湍动能k的生成项

$ {P_k} = \min \left( {{\tau _{ij}}\frac{{\partial {U_i}}}{{\partial {x_j}}},10{\beta ^ * }k\omega } \right) $

(5)

第一混合函数

$ {F_1} = \tanh \left\{ {{{\left\{ {\min \left[ {\max \left( {\frac{{\sqrt k }}{{{\beta ^ * }\omega y}},\frac{{500v}}{{{y^2}\omega }}} \right),\frac{{4{\sigma _{{\omega _2}}}k}}{{{C_{dk\omega }}{y^2}}}} \right]} \right\}}^4}} \right\} $

(6)

式中：

$ {C_{dk\omega }} = \max \left( {2\rho {\sigma _{{\omega _2}}}\frac{1}{\omega }\frac{{\partial k}}{{\partial {x_i}}}\frac{{\partial \omega }}{{\partial {x_i}}},{{10}^{ - 10}}} \right) $

(7)

F₁在近壁面处取1，在远离近壁面处取值为0，从而实现k-ε与k-ω模型的切换。相应地，输运方程中的各系数也可由混合函数F₁表示。

$ \phi = {\phi _1}{F_1} + {\phi _2}\left( {1 - {F_1}} \right) $

(8)

式中：${{a}_{1}}=\frac{5}{9},$a₂=0.44，${{\beta }_{1}}=\frac{3}{40},{{\beta }^{*}}=\frac{9}{100}$，σ_k=0.85，σ_ω=0.5，σ_ω2=0.856。

1.2 改进SST模型

已有的研究成果表明，传统的k-ε以及SST模型在模拟分离流动时，无法获得令人满意的结果，计算所得的漩涡脱落强度明显小于实验观测值^[16-17]。Medic认为，标准雷诺平均湍流模型的耗散率方程在预测分离流动时存在缺陷，将会产生较大的计算误差以至于抑制涡的脱落^[18]。在标准雷诺平均湍流模式建立中，涡一般被定义为两种尺度，一种是与时均流动发生相互作用的大尺度涡结构(large eddies)，这类涡通过与时均剪切运动的作用，从时均流动能中源源不断地提取能量以维持湍流的脉动运动。另一种尺度的涡是耗散涡结构(dissipative eddies)，这类涡尺度很小，它们通过粘性起耗散湍流脉动动能的作用^[19]。对此，Younis经过大量的理论与实验分析，认为涡在从大尺度变为小尺度的过程中，也存在能量耗散，且耗散频率与涡泄频率相同^[20]。因此，当漩涡产生时，耗散谱应当变为以下形式：

$ e\left( {\kappa ,t} \right) = \left( {{A^0} + A\left( t \right)} \right){\kappa ^s} $

(9)

式中：A⁰是常数，κ为波数，s为待定系数，A(t)是反应涡从大尺度变为小尺度的过程中，存在的能量耗散的附加项，它在平稳流动中应当趋于0。

根据湍动能k与耗散谱的关系

$ k = \int\limits_0^\infty {E\left( {\kappa ,t} \right){\rm{d}}\kappa } $

(10)

$ \frac{{{\rm{d}}\kappa }}{{{\rm{d}}t}} = - \varepsilon $

(11)

可以推出分离流动中，耗散率ε的表达式应当为

$ \frac{{{\rm{d}}\varepsilon }}{{{\rm{d}}t}} = - {C_{\varepsilon 1}}\frac{{{\varepsilon ^2}}}{k} - \frac{1}{{s + 1}}\frac{\varepsilon }{{{A^t}}}\frac{{{\rm{d}}{A^t}}}{{{\rm{d}}t}} $

(12)

根据上述ε的表达式，并通过定义A(t)的表达形式，Younis针对k-ε模型，建议对耗散率方程中的系数C_ε1进行如下修正(k-ε模型的具体形式见文献[21])：

$ c_{\varepsilon 1}^ * = {C_{\varepsilon 1}}\left( {1 + {C_t}\frac{k}{\varepsilon }\frac{1}{{Q + k}}\left| {\frac{{\partial \left( {Q + k} \right)}}{{\partial t}}} \right|} \right) $

(13)

上述湍流模型对分离流动的计算精度虽有提高，但此模型的稳定性较差，特别对于雷诺数较低的情况，有较大的误差。可将此修正思路应用于SST模型，从而提高该模型对逆压梯度的灵敏度，降低近壁区的计算难度，从而实现对分离流动更准确的预报。

本文借鉴Younis针对k-ε模型耗散率方程的修正思路，根据ε和ω的关系式ε=C_μkω，将k-ε模型耗散率方程的修正项转化为用ω表示的比耗散率的形式，并根据SST模型比耗散率方程的形式进行相应变换，得到新的修正项形式为

$ R = \beta '\frac{1}{{Q + k}}\left| {\frac{{\partial \left( {Q + k} \right)}}{{\partial t}}} \right|{P_k} $

(14)

式中：Q为时均流场中每个单位体积中的湍动能，k为湍动能，P_k为湍动能k的生成项，β′为经验常数，根据数值优化方法^[22]求得。本文经过反复调试，取β′=0.54。将R项添加到原始SST方程的比耗散率输运方程中，得到改进的SST比耗散率方程形式如下：

$ \begin{array}{l} \frac{{\partial \omega }}{{\partial t}} + {U_j}\frac{{\partial \omega }}{{\partial {x_j}}} = a{S^2} - \beta {\omega ^2} + \frac{\partial }{{\partial {x_j}}}\left[ {\left( {v + {\sigma _\omega }{v_T}} \right)\frac{{\partial \omega }}{{\partial {x_j}}}} \right] + \\ \;\;\;\;\;\;\;\;\;\;\;2\left( {1 - {F_1}} \right)\left( {{\sigma _{\omega 2}}\frac{1}{\omega }\frac{{\partial k}}{{\partial {x_i}}}\frac{{\partial \omega }}{{\partial {x_i}}} + \frac{1}{2}R} \right) \end{array} $

(15)

其余的方程与系数均与原始SST模型相同，得到改进SST湍流模型。

2 圆柱绕流数值模型 2.1 计算域及网格

本文在亚临界雷诺数下，分别用原始以及改进的SST湍流模型进行圆柱绕流数值模拟，并对模拟结果进行详细的对比分析。

本次模拟所使用的网格如图 1、图 2所示。其中，坐标系定义为：以圆柱中心点为坐标原点，来流方向为x轴方向，垂直来流方向为y轴方向，z轴方向为圆柱展向。计算域划分为：-6D＜x＜20D、-6D＜y＜6D(D表示圆柱直径)，展向厚度为一层网格。

网格划分为：圆柱周围3D的范围内采用O型网格，在其他区域使用六面体网格。网格数量根据雷诺数不同进行相应的调整，改变贴近圆柱表面的第一层网格厚度，保证y⁺＜1。边界条件为：圆柱表面采用无滑移固壁，入口采用均匀来流，出口采用压力梯度为零，其余均为对称边界。

2.2 计算方法

对SST湍流模型的改进通过开源软件OpenFOAM实现。计算过程中，各项的离散方法如下：时间项采用backward隐式欧拉格式进行离散，对流项采用带限制器的线性差分(total variation diminishing, TVD)格式进行离散，扩散项采用高斯线性守恒格式进行离散。

3 数值模拟结果分析 3.1 基本参数分析

本文选取圆柱绕流稳定阶段的数据，对其进行相应的后处理，得到各工况下，圆柱绕流的平均阻力系数、均方根升力系数以及ST数等基本参数，如表 1所示。

为了验证本文改进湍流模型的性能，分别将两种湍流模型的阻力系数、ST数数值模拟结果与文献[23-24]实验观测结果进行对比，并将数值模拟结果与对应的实验观测结果绘在同一图中，如图 3、图 4所示。

从图 3与图 4中可以看出，改进SST模型的平均阻力系数计算结果较原始SST模型略微偏大，但差异很小，可以忽略。对于ST数，改进SST模型的计算结与原始SST模型的计算结果完全一致。根据本文的理论，改进模型中所添加的能量耗散项的能量耗散频率等于涡泄频率，因此添加的R项不会改变圆柱绕流的漩涡脱落频率，因此ST数不会改变。该数值模拟的结果与理论吻合。总而言之，对于分离流动的平均阻力系数以及ST数，运用改进SST湍流模型的计算结果与原始SST模型的计算结果差异较小，且与实验值吻合度较好。

图 5为两种湍流模型的均方根升力系数数值模拟结果与文献[25]实验观测结果对比图。从图 5中可以看出，当雷诺数较小时，两种湍流模型的升力系数计算结果相近，且与实验值吻合较好，但当雷诺数较大时，运用原始SST模型计算所得的圆柱绕流均方根升力系数始终小于实验观测值，偏差较大，可见在较大雷诺数情况下原始SST模型会产生较大的误差，抑制涡的脱落，从而使得升力系数明显偏小。此时，改进SST模型的升力系数计算结果较原始SST模型的结果显著增大，并且与实验观测值的误差较小。因此，本文改进的SST模型能有效克服标准模型抑制涡脱落的缺陷，显著提升分离流动升力系数的数值模拟精度。

3.2 尾涡与流场分析

图 6与图 7分别为Re=5 000与Re=100 000时圆柱绕流的等值线涡量图，其中，L₁表示尾流区的回流区长度。从图中可以看出，两种湍流模型都能捕捉到明显的“2S模式”尾涡^[26]，并且在雷诺数较低时，尾涡的形态区别不大。然而当雷诺数较大时，改进SST模型模拟的尾流区的回流区长度明显大于标准SST模型的模拟结果。结合图 5升力系数的模拟结果可知，标准SST湍流模型在模拟较高雷诺数的分离流动时，回流区长度不足，涡脱落受到抑制，使得升力系数的计算结果明显偏小。进一步验证了改进SST湍流模型能有效克服标准SST模型在雷诺数较大时，预测分离流动漩涡脱落强度不足的缺陷，大大提升了分离流的模拟精度。

图 8为Re=10 000的情况下，沿圆柱表面周向的时均压力系数分布图。沿圆周方向均匀取48个点，记录各点的平均压力数据，以表示圆柱表面压力系数分布。图中可以看出，改进SST模型周向压力响应的最小值点约为65°处，且压力系数约为-1.2左右，与实验值基本吻合，而标准SST模型周向压力响应的最小值点约为70°处，且压力系数约为-1.3左右，其位置和大小相较于实验结果均有相对较大的偏差。可见改进SST模型对于流场的压力分布具有更高的模拟精度。

4 结论

1) 改进SST湍流模型对分离流动均方根升力系数模拟结果的准确性较原始SST模型有显著提升。

2) 当雷诺数较大时，改进SST模型模拟的尾流区的回流区长度明显大于标准SST模型的模拟结果。

3) 改进SST模型对于流场的压力分布具有更高的模拟精度。

4) 本文所改进的SST湍流模型能有效弥补原始SST模型在模拟较高雷诺数下的分离流动时，漩涡脱落强度偏小的缺陷，对于亚临界雷诺数下的分离流动的数值模拟具有良好效果。