旋转坐标系N-S方程无网格/网格混合算法

曹 骋 ，陈红全 ，张加乐 ，高煜堃

（1.常熟理工学院机械工程学院，江苏 常熟 215500；2.南京航空航天大学航空学院，南京 210016；3.安徽工业大学机械工程学院，安徽 马鞍山 243002）

0 引言

桨叶是涡桨发动机最为重要的气动部件之一，其旋转流场也是计算流体力学（Computational fluid dy⁃namics，CFD）数值模拟的难点。桨叶气动特性的优劣直接影响发动机性能。因此，发展能准确模拟旋转桨叶气动特性的CFD 算法，对涡桨发动机设计与研究具有重要的工程意义。

类似发动机桨叶这样的旋转体绕流在工程中广泛存在，如旋翼飞行器等。传统的CFD 技术在模拟这类流场时会遇到2 个难点。一个是目前的CFD 软件大多基于“网格算法”，涉及计算域网格生成。但对于存在旋转桨叶的涡桨发动机等多体复杂气动外形，要生成合适的整体网格，刻画出局部间隙等几何细节，在网格拓扑约束下并非易事，常需要用到块结构等特殊处理技术；另一个是因旋转产生的动边界问题的处理。目前的求解方法通常基于动网格技术开展，在处理动边界大位移问题时，常需要结合网格重构、网格重叠等特殊处理技术，涉及到与角速度相关的旋转特征物理量的插值运算，这会影响到算法的计算效率与精度。为了突破网格拓扑的制约，Batina等、蒲赛虎、郭曼丽等提出1 类无网格算法，计算域用无网格点离散，既可直接布点，也可用现有的结构或非结构网格点，涉及的空间导数等物理量仅依据当地点云结构逼近确定，其灵活性更适合处理旋转体等复杂的几何外形；为了准确反映旋转角速度等特征物理量，Chen等发展了一种直接求解旋转坐标系下控制方程的方法。旋转运动网格在旋转坐标系可与物体保持相对静止，无须反复重构，旋转特征物理量等也无须插值运算；像涡桨发动机的旋转桨叶、悬停旋翼等非定常绕流，在旋转坐标系下还可被看作准定常，减少了求解的复杂度。基于这些特点，该方法已在涡轮机械和旋翼绕流的模拟中获得广泛应用。但将无网格算法和旋转坐标系下控制方程结合起来应用于求解桨叶这样的旋转体绕流问题，还鲜有文献报道。

本文发展一种基于旋转坐标系的无网格/网格混合算法，用于求解桨叶等旋转物体的黏性绕流问题。基于无网格点云空间导数逼近等方法离散控制方程，并结合隐式LU-SGS 算法，发展形成求解非定常问题的双时间步隐式格式推进求解的混合算法，选取2 维振荡翼型非定常绕流、旋转圆柱黏性绕流和模拟发动机桨叶旋转运动的3 维悬停旋翼绕流进行考核运算，并与文献计算或试验结果进行比较分析。

1 控制方程

2 无网格/网格混合算法

2.1 无网格算法

2.1.1 空间导数拟合

计算域无网格布点离散后，形成点云结构，点云()由中心点及其周边卫星点组成。空间导数基于当地点云逼近确定。沿用文献[12]的做法，选用较为典型的加权最小二乘法，任一函数在点云中心点的1阶拟合导数为

式中：α、β和γ为只与点云几何位置信息相关的无网格系数；f为中心点与第个卫星点连线中点处的函数近似值。

2.1.2 空间离散

易于发现，式（1）与惯性坐标系下控制方程的主要区别在于包含矢通量、和的项及源项。对于前者，用式（3）可得

其中

分裂得到的L与在惯性坐标系下无网格算法对应的离散形式一致，可按惯性坐标系已有方法处理，最大程度地减少对已有惯性坐标系下无网格求解程序的改动。W在虚拟交界面上用中心格式确定， |中的速度分量则取为中心点和卫星点连线中点处的牵连速度分量。

源项中涉及的密度和速度等物理量直接取中心点处的值即可。式（1）中黏性矢通量部分的求解与惯性坐标系下的相同，其中的湍流黏性系数采用Spalart-Allmaras湍流模型计算，不再赘述。

2.1.3 时间离散

对式（1）进行空间离散后得到的半离散方程为

