1 叶轮设计计算
1.11 流线分点
如上图所示,分点的实质就是在流面上画特征线,组成扇形网格
Δ
s
=
Δ
u
\Delta s=\Delta u
Δs=Δu。因为流面是轴面流线绕中心线旋转一周所得到的曲面,一个流面上的全部轴面流线均相同,所以只要分相应的一条轴面流线,就等于在整个流面上绘制出了方格网。
其实,说白了,就是建立起一个空间曲面坐标系而已,以方便后续在方格网展开上进行保角变换来计算叶片各处的真实扭角。
书上给出了流线分点方法有如下几种:
(1) 逐点计算法
Δ θ \Delta\theta Δθ一般取 5 ∼ 10 5\sim10 5∼10度,从叶轮出口沿轴面流线任取 Δ s \Delta s Δs,量出 Δ s \Delta s Δs段中点的半径 R R R,然后计算其所对应的周向弧长 Δ u = Δ θ R \Delta u=\Delta \theta R Δu=ΔθR(注意这里的 Δ θ \Delta\theta Δθ应转化为弧度),若 Δ s \Delta s Δs与 Δ u \Delta u Δu相等,则分点正确。若不相等,则依据实际情况寻找新的 Δ s \Delta s Δs,再算 Δ u \Delta u Δu,直至两者相等为止。然后,从分得的点起继续向下分 2 , 3 , . . . 2,3,... 2,3,...点。
说起来容易做起来难,毋庸置疑,该方法既繁琐耗时,又容易产生累积误差,相当不精确。
(2) 作图分点法
在轴面投影图旁边,画两条夹角等于
Δ
θ
\Delta\theta
Δθ的射线,这两条射线表示夹角为
Δ
θ
\Delta\theta
Δθ的两轴面。与逐点计算分点法相同,一般取
Δ
θ
=
5
∼
10
\Delta\theta=5\sim10
Δθ=5∼10度,从出口开始,沿轴面流线尝试取长度
Δ
s
\Delta s
Δs,若
Δ
s
\Delta s
Δs的中点半径对应的两射线间的弧长
Δ
u
\Delta u
Δu,与试取的
Δ
s
\Delta s
Δs相等,则分点正确。否则另取
Δ
s
\Delta s
Δs,直到
Δ
s
=
Δ
u
\Delta s=\Delta u
Δs=Δu为止。第一个点确定以后,采用同样方法分第
2
,
3
,
.
.
.
2,3,...
2,3,...点。当流线平行于轴线时,
Δ
u
\Delta u
Δu不变,用对应的
Δ
s
\Delta s
Δs截取等长便可。特别注意的是,
Δ
s
\Delta s
Δs是轴面流线的长度,而不是射线间的垂直长度(只是在轴面流线垂直轴线时二者才相等)。由图可见,分完一条轴面流线,就等于在整个流面上画满了方格网。
本质上,方法(1)和方法(2)的是相同的,只不过一个要计算弧长,一个是直接画图量取的而已。那么方法(2)同样存在繁琐耗时和误差较大的缺陷。
事实上,自己写个小程序来实现上述分点过程就妥妥的了,既方便快捷又精度很高。在写程序之前,先看下所谓的方法3。
(3) 图解积分法
由于
Δ
u
=
Δ
θ
R
\Delta u=\Delta \theta R
Δu=ΔθR,故
Δ
θ
=
Δ
u
/
R
\Delta \theta = \Delta u / R
Δθ=Δu/R,又因为
Δ
u
=
Δ
s
\Delta u=\Delta s
Δu=Δs,则
θ
=
∫
0
θ
Δ
θ
=
∫
0
s
(
d
s
/
R
)
≈
∑
i
(
Δ
s
i
/
R
i
)
\theta = \int_0^\theta\Delta\theta=\int_0^s(ds/R)\approx\sum_i(\Delta s_i/R_i)
θ=∫0θΔθ=∫0s(ds/R)≈∑i(Δsi/Ri)
那么,以每条流线出口为起始点,选取一系列的计算点,量出各点的半径 R i R_i Ri和各间隔的轴面流线长度 Δ s i \Delta s_i Δsi。以 s s s为横坐标,作出 θ = f ( s ) \theta=f(s) θ=f(s)曲线来。然后用选定的 Δ s \Delta s Δs,从出口处截取纵坐标,过各分点做水平线和 θ = f ( s ) \theta=f(s) θ=f(s)曲线相交,从各交点向横坐标引垂线,得 1 , 2 , 3 … 1,2,3… 1,2,3…各点。按各分点 s s s长度,在轴面流线上截取相应分点即可。
这个方法倒是蛮有意思的,既简便又精度尚可,倒是可以用一用的。
不过呢,我们还是按照(1)和(2)方法的思路采用matlab编个小程序来实现就好了。
(4) 小程序精确计算法(强烈推荐)
算法思路如下:
-
周向的扇面 Δ θ \Delta\theta Δθ取 5 ° 5\degree 5°。
-
对于前盖板和后盖板,其是由直线段和圆弧段构成的,那么:
-
对于直线段,其满足 d s / d y = 1 + k 2 / k ds/dy=\sqrt{1+k^2}/k ds/dy=1+k2/k, k = tan ( α ) k=\tan(\alpha) k=tan(α)为直线的斜率, α \alpha α为直线与 x x x轴的夹角。于是 Δ s = Δ u \Delta s=\Delta u Δs=Δu等价于 d s / d y × ( y 0 − y 1 ) = ( y 0 + y 1 ) / 2 × Δ θ ds/dy\times(y_0-y_1) = (y_0+y_1)/2\times\Delta\theta ds/dy×(y0−y1)=(y0+y1)/2×Δθ,解此代数方程可以直接求得分点的坐标值 y 1 y_1 y1。进而进一步求得该直线段后续分点的坐标值 y 2 , y 3 , . . . y_2,y_3,... y2,y3,...,当然,最后该直线段会剩下一部分不满足此关系的,将其放在圆弧部分进行后续分点。
-
对于圆弧段,其左侧轴面流线上所分弧长为 Δ s = R ( θ 0 – θ 1 ) \Delta s= R(\theta_0 – \theta_1) Δs=R(θ0–θ1), θ 0 \theta_0 θ0为上一个分点的径向角度值, θ 1 \theta_1 θ1为待求点的径向角度值,其右侧周向圆弧的长度为 Δ u = ( y 0 + y 1 ) / 2 × Δ θ \Delta u=(y_0+y_1)/2\times\Delta\theta Δu=(y0+y1)/2×Δθ,其中 y 0 = O y + R c i r sin ( θ 1 ) y_0=O_y+R_{cir}\sin(\theta_1) y0=Oy+Rcirsin(θ1), y 1 = O y + R c i r sin ( θ 2 ) y_1=O_y+R_{cir}\sin(\theta_2) y1=Oy+Rcirsin(θ2)分别为上个分点和待求分点的 y y y坐标, O y O_y Oy为圆心y坐标, R c i r R_{cir} Rcir为圆弧半径。显然这里的方程是个超越方程,那么同样采用二分法或迭代法可以很容易求出分点值来。详见matlab程序。
- 接下来考虑中间流线的分点:
- 中间流线是由一系列控制点拟合而成的,在AutoCAD中采用的NURBS(非均匀有理B样条曲线),且不去管它背后的函数到底为何物。
- 我们将AutoCAD中拟合好的中间流线,拿list命令量取其长度S,定数等分为30000个点(试验发现AutoCAD最多能等分32767份,这里取个整数好计算一些),那么每段的长度为S/30000。
- 然后用AutoCAD的数据提取功能将这些点的x和y坐标输出来。
- 用matlab读取这些点坐标值,然后按从出口到进口的顺序排列好,接下来从出口点往进口点一个一个点去寻找,左侧中间流线待分曲线的累计长度为 Δ s = n S u m × ( S / 30000 ) \Delta s=nSum\times (S/30000) Δs=nSum×(S/30000),而右侧中间周向圆弧的长度为 ( y 0 + y 1 ) / 2 × Δ θ (y_0+y_1)/2\times\Delta\theta (y0+y1)/2×Δθ。不难想象,一开始左端肯定是小于右端的,而随着我们逐渐往下挨个找点,必然有一个点满足左侧突然大于等于右侧了,那么这个点就是我们想要找的分点了。
- 找好一个点后,继续下一个点的寻找工作,直到把所有点都找完了为止。
- 事实上,对于本例来讲,划分了30000段后,每一段的长度为0.0043212mm,那么这样找点的最大误差即为0.005mm,虽然还有一定的误差,但是也完全满足工程上的需要。
- 实际上,如果前盖板和后盖板懒得算的话,也可以用这里的方法来分点,只是精度上略有损失而已。
将前后盖板分点的x和y坐标直接绘制成点,点在AutoCAD中即可。也可以生成一个如下格式的文件
point x1,y1
point x2,y2
…
point xn,yn
然后把它们直接粘到AutoCAD命令行里面去,让AutoCAD一下子把这些点全部点出来。
程序跑出来的分点画在轴面投影图上是这样子的:
附上matlab程序代码:
前盖板流线分点源码
% 2016-11-06
% 离心泵设计 之 流线分点 之 前盖板流线分点
% 与中间流线分点方式不同,采用更加精准的分点方式来完成计算工作,直接算到精准等
% 与前盖板的分点方法相同
clear; clc;
DeltaTheta_Rad = 5*pi/180; % 周向弧段弧度
% A-B 上部斜线段;B-C 大圆弧
Ax = 51.007; Ay = 127.5; % 前盖板与出口边交点
Bx = 49.669; By = 112.2069; % 前盖板斜线段与圆弧段交点
Ox = -176.5668; Oy = 132; % 前盖板大圆弧段圆心
Rad = 227.1; % 大圆弧段半径
TSta = 355; TEnd = 345.090183; % 大圆弧段起始角度与终止角度,注意:顺时针排列!
O2x = 18.7289; O2y = 80; % 小圆弧段圆心,半径,起始角度与终止角度
Rad2 = 25;
T2Sta = 345.090183; T2End = 270;
% 划分斜线段
% 计算y的变化与斜线段长度的比例关系dSdY
if Ay == By % 若前盖板斜线段为90度
dSdY = 1;
dXdY = 0;
else % 若非90度
dx = abs(Ax - Bx); dy = abs(Ay - By);
dSdY = sqrt(dx^2 + dy^2) / dy; % y的变化与斜线段长度的比例关系
dXdY = (Bx - Ax) / (By - Ay);
end
idxFnd = 0; % 找到的当前分点指标
pntXY = zeros(100, 2); % 分得点坐标
pntId = zeros(100, 1); % 分得点位于哪一段上?0-斜线段,1-大圆弧,2-小圆弧,3-水平段,后面画叶片截线用
pntDs = zeros(100, 1); % 分得点距离上一个点的曲线长度DeltaS,后面画叶片截线用
% 斜线段的分点应该满足 (y0-y1)*dSdY = (y0+y1)/2*deltaTheta
y0 = Ay;
for j = 1: 100 % 默认循环100次,实际上根本没这么多点,中间break了
y1 = (dSdY - DeltaTheta_Rad/2) / (dSdY + DeltaTheta_Rad/2) * y0;
% dSLne = (y0 - y1)*dSdY;
% dSArc = (y0 + y1)/2*DeltaTheta_Rad;
% disp([dSLne, dSArc]);
if y1 < By
% 已经超出斜线段的范围了,在圆弧中继续寻找
dSPrvLft = abs(y0 - By)*dSdY; % 将剩余的斜线长度存储下来,y0值不变
break;
else
% 分得点的标识及其x和y坐标
idxFnd = idxFnd + 1;
pntXY(idxFnd, 2) = y1;
pntXY(idxFnd, 1) = Ax + dXdY * (y1 - Ay);
pntId(idxFnd) = 0;
pntDs(idxFnd) = dSdY * abs(y0 - y1);
% 更新起始点信息
y0 = y1;
end
end
% 继续划分下一段圆弧的流线分点
% 同样采用2分法来寻找分点
T0 = TSta;
TUpp = TSta; TDwn = TEnd; % [355, 345.090183]
for j = 1: 100 % 同上
if j ~= 1 % 第一个点要加上上一段直线剩余的长度
dSPrvLft = 0;
end
for idxDiv = 1: 30
TMid = (TUpp + TDwn) / 2;
dSLne = dSPrvLft + Rad * abs(T0 - TMid)*pi/180; % 流线上该段长度
y1 = Oy + Rad * sin(TMid*pi/180); % 当前点y值
dSArc = (y0 + y1)/2*DeltaTheta_Rad; % 当前右端中间圆弧长度
if dSLne < dSArc
TUpp = TMid;
else
TDwn = TMid;
end
% disp([idxDiv, TMid, dSLne, dSArc]);
end
% 分得点的x和y坐标
idxFnd = idxFnd + 1;
pntXY(idxFnd, 2) = y1;
pntXY(idxFnd, 1) = Ox + Rad * cos(TMid*pi/180);
pntId(idxFnd) = 1;
pntDs(idxFnd) = dSLne;
% 更新起始点信息
T0 = TMid;
y0 = y1;
TUpp = TMid;
TDwn = TEnd;
% 若分点到达了大圆弧的末端,则直接跳出,进行后续小圆弧的分点工作
if ( abs(TMid-TEnd)<1e-6 )
if (abs(dSLne-dSArc)>1e-6) % 左右两段未相等呢
dSPrvLft = dSLne; % 存储末端剩余长度
pntXY(idxFnd, :) = [0, 0]; % 最后点洗掉
pntId(idxFnd) = 0;
pntDs(idxFnd) = 0;
idxFnd = idxFnd - 1; % 指标-1
else
dSPrvLft = 0;
end
break;
end
end
% 继续划分下一段小圆弧的流线分点
% 同样采用2分法来寻找分点
T0 = T2Sta;
TUpp = T2Sta; TDwn = T2End; % [345.090183, 270]
for j = 1: 100 % 同上
if j ~= 1 % 第一个点要加上上一段直线剩余的长度
dSPrvLft = 0;
end
for idxDiv = 1: 30
TMid = (TUpp + TDwn) / 2;
dSLne = dSPrvLft + Rad2 * abs(T0 - TMid)*pi/180; % 流线上该段长度
y1 = O2y + Rad2 * sin(TMid*pi/180); % 当前点y值
dSArc = (y0 + y1)/2*DeltaTheta_Rad; % 当前右端中间圆弧长度
if dSLne < dSArc
TUpp = TMid;
else
TDwn = TMid;
end
% disp([idxDiv, TMid, dSLne, dSArc]);
end
% 分得点的x和y坐标
idxFnd = idxFnd + 1;
pntXY(idxFnd, 2) = y1;
pntXY(idxFnd, 1) = O2x + Rad2 * cos(TMid*pi/180);
pntId(idxFnd) = 2;
pntDs(idxFnd) = dSLne;
% 更新起始点信息
T0 = TMid;
y0 = y1;
TUpp = TMid;
TDwn = T2End;
% 若分点到达了小圆弧的末端,则直接跳出,进行后续小圆弧的分点工作
if ( abs(TMid-T2End)<1e-6 )
if (abs(dSLne-dSArc)>1e-6) % 左右两段未相等呢
dSPrvLft = dSLne; % 存储末端剩余长度
pntXY(idxFnd, :) = [0, 0]; % 最后点洗掉
pntId(idxFnd) = 0;
pntDs(idxFnd) = 0;
idxFnd = idxFnd - 1; % 指标-1
else
dSPrvLft = 0;
end
break;
end
end
% 继续最后水平直线段的分点
Dx = 18.7289; Dy = 55;
Ex = 0; Ey = 55;
y1 = Dy;
x1 = Dx - ((y0 + y1)/2*DeltaTheta_Rad - dSPrvLft);
if x1 > Ex; % 未到达末端
idxFnd = idxFnd + 1;
pntXY(idxFnd, :) = [x1, y1];
pntId(idxFnd) = 3;
pntDs(idxFnd) = (y0 + y1)/2*DeltaTheta_Rad;
% 继续分点,全部为水平段,直接计算,逐个减去DeltaX即可
deltaX = y1 * DeltaTheta_Rad;
x2 = x1 - deltaX;
while(x2 > Ex)
idxFnd = idxFnd + 1;
pntXY(idxFnd, :) = [x2, y1];
pntId(idxFnd) = 3;
pntDs(idxFnd) = deltaX;
x2 = x2 - deltaX;
end
end
pntXY = pntXY(1: idxFnd, :);
pntId = pntId(1: idxFnd);
pntDs = pntDs(1: idxFnd);
pntXYIdDs = [ pntXY, pntId, pntDs ];
% pntXYIdDs = [pntXY(1: idxFnd, :), pntId(1: idxFnd), pntDs(1: idxFnd)];
save LiuXianFenDian_QianGaiBanLiuXian_FenDeDian.txt pntXYIdDs -ascii;
后盖板流线分点源码
% 2016-11-05
% 离心泵设计 之 流线分点 之 后盖板流线分点
% 与中间流线分点方式不同,采用更加精准的分点方式来完成计算工作,直接算到精准等
clear; clc;
DeltaTheta_Rad = 5*pi/180; % 周向弧段弧度
% A-B 上部斜线段;B-C 大圆弧
Ax = 69.007; Ay = 127.5; % 后盖板与出口边交点
Bx = 69.9893; By = 71.2217; % 后盖板斜线段与圆弧段交点
Ox = 0; Oy = 70; % 后盖板圆弧段圆心
Rad = 70; % 圆弧段半径
TSta = 1; TEnd = 270; % 圆弧段起始角度与终止角度,注意:顺时针排列!
% 划分斜线段
% 计算y的变化与斜线段长度的比例关系dSdY
if Ay == By % 若后盖板斜线段为90度
dSdY = 1;
dXdY = 0;
else % 若非90度
dx = abs(Ax - Bx); dy = abs(Ay - By);
dSdY = sqrt(dx^2 + dy^2) / dy; % y的变化与斜线段长度的比例关系
dXdY = (Bx - Ax) / (By - Ay);
end
idxFnd = 0; % 找到的当前分点指标
pntXY = zeros(100, 2); % 分得点坐标
pntId = zeros(100, 1); % 分得点位于哪一段上?0-斜线段,1-圆弧段,后面画叶片截线用
pntDs = zeros(100, 1); % 分得点距离上一个点的曲线长度DeltaS,后面画叶片截线用
% 斜线段的分点应该满足 (y0-y1)*dSdY = (y0+y1)/2*deltaTheta
y0 = Ay;
for j = 1: 100 % 默认循环100次,实际上根本没这么多点,中间break了
y1 = (dSdY - DeltaTheta_Rad/2) / (dSdY + DeltaTheta_Rad/2) * y0;
% dSLne = (y0 - y1)*dSdY;
% dSArc = (y0 + y1)/2*DeltaTheta_Rad;
% disp([dSLne, dSArc]);
if y1 < By
% 已经超出斜线段的范围了,在圆弧中继续寻找
dSPrvLft = abs(y0 - By)*dSdY; % 将剩余的斜线长度存储下来,y0值不变
break;
else
% 分得点的标识及其x和y坐标
idxFnd = idxFnd + 1;
pntXY(idxFnd, 2) = y1;
pntXY(idxFnd, 1) = Ax + dXdY * (y1 - Ay);
pntId(idxFnd) = 0;
pntDs(idxFnd) = dSdY * abs(y0 - y1);
% 更新起始点信息
y0 = y1;
end
end
% 继续划分下一段圆弧的流线分点
% 同样采用2分法来寻找分点
T0 = TSta;
TUpp = TSta; TDwn = TEnd - 360; % 1: -1: -90
for j = 1: 100 % 同上
if j ~= 1 % 第一个点要加上上一段直线剩余的长度
dSPrvLft = 0;
end
for idxDiv = 1: 30
TMid = (TUpp + TDwn) / 2;
dSLne = dSPrvLft + Rad * abs(T0 - TMid)*pi/180; % 流线上该段长度
y1 = Oy + Rad * sin(TMid*pi/180); % 当前点y值
dSArc = (y0 + y1)/2*DeltaTheta_Rad; % 当前右端中间圆弧长度
if dSLne < dSArc
TUpp = TMid;
else
TDwn = TMid;
end
% disp([idxDiv, TMid, dSLne, dSArc]);
end
% 分得点的x和y坐标
idxFnd = idxFnd + 1;
pntXY(idxFnd, 2) = y1;
pntXY(idxFnd, 1) = Ox + Rad * cos(TMid*pi/180);
pntId(idxFnd) = 1;
pntDs(idxFnd) = dSLne;
% 更新起始点信息
T0 = TMid;
y0 = y1;
TUpp = TMid;
TDwn = TEnd - 360;
end
pntXY = pntXY(1: idxFnd, :);
pntId = pntId(1: idxFnd);
pntDs = pntDs(1: idxFnd);
pntXYIdDs = [ pntXY, pntId, pntDs ];
save LiuXianFenDian_HouGaiBanLiuXian_FenDeDian.txt pntXYIdDs -ascii;
中间流线分点源码
% 2016-11-05
% 离心泵设计 之 流线分点 之 中间流线分点
% 在AutoCAD中将中间流线等距离划分为上万个节点,然后用数据提取功能将这些点的坐标
% 提取出来;用matlab读入这些点的坐标,然后从出口点开始逐个向下寻找匹配,直至找到
% 划分点为止,该划分点应满足累加长度与右侧圆弧段的长度大致相等(首个累加长度大于
% 圆弧段长度的点即可)。
% 注意:轴线y坐标为0!!
clear; clc;
DeltaTheta_Deg = 5; % 周向弧段弧度
DeltaTheta_Rad = DeltaTheta_Deg*pi/180;
A = load('ZhongJianLiuXianLiSanDian_SanWanGe.txt'); % 读入AutoCAD导出的流道中线大批离散点数据
A = A(end: -1: 1, :); % 逆序,从出口到进口
LghAll = 129.6360; % 流道中线总体的长度,AutoCAD中用list命令量取
DivNum = 30000; % 流道中线所划分的段数,AutoCAD中用定数等分划分流道中线
DisEch = LghAll / DivNum; % 流道中线每个小微段的曲线长度
% 从出口到进口逐渐寻找分点
y0 = A(1, 2); % 右侧大圆弧半径
idxSum = 0; % 累加等长小微段的数目
pntIdx = zeros(100, 1); % 假设最多找到了100个分点,其标识,多的抛掉
pntXY = zeros(100, 2); % 分得点坐标
idxFnd = 0; % 找到的当前分点指标
disp(' 轴向长度 周向长度 误差(%) ');
for j = 2: size(A, 1)
idxSum = idxSum + 1; % 累加段数 + 1
DisNow = idxSum * DisEch; % 左侧累加曲线长度
DisArc = (y0 + A(j, 2)) / 2 * DeltaTheta_Rad; % 右侧轴向中间圆弧长度
if (DisNow >= DisArc) % 首个流道中线分点长度大于右侧中间圆弧长度的点即为找到的划分点
idxFnd = idxFnd + 1;
pntIdx(idxFnd) = j; % 将分得点的j指标记录下来
pntXY(idxFnd, :) = A(j, :); % 分得点的坐标
disp([DisNow, DisArc, (DisNow-DisArc)/DisNow*100]); % Show Something
% 准备下一分点
y0 = A(j, 2);
idxSum = 0;
end
end
% 释放多余空间
pntIdx = pntIdx(1: idxFnd);
pntXY = pntXY(1: idxFnd, :);
pntSav = [pntIdx, pntXY];
save LiuXianFenDian_ZhongJianLiuXian_FenDeDian.txt pntSav -ascii;