叶轮机叶片参数化设计与网格生成方法*

诸永定 王 威 王 军 李 斌 尹国庆

（1.宁波方太厨具有限公司；2.华中科技大学能源与动力工程学院）

0 引言

在叶轮模型的几何表达时，通过参数化造型及拟合，可以将一个复杂的叶轮几何模型用若干个简单的控制参数来表达，这样避免了CAD 模型中复杂的三维曲面和离散点模型中庞大的离散点数据，同样也便于改变控制参数来实现叶轮的改型。

CFD 数值分析方法手段是现代叶轮机械设计的特征，其分析精度与效果评价与其内流流道网格划分密切相关。

由于叶轮机流道及叶片形状具有高度三维性以及内部流动的复杂性，网格质量的正交性以及网格在流场中大梯度区域的密度对数值计算结果影响很大，因此高质量的网格生成技术往往是提高求解效率并正确模拟复杂流动的关键。

本文选取NASA Rotor67 作为研究对象，分别采用CST 和3 次B 样条参数化方法来拟合Rotor67 二维叶型型线、积叠线以及子午通道型线。通过对控制点的反算求解，最终实现对Rotor67 的参数化拟合。对于参数化拟合后的叶轮机叶片，采用代数方法生成H型网格，通过给定的边界上的节点分布，插值得到计算区域内网格节点分布，并采用雅各比比率（Jacoby ratio）对生成的网格质量进行评价。为了验证参数化拟合的精确性以及网格的质量，本文计算了设计工况下Rotor67 特性曲线并与实验数据进行比对，论证了方法的可行性。

1 叶轮的参数化造型

通常原始叶片型面是由多个截面二维叶型型线按一定的规律积叠而成的空间曲面[1]，因此本文参数化拟合主要思路是先参数化二维叶型，再通过积叠线构成一个三维叶片方式完成。

1.1 叶片成型的积叠线的参数化拟合

积叠线的参数化是为了更好的控制弯、掠规律，并探索新型的三维复杂造型。为计算方便，本文选择二维叶型的前缘点作为积叠点，叶片的弯掠规律的控制便是通过控制积叠点形成的积叠线的形状来实现的。叶片积叠线属于三维空间曲线，将其分为切向弯和前后掠两条二维曲线分别进行参数化拟合，同样采用3 次B 样条进行参数化拟合。本文采用14 个型值点拟合出Rotor67 的参数化弯、掠积叠线，如图1 所示。

图1 B样条参数化弯、掠积叠线Fig.1 Parameterised shape of stacking line by B-Spline method

1.2 二维叶型型线的参数化拟合

叶型的位置和范围确定后，即可以构造叶型型线并确定其在流面上的位置。对于一般的二维叶型，应根据所给离散点，拟合出二维叶型型线。叶型型线的拟合，越接近实际叶型形状，对后面的流场计算产生的误差就越小，同时为了减少气流经过叶栅时的损失，要求叶栅型面光滑。由于CST 参数化方法具有设计变量少，可调节，设计空间广等优点[2]，本文采用CST 参数化方法对二维叶型进行拟合。

CST参数化下的通用叶型可以表示如下[3]：

N1和N2定义了所表示几何外形具有的基本几何特性，对于双圆弧翼型取N1=1，N2=1。形状函数S(ψ)通常采用n阶Bernstein 多项式的加权和作为表达式，即：

式中，n为多项式中的指数；i为多项式阶数；bi(i=0,1,2,…n)为Bernstein系数。

对于三维叶片可由多个叶型型面组成，如图2 所示。因此可定义三维叶片的方程为：

图2 CST 叶型截面定义Fig.2 Blade sectional definition by CST method

由以上分析可知，叶型型线的CST 参数化方法的控制参数为类函数中的N1,N2，形状函数中的Bernstein系数b0,b1…bn以及翼型后缘厚度ζT。为了计算Bernstein 系数，需要得到一组控制点，由于Bernstein多项式的极值点均匀分布在翼型的弦上：

代入式（4）即是所需要的一组控制点[ψi,ζ(ψi)]，因此有：

从而可以求解出一组CST参数b0,b1…bn。

图3展示了Rotor67不同截面位置叶型的CST翼型拟合效果。

