基于奇点分布法的轴流叶轮翼型设计方法

摘要

本发明公开了一种基于奇点分布法的轴流叶轮翼型设计方法，它包括以下步骤：S1：计算总体及各流面的基本几何和流动参数；S2：计算各流面基本翼型骨线节点上各分速度及骨线节点位置；S3：循环完成所有流面上的骨线运算，程序结束。翼型骨线上的旋涡密度函数为：

说明书

技术领域

本发明涉及一种轴流叶轮翼型设计方法，特别是涉及一种基于奇点分布法的轴流叶轮翼型设计方法。

背景技术

高扬程、大流量的轴流泵由于其特有的外特性决定了这一泵型在国民经济众多部门有不可代替的作用，属于化工、石油、矿山、冶金等行业生产线上的关键设备。近年里，在市场需求的强势推动下，国内轴流泵制造业发展迅速，产品结构日趋多样化，适应的流量、扬程范围不断延伸，国内轴流泵专兼业的生产厂家超过数百家，产量达数万台。

设计是优秀产品的基本保证，产品性能很大程度上取决于设计质量。目前国内现状是：很多企业，甚至包括一些大型企业，在与客户签订了供货合同后，仍以相似计算法获取产品图样，即寻找一个比转速与待设计产品接近的已有叶轮，根据要求的扬程、流量等给定参数修正原产品的几何尺寸以完成设计。尽管这种做法在目前十分普遍，但其局限性非常明显：由此形成的产品水力性能，包括效率等关键指标，很难超越模拟对象，设计过程中的技术创新很少，难以提高企业整体技术水平，有时要找到一个比转速和外特性都符合要求的已有优秀模型并不容易，等等。

奇点分布法是一种理论严密，理论层次较高，计算难度较大，所含经验与假设成份较少的轴流叶轮设计方法。这一方法较早用于轴流式水轮机转轮的设计，并取得了较为理想的效果。这一事实是奇点分布法的先进性与可行性的直接反映。在长期的设计实践过程中，我们逐步认识到现有设计理论与计算方法仍有不足与缺陷。

根据用户提供的泵的流量、扬程、转速等性能参数设计叶轮叶片的第一步，是要在若干个（一般为5个）与叶轮同轴的圆柱形流面的展开平面上形成翼型。在奇点分布法中，平面翼型是以翼型骨线为依据按一定规律加厚的结果。因此，在展开流面上获得正确的翼型骨线是实现先进、高效产品的基础，大量的实验结果表明叶轮性能主要是由翼型骨线的形态决定的，与翼型的最大厚度和厚度分布关系不大。为得到平面上的翼型骨线，奇点分布法的基本原则是：在各展开为平面的流面上，均匀流中排列一无穷直列叶栅，栅中各翼型是几何全等的。现设想将这一叶栅从平面均匀流流场中抽去，而在原所有翼型骨线上布置相同规律的连续涡，这些涡层产生的平面诱导速度场叠加均匀流场后的合成流场应与原固态所有真实翼型产生的平面流场等同。

泵在设计点的扬程、流量等参数是用户给定的，在一计算平面流面上，为实现这些参数和尽可能小的水力损失，平面上速度分布就必须是特定的。平面流场又是由给定均匀流场与涡层的诱导速度合成的，涡层产生的各点诱导速度是由各翼型骨线布置的相同点涡和旋涡层决定的，因此，涡层沿骨线的强度分布，或者说漩涡密度，必须能正确反映平面流速场的要求。

如前所述，奇点分布法最初成功用于轴流式转轮翼型的设计。也许就是这个原因，目前所有文献中介绍的翼型骨线旋涡密度都是转轮翼型设计经验的提炼，也只能用于转轮翼型的设计。而在一给定平面均匀流场中，在同一坐标系下，叶轮和转轮流面上均匀流流向平面叶栅的方向相反，因而栅中翼型前，后缘点及奇点、驻点位置正好颠倒。由于这两点处的边界条件是不同的，不可能存在同一旋涡分布规律能适应两种边界条件。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于奇点分布法的轴流叶轮翼型设计方法，给出了一种适合轴流叶轮翼型骨线设计的与传统形式不同的新的旋涡密度函数，计算步骤相对简单、计算量相对较小；并通过软件算法以实现奇点分布法巨大的计算量，设计方便、速度快，且可有效提高叶轮的设计及制造精度，适用于工业化设计和生产。

本发明的目的是通过以下技术方案来实现的：基于奇点分布法的轴流叶轮翼型设计方法，它包括以下步骤：

S1：计算总体及各流面的基本几何和流动参数；

S2：计算各流面基本翼型骨线节点上各分速度及骨线节点位置；

S3：循环完成所有流面上的骨线运算，程序结束。

进一步地，步骤S1包括以下子步骤：

S101：确定计算流面个数；

S102：循环计算各流面半径，计算各流面叶栅稠密度l/t、栅距t及弦长l；

S103：根据用户给定的泵在设计点的性能参数，循环计算各计算流面上均匀流大小W∞和流动方向角β∞，进一步计算各流面上均匀流的水平及铅垂分量W∞u、W∞z，

S104：确定各流面的冲角α后，循环计算各流面上翼型安放角β_e；

S105：循环计算各流面上绕单个翼型的环量Γ_B；

S106：以插值方法获取反映孤立翼型和栅中翼型性能差异的系数l_n的插值方程；

S107：根据各流面上已知的翼型安放角β_e和叶栅稠密度l/t确定流面的A₀值，并由此及Γ_B值循环计算各流面的系数A₁。

进一步地，步骤S2包括一个基本翼型骨线节点上各分速度的计算步骤和一个基本翼型骨线节点位置计算步骤：

其中，基本翼型骨线节点上各分速度的计算步骤包括：

S211：将基本翼型骨线节点各点的速度分解为三项：

（1）由基本翼型骨线上的连续分布涡产生的诱导速度V_1u,V_1z；

（2）由除开基本翼型的其余翼型的连续分布涡产生的诱导速度V_2u,V_2z；

（3）均匀流的速度分量W∞u，W∞z；

S212：将一流面的长l翼型骨线等分成6段，在各节点上，由于基本翼型上的连续分布涡产生的诱导速度的两个分量V_1u,V_1z，由下式计算，$>\left(\begin{array}{c}{V}_{1u}=-V_y\sin \mathrm{\beta e}=-(\frac{A_0}{2}+A_1{s}_0/l)\sin \mathrm{\beta e}\\ V_{1z}=V_y\cos \mathrm{\beta e}=(\frac{A_0}{2}+A_1{s}_0/l)\cos \mathrm{\beta e}\end{array}\right);$>

S213：计算各分点S₀上的a(s_o,-3l/6),a(s_o,-2l/6),......,a(s_o,3l/6)值及b值，计算流面翼型骨线各分点上的V_2u,V_2z：

将一流面上的7个分点从下向上编号，编号分别为0，1，……6，在程序中形成一个7×7的方矩数组a(7,7)，数组元素a(m,k)(m,k＝0,1,2,......6)表示编号为m的点涡在编号为k的计算点产生的诱导速度的水平分量；

该方阵关于主对角线反对称，且平行于主对线的直线上的元素相等；

正弦函数和双曲正弦函数都是奇函数，余弦函数和双曲余弦函数都是偶函数，由函数的奇偶性，可以看出a(i,j)＝-a(j,i)，其中i∈[0,6]，j∈[0,6]，每对关于主对角线对称的元素都有这一特性；