式中：R()为残值，包含了对流矢通量、黏性矢通量和源项经空间离散后的部分。

对于非定常流动，式（8）可通过引入虚拟时间和非定常残值，利用显式双时间步在虚拟时间层上推进，以求得下一物理时间步的解W。本文结合LU-SGS 格式形成隐式双时间步推进算法。沿用Yoon等的思路，最终得到的2步推进格式为

式中：为对下一物理时间步守恒变量值W的逼近，即当虚拟时间层→时，有→W；()和()为点云中心点的卫星点集分割出来的2个子集。

需要注意的是式（12）中包含了源项产生的部分

上述算法将与网格算法集成，结合第2.3节中的跨区交界面处理技术，发展形成基于LU-SGS 隐式推进的无网格/网格混合算法。

2.2 边界条件

以上的速度矢量均为绝对速度矢量。镜像点处的压强和密度取为对应原像点处的值。

远场则基于带源项的特征分析无反射边界条件确定。

2.3 无网格/网格混合算法实施

计算域采用整体网格和物面附近局部无网格离散形成混合算法，实施的关键是无网格和网格分区交界面的处理。本文通过引入无网格和网格对偶点，实现了混合算法绕流信息的跨区交换。为了图示清晰，以2维为例加以说明。对偶点一般取为交界面处2层网格的格点，如图1所示。

第1 层称为无网格对偶点（图1 空心圆点），第2层称为网格对偶点（图1方点）。这些对偶点既被当做点云中的无网格点，也被当作网格单元的格点。作为前者时，纯无网格点（图1 实心圆点）和无网格对偶点周围点云模板完整（图1 虚线连接的点云结构），可用无网格算法计算；作为后者时，纯网格节点和网格对偶点周围网格模板完整（格点格式下），可用格点格式有限体积法求解。2 种算法均将物理量存储于对偶点处，相互间可直接取用以实现跨区信息交换，无须插值。

由上述可知，采用原有方法挑选的对偶点，网格部分需要采用格点格式求解。而在求解非定常流动时需要额外的计算修正，涉及到所谓的质量矩阵运算。为了避免该问题，本文利用无网格取点的灵活性调整了对偶点的选取方法。将对偶点选取在网格区域靠近无网格区的若干层单元的中心，如图2 所示（图中各点含义同图1）。

图1 调整前选取的对偶点

图2 调整后选取的对偶点

调整位置后，存储于对偶点处的物理量亦可认为存储于相应单元的中心，在计算时可被当作相应单元中心处的值直接取用，因此仍然无须插值便可实现跨区信息交换。但同时，除无网格对偶点所对应的单元外（网格模板不完整，但点云模板完整，使用无网格算法计算），其余网格均可统一采用格心格式的有限体积法统一求解，进而避免质量矩阵的计算。数值计算表明，上述处理方法能满足跨区交界面数值通量传递的要求，表现为跨区等值线的连续（见第3.1节图7）。

3 结果分析

本文对上述发展的混合算法编制了对应的计算程序，并进行了多个典型算例的考核运算。这里先给出NACA 0012翼型的俯仰振荡非定常动边界绕流和旋转圆柱黏性绕流算例，再给出旋翼悬停飞行绕流算例。

3.1 NACA 0012翼型俯仰振荡绕流

NACA 0012翼型俯仰振荡绕流算例源自AGARD报告，其中翼型的俯仰可视为旋转运动，已被许多文献用于考核求解旋转体绕流问题。初始时翼型位于平均攻角0.016°处，而后绕25%弦长点作俯仰振荡，最大振幅对应攻角= 2.51°，来流马赫数= 0.755，雷诺数= 5.5×10，试验振荡频率= 62.5 Hz，攻角()随时间变化为

计算采用的局部无网格区和网格区如图3 所示。远场边界落在离物体10倍弦长处，涉及16320个无网格点和4331 个网格单元。为准确捕捉黏性流动的特征，物面附近的无网格区采用了沿法向的布点方式，第1层无网格点距离物面高度为1×10。

图3 NACA 0012翼型周围的无网格区和网格区

