【www.zhangdahai.com--其他范文】
叶长亮 王福军 汤 远 陈 俊 郑 源
(1.河海大学能源与电气学院，南京 211100；
2.中国农业大学水利与土木工程学院，北京 100083；
3.上海交通大学船舶海洋与建筑工程学院，上海 200240)
轴流泵具有大流量、低扬程的特点，在农业灌溉领域应用广泛[1-2]。水翼作为轴流泵叶片的简化模型，其绕流特性对于研究轴流泵叶片的绕流分布具有参考意义。文献[3-5]通过研究水翼分离流动特性为轴流泵内部流动特性提供了借鉴，然而，目前对轴流泵叶片附近的转捩流动仍缺乏认识。
现有的针对水力机械的转捩研究，大多采用数值模拟方法[6-7]。文献[8]在研究喷水推进轴流泵特性时，发现采用转捩模型不仅提高了外特性的预测精度，在叶片表面的压力分布和摩阻系数分布上也与试验值误差较小。文献[9-10]在船用导管螺旋桨研究中发现，采用转捩模型能成功地捕捉到层流状态到湍流状态的转换。文献[11]也采用相同的计算策略研究了导管螺旋桨上的转捩流动，得到的结论也验证了考虑转捩的重要性。由于转捩流动的复杂性，针对不同的问题，往往需要结合流动特征对转捩预测模型进行修正。例如文献[12]采用新的关联函数直接建立转捩区域长度参数、临界动量雷诺数与当地流场的关系，修正了SSTγ-Reθt转捩模型，提高了水翼绕流计算精度；
文献[13]在关联函数中引入了雷诺数和马赫数，以体现高超声速流动的可压缩性；
文献[14]建立了横流转捩关联函数，使得γ-Reθt转捩模型具备预测横流转捩的能力；
为了提高对大曲率水翼绕流边界层转捩的预测能力，文献[15-16]提出的转捩起始位置关联函数应用于γ-Reθt转捩模型并取得了较好的效果。
综合分析发现，由γ-Reθt转捩模型和SSTk-ω湍流模型耦合的计算模型(本文称作SSTγ-Reθt转捩模型)广泛应用于边界层转捩预测中，具有计算效率高、适用转捩类型多的特点，本文采用SSTγ-Reθt转捩模型为核心数值模型预测边界层转捩流动。对于轴流泵而言，其内部流动雷诺数高，边界层转捩流动特征复杂，需结合流场特性关联SSTγ-Reθt转捩模型的预测性能。采用广泛用于类比轴流泵叶片的NACA0009、Donaldson修型尾缘水翼和NACA0016水翼等，结合高雷诺数绕流边界层转捩特性开展SSTγ-Reθt转捩模型预测性能研究，并基于流场特性开展转捩模型的修正，以期为高雷诺数叶片绕流边界层转捩预测提供参考。
1.1 模型建立与网格划分
文献[17]对NACA0009钝尾缘水翼绕流特性进行了试验研究，并获得了边界层流场特性。为避免大范围的水流旋转，降低来流的紊流程度，在试验段上游增加了特殊装置使来流更加稳定；
试验段测试的水翼攻角为0°；
入口距水翼前缘0.5L(L表示水翼弦长)，出口距水翼后缘6L。试验结果表明，由于水翼的影响，进口各位置的流线方向平均速度与整个截面的流线方向平均速度存在一定的差异(在2%以内)。不同雷诺数下，试验段进口湍流度约为1%，水翼边界层具有明显的层流边界层区同时观测到自然转捩过程。试验测量的水翼边界层流场结果可为数值计算结果提供参考，从而分析不同湍流模型对水翼边界层转捩预测的影响。
1.2 网格无关性验证
文献[18]开展了SSTγ-Reθt在翼型流场应用中网格的适用性研究，就水翼壁面单元而言，网格在流向上的分辨率为50～150，壁面法向上的分辨率为0.77～0.97，展向上的分辨率为10～40。为了更好地在靠近壁面的边界层区域生成密集的网格，使用了O形网格，图1为计算域离散网格。为了满足SSTγ-Reθt方法的要求，提高计算精度，捕捉流场的基本特征，对近壁区域进行了局部网格细化。X、Y和Z方向的网格分辨率分别为98.6、1.03和38.4。
图1 NACA0009钝尾缘水翼计算域网格
本文对SSTγ-Reθt转捩模型进行了网格无关性研究。4组网格中每组水翼壁面网格的最大y+值都接近1。网格M4是由2 650 156节点构成的粗网格。网格M3在水翼节点数上与M4相同，但在垂直于水翼壁面的方向以及翼型上下游区域被细化，因此节点数达到3 985 230。与网格M3相比，网格M2水翼表面包括翼型前缘和尾缘的节点数量都增加了1倍。网格M1除了沿水翼壁面网格节点数是网格M2的两倍外，其他区域的节点数与网格M2相同。选取水翼尾涡脱落主频作为水翼绕流中关键参数，对比4组网格数量下尾涡脱落主频，如表1所示。尾涡脱落频率比较表明，网格M1和网格M2之间的差异非常小，与试验值的误差均在0.5%以内，考虑到计算效率，因此网格M2将用于进一步的计算。
表1 不同网格数量下翼型尾涡脱落主频对比
为了精确获得水翼绕流的边界层分离和尾迹脱落频率等流动特性，需要足够小的网格单元。本文采用ICEM CFD对计算区域划分六面体网格。控制方程基于有限体积法在同位网格上进行离散，对流项离散后的界面值采用混合有界迎风格式进行插值，最高具有二阶精度；
扩散项离散后的界面值采用中心差分格式进行插值；
对于非定常雷诺平均模拟(URANS)，瞬态项离散采用具有二阶精度的二阶向后欧拉格式；
采用耦合式解法联立求解N-S离散方程组变量，采用多重网格技术以加快迭代收敛速度。时间步长为5×10-5s，平均库朗数小于1。该部分计算工作在中国农业大学的高性能计算中心运行，使用256 GB内存和128个核进行并行计算。为了获得水翼尾迹区脱落涡的非定常特性，在尾迹区设置了监测点。
1.3 控制方程
γ-Reθt转捩模型由转捩起始动量雷诺数Reθt输运方程以及间歇因子γ输运方程构成，其中，Reθt是判断转捩起始的条件，γ用于表征流动状态，公式为
(1)
(2)
式中t——时间ρ——密度
∂Uj/∂xj——j方向的速度偏导
Pγ——间歇因子生成项
Eγ——间歇因子破坏项
μt——涡粘系数μ——粘度
σf——常数，取1.0
Pθt——转捩源项
σθt——常数，取2.0
当转捩发生在层流分离点附近时，由于分离剪切层的湍动能较小，预测湍动能发展至湍流再附着的位置与试验获得的实际湍流再附着位置相差较大。为了解决这一问题，γ-Reθt转捩模型对间歇因子γ进行了修正，修正的主要思路是，一旦层流边界层分离时，允许局部间歇因子γ超过1。这将导致k的增加，因此将会使再附着点的位置更早发生[18]，具体修正公式为
(3)
γeff=max(γ,γsep)
(4)
式中γsep——分离流动计算时采用的间歇因子
γeff——修正后有效间歇因子
Rev——涡量雷诺数
s1——常数，取8.0
Reθc——转捩临界动量厚度雷诺数
Freattach、Fθt——内部关联函数
图2为SSTγ-Reθt转捩模型与SSTk-ω湍流模型计算得出的水翼边界层切向时均速度分布的结果。图中，y为垂直于水翼壁面的距离，h为翼型尾部厚度，Uxave为平行于翼型表面的速度分量，Ute为边界层外缘时均速度。图中表示了相对弦长x/L(x方向相对水翼弦长L的距离)分别为0.1、0.3、0.6、0.8、0.9以及0.99处边界层速度分布情况。为了分析不同位置的速度分布，位置x/L为0.3、0.6、0.8、0.9和0.99处的速度分布曲线沿横坐标分别移动了1、2、3、4、5个单位坐标。可以看出，在水翼前缘，即0.1倍翼型弦长处，两个模型预测的边界层速度分布与试验值较为吻合，但是沿着流动方向，SSTk-ω湍流模型计算的边界层速度分布误差越来越大；
SSTγ-Reθt转捩模型只是在靠近水翼尾部区域，即0.8倍翼型弦长处计算的误差较大。
图2 不同湍流模型下水翼边界层切向时均速度分布
本文边界层厚度为水翼边界层内从壁面(速度为零)开始，沿法线方向当x方向的速度Ux至外部速度Ue的99%位置之间的距离，记δ，即δ=y|Ux=0.99Ue。不同模型计算的无量纲边界层相对厚度δ/h沿水翼翼弦长的分布如图3所示。随着边界层的发展，边界层厚度随着与水翼前缘距离的增加而增加。两种模型计算得到的边界层厚度的差异较大，SSTk-ω模型计算的厚度远大于试验结果，SSTγ-Reθt转捩模型计算的水翼前缘的边界层厚度与试验结果吻合较好。
图3 不同湍流模型下水翼边界层相对厚度δ/h分布
图4为两种模型下水翼边界层形状因子H12分布，形状因子H12定义为边界层位移厚度δ1与其动量厚度δ2之比。一般认为，当H12超过2.6时，边界层为层流；
当H12小于1.5时，边界层为湍流；
对于H12在1.5～2.6范围时，边界层处于转捩状态。可以看出，SSTγ-Reθt转捩模型预测的形状因子在x=0.2L后偏离试验值，并在0.65L附近下降至1.5，这说明该位置处已完成转捩，比试验中观察到的0.85L转捩结束较为提前。
图4 不同湍流模型下水翼边界层形状因子H12分布
进一步分析SSTγ-Reθt转捩模型对不同雷诺数下边界层转捩的预测效果，图5给出了边界层完全转捩位置xt随雷诺数ReL的变化。可以看出，在雷诺数低于1.6×106时，边界层转捩位置与实测值比较一致，而在雷诺数高于1.6×106时，随着雷诺数的增大，边界层转捩位置逐渐偏离实测值，说明SSTγ-Reθt转捩模型对高雷诺数水翼边界层转捩预测效果不佳。
图5 不同雷诺数下边界层完全转捩位置
为进一步探究SSTγ-Reθt转捩模型在高雷诺数下预测效果不佳的原因，给出图6所示数值计算得到的雷诺数为2.0×106时，湍流度Tu沿流向的发展，其中负值表示水翼头部的上游位置。可以看出，沿着流动方向，湍流度Tu快速衰减，当到达水翼前缘时，湍流度未达到1%，这与试验所测不相符。因此，高雷诺数下计算域进口必需设置较大的Tu，但过大的Tu会加速间歇因子γ预测值的增加，导致边界层转捩过程发展过快，文献[18]也指出，过大的Tu会导致预测的壁面摩阻系数偏离层流值。因此高雷诺数下自由流场中湍流度衰减过快是导致转捩位置预测不准确的主要原因。
图6 湍流度沿流向发展
为解决高雷诺数条件下，SSTγ-Reθt转捩模型(后续称为原转捩模型)预测边界层转捩提前的问题，对原转捩模型进行基于流动特征的修正，以提高在高雷诺数下边界层转捩预测的精度。
在原转捩模型中，Reθt方程的计算与自由流的湍流度有关。当地自由流湍流度预测的精度直接影响了转捩位置的预测精度。为了解决自由流动中湍流度衰减预测误差大而导致壁面附近湍流度预测精度较低的问题，文献[7]提出了环境源项法来控制自由流动中湍流度的衰减，即在指定区域保持环境湍流度和湍流耗散率不变，在其余区域按给定的衰减率自由衰减。文献[19]根据这一思想，对SSTk-ω湍流模型输运方程源项进行修改，公式为
(5)
(6)
式中ka、ωa——指定区域的环境湍动能、环境湍流比耗散率
τ——应力
β*、β3——湍流模型系数
γ3、σk3、σω3——湍流模型常数
νt——运动涡粘系数
F1——湍流模型的函数
ui、uj——i、j方向的速度分量
ω——湍流比耗散率
σω2——常数，取0.856
根据源项均衡态假设[20]有
(7)
为此，根据NACA0009钝尾缘水翼试验，选取试验中来流湍流度1%为参考对ωa取值进行校准，并定义涡粘系数比RT=ka/(νωa)对ωa进行了无量纲化，ν表示运动粘度系数，具体取值如表2所示。改变ωa，仍保证收敛精度高于10-4，从而保证计算精度。
表2 湍流度1%下环境湍流比耗散率标定设置
选取形状因子这一参数，对不同ωa取值下的结果进行分析，由图7可以看出，与试验值相比，改变ωa取值，在0.2L处H12取值变化不大，而后同一位置处H12的取值随着ωa的增加而减小；
可以看出当ωa为4 000，即涡粘系数比为15时，H12整体分布与试验较为一致。因此，对于同一湍流度下不同雷诺数流动，只需根据RT为常数即可确定ωa。
图7 1%来流湍流度下不同环境湍流比耗散率预测的边界层形状因子分布
间歇因子γ是表征流动的一个重要参数。在边界层层流区γ=0，而在边界层外γ=1。在转捩区域，随着湍流强度的变化，γ介于0和1之间。图8为采用不同涡粘系数比RT计算得出的水翼边界层间歇因子γ的分布。当ωa大于12 000 s-1，即涡粘系数比RT小于5时，边界层间歇因子γ在0.4L位置处开始后迅速增加，转捩开始发生，随后在0.6L发展为全湍流边界层。当涡粘系数比为15时，间歇因子在0.8L处增加至1，即在该位置处完成转捩，与试验值相符。继续增加涡粘系数比，预测的转捩位置逐渐向后移动。当涡粘系数比大于20时，转捩结束位置的预测值稳定在0.9L，预测的转捩区长度大于试验观测值。
图8 不同涡粘系数比下间歇因子γ分布
在弦长雷诺数ReL=2×106条件下，验证了基于环境源项法修正的转捩模型能够有效预测转捩位置。进一步针对不同弦长雷诺数，采用上述方法，获得不同湍流比耗散率以及涡粘系数比，如表3所示。可以看出，随着雷诺数的增加，对应的涡粘系数比也随之增加。
表3 不同雷诺数下指定区域环境湍动能ka和湍流比耗散率ωa设置
图9为不同雷诺数下边界层转捩位置，通过与试验结果的对比，虽然在高雷诺数区间预测的转捩位置较实测值提前，而在低雷诺数区间预测的转捩位置则稍有延后，但整体来看，修正模型在很大程度上提高了对转捩位置预测的准确性。
图9 不同雷诺数下边界层转捩位置
因此，在运用该修正转捩模型时，首先确定来流湍流强度为1%。然后根据来流速度或者弦长雷诺数计算ka，再根据前文给定的推荐值选取ωa用于后续计算。以ReL=2.6×106为例，计算得到的ωa为5 200 s-1，涡粘系数比为19.5，因此，选择该值进行转捩数值计算，图10为边界层形状因子分布，可以看到修正模型预测的从0.7L处迅速减小，边界层迅速从层流状态过渡到湍流状态，这与试验测量结果是比较吻合的。
图10 边界层形状因子分布(ReL=2.6×106)
4.1 Donaldson修型尾缘水翼
Donaldson修型尾缘水翼在钝尾缘水翼基础上，采用一条45°直线和一条三次多项式曲线对尾缘形状进行了修改[17]，由于试验均在同一装置下进行，本算例的计算域及边界条件与NACA0009钝尾缘水翼算例基本一致，计算域网格如图11所示。
图11 Donaldson修型尾缘水翼计算域网格
文献[17]对该水翼在不同雷诺数下的流动进行了试验研究，试验测量结果可为本算例验证提供参考。将基于环境源项法修正模型应用于Donaldson修型尾缘水翼绕流计算中，不同雷诺数下环境湍动能和湍流比耗散率的设置如表4所示。
表4 不同雷诺数下环境湍动能和湍流比耗散率设置
图12为翼型弦长雷诺数为2×106时吸力面边界层x方向时均速度分布，可以看出，修正模型预测的边界层厚度略高于试验值，但相对原始模型已有一定的提高。
图12 吸力面边界层时均速度分布(ReL=2.0×106)
图13给出了不同雷诺数下尾迹区尾涡脱落频率，原转捩模型和修正转捩模型预测结果均与试验结果有较大偏差，其原因可能与尾缘修型对尾涡脱落的影响机理有关，试验研究发现，尾缘修型后，原本交替排列的脱落涡会产生错位，导致尾缘上下脱落涡产生相互作用，影响尾涡脱落频率。但整体来看修正模型提高了对尾涡脱落频率的预测精度，相比原转捩模型最大提高约8%。
图13 不同雷诺数下尾迹尾涡脱落频率
4.2 NACA0016水翼
文献[21]在美国海军大型空化水洞中进行了NACA0016的水翼试验。水翼弦长L为2.133 6 m，顶面呈拱形，底面平行于试验段两侧，试验段宽度为1.43L。本算例采用二维计算，为了提高近壁网格的正交性，在水翼周围采用O形网格块来细化网格，壁面网格法向扩展比为1.1，得到水翼表面的最大y+约为1，最终划分的网格总数约1.7×105，计算域网格如图14所示。
图14 NACA0016水翼计算域网格
将基于环境源项法的修正模型应用于NACA0016水翼绕流流动，对弦长雷诺数ReL=1.7×107的流动进行计算，计算得出ka=0.001 35、ωa=2 249.7 s-1。图15为距水翼前缘0.93L处，原转捩模型、修正转捩模型计算边界层速度值与试验值对比图，修正转捩模型由于准确预测了此处的层流边界层，计算的边界层速度分布以及边界层厚度均与试验结果比较一致。
图15 x=0.93L处边界层速度分布
图16为数值计算得到的雷诺数为1.7×107时，湍流度Tu沿流向的发展。可以看出，沿着流动方向，采用修正转捩模型，湍流度Tu快速衰减，当到达水翼前缘时，湍流度约为1%，与试验所测相符，修正后的转捩模型较好地保证了水翼前缘的湍流度，进而较为准确地预测出转捩区边界层参数。
图16 模型修正前后计算湍流度沿流向发展
图17为NACA0016水翼吸力面边界层相对厚度δ/L以及摩阻系数Cf的试验值和模拟值分布图。可以看出，在靠近水翼中部附近，原转捩模型与修正转捩模型计算的边界层相对厚度分别为0.001 53和0.000 91，相对误差分别为112.3%和27.57%；
原转捩模型与修正转捩模型计算的摩阻系数分别为0.005 22和0.000 71，相对误差分别为1 100%和74.0%。在水翼尾缘附近，原转捩模型和修正转捩模型预测的相对厚度分别为0.012 03和0.010 80，相对误差分别为16.4%和4.85%；
原转捩模型与修正转捩模型计算的摩阻系数分别为0.002 6和0.002 5，相对误差分别为30%和25%。总体来看，修正转捩模型在边界层相对厚度、摩阻系数的预测精度上比原转捩模型均提高3倍以上。
图17 水翼边界层相对厚度和壁面摩阻系数Cf分布
(1)相比于SSTk-ω模型，采用SSTγ-Reθt转捩预测模型能够显著提高边界层流场分布计算精度。即相同网格条件下，在近壁区边界层流场的预测上，SSTγ-Reθt转捩模型的预测精度高于SSTk-ω湍流模型，尤其在低雷诺数条件下，SSTγ-Reθt转捩模型的预测精度与试验值更为接近。
(2)对于高雷诺数边界层转捩流动，壁面附近湍流度预测精度不够，导致SSTγ-Reθt转捩模型预测精度有所下降，通过引入环境湍动能和环境湍流比耗散率参数，建立湍流比耗散率与雷诺数的关系，修正转捩模型的输运方程环境源项，能够有效提高其对高雷诺数边界层转捩流动的预测精度。
(3)基于SSTγ-Reθt转捩修正模型，针对Donaldson修型尾缘水翼高雷诺数条件下的尾涡脱落频率等典型流动特征预测精度相比原转捩模型最大提高约8%；
针对NACA0016水翼中部转捩区域的边界层相对厚度、摩阻系数的预测精度相比原转捩模型均提高3倍以上，水翼尾缘局部区域的边界层相对厚度、摩阻系数预测精度也有显著提升。
