螺旋桨设计之Bseries使用AU系列的方法设计（附加代码，附加学习资料)

《船舶原理》P105-P159中对AU系列的图谱法设计进行了非常详细的介绍，包括怎样使用图谱，中间计算公式，得到型面数据如何绘制二维三维图像进行CAD建模。然而，AU系列螺旋桨的图谱非常难找，反而B系列的图谱要好找得多。例如：《荷兰船模水池B系列螺旋桨新旧图谱的比较-孙永华》就给出了一个有雷诺数修正之后的图谱的参数表。然而，B系列螺旋桨图谱的坐标定义和AU系列有差别，所以要使用AU系列的设计方法设计B系列螺旋桨，就需要对B系列螺旋桨图谱进行处理，使得其符合AU系列的设计要求。本人给出转换代码，并且给出使用AU系列设计方法设计盘面比0.65的B系列螺旋桨设计的初始设计和终结设计的代码。供大家参考。
gendata.m

clear;clc;
data1 = [ 0.00880496 0 0 0 0 0.00379368 0 0 0 0
 -0.204554 1 0 0 0 0.00886523 2 0 0 0
0.166351 0 1 0 0  -0.032241 1 1 0 0
0.158114 0 2 0 0 0.00344778 0 2 0 0
 -0.147581 2 0 1 0  -0.0408811 0 1 1 0
 -0.481497 1 1 1 0  -0.108009 1 1 1 0
0.415437 0 2 1 0  -0.0885381 2 1 1 0
0.0144043 0 0 0 1 0.188561 0 2 1 0
 -0.0530054 2 0 0 1  -0.00370871 1 0 0 1
0.0143481 0 1 0 1 0.00513696 0 1 0 1
0.0606826 1 1 0 1 0.0209449 1 1 0 1
 -0.0125894 0 0 1 1 0.00474319 2 1 0 1
0.0109689 1 0 1 1  -0.00723408 2 0 1 1
 -0.133698 0 3 0 0 0.00438388 1 1 1 1
0.00638407 0 6 0 0  -0.0269403 0 2 1 1
 -0.00132718 2 6 0 0 0.0558082 3 0 1 0
0.168496 3 0 1 0 0.0161886 0 3 1 0
 -0.0507214 0 0 2 0 0.00318086 1 3 1 0
0.0854559 2 0 2 0 0.015896 0 0 2 0
 -0.0504475 3 0 2 0 0.0471729 1 0 2 0
0.010465 1 6 2 0 0.0196283 3 0 2 0
 -0.00648272 2 6 2 0  -0.0502782 0 1 2 0
 -0.00841728 0 3 0 1  -0.030055 3 1 2 0
0.0168424 1 3 0 1 0.0417122 2 2 2 0
 -0.00102296 3 3 0 1  -0.0397722 0 3 2 0
 -0.0317791 0 3 1 1  -0.00350024 0 6 2 0
0.018604 1 0 2 1  -0.0106854 3 0 0 1
 -0.00410798 0 2 2 1  0.00110903 3 3 0 1
 -0.000606848 0 0 0 2  -0.000313912 0 6 0 1
 -0.0049819 1 0 0 2 0.0035985 3 0 1 1
0.0025983 2 0 0 2  -0.00142121 0 6 1 1
 -0.000560528 3 0 0 2  -0.00383637 1 0 2 1
 -0.00163652 1 2 0 2 0.0126803 0 2 2 1
 -0.000328787 1 6 0 2  -0.00318278 2 3 2 1
0.000116502 2 6 0 2 0.00334268 0 6 2 1
0.000690904 0 0 1 2  -0.00183491 1 1 0 2
0.00421749 0 3 1 2 0.000112451 3 2 0 2
5.65229E-05 3 6 1 2  -2.97228E-05 3 6 0 2
 -0.00146564 0 3 2 2 0.000269551 1 0 1 2];
