基于三维动网格的旋转叶轮机械非定常流场的计算方法

摘要

本发明提供了一种基于三维动网格的旋转叶轮机械非定常流场的计算方法，应用于旋转叶轮机械技术领域，该方法基于三维动网格技术，实现了旋转叶轮机械的非定常流场计算，并且具有较传统滑移网格法更快的迭代速度。

说明书

技术领域

本发明涉及旋转叶轮机械技术领域，更具体地说，涉及一种用于旋转叶轮 机械的非定常流场的计算方法。

背景技术

现阶段常规的旋转叶轮机械的非定常流动计算方法主要有滑移网格法，该 方法使用多参考系（Multiple Reference Frame，简称MRF），将叶轮计算域作为 一个滑移子域设置在旋转参考系（非惯性系）中，滑移网格随参考系一起转动 无须重新生成，并可保持初始网格的质量，其余计算域设在惯性系。两个参考 系之间利用滑移界面进行数据对接，从而实现整体流场的计算。但是滑移网格 法使用多参考系之间的数据对接，影响了时间上的连贯性，所以滑移网格法的 迭代速度慢。

解决非定常流动问题的常用办法还有动网格技术，动网格技术主要用来解 决流场形状因边界运动而随时间改变的问题，在工业中的应用主要有偏心泵或 齿轮泵的转动、阀门的开关过程、活塞在缸中的往复运动等等。但动网格方法 在应用到三维情形时，复杂性陡增，变形后的网格往往出现负体积而导致计算 出错终止，因此动网格技术应用目前仍局限于二维或准三维（二维域在法向拉 伸而成）的案例。

鉴于滑移网格法迭代速度慢，现有的动网格技术不适用于旋转叶轮机械的 三维非定常流场计算中，因此有必要对旋转叶轮机械的非定常流场计算方法进 行改进。非定常流动分析对于了解旋转叶轮机械的动态特性、提高其性能和可 靠性，具有重要的科学和工程应用价值。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于三维动网格的旋转 叶轮机械非定常流场的计算方法，该方法实现了旋转叶轮机械的非定常流场计 算，与传统的滑移网格法相比，本发明方法具有更快的迭代速度。

为了达到上述目的，本发明采用如下的技术方案：基于三维动网格的旋转 叶轮机械非定常流场的计算方法，旋转叶轮机械是指通过叶轮旋转所产生的离 心力进行增压的机械；其特征在于，包括如下步骤：

第一步，将旋转叶轮机械中设有叶轮的旋转工作区域以及与旋转工作区域 直接连接的区域定义为计算域，使用机械制图软件构造计算域的三维实体，形 成三维实体文件；

第二步，使用网格划分软件读取三维实体文件，对计算域进行网格划分， 得到三维的初始网格，形成网格文件；

第三步，使用计算流体动力学软件读取网格文件，并进行物性参数设置； 所述物性参数设置包括运动边界的设置；所述运动边界的设置是指，运动边界 的界面设置为叶轮计算域的边界面，运动规律设置为绕轴转动，运动方式设置 为Profile方式；

第四步，使用计算流体动力学软件，依次对计算周期中各个时间步的数值 进行计算；每个时间步的数值计算在计算所得的数值收敛后完成；当一个时间 步的数值计算完成后，通过更新叶轮计算域网格和临近叶轮计算域部分的网格 节点，重构下一时间步网格，重构下一时间步网格采用两种方案：（1）对网格 尺寸和畸变率在容许范围内的区域，采用弹簧光滑法结合动态分层法得到下一 时间步的网格节点，然后通过守恒定律和插值运算从现有网格得到下一时间步 网格的物理量，从而重构下一时间步的网格；（2）对网格尺寸和畸变率超出容 许范围的区域，重新划分网格；下一时间步网格重构后，计算下一时间步的数 值，直到完成最后时间步的计算。

采用本发明方法可实现旋转叶片泵的全三维非定常流场数值计算。使用本 发明方法，动网格计算仅在一个参考系中进行，新旧网格节点的拓扑关系保证 了良好的计算精度和时间上的连惯性；而传统的滑移网格计算因多参考系之间 的数据对接影响了时间上的连贯性、降低了迭代速度。在同样的计算模型、计 算网格、初始条件、边界条件及软件设置条件下，在经过一段启动时间后，本 发明方法与滑移网格的计算结果趋近一致，但本发明方法具有较快的迭代速度。