首先计算主对角线右上方的全部元素，以u1,z1表示流面骨线上任意两个相邻计算点的横坐标与纵坐标之差，在同一流面上这是两个常数；在变量p一定时（p=0，1，2，……6），q从p+1到6循环，计算p，q代表的两点的横、纵坐标之差2π(u₀-u)/t,2π(z₀-z)/t，进一步计算sh[2π(z₀-z)/t],ch[2π(z₀-z)/t],cos[2π(u₀-u)/t]，将这些值代入式$>a(s_0,s)=-\frac{1}{2}·\frac{\mathrm{sh}\left[\frac{2\pi }{t}(z_0-z)\right]}{\mathrm{ch}\left(\frac{2\pi }{t}(z_0-z)\right)-\cos \left(\frac{2\pi }{t}(u_0-u)\right)}+\frac{t}{2\pi }\frac{z_0-z}{{(z_0-z)}^2+{(u_0-u)}^2}$>中，即可获得方阵右上角全部a(s₀.s)值；由方阵元素的负对称性，形成主对角线元素左下方全部元素；主对线上n+1个元素全部置0；

在利用式：

$>V_{2U}=\frac{1}{2560}·\frac{\mathrm{\pi l}}{t}\left(\begin{array}{c}334a(x_0,-\frac{3l}{6})A_0+630a(x_0,-\frac{2l}{6})A_0+210a(x_0,-\frac{2l}{6})A_1\\ -180a(x_0,-\frac{l}{6})A_0-120a(x_0,-\frac{l}{6})A_1+460a(x_0,0)A_0\\ +460a(x_0,0)A_1-90a(x_0,\frac{l}{6})A_0-120a(x_0,\frac{l}{6})A_0\\ +126a(x_0,\frac{2l}{6})A_0+210a(x_0,\frac{2l}{6})A_1\end{array}\right)$>计算一个流面骨线上的V_2u时，调用方阵元素即可得到a(s_o,-3l/6),a(s_o,-2l/6),......,a(s_o,3l/6)值，以同样方法形成b方程的全部元素；

基本翼型骨线节点位置计算步骤为：在确定了骨线上各分点的全部速度分量后，以公式确定各分点上的β角，并得到最终骨线上各分点位置的坐标。

翼型骨线上的旋涡密度函数为：或$>r\left(s\right)=A_0\sqrt{\frac{1-s/(l/2)}{1+s/(l/2)}}+A_1\sqrt{1-{[s/(l/2)]}^2};$>附加项为：$>f\left(\theta \right)=A_0\mathrm{tg}\frac{\theta }{2},$>或$>f\left(s\right)=A_0\sqrt{\frac{1-s/(l/2)}{1+s/(l/2)}}.$>其中，l表示计算翼型的弦长，s表示弦上计算节点的曲线坐标，A₀与A₁是密度函数r(s)的两个系数，变换量θ由cosθ＝s/(l/2)定义。

这一旋涡密度函数与现有文献给出的密度函数完全不同。

本发明的有益效果是：

1）提供了一种适用于轴流叶轮的基于奇点分布法的翼型设计方法，给出了一种合理的适合轴流叶轮翼型骨线设计的新的旋涡密度函数和附加项，计算步骤相对简单、计算量相对较小；附加项能够形成几何形态和水力性能均令人满意的翼型骨线，旋涡密度函数能够保证置于均匀流的翼型满足应有的边界条件；

2）通过软件算法实现了奇点分布法巨大的计算量，设计方便、速度快，且可有效提高叶轮的设计及制造精度，可适用于工业化设计和生产。

附图说明

图1为平面叶栅结构示意图；

图2为绕翼型骨线环量示意图；

图3为均匀流绕流翼型示意图；

图4为直线坐标系s和直角坐标系u-o-z结构示意图；

图5为基本翼型上点涡产生的诱导速度示意图。

具体实施方式

下面结合附图进一步详细描述本发明的技术方案，但本发明的保护范围不局限于以下所述。

基于奇点分布法的轴流叶轮翼型设计方法，它包括以下步骤：

S1：计算总体及各流面的基本几何和流动参数；

S2：计算各流面基本翼型骨线节点上各分速度及骨线节点位置；

S3：循环完成所有流面上的骨线运算，程序结束。

进一步地，步骤S1包括以下子步骤：

S101：确定计算流面个数；

S102：循环计算各流面半径，计算各流面叶栅稠密度l/t、栅距t及弦长l；

S103：根据用户给定的泵在设计点的性能参数，循环计算各计算流面上均匀流大小W∞和流动方向角β∞，进一步计算各流面上均匀流的水平及铅垂分量W∞u、W∞z，

S104：确定各流面的冲角α后，循环计算各流面上翼型安放角β_e；

S105：循环计算各流面上绕单个翼型的环量Γ_B；

S106：以插值方法获取反映孤立翼型和栅中翼型性能差异的系数l_n的插值方程；

S107：根据各流面上已知的翼型安放角β_e和叶栅稠密度l/t确定流面的A₀值，并由此及Γ_B值循环计算各流面的系数A₁。

进一步地，步骤S2包括一个基本翼型骨线节点上各分速度的计算步骤和一个基本翼型骨线节点位置计算步骤：

其中，基本翼型骨线节点上各分速度的计算步骤包括：

S211：将基本翼型骨线节点各点的速度分解为三项：

（1）由基本翼型骨线上的连续分布涡产生的诱导速度V_1u,V_1z；

（2）由除开基本翼型的其余翼型的连续分布涡产生的诱导速度V_2u,V_2z；

（3）均匀流的速度分量W∞u，W∞z；

S212：将一流面的长l翼型骨线等分成6段，在各节点上，由于基本翼型上的连续分布涡产生的诱导速度的两个分量V_1u,V_1z，由下式计算，$>\left(\begin{array}{c}{V}_{1u}=-V_y\sin \mathrm{\beta e}=-(\frac{A_0}{2}+A_1{s}_0/l)\sin \mathrm{\beta e}\\ V_{1z}=V_y\cos \mathrm{\beta e}=(\frac{A_0}{2}+A_1{s}_0/l)\cos \mathrm{\beta e}\end{array}\right).$>

S213：计算各分点S₀上的a(s_o,-3l/6),a(s_o,-2l/6),......,a(s_o,3l/6)值及b值，计算流面翼型骨线各分点上的V_2u,V_2z：

将一流面上的7个分点从下向上编号，编号分别为0，1，……6，在程序中形成一个7×7的方矩数组a(7,7)，数组元素a(m,k)(m,k＝0,1,2,......6)表示编号为m的点涡在编号为k的计算点产生的诱导速度的水平分量；

该方阵关于主对角线反对称，且平行于主对线的直线上的元素相等；

正弦函数和双曲正弦函数都是奇函数，余弦函数和双曲余弦函数都是偶函数，由函数的奇偶性，可以看出a(i,j)＝-a(j,i)，其中i∈[0,6]，j∈[0,6]，每对关于主对角线对称的元素都有这一特性；