data2 = [0.00083265 2 0 1 2
0.00155334 0 2 1 2
0.000302683 0 6 1 2
 -0.0001843 0 0 2 2
 -0.000425399 0 3 2 2
8.69243E-05 3 3 2 2
 -0.0004659 0 6 2 2
5.54194E-05 1 6 2 2];
data3 = [0.000353485 0 0 0 0 0 -0.000591412 0 0 0 0 0
-0.00333758 0 1 0 2 0 0.00696898 0 0 1 0 0
-0.00478125 0 1 1 1 0 -6.67E-05 0 0 6 0 1
0.000257792 2 1 0 2 0 0.0160818 0 2 0 0 0
6.43E-05 1 0 6 2 0 -0.000938091 1 0 1 0 0
-1.11E-05 2 0 6 2 0 -0.00059593 1 0 2 0 0
-2.76E-05 2 1 0 2 1 7.82E-05 2 0 2 0 0
0.0000954 1 1 1 1 1 5.22E-06 1 1 0 2 1
3.20E-06 1 1 3 1 2 -8.85E-07 2 1 1 1 1];
data4 = [2.30E-05 1 0 6 0 1
-1.84E-06 2 0 6 0 1
-0.00400252 1 2 0 0 0
0.000220915 2 2 0 0 0];
Ctn = data1(:,1);tsn = data1(:,2);ttn = data1(:,3);tun = data1(:,4);tvn = data1(:,5);
Cqn = [data1(:,6);data2(:,1)];qsn = [data1(:,7);data2(:,2)];qtn = [data1(:,8);data2(:,3)];
qun = [data1(:,9);data2(:,4)];qvn = [data1(:,10);data2(:,5)];
Cdtn = data3(:,1);tdsn = data3(:,2);tdtn = data3(:,3);tdun = data3(:,4);tdvn = data3(:,5);tdwn = data3(:,6);
Cdqn = [data3(:,7);data4(:,1)];qdsn = [data3(:,8);data4(:,2)];qdtn = [data3(:,9);data4(:,3)];
qdun = [data3(:,10);data4(:,4)];qdvn = [data3(:,11);data4(:,5)];qdwn = [data3(:,12);data4(:,6)];
clear data1 data2 data3 data4
%% 图谱计算
J = 0:0.0005:1.5;%进速系数
Z = 3;%叶数
P_D = 0.25:0.0002:1.2;%螺距比
pmb = 0.65;%盘面比
len = length(P_D);
Kt = zeros(len,length(J));
Kq = zeros(len,length(J));
eta = zeros(len,length(J));
for i = 1:len
    for j = 1:length(J)
        Kt(i,j) = sum(Ctn.*J(j).^tsn.*P_D(i).^ttn.*pmb.^tun.*Z.^tvn);
        if Kt(i,j)<=0 continue;end
        Kq(i,j) = sum(Cqn.*J(j).^qsn.*P_D(i).^qtn.*pmb.^qun.*Z.^qvn);
        eta(i,j) = Kt(i,j).*J(j)/(2.*pi.*Kq(i,j));
    end