进一步的方案是：所述第四步中的采用弹簧光滑法结合动态分层法得到下 一时间步的网格节点是指，对网格尺寸和畸变率在容许范围内的区域，采取弹 簧光滑法，将该区域的网格边视为网格节点连接的弹簧，以边界位移作为弹簧 的边界条件，通过求解弹簧的力平衡方程得到下一时间步的网格节点的位置， 同时通过动态分层法，根据边界的位移量动态地增加或减少边界上的网格层， 使计算域网格层保持密度。

所述第四步中的计算周期是指，叶轮计算域旋转360°所需的时间。

所述第四步中的数值收敛的判断方法有两种：（一）以残差值的变化判断； （二）编写程序对数值进行监测。

更进一步的方案是：所述第三步中的物性参数设置还包括：设置计算域入 口为压力边界条件；设置计算域出口为质量流量边界条件；设置非定常流动计 算采用标准k-ε湍流模型；设置时间步长Δt；设置非定常流动计算的初始条件， 所述初始条件采用定常流动的收敛解。

所述第三步中的时间步长Δt的取值范围为：其中，n为叶轮的转 速值，Z为叶片数，k为＞1的整数。

所述第三步中的初始条件采用的定常流动的收敛解是通过转子冻结法计算 得出。

所述第二步中的网格划分软件采用ICEM软件。

所述第三步和第四步中的计算流体动力学软件采用Ansys-Fluent软件。

本发明相对于现有技术具备如下的突出优点和效果：

采用本发明方法可实现旋转叶片泵的全三维非定常流场数值计算。使用本 发明方法，动网格计算仅在一个参考系中进行，新旧网格节点的拓扑关系保证 了良好的计算精度和时间上的连惯性；而传统的滑移网格计算因多参考系之间 的数据对接影响了时间上的连贯性、降低了迭代速度。在同样的计算模型、计 算网格、初始条件、边界条件及软件设置条件下，在经过一段启动时间后，本 发明方法与滑移网格的计算结果趋近一致，但本发明方法具有较快的迭代速度。

附图说明

图1是本发明方法的流程图；

图2是离心泵的计算域的网格示意图；

图3（a）是t=0时计算域表面的局部网格示意图；

图3（b）是t=1Δt时计算域表面的局部网格示意图；

图4（a）是与图3（a）对应的计算域中心截面的局部网格示意图；

图4（b）是与图3（b）对应的计算域中心截面的局部网格示意图；

图5（a）是图4（a）中叶轮出口的的局部放大图；

图5（b）是图4（b）中叶轮出口的的局部放大图；

图6是利用本发明方法迭代计算的残差记录图；

图7是利用滑移网格法迭代计算的残差记录图；

图8是叶轮无量纲径向力F_x′图；

图9是叶轮无量纲径向力F_y′图。

具体实施方式

下面结合实施例对本发明作进一步详细的描述，但本发明的实施方式不限 于此。

实施例

本发明非定常流场的计算方法，应用于旋转叶轮机械，旋转叶轮机械是指 通过叶轮旋转所产生的离心力进行增压的机械，例如离心泵、离心压缩机、混 流泵、轴流泵、轴流压缩机、离心风机、轴流风机等；本发明方法的流程图见 图1，包括如下步骤：

第一步，将旋转叶轮机械中设有叶轮的旋转工作区域以及与旋转工作区域 直接连接的区域定义为计算域，使用机械制图软件（如Pro/E软件）构造计算 域的三维实体，形成三维实体文件；

第二步，使用网格划分软件（如ICEM软件）读取三维实体文件，对计算域 进行网格划分，得到结构/非结构混合网格初始单元，形成网格文件；

第三步，使用计算流体动力学软件（如Ansys-Fluent）读取网格文件，并 进行物性参数设置；

物性参数设置包括设置运动边界，运动边界的设置是指，运动边界的界面 设置为叶轮计算域的边界面，运动规律设置为绕轴转动，运动方式设置为 Profile方式；运动边界的设置和控制是动网格计算的重点；对于已知运动规律 的运动边界，需要定义运动边界的运动方式，一般可采用Profile（动边界文件） 和UDF（用户自定义函数）两种方式来控制；Profile方式适合较为简单的运动 规律，如平移和转动；对于比较复杂的动边界运动，则需要采用UDF方式进行 控制；虽然本方法的运动边界是叶轮计算域的边界面，曲面形状比较复杂，但 是运动规律却是简单的绕轴转动，因此本方法采用Profile方式定义运动边界 的运动方式；