首先计算主对角线右上方的全部元素，以u1,z1表示流面骨线上任意两个相邻计算点的横坐标与纵坐标之差，在同一流面上这是两个常数；在变量p一定时（p=0，1，2，……6），q从p+1到6循环，计算p，q代表的两点的横、纵坐标之差2π(u₀-u)/t,2π(z₀-z)/t，进一步计算sh[2π(z₀-z)/t],ch[2π(z₀-z)/t],cos[2π(u₀-u)/t]，将这些值代入式$>a(s_0,s)=-\frac{1}{2}·\frac{\mathrm{sh}\left[\frac{2\pi }{t}(z_0-z)\right]}{\mathrm{ch}\left(\frac{2\pi }{t}(z_0-z)\right)-\cos \left(\frac{2\pi }{t}(u_0-u)\right)}+\frac{t}{2\pi }\frac{z_0-z}{{(z_0-z)}^2+{(u_0-u)}^2}$>中，即可获得方阵右上角全部a(s₀.s)值；由方阵元素的负对称性，形成主对角线元素左下方全部元素；主对线上n+1个元素全部置0；

在利用式：

$>V_{2U}=\frac{1}{2560}·\frac{\mathrm{\pi l}}{t}\left(\begin{array}{c}334a(x_0,-\frac{3l}{6})A_0+630a(x_0,-\frac{2l}{6})A_0+210a(x_0,-\frac{2l}{6})A_1\\ -180a(x_0,-\frac{l}{6})A_0-120a(x_0,-\frac{l}{6})A_1+460a(x_0,0)A_0\\ +460a(x_0,0)A_1-90a(x_0,\frac{l}{6})A_0-120a(x_0,\frac{l}{6})A_0\\ +126a(x_0,\frac{2l}{6})A_0+210a(x_0,\frac{2l}{6})A_1\end{array}\right)$>计算一个流面骨线上的V_2u时，调用方阵元素即可得到a(s_o,-3l/6),a(s_o,-2l/6),......,a(s_o,3l/6)值，以同样方法形成b方程的全部元素；

基本翼型骨线节点位置计算步骤为：在确定了骨线上各分点的全部速度分量后，以公式确定各分点上的β角，并得到最终骨线上各分点位置的坐标。

翼型骨线上的旋涡密度函数为：或$>r\left(s\right)=A_0\sqrt{\frac{1-s/(l/2)}{1+s/(l/2)}}+A_1\sqrt{1-{[s/(l/2)]}^2};$>附加项为：$>f\left(\theta \right)=A_0\mathrm{tg}\frac{\theta }{2},$>或$>f\left(s\right)=A_0\sqrt{\frac{1-s/(l/2)}{1+s/(l/2)}}.$>这里l表示计算翼型的弦长，s表示弦上计算节点的曲线坐标，A₀与A₁是密度函数r(s)的两个系数，变换量θ由cosθ＝s/(l/2)定义。

如图1所示为一在平面均匀流场中排列的栅距为t，弦长为l的无限平面直列叶栅，这一流面实际是一圆柱形流面沿其一母线剖分展开成的平面，栅中平板实际是待定翼型的弦。

由于最终所获的翼型骨线，与所有实际应用翼型骨线一样，都是微弯的，因而骨线与翼型弦线差别不大，布置在骨线上的连续分布旋涡可以用连续分布在弦上的涡列代替。任选一平板作为基本翼型，并以其中点为原点建立与平板重合的直线坐标轴s，正面向上。

在奇点分布法中，要求设置在翼型骨线或弦上的连续分布涡产生与固态翼型形成相同的平面流场。现假设布置在弦线上的逆时针连续涡的旋涡密度函数r(s)有如下形式，这也是所有栅中翼型弦线上分布涡共同的旋涡强度分布规律。

$>\left(\begin{array}{c}r\left(s\right)=f\left(s\right)+A_1\sqrt{1-{(S/\frac{l}{2})}^2}+{2A}_2(S/\frac{l}{2})\sqrt{1-{(S/\frac{l}{2})}^2}+{3A}_3(4{(S/\frac{l}{2})}^2-1)\sqrt{1-{(S/\frac{l}{2})}^2}\\ +{4A}_4(S/\frac{l}{2})(2{(S/\frac{l}{2})}^2-1)\sqrt{1-{(S/\frac{l}{2})}^2}+......\end{array}\right)---\left(1\right)$>

令式（1）中可以将r(s)、f(s)这一s的函数写成以正弦函数为基的广义富里叶级数r(θ)：

$>r\left(\theta \right)=f\left(\theta \right)+\underset{k=1}{\overset{\infty }{\Sigma }}{A}_k\mathrm{Sink\theta }---\left(2\right)$>

式（1）或（2）给出的密度函数中含有无穷多项，事实上这是不必要的。为此，计算涡层产生的诱导速度绕流一翼型骨线的环量Γ_B。

应该指出，沿围绕一平板的闭曲线由诱导速度产生的环量只能由本平板上的分布连续涡生成，栅中其它翼型弦线上连续涡对此并无贡献，对此简要说明如下。

如图2所示，一弦线上一微段长ds，微段上一点的直线坐标为s，则微段上各点的旋涡密度为r(s)。在点涡为逆时针条件下，由微段上的连续涡产生微段两侧的切向诱导速度Vτ等值反向，方向见图2。诱导速度沿包围微段的闭曲线ABCDA的环量Γ_B显然为2Vτds。本翼型，及其它所有翼型上的点涡在ds两侧产生的诱导速度等值同向，沿闭曲线ABCDA的环量为0。

由斯托克斯定理，dΓ_B也等于r(s)ds，于是，一翼型骨线或弦线上连续涡的诱导速度绕一翼型的环量Γ_B为：

$>\left(\begin{array}{c}{\Gamma }_B={\int }_{-l12}^{l12}r\left(s\right)\mathrm{ds}={\int }_{\pi }^{0}r\left(\theta \right)(-l/2)\sin \mathrm{\theta d\theta }\\ ={\int }_0^{\pi }(f\left(\theta \right)+\underset{k=1}{\overset{\infty }{\Sigma }}{A}_k\sin \mathrm{k\theta })(l/2)\sin \mathrm{\theta d\theta }\end{array}\right)$>

利用-2sinkθsinθ＝cos(k+1)θ-cos(k-1)θ，易于证明上面定积分中含A₂、A₃、A₄......的各项均为0，这些项对环量Γ_B无贡献，因而密度函数Υ(θ)简化为：

Υ(θ)=f(θ)+A₁sinθ    （3）

目前，几乎所有的参考资料中都假定附加项f(θ)有如下形式：

$>f\left(\theta \right)=A_0\mathrm{ctg}\frac{\theta }{2}---\left(4\right)$>

或$>f\left(s\right)=A_0\sqrt{\frac{1+s/(l/2)}{1-s/(l/2)}}---\left(5\right)$>

经反复计算分析，式（4）、（5）定义的附加项能满足轴流式转轮流面翼型的边界条件，但不适合轴流式叶轮的设计。在后来的设计实践中，采用如下定义的附加项f(θ),f(s)：

$>f\left(\theta \right)=A_0\mathrm{tg}\frac{\theta }{2}---\left(6\right)$>

或$>f\left(s\right)=A_0\sqrt{\frac{1-s/(l/2)}{1+s/(l/2)}}---\left(7\right)$>

翼型骨线上的旋涡密度函数成为：

$>r\left(\theta \right)=A_0\mathrm{tg}\frac{\theta }{2}+A_1\mathrm{Sin\theta }---\left(8\right)$>