end
%% Bp-delta图
%得到P_D-bp-delta-eta的数据组
rawdata = [];
for i = 1:len
    long = sum(Kq(i,:)~=0)+1;
    sqrtbp = sqrt(33.3.*Kq(i,1:long).^0.5./J(1:long).^2.5);
    delta = 33.86./J(1:long);
    rawdata(end+1:end+long,:) = [P_D(i).*ones(long,1),sqrtbp',delta',eta(i,1:long)'];
end
rawdata = rawdata(~ismember(rawdata(:,2),inf),:);

%delta等值线
dta = 40:10:240;
sigma_delta = 0.5;
figure(2);hold on; box on;%delta等值线
for i = 1:length(dta)
        tempdelta = dta(i);
        tempdata = rawdata(rawdata(:,3)<=tempdelta+sigma_delta & rawdata(:,3)>=tempdelta-sigma_delta,:);%挑选出delta在误差范围内的数据
        tempdata(tempdata(:,2)==0,:) = [];
        %挑选出每个P_D最接近的delta所对应的bp
        member = unique(tempdata(:,1));
        pour = [];
        for j = 1:length(member)
            temp = tempdata(tempdata(:,1)==member(j),:);
            temp(:,3) = abs(temp(:,3)-dta(i));
            [~,idx]=sort(temp(:,3),'ascend');
            temp=temp(idx,:);
            pour(end+1,:) = temp(1,:);
        end
        [~,idx]=sort(pour(:,2));
        pour=pour(idx,:);
        if size(pour,1)~=1
            plot(pour(10:end,2),pour(10:end,1),LineWidth=3,Color='k')
        else
            disp(i)
        end
end

%eta等值线
eta = 0.3:0.1:0.8;
sigma_eta = 0.0001;
figure(3);hold on; box on;%eta等值线
xiaolvdata = [];
for i = 1:length(eta)
        tempeta = eta(i);
        tempdata = rawdata(rawdata(:,4)<=tempeta+sigma_eta & rawdata(:,4)>=tempeta-sigma_eta,:);%挑选出delta在误差范围内的数据
        [~,idx]=sort(tempdata(:,4));
        tempdata=tempdata(idx,:);
        %提取效率大于给定效率中效率差最小的一个数和效率大小于给定效率中效率差最小的一个数
        member = unique(tempdata(:,1));
        pour = [];
        for j = 1:length(member)
             idx = find(tempdata(:,1)==member(j));
             pour1 = tempdata(idx,:);
             pour2 = tempdata(idx,:);
            if sum(pour1(:,4)<=eta(i))
            pour1 = pour1(pour1(:,4)<=eta(i),:);
            [~,idx]=sort(pour1(:,4),'descend');
            pour1=pour1(idx,:);
            pour1 = pour1(1,:);
            pour(end+1,:) = pour1;
            end
            if sum(pour2(:,4)>eta(i))
            pour2 = pour2(pour2(:,4)>eta(i),:);
            [~,idx]=sort(pour2(:,4),'ascend');
            pour2=pour2(idx,:);
            pour2 = pour2(1,:);
            pour(end+1,:) = pour2;
            end
        end
        if sum(sum(pour))
            scatter(pour(:,2),pour(:,1),MarkerEdgeColor='k',SizeVariable=1)
            xiaolvdata = [xiaolvdata;pour];
        end
end
title('$$\sqrt{\mathrm{B}_{\mathrm{P}}}-\mathrm{\delta}$$','Interpreter','latex')
xlabel("$\sqrt{\mathrm{B}_{\mathrm{p}}}$",'Interpreter','latex')
ylabel("$$\frac{\mathrm{P}}{\mathrm{D}}$$",'Interpreter','latex')
legend(['$$\mathrm{\eta}$$=0.3'],['$$\mathrm{\eta}$$=0.4'],['$$\mathrm{\eta}$$=0.5'],['$$\mathrm{\eta}$$=0.6'],['$$\mathrm{\eta}$$=0.7'],['$$\mathrm{\eta}$$=0.8'],'Interpreter','latex')

%寻找最佳效率线
bp = 2:0.1:18;
sigma_bp = 0.01;
pourfinal = [];
for i = 1:length(bp)
    tempbp = bp(i);
    tempdata = rawdata(rawdata(:,2)<=tempbp+sigma_bp & rawdata(:,2)>=tempbp-sigma_bp,:);%挑选出delta在误差范围内的数据
    [~,idx]=sort(tempdata(:,4),'descend');
    tempdata=tempdata(idx,:);
    pourfinal = [pourfinal;tempdata(1,:)];
end
[~,idx]=sort(pourfinal(:,2),'descend');
pourfinal=pourfinal(idx,:);
figure(4);hold on;box on;
huiguix = pourfinal(:,2);
huiguiy = pourfinal(:,1);
[fitresult, ~] = fit( huiguix, huiguiy, 'power2');
plot(fitresult,huiguix,huiguiy)
title("回归得到的最佳效率曲线")
xlabel("$\sqrt{\mathrm{B}_{\mathrm{p}}}$",'Interpreter','latex')
ylabel("$$\frac{\mathrm{P}}{\mathrm{D}}$$",'Interpreter','latex')
cof = coeffvalues(fitresult);%预测行驶： cof(1).*bp.^cof(2)+cof(3)

enddesign.m

clear;clc;
load pmb0.65.mat
D = 2;%螺旋桨直径
w = 0.12;%伴流分数
t = 0.2;%推力减额
eta_h = (1-t)/(1-w);
eta_s = 0.98;
eta_e = 0.98;
Vs = 20;
Va = Vs*(1-w);%螺旋桨进速
Pe = 26*1.36*Vs*0.514;%有效马力
N = 290:5:320;%给定转速
zhijingxishu = N.*D./Va;%直径系数
sigma_zjxs = 0.01;
chubusheji = zeros(length(zhijingxishu),4);
for i = 1:length(zhijingxishu)
    tempzjxs = zhijingxishu(i);
    tempdata = rawdata(rawdata(:,3)<=tempzjxs+sigma_zjxs & rawdata(:,3)>=tempzjxs-sigma_zjxs,:);%挑选出delta在误差范围内的数据
    tempdata(tempdata(:,2)==0,:) = [];
% scatter(tempdata(:,2),tempdata(:,1))
    bijiaodata = abs(cof(1).*tempdata(:,2).^cof(2)+cof(3)-(tempdata(:,1)));
    idx = find(bijiaodata==min(bijiaodata));
    if isempty(idx)
        continue
    else
    chubusheji(i,:) = tempdata(idx(1),:);
    end
end
Pd = chubusheji(:,2).^4.*Va.^5./N'.^2;
Ps = Pd./eta_s./eta_e;
Pte = Pd.*chubusheji(:,4).*eta_h;
figure(5);hold on; box on;
[AX1,H1,H2]=plotyy(N,[Pte,Pe.*ones(size(N,2),size(N,1))],N,[chubusheji(:,1),chubusheji(:,4)]);
set(AX1(1),'XColor','k','YColor','k'); %X轴和第一个Y轴的颜色
set(AX1(2),'XColor','k','YColor','k'); %X轴和第二个Y轴的颜色
HH1=get(AX1(1),'Ylabel');
set(HH1,'String','Pe、Pte','color','k'); %第一个Y轴标签的的颜色
HH2=get(AX1(2),'Ylabel');
set(HH2,'String','P/D、η0','color','k'); %第二个Y轴标签的的颜色
set(H1(1),'Linestyle','-','color','k','Linewidth',3,'LineStyle','-','Marker','d');
set(H1(2),'Linestyle','-','color','k','Linewidth',3,'LineStyle','-','Marker','o');
set(H2(1),'Linestyle','-','color',zchcolor(5),'Linewidth',3,'LineStyle','-','Marker','p');
set(H2(2),'Linestyle','-','color',zchcolor(5),'Linewidth',3,'LineStyle','-','Marker','v');
xlabel('螺旋桨转速r/min');
legend([H1(1),H1(2),H2(1),H2(2)],{'Pte';'Pe';'P_D';'eta0'});

%求Pe和Pte的交点处的转速，马力，P_D，eta_0
% [fitresult, ~] = fit( N', Pte, 'poly2');
% p = coeffvalues(fitresult);%预测行驶： p(1).*n.^2+n.*p(2)+p(3)
% syms x
% N_best = solve(p(1).*x.^2+x.*p(2)+p(3)==Pe,x);N_best = double(N_best);N_best = round(N_best(N_best>0));
bijiaopte = abs(Pte-Pe);
idx = find(bijiaopte==min(bijiaopte));
disp("************初步设计************")
N_best = N(idx);disp(['最佳转速为:',num2str(N_best),'r/min'])
Pzhuji = Ps(idx);disp(['所需主机马力为：',num2str(Pzhuji),'hp'])
P_D_best = chubusheji(idx,1);disp(['P/D为：',num2str(P_D_best)])
eta_0_best = chubusheji(idx,4);disp(['η0为：',num2str(eta_0_best)])
disp("****************************")
%% 终结设计
% clear;clc;
% D = 1.15;%螺旋桨直径
% w = 0.12;%伴流分数
% t = 0.14;%推力减额
% eta_h = (1-t)/(1-w);
% eta_s = 0.98;
% eta_e = 0.98;
% Vs = 20;
% Va = Vs*(1-w);%螺旋桨进速
% Pe = 9.3*1.36*Vs*0.514;%有效马力
% load rawdata2.mat
V = [Vs-2,Vs-1,Vs,Vs+1,Vs+2];
Va = V.*(1-w);
Pd = eta_s*eta_e*Pzhuji;
Bp = (N_best.*Pd.^0.5./Va.^2.5).^0.5;
% 找根号bp和最佳效率曲线的交点
sigma_Bp = 0.00005;
zhongjiesheji = zeros(length(Bp),4);
for i = 1:length(Bp)
    tempbp = Bp(i);
    tempdata = rawdata(rawdata(:,2)<=tempbp+sigma_Bp & rawdata(:,2)>=tempbp-sigma_Bp,:);%挑选出delta在误差范围内的数据
    tempdata(tempdata(:,2)==0,:) = [];
    bijiaodata = abs(cof(1).*tempdata(:,2).^cof(2)+cof(3)-(tempdata(:,1)));
    idx = find(bijiaodata==min(bijiaodata));
    if length(idx)>1
        disp("1")
    end
    zhongjiesheji(i,:) = tempdata(idx(1),:);
end
Pte = Pd.*zhongjiesheji(:,4).*eta_h;
drag = [21407,23651,26001,28457,31017]./1000;
Pee = drag.*V.*1.36.*0.514;
%绘制有效马力曲线
% figure(100);plot(V,Pee,'LineWidth',3);xlabel("航速/kn");ylabel("马力/hp");title("水下航行器有效马力曲线");


figure(6);hold on;box on;
[AX2,H3,H4] = plotyy(V,[Pte,Pee'],V,[zhongjiesheji(:,1),zhongjiesheji(:,4),zhongjiesheji(:,3).*Va'./N_best]);
set(AX2(1),'XColor','k','YColor','k'); %X轴和第一个Y轴的颜色
set(AX2(2),'XColor','k','YColor',zchcolor(5)); %X轴和第二个Y轴的颜色
HH1=get(AX2(1),'Ylabel');
set(HH1,'String','Pte、Pee','color','k'); %第一个Y轴标签的的颜色
HH2=get(AX2(2),'Ylabel');
set(HH2,'String','P/D、η0、D','color','k'); %第二个Y轴标签的的颜色
set(H3(1),'Linestyle','-','color','k','Linewidth',3,'LineStyle','-','Marker','o');
set(H3(2),'Linestyle','-','color','k','Linewidth',3,'LineStyle','-','Marker','d');
set(H4(1),'Linestyle','-','color',zchcolor(5),'Linewidth',3,'LineStyle','-','Marker','p');
set(H4(2),'Linestyle','-','color',zchcolor(5),'Linewidth',3,'LineStyle','-','Marker','v');
set(H4(3),'Linestyle','-','color',zchcolor(5),'Linewidth',3,'LineStyle','-','Marker','+');
xlabel('航速/kn');
lgd = legend([H3(1),H3(2),H4(1),H4(2),H4(3)],{'Pte';'Pe';'P_D';'eta0';'D'});
lgd.NumColumns = 2;
leg.Orientation = 'vertical';
%求Pe和Pte的交点处的转速，马力，P_D，eta_0
disp("************终结设计************")
[fitresult, ~] = fit( V', Pte, 'poly2');%把Pet拟合
p1 = coeffvalues(fitresult);%预测行驶： p(1).*n.^2+n.*p(2)+p(3)
[fitresult, ~] = fit( V', Pee', 'poly2');%把Pe拟合
p2 = coeffvalues(fitresult);%预测行驶： p(1).*n.^2+n.*p(2)+p(3)
syms x%求解可达航速
V_available = solve(p1(1).*x.^2+x.*p1(2)+p1(3)==p2(1).*x.^2+x.*p2(2)+p2(3),x);
V_available = double(V_available);V_available = max(V_available);
%拟合螺旋桨直径
D_available = zhongjiesheji(:,3).*Va'./N_best;
[fitresult, ~] = fit( V', zhongjiesheji(:,3).*Va'./N_best, 'poly2');Dp1 = coeffvalues(fitresult);
D_best = Dp1(1).*V_available.^2+V_available.*Dp1(2)+Dp1(3);
%拟合螺距比
P_D_available = zhongjiesheji(:,1);
[fitresult, ~] = fit( V', zhongjiesheji(:,1), 'poly2');P_D_p1 = coeffvalues(fitresult);
P_D_best = P_D_p1(1).*V_available.^2+V_available.*P_D_p1(2)+P_D_p1(3);
%拟合螺旋桨效率
eta_available = zhongjiesheji(:,4);
[fitresult, ~] = fit( V', zhongjiesheji(:,4), 'poly2');etap1 = coeffvalues(fitresult);
eta_best = etap1(1).*V_available.^2+V_available.*etap1(2)+etap1(3);
%拟合有效推马力
[fitresult, ~] = fit( V', Pte, 'poly2');Ptep1 = coeffvalues(fitresult);
Pte_best = Ptep1(1).*V_available.^2+V_available.*Ptep1(2)+Ptep1(3);

disp(['可达航速为:',num2str(V_available),'kn'])
disp(['螺旋桨直径为：',num2str(D_best),'m'])
disp(['P/D为：',num2str(P_D_best)])
disp(['η0为：',num2str(eta_0_best)])
disp(['Pte为：',num2str(Pte_best),'hp'])
disp("****************************")

zchcolor.m

function map = zchcolor(N)
color = [102,103,171;...%长春花蓝
    147,149,151;...%极致灰
    245,223,77;...%靓丽黄
    15,76,129;...%经典蓝
    225,111,97;...%活珊瑚橘
    95,75,139;...% 紫外光
    136,176,75;...% 草木绿
    247,202,201;...% 水晶粉
    146,168,209;...% 宁静蓝
    150,79,76;...% 玛萨拉酒红
    177,99,163;...% 璀璨紫兰花
    0,152,116;...% 翡翠绿
    226,73,47;...% 探戈橘
    203,101,134;...% 忍冬红
    69,181,170;...% 松石绿
    240,192,90;...% 含羞草黄
    90,91,159;...% 鸢尾蓝
    230,0,0;...% 中国红
    ];
color = color./255;
map = color(N,:);

使用说明：将三个文件放在一个文件夹内，首先运行gendata.m，然后将工作区另存为pmb0.65.mat，然后即可运行enddesign.m，如果要使用不同盘面比的图谱，可以在gendata.m中修改参数pmd参数即可。
学习资料上传百度网盘了，包括PropExpert和一些船舶螺旋桨方面的书籍《船舶原理（上下册）上海交通大学出版社》《船舶螺旋桨理论与应用》《船用螺旋桨技术研究及系列图谱》《现代泵理论与设计》《鱼雷电动力技术》还有一些学术论文等等。
链接：https://pan.baidu.com/s/1sB2E6AziS1Cwlcu5CuIdSQ?pwd=1234
提取码：1234