物性参数设置还包括，设置计算域入口为压力边界条件；设置计算域出口 为质量流量边界条件，可根据工况点流量设置质量流量边界条件；设置非定常 流动计算采用标准k-ε湍流模型；设置时间步长Δt，Δt不大于相邻的叶轮叶 片掠过同一位置的时间差，时间步长Δt根据叶轮的转速值n和叶片数Z确定， 其中，n为叶轮的转速值，Z为叶片数，k为＞1的整数；设置非定常流 动计算的初始条件，初始条件的值采用旋转叶轮机械内定常流动的收敛解，定 常流动的收敛解可通过转子冻结法计算得出；

第四步，使用计算流体动力学软件（如Ansys-Fluent软件），以叶轮计算 域旋转360°为一个计算周期，依次对计算周期中各个时间步的数值进行计算， 每个时间步的数值计算在计算所得的数值收敛后完成；当一个时间步的数值计 算完成后，通过更新叶轮计算域和网格节点，重建三维实体，从而重建下一时 间步的网格，计算下一时间步的数值，直到完成最后时间步的计算；数值收敛 的判断方法有两种：（一）以残差值的变化判断；（二）编写程序对数值进行监 测；在不同的时间步中，因运动边界的运动，计算域发生了变形，因此在计算 前要获取下一时间步的网格，通过更新叶轮计算域网格和临近叶轮计算域部分 的网格节点，重构下一时间步网格，重构下一时间步网格采用两种方案：（1） 对网格尺寸和畸变率在容许范围内的区域，采用弹簧光滑法结合动态分层法得 到下一时间步的网格节点，然后通过守恒定律和插值运算从现有网格得到下一 时间步网格的物理量，从而重构下一时间步的网格；（2）对网格尺寸和畸变率 超出容许范围的区域，重新划分网格；下一时间步网格重构后，计算下一时间 步的数值，直到完成最后时间步的计算；

采用弹簧光滑法结合动态分层法得到下一时间步的网格节点是指，对网格 尺寸和畸变率在容许范围内的区域，采取弹簧光滑法，将该区域的网格边视为 网格节点连接的弹簧，以边界位移作为弹簧的边界条件，通过求解弹簧的力平 衡方程得到下一时间步的网格节点的位置，同时通过动态分层法，根据边界的 位移量动态地增加或减少边界上的网格层，使计算域网格层保持密度。

在第四步之后还可以包括第五步，计算后处理；计算后处理包括，显示非 定常流场计算结果，预测稳态的特性曲线，得到压力的脉动值（包括脉动频率 及脉动幅值）等非定常结果。

本发明的突出优点是：采用本发明方法可实现旋转叶片泵的全三维非定常 流场数值计算。使用本发明方法，动网格计算仅在一个参考系中进行，新旧网 格节点的拓扑关系保证了良好的计算精度和时间上的连惯性；而传统的滑移网 格计算因多参考系之间的数据对接影响了时间上的连贯性、降低了迭代速度。 在同样的计算模型、计算网格、初始条件、边界条件及软件设置条件下，在经 过一段启动时间后，本发明方法与滑移网格的计算结果趋近一致，但本发明方 法具有较快的迭代速度。

为验证本发明方法计算旋转叶轮机械的非定常流场是可行性的，采用离心 泵作为实验对象进行计算及分析。选取单级单吸IS型管道离心泵，泵设计工况 参数为：转速n=2900r/min，流量Q=155m3/h，扬程H=64m。工作介质为水， 密度ρ=998.2kg/m3，动力粘度μ=1.003×10-3Pa·s。计算域由吸入管、叶轮及泵 壳组成，叶轮的叶片数Z=5。应用Pro-E建立水泵的三维计算域，使用Gambit 进行计算域网格的划分，得到如图2所示的非结构性网格单元。其中吸入管77200 单元、叶轮195008单元、蜗壳105176单元，网格单元总数为377384，网格节 点总数为83648。