或$>r\left(s\right)=A_0\sqrt{\frac{1-s/(l/2)}{1+s/(l/2)}}+A_1\sqrt{1-{[s/(l/2)]}^2}---\left(9\right)$>

现利用式（8）、（9）计算绕流一个翼型的环量Γ_B：

绕翼型的速度环量是叶轮产生扬程的基础，它们之间有熟悉的关系：

$>{\Gamma }_B=\frac{2\pi }{\mathrm{\omega z}}·\frac{\mathrm{gH}}{\mathrm{\eta h}}---\left(11\right)$>

（11）式右端的叶轮叶片数Ζ，叶轮旋转角速度ω，扬程H及叶轮水力效率η_h都是给定或前期已确定的物理量。于是，在计算流面上Γ_Β是一确定量。式（10）表明，A₀、A₁两项对形成要求的翼型环量都有作用。从后面分析中可以看到，这一事实对形成合理的翼型骨线有重要意义。

现假设一平面流面上均匀流大小为W∞，与平板翼型夹角即冲角为α，如图3所示。对运动叶栅，均匀流实际是指叶栅前后相对速度矢量的几何平均值，其大小W∞及其与叶栅列线夹角β∞都是可计算的已知量，冲角α则由设计人员给出。

下面将证明，在Α₀＝2αW∞的条件下，本文定义的旋涡密度公式（8）、（9）能获取较为理想的翼型骨线形态。为此，首先计算骨线或弦上的分布连续涡在骨线或弦上一点S₀处产生的诱导速度V_y。

在图3中，S、S₀两点直线坐标分别为s、s₀，这里S₀点位置视为固定的，因而s₀在计算中为一常数，S点位置则在区间[-l/2，l/2]内变化，连续涡在S₀处产生的诱导速度V_y为：

$>\left(\begin{array}{c}{V}_y=\frac{1}{2\pi }{\int }_{-l/2}^{l/2}\frac{r\left(s\right)\mathrm{ds}}{s_0-s}\\ =\frac{1}{2\pi }{\int }_{-l/2}^{l/2}\left\{\right[A_0\sqrt{\frac{1-s/(l/2)}{1+s/(l/2)}}+A_1\sqrt{1-{[s/(l/2)]}^2}]/(s_0-s)\}\mathrm{ds}\end{array}\right)$>

上面被积函数中，s₀＞s，因而被积函数的分母及分式值都正的。在图3中，因为S处点涡为逆时针方向，S₀处的诱导速度方向将垂直于数轴s且指向左上。因此，上面定积分值为正时，S₀处的诱导速度方向应如图3所示。如果把诱导速度视为矢量，那么y轴正向应指向左上，这样才能保证速度矢量在y轴投影为正值。y轴与数轴的交点是任意的，为方便，现置于区间（-l/2,l/2)的中点。

除令s/(l/2)＝cosθ外，再令s₀/(l/2)＝cosθ₀，在计算过程中，s₀和θ₀都视为常量。于是有：

$>\left(\begin{array}{c}{V}_Y=\frac{1}{2\pi }{\int }_{\pi }^0\frac{A_0\mathrm{tg}\frac{\theta }{2}+A_1\sin \theta }{\frac{l}{2}(\cos {\theta }_0-\cos )}·(-\frac{l}{2})\sin \mathrm{\theta d\theta }\\ =\frac{1}{2\pi }{\int }_0^{\pi }\frac{A_0(1-\cos )+A_1(1-\cos 2\theta )/2}{\cos {\theta }_0-\cos }\mathrm{d\theta }\end{array}\right)$>

上式出现了两个收敛的广义积分，利用格劳尔公式：

$>{\int }_0^{\pi }\frac{\cos \mathrm{k\theta d\theta }}{\cos {\theta }_0-\cos \theta }=-\pi \frac{\sin k{\theta }_0}{\sin {\theta }_0}$>

并将1视为cos0·θ，可以得到：

$>V_y=\frac{1}{2\pi }(\pi A_0+\frac{A_1}{2}·2\pi \cos {\theta }_0)$>

$>=\frac{1}{2}(A_O+\frac{A_1{s}_0}{l/2})$>

由于S₀在区间[-l/2,l/2]内的位置事实上是任意的，因而有：

$>V_y=\frac{1}{2}(A_0+\frac{A_{1}s}{l/2})---\left(12\right)$>

骨线或弦上任一点处的速度应由骨线上涡层在此产生的诱导速度和均匀流速度的叠加。由图3，均匀流在y轴，s轴的投影分别为-W∞sinα≈-W∞α，W∞cosα，于是，骨线上各点处的合成速度在s轴，y轴两个投影W_s，W_y为：

$>W_y=-W\infty \alpha +v_y=-W\infty \alpha +\frac{A_0}{2}+\frac{A_{1}s}{l}---\left(13\right)$>

W_s＝W∞cosα    （14）

分布在弦上的涡层事实上应分布在涡层代替的固态骨线曲线上，设骨线方程为y＝f(s)，骨线上各点处的斜率为由于固态骨线的不可穿越性，骨线上各点处的合成速度必须与骨线相切。在骨线微弯的条件下，骨线上各点合成速度在y轴，s轴的投影也由式（13）、（14）表达，合成速度与s轴夹角正切为W_y/W_S，于是得到：

$>\frac{\mathrm{dy}}{\mathrm{ds}}=\frac{W_y}{W_s}=\frac{-W\infty \alpha +\frac{A_0}{2}+\frac{A_{1}s}{l}}{W\infty \cos \alpha }$>

积分上式，得到：

$>y=\frac{A_1}{\mathrm{lW}\infty \cos \alpha }·\frac{s^2}{2}+(-\frac{\alpha }{\cos \alpha }+\frac{A_0}{2W\infty \cos \alpha })s+C---\left(15\right)$>

积分常数C可由骨线应满足的边界条件确定：当时，y=0，如果A₀＝2αW∞，骨线方程（15）成为：

$>y=\frac{A_1}{\mathrm{lW}\infty \cos \alpha }·\frac{s^2}{2}+C$>

积分常数$>C=-\frac{A_1}{\mathrm{lW}\infty \cos \alpha }·\frac{1}{2}·{\left(\frac{l}{2}\right)}^2$>

骨线方程成为：

$>y=\frac{A_1}{\mathrm{lW}\infty \cos \alpha }·\frac{1}{2}(s^2-{\left(\frac{l}{2}\right)}^2)$>

这是一条对称于y轴的抛物线，抛物线的最高点对应骨线的拱度f，应等于s＝0时的|y|值，于是有：

$>f=\frac{A_1}{\mathrm{lW}\infty \cos \alpha }·\frac{1}{2}·{\left(\frac{l}{2}\right)}^2---\left(16\right)$>

如果不设置附加项，即令A₀＝0，或者A₀点虽然不为0，但不等于2αW∞，都不存在能满足两端点边界条件的积分常数C值，即不存在能产生给定流场的骨线。

式（10）表明，在绕流翼型的环量Γ_B一定时，由于A₀的存在减小了A₁值，由式（14）可看出，这时减小了f值。实际水流都是有粘性的，较小的f值有利于减小翼型的阻力系数，也减小了翼型表面脱流的可能。事实上，所有优秀翼型都有拱度不大的微弯持点。