升力系数和力矩系数随瞬时攻角的变化如图4、5所示。从图中可见迟滞环现象。图4中的升力系数相对试验值整体偏低，但与文献中的计算值吻合得较好；对于力矩系数，所给结果只在个别试验值处有较大偏差，但整体与试验值和计算值都符合较好。4个典型时刻翼型表面的压强系数分布如图6所示，相应的试验值也在图中一并给出，供比较。不同时刻流场的压强系数等值线如图7 所示。展示出翼型俯仰振荡时，激波等流场特征在翼型上下表面交替演化的情况。

图4 升力系数随攻角变化

图5 力矩系数随攻角变化

图6 不同振荡位置翼型表面对应压强系数

图7 不同时刻流场的压强系数等值线

3.2 旋转圆柱黏性绕流

旋转圆柱黏性绕流在增升减阻、改变旋涡释放的规律、降低流动引起的振动等方面具有重要理论价值，历来为流体力学学者所重视。Coutancea 等也对该流动现象进行了大量观察和测量，为数值模拟提供了丰富的参考依据。这里给出本算例发展的算法获得的运算结果。

设无穷远处来流速度为，方向沿轴正向，圆柱直径为，从静止开始突然以角速度按逆时针方向旋转，圆柱旋转的角速度= 2。其中为转速比，本文沿用文献[21]取=0.5。来流条件为=0.1，=200。设无量纲时间=(，本文计算了∈[0,4]范围内的流场情况。计算采用的无网格区和网格区如图8所示。

图8 圆柱周围的无网格区和网格区

由于外形规则，网格区使用了4500 个结构网格单元，无网格区则填充了11476 个无网格点，远场边界落在离物体10 倍圆柱直径处。物理时间步长取为总物理时长的1/800，即Δ=/800=/200。

不同无量纲时刻的计算流线和试验流场照片如图9 所示。从图中可见，本文发展算法所得结果与试验结果基本一致，表现在捕捉的旋涡在强度、位置等特征上与对应时刻试验中观察到的现象基本吻合。3 个无量纲时刻背风区上的方向速度分布如图10 所示。从图中可见，各时刻的速度分布曲线与试验值均较为接近。

图9 各时刻流线与试验流场比较

图10 背风区y = 0 上x 方向速度

3.3 悬停旋翼黏性绕流

悬停旋翼因绕固定轴旋转，相当于发动机旋转桨叶，因此，给出更具挑战的3 维Caradonna 旋翼模型悬停飞行算例来模拟发动机旋转桨叶绕流。该模型由2 片桨叶构成，桨叶平面呈长方形，展弦比为6。桨叶剖面为NA⁃CA0012翼型，沿展向无扭转。计算状态取为：桨尖马赫数= 0.800（约2268 r/min），桨距角θ= 5°，雷诺数= 3.58×10。该旋转引起的非定常流动在本文采用的旋转坐标系下可视为准定常流动，因此，计算可按准定常进行。采用网格和无网格混合算法，计算区域离散整体用了235379 个网格单元和物面附近1200192 无 网格点，如图11 所示。从图中可见3 维Caradonna 旋翼模型和典型截面无网格和网格离散情况。为了准确捕捉黏性流动特性，法向第1 层无网格布点与物面距离为1×10，远场落在离物体15倍的展向翼型弦长处。

图11 旋翼表面及典型截面无网格和网格离散

在不同展向位置=0.50、0.68、0.80、0.96 的表面压强系数计算值和试验值比较如图12 所示。从图中可见，除了在个别点处与试验值有所偏差外，计算结果与试验值总体吻合，分布趋势一致。

图12 旋翼不同展向位置表面压强系数

4 结论

（1）发展的混合算法是直接基于3 维旋转坐标系下的N-S控制方程，相较于传统的基于固定坐标系下的动网格算法，网格和局部无网格点云随物体旋转而无须重构，旋转角速度等重要的特征物理量也无须插值运算，有效避免了因重构、插值等操作引起的误差；

（2）混合算法绕流信息的跨区传递可利用无网格区点云取点的灵活性，借助于本文引入的无网格和网格对偶点实现。算例图示等值线能跨区光滑过渡，证明本文发展的对偶点调整选取方法切实可行；

（3）发展的算法成功结合了LU-SGS 时间推进，并通过了模拟发动机桨叶旋转运动的3 维悬停旋翼等典型复杂非定常绕流算例的考核运算，数值模拟结果与文献计算或试验结果吻合，展示出算法在模拟发动机桨叶等旋转部件绕流问题方面具有重要的工程应用前景。