计算使用Fluent流动软件，选取标准k-e湍流模型。采用如下边界条件：1） 进出口条件按流量值给定；2）壁面采用无滑移固壁条件并由标准壁面函数确定 固壁附近流动。由泵转速与叶轮叶片数计算得到叶轮的旋转周期为2.069×10^-2s， 叶片掠过周期为4.138×10^-3s，因此选取计算时间步长Δt=4.138×10^-4s，即一个 叶片掠过周期使用10个时间步进行计算。

对于传统的滑移网格计算，按常规的方法将叶轮域设为旋转的滑移子域 （Moving Mesh）并给定转向和转速，其余计算域设为静止域。

对于动网格计算，使用Profile文件定义计算域边界面的转向和转速。为简 化计算，变形网格仅限于叶轮计算域，将所有计算域在惯性系中设为静止域。 采用以下三种方式实现网格变形：弹簧光滑法（Spring-based smoothing）、动态 分层法（Dynamic layering）以及局部网格重构法（Local remeshing）。弹簧光滑 法将网格边视为网格节点连接的弹簧，以边界位移作为弹簧的边界条件，通过 求解弹簧的力平衡方程得到下一时间步的网格节点的位置，但弹簧光滑法一般 只能处理小变形问题。动态分层法是根据边界的位移量动态地增加或减少边界 上的网格层，使计算域网格层保持一定的密度。网格重构法是对弹性光滑法的 补充，以网格尺寸和畸变率作为评判标准，对变形过大、局部网格发生严重畸 变的区域，则重新划分网格。下一时间步的网格上的物理量通过守恒定律和插 值从现有网格中获得。

图3（a）和图3（b）分别是t=0和t=1Δt时计算域表面的局部网格示意图。 因叶轮计算域表面是流场与旋转固体的交界面，因此叶轮计算域表面的网格除 了随体转动，没有看出明显的变形或局部重构。图4（a）和图4（b）分别是与 图3（a）和图3（b）对应的计算域中心截面的局部网格示意图。为便于观察， 图5（a）和图5（b）分别给出图4（a）和图4（b）中叶轮出口叶片附近的局部 放大效果。由图5（a）和图5（b）可见，计算域内部网格除了随体转动，还出 现了不同程度的变形和局部重构（见方框标示部分）。

图6和图7分别给出动网格和滑移网格迭代计算离心泵非定常流场的残差 记录，图中横坐标表示迭代步数，纵坐标表示方程迭代计算的残差值，6条线分 别是流动连续方程、动量方程（三个分量）、k方程及e方程残差值随迭代步数 的变化。曲线的每一次脉动代表某一时间步迭代收敛并开始进入下一时间步的 计算。由图可见，与滑移网格相比，动网格具有较快的收敛速度。对同样的工 作时间段（t≈2.0s），动网格的残差上限在10^-1以下（见图6和图7纵坐标上限）， 滑移网格的残差上限则在1.0左右，相差了一个数量级。动网格迭代总步数为 15500（见图6和图7横坐标上限），而滑移网格的迭代总步数则需要47000，换 句话说，在这个算例中动网格的迭代速度几乎是滑移网格迭代速度的3倍。究 其原因，是因为动网格计算仅在一个惯性系中进行，迭代过程中尽管有网格的 变形和重构等不利因素，但新旧网格节点的拓扑关系保证了良好的计算精度和 时间上的连惯性。滑移网格计算则是在多参考系MRF下进行，多参考系之间的 数据对接影响了时间上的连贯性而导致迭代速度的下降。

图8和图9分别给出了按动网格和常规滑移网格方法计算得到的叶轮所受 的无量纲径向力F_x′和F_y′随时间的变化曲线。无量纲径向力的定义为：

$F_i^{\prime }=F_i/\frac{1}{2}\mathrm{\rho A}{U}_{\mathrm{in}}^2,i=x,y$

其中A是叶轮受力总面积，F_x和F_y分别为叶轮所受x、y方向径向力，U_in为 泵入口平均流速。由图8和图9可见，在经历了一段叶轮启动时间（t≈0.02s， 约1个旋转周期）后，径向力值随时间作规则的周期脉动。在任一个叶轮旋转 周期内，径向力出现5次脉动信号，脉动频率与叶轮的叶片数Z相对应。在经 历约5个叶轮旋转周期后（t≈0.1s），两种网格计算数值结果逐渐趋近一致，径 向力脉动进入了相对稳定的阶段，因此可认为从此往后的计算结果基本接近真 实情况。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实 施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、 替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。