上面的定量计算表明，本申请提出的附加项式（6）、（7）式能够形成几何形态和水力性能均令人满意的翼型骨线。类似的分析计算表明，适合转轮翼型骨线设计的传统的附加项式（4）、（5）也有这样的功能，但两种附加项有不同的函数形式、应用范围和计算结果。

本申请定义的旋涡密度表达式$>2\mathrm{\alpha W}\infty \sqrt{\frac{1-s/(l/2)}{1+s/(l/2)}}+A_1\sqrt{1-{[s/(l/2)]}^2},$>能保证置于均匀流的翼型满足应有的边界条件，包括：

1）骨线表面各点由均匀流与本骨线连续旋涡在骨线上的诱导速度叠加后的合成速度与骨线相切，无骨线的法向分量，这是由骨线的不可穿透性决定的。

2）满足库塔-恰布雷金定理，即均匀流绕流带尖锐前、后缘物体时，前、后缘分别是奇点和驻点，或速度是有限值。

骨线本身就是在骨线上各点切线斜率等于合成流动方向条件下推出的，边界条件1）自然满足。

由式（13）、（14），翼型骨线前缘、后缘两点s坐标为-l/2，l/2，在A₀＝2αW∞时，前缘两个速度分量$>W_s=W\infty \cos \alpha ,W_y=-\frac{A_1}{2},$>后缘两个速度分量$>W_s=W\infty \cos \alpha ,W_y=\frac{A_1}{2},$>因而前后缘的速度都是有限值。现考查由本申请定义的旋涡密度函数在骨线表面诱导的切向速度V_τ的分布，由2V_τds＝r(s)ds，有：

$>V_{\tau }=\pm \frac{1}{2}(A_0\sqrt{\frac{1-s/(l/2)}{1+s/(l/2)}}+A_1\sqrt{1-{[s/(l/2)]}^2})$>

代入骨线前、后缘点的s生标-l/2，l/2，得到前、后缘点切向速度分别为无穷大与0，叠加有限值的速度值后，前后缘速度分别成为无穷大与有限值，满足库塔-恰布雷金定理。

传统的旋涡密度也能产生合理的翼型骨线，在转轮叶栅中，均匀流基本从翼型的s轴的正向流向负向，奇点与驻点与叶轮翼型是颠倒的，如果将这种旋涡密度借用到叶轮翼型，则不能满足库塔-恰布雷金假定。由此可见，由于展开流面上叶栅处于不同的均匀流场中，因而不存在统一的边界条件，也不存在同时能适应不同的边界条件的统一旋涡密度函数，后继导出的若干公式也不能互相借用。

以奇点分布法确定各计算流面上翼型骨线，实际上是要确定骨线上若干节点的坐标，通过这些节点的插值曲线即为翼型骨线。由于这一方法理论深刻，计算步骤繁杂，计算量巨大，以人工计算几乎不可能实现任何一次设计任务。这是这一方法难以为企业接受，很难进入设计实践领域的重要原因。以算法语言完成计算过程，快速，准确地获取各翼型骨线上若干个节点位置的解析解是将先进理论转化为现实生产力，提高企业自主创新能力的需要。

设在平面上有一栅距为t，弦长为l的无限平面直列平板叶栅，现任选其中一个翼型为基本翼型，其它翼型统称其它翼型。

建立两个坐标系，除前文所述的直线坐标s轴外，另一个平面直角坐标系u-o-z的原点也在基本翼型平板的中点，坐标系的横轴u和纵轴z分别平行直列叶栅的列线和栅轴，如图4所示。

设想在栅中所有翼型上各取一点，通过它们的连线平行于u轴，其中一个点位于基本翼型骨线上，其直角坐标为（u，z），其余各点的直角坐标显然为u+nt（n=±1，±2，±3，……不取0），z。在所有点上设置一强度同为Γ的逆时针点涡，所有这些点涡在基本翼型上直线坐标为s₀，直角坐标为(u₀,z₀)的给定点S₀处产生的诱导速度的两个直角坐标分量V_u,V_z，由平面势流理论应为：

$>V_u=-\frac{\Gamma }{2t}·\frac{\mathrm{sh}\left[\frac{2\pi }{t}(z_0-z)\right]}{\mathrm{ch}\left[\frac{2\pi }{t}(z_0-z)\right]-\cos \left[\frac{2\pi }{t}(u_0-u)\right]}$>

$>V_Z=\frac{\Gamma }{2t}·\frac{\sin \left[\frac{2\pi }{t}(u_0-u)\right]}{\mathrm{ch}\left[\frac{2\pi }{t}(z_0-z)\right]-\cos \left[\frac{2\pi }{t}(u_0-u)\right]}$>

事实上各翼型上分布同一规律的，逆时针方向的连续涡，旋涡密度是点的s坐标的函数r(s)。前文已说明，传统的旋涡密度函数不适应轴流泵叶栅的翼型边界条件。本申请中的旋涡密度r(s)有如下的创新形式：

$>r\left(s\right)=A_0\sqrt{\frac{1-s/(l/2)}{1+s/(l/2)}}+A_1\sqrt{1-{[s/(l/2)]}^2}$>

以这一旋涡密度函数代替上式中的Γ，得到各翼型上直线坐标为s，直角坐标为u+nt＝（n=0，±1，±2，……），z的所有逆时针点涡在基本翼型S₀(u₀,z₀)处产生的诱导速度的两个分量dvz和dvz：

$>\mathrm{dvu}=\frac{-(A_0\sqrt{\frac{1-s/(l/2)}{1+s/(l/2)}}+A_1\sqrt{1-{(s/(l/2))}^2})}{2t}\frac{\mathrm{sh}\left[\frac{2\pi }{t}(z_0-z)\right]}{\mathrm{ch}\left[\frac{2\pi }{t}(z_0-z)\right]-\cos \left[\frac{2\pi }{t}(u_0-u)\right]}$>

$>\mathrm{dvz}=\frac{(A_0\sqrt{\frac{1-s/(l/2)}{1+s/(l/2)}}+A_1\sqrt{1-{[s/(l/2)]}^2})}{2t}\frac{\sin \left[\frac{2\pi }{t}(u_0-u)\right]}{\mathrm{ch}\left[\frac{2\pi }{t}(z_0-z)\right]-\cos \left[\frac{2\pi }{t}(u_0-u)\right]}---\left(14\right)$>

在区间[-l/2,l/2]积分式（14），最终得到栅中所有翼型的连续涡在基本翼型S₀(u₀,z₀)处产生的诱导速度的两个分量V_u，V_Z：

$>V_u=\frac{{\int }_{-l/2}^{l/2}-(A_0\sqrt{\frac{1-s/(l/2)}{1+s/(l/2)}}+A_1\sqrt{1-{(s/(l/2))}^2})}{2t}\frac{\mathrm{sh}\left[\frac{2\pi }{t}(z_0-z)\right]}{\mathrm{ch}\left[\frac{2\pi }{t}(z_0-z)\right]-\cos \left[\frac{2\pi }{t}(u_0-u)\right]}\mathrm{ds}$>

$>V_Z=\frac{{\int }_{-l/2}^{l/2}(A_0\sqrt{\frac{1-s/(l/2)}{1+s/(l/2)}}+A_1\sqrt{1-{(s/(l/2))}^2})}{2t}\frac{\sin \left[\frac{2\pi }{t}(u_o-u)\right]}{\mathrm{ch}\left[\frac{2\pi }{t}(z_0-z)\right]-\cos \left[\frac{2\pi }{t}(u_0-u)\right]}\mathrm{ds}$>

