高速旋转流场中的网格效应

盛 誉 路 昕 凌乃阳 周陈龙（核工业理化工程研究院，天津300180）

0 引言

近年来，随着计算机科学技术的迅速发展，计算流体动力学（CFD）[1]得到了越来越广泛的应用，已成为许多工程领域流动问题研究、设计的主要手段之一。但是，随着发展的逐渐深入，CFD也面临着越来越多的困难，特别是CFD中的网格效应问题，迄今为止很难给出一个明确的答案，常常只能通过经验解决。

CFD是通过数值方法求解在空间域上离散得到的控制方程，而在空间域上离散控制方程必须使用网格[4]。可见，网格分布将直接影响数值计算的精度和稳定性。对于旋转流体中的一般工程问题来说，雷诺数不大，网格效应对计算结果的影响不会太明显，在某些条件下甚至可以得到精确解。但是，对于高速旋转流场，径向压强和密度成指数分布。为了确定流场，数值求解所需网格的分布应与流场的物理规律相协调。在为数不多的关于数值计算网格效应的研究论文中，至今尚未发现有针对这种高雷诺数下旋转流场网格效应问题进行的研究，而如何确定合适的网格、近壁网格尺寸多少，以便使计算的精度高、可靠性好，也未有明确的答案。

为解决以上问题，文章采用有限体积法，通过求解柱坐标系下的流体动力学方程[2]，对高速旋转状态下圆筒空腔内部三维流场的网格尺度效应问题进行数值研究。

1 计算域模型与网格划分

本文研究的是空筒的旋转问题，采用的模型如图1所示。计算模型为半径ra=0.06 m、高h=0.01 m。针对其内部高速旋转的流场研究侧壁附近网格对流场计算准确性的影响问题。

图1 流场计算域模型

网格分布直接影响到数值计算的精度和稳定性，采用了结构化的六面体网格，如图2和图3所示，根据引言对流场的分析可知，流动参数在轴向、角向梯度不大，径向梯度很大，且大部分气体都集中在侧壁附近。因此，网格划分的原则是在近侧壁10%半径的区域内单独分区，且径向从中心轴线到侧壁逐渐加密。

图2 流场计算域网格

图3 流场计算域对称面网格

对近壁区的流动，采用壁面函数[3]的处理方法。壁面函数法依赖于近壁第一层网格，在划分网格时，需要把第一个网格节点布置在对数分布规律成立的范围内，通常通过调节无量纲长度参数y+保证计算结果的准确度[5]。近壁面网格无量纲长度参数y+的定义式为：

其中：Δy为近壁第一个网格节点到壁面的距离；ρ为流体密度；uτ为壁面摩擦速度；τw为壁面切应力；up、uw分别为近壁第一个网格节点处的速度、壁面速度。

表1给出了不同的近壁面Δy所对应的y+值。为了保证不同方案间计算结果的可比性，在改变近壁网格距离时仅对近壁面单独分区的网格密度进行调整[6]。为了保证数值计算的收敛性，一般应保证网格的纵横比（最大边/最小边长度的比率）小于1000，在求解高雷诺数的边界层问题时允许达到105～106。这里流场的雷诺数达到1.2×106，考虑到计算周期，采用的计算网格最大纵横比为2.8×105，网格节点总数均为200000。

表1 不同近壁网格距离的y+值

2 计算方法

2.1 流体动力学基本方程

在高速旋转的空筒内部，气体流动所遵守的流体动力学基本方程可以写成以下的通用形式：

式中，φ为通用变量，代表u、v、w、T等求解变量；Γ为广义扩散系数；S为广义源项。对于不同的方程，Γ、S具有相应的形式，表2给出了这些符号与方程的对应关系，其中，cp为比热容；k为导热系数；S为黏性耗散项。

表2 通用控制方程中符号对应的形式

2.2 初始条件与边界条件

在建立了流体动力学运动基本方程以后，需要确定运动的初始状态和边界状态才具有唯一的形态。初始条件为气体旋转运动的等温刚体条件，即对于空筒内部的气体，气体温度T=T0=300 K，旋转速度ω=ω0=3 000 rad/s。空筒的上下端盖和侧壁定义为壁面边界，壁面旋转速度和温度与空筒内部气体的初始状态保持一致。

采用Fluent16.0中基于密度的求解器中的隐式方法对流场求解。根据流场特性，雷诺数达到106量级，因此，需要采用湍流模型进行求解，根据经验选取Transition SST模型。

3 网格尺度对计算结果的影响

针对近壁面第一层网格在壁面法线方向长度分别为1.0×10-4m（方案1）、1.0×10-5m（方案2）、1.0×10-6m（方案3）、1.0×10-7m（方案4）时的旋转流场进行了数值计算。由于最小网格尺度不同，计算所需的时间步长也有所变化，网格尺度越小[7]，计算时间步长也越小。对于尺度很小的网格，若强制采用大的时间步长，计算容易发散，甚至得到非物理解。方案1和方案2的计算在时间步长Δt=1.0×10-7s的情况下很快收敛，方案3和方案4的计算相当缓慢，方案4的Δt甚至小于了1.0×10-8s，经过长时间的计算得到了收敛结果。

图4显示了四种方案所得到的气体旋转速度沿径向的分布情况。在稳定状态下，圆筒内的气体运动理论上遵循刚体运动的原则，即气体的旋转速度沿径向呈线性分布。从图中可以看出，方案1、2、3所计算出的气体旋转速度并不完全遵守气体沿径向呈线性分布的规律，大部分气体的旋转速度均滞后于相应的理论值。但当近壁面第一层网格法向长度减小至1.0×10-7m（即方案4）时，气体的旋转速度曲线已基本与理论值吻合。图5显示了通过方案4所计算出的气体旋转速度云图。

图4 旋转速度沿径向分布

图5 旋转速度分布云图

在高速旋转的圆筒内部，气体压强沿径向呈指数分布，越靠近外壁，气体压强分布越陡。图6显示了四种方案所得到的气体压强沿径向的分布情况。方案1,2,3所得到的气体压强分布差别不大，大部分区域的气体基本遵守径向指数分布的原则，而在径向上0.6倍半径的区域附近，气体压强明显低于理论值。对比之下，方案4所得到的结果虽与理论值稍有偏差，但其分布曲线已与理论值曲线十分接近。图7显示了通过方案4计算出的气体压强云图。

图6 压强沿径向分布

图7 压强分布云图

表3给出了不同方案间的转子壁面功耗比较，近壁面网格尺度从10-4递减至10-7的过程中，计算出的功耗逐步减小。根据资料，相同转速与气体介质条件下，半径0.06 m，高0.1 m长转子负载运行时，转子壁面消耗的能量约为0.52 W。本文的研究对象可视为长转子轴向的一部分。由于这里圆筒壁面的轴向高度仅为长转子壁面轴向高度的约10%，可知方案1、2、3所计算出的壁面功耗偏大，方案4所计算出的壁面功耗接近实验值。

表3 不同方案间的功耗比较

4 结论

本文采用有限体积法，求解了柱坐标系下的N-S方程，对雷诺数达到106量级的高速旋转流场进行了数值研究。针对不同网格尺度的流场计算结果，采用数值分析、理论、实验相结合的方法进行了对比。得到如下结论：

（1）近壁面网格的法向尺度达到10-7时，计算结果最接近实际流动状态。气体的旋转速度和气体压力在径向上分别接近线性分布和指数分布，转子壁面功耗最接近实验值。

（2）在高雷诺数的流场计算时，应根据实际流动条件来选择合适的近壁面网格尺度。