图3 Rotor 67叶型的CST拟合效果Fig.3 Rotor 67 profile parameterised CST method

1.3 子午通道型线参数化拟合

子午通道型线包括轮毂线和机匣线，其形状对叶轮机的工作性能有重要影响，合理的子午面流道变化不仅能有效改变叶轮机的流量与压比，还能在一定程度上抑制泄漏涡发展，减小其造成的损失，使得总体性能获得提升。

图4 B样条参数化子午通道型线Fig.4 Parameterised shape of meridional plane by BSpline method

NASA Rotor 67 的子午通道型线如图4 所示，可以看出轮毂沿流向逐渐增大，机匣沿流向逐渐减小，由于影响叶轮机械性能的主要因素在轮毂和机匣的下压与上抬[4]，因此将控制参数选择在径向变化量较大的地方，其它数据点根据插值得到。为了比较好的拟合端壁型线，并且比较少地采用拟合参数，本文选择5 个控制点的三次B样条曲线分别拟合轮毂线和机匣线。图4为Rotor 67子午通道型线参数化拟合效果。

通过分别对叶轮叶片和子午流道型线的参数化拟合，最终实现对Rotor 67的参数化拟合，如图5所示。

图5 Rotor 67参数化拟合效果Fig.5 Parameterised 3D model of Rotor 67

2 叶片网格生成

本文采用代数方法生成叶片计算域的H 型结构网格，根据边界节点通过直接的代数计算与插值得到通道中计算区域内结构网格节点分布，生成分块结构化网格。

先计算叶片区子午面上网格点坐标()z,r。通过上节参数化得到的叶型截面数据和轮毂、机匣的数据，用代数方法求得叶片前、尾缘与轮毂、机匣的交点，然后根据预先设置的叶片区轴向与径向的节点分布设置，可以得到叶片区子午面上边界的节点分布，进而采用代数插值得到叶片区子午面上网格分布[5]，如图6所示。同样的方法可以得到叶片区上游和下游子午面上的网格。

图6 叶片区子午面网格分布Fig.6 Grid in the meridional plane

将叶片区子午面上的节点信息与叶型截面数据联立求解，可以分别得到叶片吸力面与压力面的网格节点分布，如图7所示：

图7 叶片表面网格分布Fig.7 Grid in the blade surface

将生成的子午流面上的节点根据设定的周向节点分布进行周向旋转，可以得到整个流道网格点的轴向坐标、周向坐标和径向坐标()z,r,x，图8 给出了三个流面上的网格生成情况。

图8 S1，S2，S3流面的生成网格Fig.8 Grid in the S1,S2,S3 stream surface

在流动计算的网格生成中，在保持某一方向上网格总数不变的前提下，为了捕捉局部地区高梯度的物理量变化，须要加密局部网格，同时为了保证计算的精度，还须使网格的跨度光滑过渡。对于周向网格和径向网格节点分布要求两端密而中间稀，本文采用如下双曲正切关系式来控制栅距内周向网格及轮毂和机匣的附面层径向网格的疏密，即

式中，M为节点个数；s1，sM分别为起始与终止点；α，β分别为阻尼和伸展因子；当α=1,β=1 时，网格节点均匀分布；当α ＞1,β ＜1 时，网格节点在两端对称加密；α越大，β越小，则两端节点密度越大。

同样，对于叶轮进口截面至叶栅通道进口和叶栅通道出口至叶轮出口截面的网格线上的节点分布要求从一段向另一端节点分布由密到疏，本文采用下列形式来进行节点分布，即

为了保证连接处网格的衔接质量，让控制参数α满足给定第一个网格间距的已知条件：

其中，Δs=s2-s1；Δq=

图9 为加密后叶片流道的H 型网格。采用网格自动生成算法生成网格后，需要对所生成的网格进行质量检测，评价其有效性。本文选用雅各比比率，即最小雅各比矩阵与最大雅各比矩阵行列式的比值作为六面体网格单元的质量评价标准[6]。单元的雅各比矩阵反映了网格的基本质量指标，包括角度、长度和面积等。