（15）

式（15）给出了计算基本翼型骨线上若干个给定计算点处的诱导速度的两个分量的表达式。

基本翼型上S₀(u₀,z₀)是一给定点，但基本翼型上点涡坐标（s，z）则是积分变量，当u＝u₀,z＝z₀时，两被积函数均为0/0，这一不确定性使积分不能进行。

为此，将基本翼型上计算点S₀(u₀,z₀)处的诱导速度分解为两个部分：

基本翼型骨线上的连续涡在S₀点产生的诱导速度V_1u，V_1z，除开基本翼型的其余所有翼型的连续涡在S₀处产生的诱导速度V_2u，V_2z，从而有：

V_u＝V_1u+V_1z

V_z＝V_2u+V_2z

设平面无穷直列叶栅所在的平面均匀流的大小为W∞，流速矢量与叶栅列线夹角为β∞，那么均匀流的两个直角坐标分量W∞u和W∞z以式（16）计算，在同一流面各点处是两个常量：

$>\left(\begin{array}{c}W\infty u=W\infty \cos \beta \infty \\ W\infty z=W\infty \sin \beta \infty \end{array}\right)---\left(16\right)$>

于是，基本翼型S₀处合成的相对速度的两个分量为：

$>\left(\begin{array}{c}{W}_u=W\infty u+V_u=W\infty u+V_{1u}+V_{2u}\\ W_z=W\infty z+V_z=W\infty z+V_{1z}+V_{2z}\end{array}\right)---\left(17\right)$>

在一展开为平面的计算流面上，旋涡密度r(s)中的两个系数不能任意给定。对平面孤立翼型，前文已导出A₀＝2αW∞，这里冲角α指平面上均匀流与翼型骨线的夹角，其值应不大于1°。程序中，将轮毂，轮缘处的α值分别确定为1°和0.25°，中间流面上这一值接线性变化规律确定。

栅中翼型的动力特性与弧点翼型有区别，对平板直列叶栅，反映它们不同的A_O值的系数l_n已由实验完成，且将l_n值与叶栅稠密度lt和翼型安放角β_e的函数关系以若干条曲线表示。程序中以拉格朗日插值方法求解了这些曲线的回归方程，共10个，供计算中调用。由此得到栅中翼型A_O的最终表达式：A₀＝2αW∞ln。

各计算流面上的翼型安放角βe＝β∞+α，已计算了各流流面的β∞和α后，也要事先算出。

前文式（10）中已导出了A₀,A₁与Γ_B的关系。由式（11），Γ_B可以由事先确定的叶轮叶片数z，叶轮旋转角速式ω，扬程H及叶轮水力效率η_h计算，于是可以由式（10）获得各计算流面的A₁值。这样本申请定义的旋涡密度函数式（9）在各计算流面上都有了确定的形式。

下面讨论V_1u和V_1z的计算：

基本翼型上一给定点S₀（直线坐标为s₀，直角坐标为(u₀,z₀)）上由本翼型的逆时针连续涡产生的诱导速度与翼型骨线垂直，其大小V_y已在前文中导出：由于连续涡为逆时针，V_y的两个直角坐标分量为：

$>\left(\begin{array}{c}{V}_{1u}=-V_y\sin \mathrm{\beta e}=-(\frac{A_0}{2}+A_1{s}_0/l)\sin \mathrm{\beta e}\\ V_{1z}=V_y\cos \mathrm{\beta e}=(\frac{A_0}{2}+A_1{s}_0/l)\cos \mathrm{\beta e}\end{array}\right)---\left(18\right)$>

V_1u和V_1z除开可用点的直线坐标s₀表达外，还可以用直角坐标表示。

基本翼型骨线上一点(u,z)处的逆时针旋涡强度为r(s)，它在计算点(u₀,z₀)处产生的诱导速度大小V_y为如图5。V_y的两个分量dv_1u和dv_1z分别为：

$>\left(\begin{array}{c}{\mathrm{dV}}_{1u}=-V_y\sin \mathrm{\beta e}=-\frac{1}{2\pi }·\frac{r\left(s\right)\mathrm{ds}}{s_0-s}·\sin \mathrm{\beta e}\\ =-\frac{1}{2\pi }·\frac{r\left(s\right)\mathrm{ds}}{s_0-s}·\frac{z_0-z}{s_0-s}=-\frac{1}{2\pi }·r\left(s\right)\mathrm{ds}·\frac{z_0-z}{{(z_0-z)}^2+{(u_0-u)}^2}\end{array}\right)$>

$>\left(\begin{array}{c}{\mathrm{dV}}_{1z}=V_y\cos \mathrm{\beta e}=\frac{1}{2\pi }·\frac{r\left(s\right)\mathrm{ds}}{s_0-s}·\cos \mathrm{\beta e}\\ =\frac{1}{2\pi }·\frac{r\left(s\right)\mathrm{ds}}{s_0-s}·\frac{u_0-u}{s_0-s}=\frac{1}{2\pi }·r\left(s\right)\mathrm{ds}·\frac{u_0-u}{{(z_0-z)}^2+{(u_0-u)}^2}\end{array}\right)$>

于是有：

$>\left(\begin{array}{c}{V}_{1u}={\int }_{-l/2}^{l/2}-\frac{1}{2\pi }·\frac{r\left(s\right)\mathrm{ds}(z_0-z)}{{(z_0-z)}^2+{(u_0-u)}^2}\\ V_{1Z}={\int }_{-l/2}^{l/2}\frac{1}{2\pi }·\frac{r\left(s\right)\mathrm{ds}(u_0-u)}{{(z_0-z)}^2+{(u_0-u)}^2}\end{array}\right)---\left(19\right)$>

（19）式中的r(s)如前所述，指本申请定义的旋涡密度表达式，即式（9）。

下面讨论V_2u和V_2Z的计算：

由于V_2u和V_2Z分别指基本翼型点S₀(u₀,z₀)处由除开基本翼型其它所有翼型上的连续涡在此产生的诱导速度的两个分量，由这一定义，有：

$>\left(\begin{array}{c}{V}_{2u}=V_u-V_{1u}={\int }_{-l/2}^{l/2}\frac{r\left(s\right)\mathrm{ds}}{t}(-\frac{1}{2}\frac{\mathrm{sh}\left(\frac{2\pi }{t}(z_0-z)\right)}{\mathrm{ch}\left(\frac{2\pi }{t}(z_0-z)\right)-\cos \left(\frac{2\pi }{t}(u_0-u)\right)}+\frac{t}{2\pi }·\frac{z_0-z}{{(z_0-z)}^2+{(u_0-u)}^2})\\ V_{2Z}=V_z-V_{1z}={\int }_{-l/2}^{l/2}\frac{r\left(s\right)\mathrm{ds}}{t}\left(\frac{1}{2}\frac{\sin \left(\frac{2\pi }{t}(u_0-u)\right)}{\mathrm{ch}\left(\frac{2\pi }{t}(z_0-z)\right)-\cos \left(\frac{2\pi }{t}(u_0-u)\right)}\right)-\frac{t}{2\pi }·\frac{u_0-u}{{(z_0-z)}^2+{(u_0-u)}^2}\end{array}\right)---\left(20\right)$>

为方便，令：

$>a(s_0,s)=-\frac{1}{2}·\frac{\mathrm{sh}\left[\frac{2\pi }{t}(z_0-z)\right]}{\mathrm{ch}\left(\frac{2\pi }{t}(z_0-z)\right)-\cos \left(\frac{2\pi }{t}(u_0-u)\right)}+\frac{t}{2\pi }·\frac{z_0-z}{{(z_0-z)}^2+{(u_0-u)}^2}$>

$>b(s_o-s)=\frac{1}{2}·\frac{\sin \left[\frac{2\pi }{t}(u_0-u)\right]}{\mathrm{ch}\left(\frac{2\pi }{t}(z_0-z)\right)-\cos \left(\frac{2\pi }{t}(u_0-u)\right)}-\frac{t}{2\pi }·\frac{u_0-u}{{(z_0-z)}^2+{(u_0-u)}^2}---\left(21\right)$>

a(s₀,s)和b(s_o,s)均无量纲，其物理意义是：除开基本翼型，其它翼型上所有位于（u+nt(n＝±1,±2,±3,......不取0)，z），强度为1，相距为t的逆时针点涡在基本翼型上(u₀,z₀)处产生的诱导速度乘以栅距t后的水平和铅垂分量，式中的(u,z)是基本翼型上一点涡的坐标。

在定积分式（20）中，当u、z趋近于u₀、z₀时，式中的4项都成为0/0，但积奇函数的差值有确定值0。对此证明如下：

将e^x,e^-x展开为麦克劳林级数，有：

e^x＝1+x+x²/2!+x³/3!+......

e^-x＝1-x+x²/2!-x³/3!+......

从而：

chx＝(e^x+e^-x)/2＝1+x²/2!(在x为微量时)

shx＝(e^x+e^-x)/2＝x(在x为微量时)

将sinx,cosx展开为麦克劳林级数，有：

sinx＝x-x³/3!+x⁵/5!-x⁷/7!+......＝x(在x为微量时)

cosx＝1-x²/2!+x⁴/4!-x⁶/6!+......＝1-x²/2!(在x为微量时)

当u、z充分接近u₀、z₀时，都是微量，因此：

$>\left(\begin{array}{c}\mathrm{ch}\left[\frac{2\pi }{t}(z_0-z)\right]-\cos \left[\frac{2\pi }{t}(u_0-u)\right]\\ =1+{\left[\frac{2\pi }{t}(z_0-z)\right]}^2/2!-\{1-{\left[\frac{2\pi }{t}(u_0-u)\right]}^2\}/2!\\ =\frac{1}{2}{\left(\frac{2\pi }{t}\right)}^2{(z_0-z)}^2+\frac{1}{2}{\left(\frac{2\pi }{t}\right)}^2{(u_0-u)}^2\end{array}\right)$>

同时：

$>\mathrm{sh}\left[\frac{2\pi }{t}(z_0-z)\right]=\frac{2\pi }{t}(z_0-z)$>

这样，在u、z充分接近u₀、z₀时，

$>a(s_0,s)=-\frac{1}{2}·\frac{2\pi (z_0-z)/t}{\frac{1}{2}{\left(\frac{2\pi }{t}\right)}^2{(z_0-z)}^2+\frac{1}{2}{\left(\frac{2\pi }{t}\right)}^2{(u_0-u)}^2}+\frac{t}{2\pi }·\frac{z_0-z}{{(z_0-z)}^2+{(u_0-u)}^2}=0$>

同样可以证明，在相同的条件下，b(s₀,s)＝0。这表明，式（20）中的两个被积函数是没有间断点的连续函数，积分可进行。

把由式（21）定义的a(s_o,s)代入式（20），得到V_2u的表达式：

$>V_{2u}=\frac{1}{t}{\int }_{-l/2}^{l/2}a(s_0·s)\{A_0\sqrt{\frac{1-s/(l/2)}{1+s/(l/2)}}+A_1\sqrt{1-{(s/(l/2))}^2}\}\mathrm{ds}---\left(22\right)$>

为积分方便，作换元s/(l/2)＝x,ds＝(l/2)dx,代入（22），得到：

$>\left(\begin{array}{c}{V}_{2u}=\frac{1}{t}{\int }_{-1}^{1}a(x_0,x)(A_0\sqrt{\frac{1-x}{1+x}}+A_1\sqrt{1-x^2})\frac{l}{2}\mathrm{dx}\\ =\frac{l}{2t}{\int }_{-1}^{1}a(x_0,x)(A_0\sqrt{\frac{1-x}{1+x}}·\sqrt{\frac{1-x}{1-x}}+A_1\sqrt{1-x^2}·\frac{\sqrt{1-x^2}}{\sqrt{1-x^2}})\mathrm{dx}\\ =\frac{l}{2t}{\int }_{-1}^{1}a(x_0,x)\left[\frac{A_0(1-x)+A_1(1-x^2)}{\sqrt{1-x^2}}\right]\mathrm{dx}\end{array}\right)---\left(23\right)$>

式（23）代表满足x＝s/(l/2)的s值所在点涡，对满足x₀＝s₀/(l/2)的s_ο值所在点涡的作用，s，s_ο值角坐标分别为(u,z)，(u₀,z₀)，这一作用值即式（21）第一式。

为求解定积分式（23），要用到涅伯米宁公式：

$>{\int }_{-1}^1\frac{f\left(x\right)\mathrm{dx}}{\sqrt{1-x^2}}=\frac{\pi }{1280}[167f(-1)+378f(-\frac{2}{3})-135f(-\frac{1}{3})+460f\left(0\right)-135\left(f\right)\left(\frac{1}{3}\right)+378f\left(\frac{2}{3}\right)+167f\left(1\right)]$>

现式（23）中f(x)＝a(x₀,x)[A₀(1-x)+A₁(1-x²)]，因而

$>167f(-1)=a(x_0,-\frac{3l}{6})\left({2A}_0\right)167=a(x_0,-\frac{3l}{6})·334A_o$>

$>\left(\begin{array}{c}378f(-\frac{2}{3})=a(x_0,-\frac{2l}{6})(\frac{5}{3}{A}_0+\frac{5}{9}{A}_1)378=a(x_0,-\frac{2l}{6})(630A_0+210A_1)\\ -135f(-\frac{1}{3})=a(x_0,-\frac{l}{6})(\frac{4}{3}{A}_0+\frac{8}{9}{A}_1)(-135)=a(x_0,-\frac{l}{6})(-180A_0-120A_1)\end{array}\right)$>

$>\left(\begin{array}{c}460f\left(0\right)=a(x_0,0)(A_0+A_1)460=a(x_0,0)(460A_0+460A_1)\\ -135f\left(\frac{1}{3}\right)=a(x_0,\frac{l}{6})(\frac{2}{3}{A}_0+\frac{8}{9}{A}_1)(-135)=a(x_0,\frac{l}{6})(-90A_0-120A_1)\end{array}\right)$>

$>378f\left(\frac{2}{3}\right)=a(x_0,\frac{2l}{6})(\frac{1}{3}{A}_0+\frac{5}{9}{A}_1)378=a(x_0,\frac{2l}{6})(126A_0+210A_1)$>

$>167f\left(1\right)=a(x_0,\frac{3l}{6})(0+0)=0$>