假设xk∈ℜ3为任意六面体单元的第k个节点的方向矢量且xk=[xk,yk,zk]T,k=1,2,…8。xk,i∈ℜ3(i=1,2,3) 表示与该节点相邻的所有节点。ek,i=xk,i-xk(i=1,2,3) 为节点k和其相邻节点的方向矢量按顺序构成三个边矢。节点k的雅各比矩阵行列式的值Jk=det[ek,1ek,2ek,3]，将六面体单元上8 个节点的单位雅各比矩阵行列式的值Jk的最小值和最大值的比重定义为雅各比比率JR，即

其中，Jmin=min(Jk),Jmax=max(Jk)，如果一个单元Jmin≤0 则该单元不满足计算要求，在0～+1 之间，Jmin越大，说明单元的质量越好。

图9 Rotor 67计算网格Fig.9 Computional grid of Rotor 67

为了便于气动热力计算，需要将网格数据转化为计算网格文件，对于一个网格文件，须包含网格节点、网格单元、网格域、网格边界等网格信息信息。ANSYS求解器可以读取多种形式的网格文件，本文采用扩展名为cfx5 格式为网格文件，*.cfx5 网格文件是按照非结构化网格的数据形式存储计算网格的，即网格文件仅存储一系列离散的网格单元，不直接存储网格邻接关系。它可以直接采用ASC-II字码存储网格，方便阅读，同时也易于实现程序输出。

3 计算结果与讨论

为了验证参数化拟合精度以及网格生成质量，生成的网格数据文件导入Fluent 中进行数值计算。Strazisar[7]等人对Rotor 67 进行了详细、高精度的热态几何、流场结构和整机性能试验测试，因而作为流动机理研究和叶轮机CFD校核的标准案例而广为使用。Rotor 67的基本设计参数如表1所示。

表1 Rotor 67基本设计参数Tab.1 Rotor 67 design parameters

将采用上述方法生成的网格导入Fluent 中进行数值计算。本文设置单流道网格节点数为50×50×150（切向×径向×轴向），表2为设置网格节点分布规律的系数值。α,β分别为式(7)中的阻尼和伸展因子，而进出口段流向节点分布规律的系数α0值根据第一个网格间距迭代求得。

表2 网格节点分布系数值Tab.2 Coefficients of node distribution

网格质量如图10所示，由图可知，计算域的网格雅各比比率在0.3 以上，70%以上的网格其雅各比比率均在0.7，网格正交性满足数值计算要求。

图10 网格质量检查Fig.10 Grid quality checking

采用密度基隐式求解，二阶离散格式。湍流模型采用标准SSTk-ω模型，壁面为无滑移边界条件。计算中，与试验条件一致，设置进口总压为101325Pa、进口总温为288.15K。

Rotor 67 在设计转速下通过实验测得的堵塞流量为34.96kg/s，本文数值计算得到的堵塞流量为34.48kg/s，在实验测量结果的误差范围之内。图11是计算得到的Rotor 67 在设计转速下的特性曲线与实验值的对比图。从图中可以看出，数值计算得到的压气机设计转速下的效率与实验值相比具有较为一致的趋势，但数值偏低，偏差基本在5%以内；数值计算得到的不同流量下的总压比和实验值相接近，且变化趋势与实验相一致，计算具有满意的可信度。

图11 Rotor 67设计转速下的特性曲线与实验对比Fig.11 Comparison of simulated and measured overall performance curves for Rotor 67

图12给出了最高效率点附近不同叶高截面相对马赫数的等值线图以及相应的实验结果。经过对比发现，流场中的相对马赫数与实验结果符合较好，能够很好地描述流场细节。

图12 最高效率点附近不同叶高处相对马赫数分布等值线图Fig.12 Relative Mach number contours line of different blade height near peak efficiency

4 结论

1）建立叶轮叶片三维实体的参数化造型一种方法，采用CST和3次B样条参数化方法完成叶片的参数化拟合，运用代数方法建立叶轮叶片六面体网格生成。

2）开发前处理计算程序，可以高效、快速地对叶片进行参数化拟合并对流道进行网格划分，通过边界曲线和边界曲面上的网格节点分布控制整各流道网格，输出网格数据能够被调用，从而能快速完成数值模拟计算，有效缩短CFD前处理时间，为叶片的参数优化设计打下基础。

3）用Rotor 67 转子实例的参数化拟合以及网格划分的计算结果验证本文方法的可行性，拟合精度以及网格生成质量满足计算分析要求。