将以上结果代回式（23），得到计算V_2U的表达式：

$>V_{2U}=\frac{1}{2560}·\frac{\mathrm{\pi l}}{t}\left(\begin{array}{c}334a(x_0,-\frac{3l}{6})A_0+630a(x_0,-\frac{2l}{6})A_0+210a(x_0,-\frac{2l}{6})A_1\\ -180a(x_0,-\frac{l}{6})A_0-120a(x_0,-\frac{l}{6})A_1+460a(x_0,0)A_0\\ +460a(x_0,0)A_1-90a(x_0,\frac{l}{6})A_0-120a(x_0,\frac{l}{6})A_0\\ +126a(x_0,\frac{2l}{6})A_0+210a(x_0,\frac{2l}{6})A_1\end{array}\right)---\left(24\right)$>

式（24）中，表示基本翼型骨线上直线坐标为直角坐标为(u,z)的逆时针点涡（在s等于确定值l/6的解条件下，这点的直角坐标(u,z)是确定易于计算），在s₀＝(l/2)x₀处产生的对应值，由于s₀处的直角坐标(u₀,z₀)也易于计算，因而a(s₀,l/6)可由式（21）第一式获得。

其余a值可以同样说明。

当(u,z)等于(u₀,z₀)时，a值为0。

在计算式（20）的V_2z时，只要把系数a换成对应b值即可。

将展开为平面的流面上翼型骨线等分成6等份，得到含首尾在内的7个分点。利用式（24）分别计算各分点V_2U值，同时计算7个V_2z的值。

到此，7个分点上V_1u,V_1z,V_2u,V_2z,W∞u,W∞z都已得到，于是可以计算7个分点上合成的相对速度矢量与列线的夹角β₀，β₁，β₂，......β₆，

$>\beta =\arctan \frac{W\infty z+V_{1z}+V_{2z}}{W\infty u+V_{1u}+V_{2u}}---\left(25\right)$>

过起点作线段l₀₁＝l/6与列线夹角等于(β₀+β₁)/2得1点，过1点作线段l₁₂＝l/6与列线夹角等于(β₁+β₂)/2得到2点，等等，直到确定平面点6的坐标。通过0，1，2，…..，6点的样条曲线即为这一流面上的翼型骨线。

大量实验证明，决定翼型水力性能的是翼型骨线形态。在获得各流面上骨线后加厚叶片，得到完整翼型。

对应用算法语言程序的说明：

可以看出，奇点分布法计算工作量巨大，难以人工计算得到结果。以QB高级语言编写了实现全部计算过程的计算软件，使这一方法应用于设计实践成为可能。

第一部分：计算总体及各流面的基本几何和流动参数。

确定计算流面个数，流面个数可确定为任意值，通常流面个数取5个。

循环计算各流面半径，计算根据最新资料推荐的各流面叶栅稠密度l/t，栅距t及弦长l。

根据用户给定的泵在设计点的流量，扬程，转速等性能参数，循环计算各计算流面上均匀流大小W∞和流动方向角β∞，进一步计算各流面上均匀流的水平及铅垂分量W∞u,W∞z。

确定各流面的冲角α后，循环计算各流面上翼型安放角β_e。

循环计算各流面上绕单个翼型的环量Γ_B。

以插值方法获取反映孤立翼型和栅中翼型性能差异的系数l_n的插值方程，共10个。根据各流面上已知的翼型安放角β_e和叶栅稠密度l/t确定流面的A₀值，并由此及Γ_B值循环计算各流面的系数A₁。

到此，各流面上基本翼型几何参数和基本流动参数已确定，这是确定各翼型骨线上各节点位置的基础。

第二部分：各流面基本翼型骨线节点上各分速度的计算及骨线节点位置的计算。

将一流面的长l翼型骨线等分成6段，得到7个分点，它们的直线坐标分别为-3l/6,-2l/6,......,3l/6,由于各流面上翼型弦长及翼型安放角β_e在第一部分计算，各流面的系数A₀和A₁也在第一部分计算，由式（18）即可分别求出7个分点的V_1u,V_1z。

各流面上均匀流的两个速度分量V∞u,V∞z也在第一部分计算，它们在同一流面各点是两个常数。

在计算流面翼型骨线各分点上的V_2u,V_2z时，必须计算各分点S₀上的a(s_o,-3l/6),a(s_o,-2l/6),...a.(.s_o,,3l/6)值及b值。

将一流面上的7个分点从下向上编号，编号分别为0，1，……6。在程序中形成一个7×7的方矩数组a(7,7)。数组元素a(m,k)(m,k＝0,1,2,......6)表示编号为m的点涡在编号为k的计算点产生的诱导速度的水平分量。

这是一个关于主对角线反对称的方阵，且平行于主对线的直线上的元素相等。因为该方阵关于主对角线反对称，因此可以简化方阵元素计算工作量，现说明如下。

正弦函数和双曲正弦函数都是奇函数，余弦函数和双曲余弦函数都是偶函数。在a(2,4)的表达式中，z₀-z＝(-l/6-l/6)sinβ_e,u₀-u＝(-l/6-l/6)cosβ_e，在a(4,2)的表达式中，z₀-z＝(l/6-(-l/6))sinβ_e,u₀-u＝(l/6-(-l/6))cosβ_e，由函数的奇偶性，可以看出a(2,4)＝-a(4,2)，其它每对关于主对角线对称的元素也有这一特性。

首先计算主对角线右上方的全部元素。以u1,z1表示流面骨线上任意两个相邻计算点的横坐标与纵坐标之差，在同一流面上这是两个常数。在变量p一定时（p=0，1，2，……6），q从p+1到6循环，计算p，q代表的两点的横、纵坐标之差2π(u₀-u)/t,2π(z₀-z)/t，进一步计算sh[2π(z₀-z)/t],ch[2π(z₀-z)/t],cos[2π(u₀-u)/t]将这些值代入式（21）即可获得方阵右上角全部a(s₀.s)值。

上面计算简化了计算过程：并不需要确定两个点的横纵坐标的实际值，只需求出它们的差值即可。

由方阵元素的负对称性，形成主对角线元素左下方全部元素。

主对线上7个元素全部置0。

在利用式（24）计算一个流面骨线上的V_2u时，只需调用方阵元素即可得到a(x₀,-2l/6)等等值。

以同样的原则和方法形成b方程的全部元素。

在确定了骨线上7个分点的全部速度分量后，以式（25）确定7个分点上的β角，并可利用前面所述方法得到最终骨线上7个点位置的坐标。

轴流泵有大功率大尺寸的特点，在翼型骨线较长时，只由7个点难以确定骨线曲线的准确形态。为此，程序添加了一部分内容。以5段相切圆弧通过7个已确定点，其中第一段圆弧由7个分点的前3个决定。然后在圆弧上插入8个点，程序最终将输出骨线曲线上15个点而不仅是7个点的坐标。

按上述原则循环完成所有流面上的骨线运算，程序结束。

在各流面上据已有翼型骨线加厚叶片，加厚时，以型线为工作面向背面加厚。翼型厚度变化规律如下表所示：

x/l00.050.0750.10.20.30.40.50.60.70.80.90.951.0S/Smax00.2960.4050.4890.7780.920.9781.00.8830.7650.5440.3560.20

表中S指翼型厚度，Smax指翼型给定的翼型最大厚度。