主程序

clear;

global r;

global R;

global phi;

global lamb;

V1=7;P=400;lamb0=6;B=3;

D=2.5;

R=D/2;

N=60*lamb0*V1/pi/D;

omiga=V1*lamb0/R;

r=0.1*R;                                                                                   %当前只进行0.1R处翼型截面的气动力计算试一试

a0=lamb0^2/2000+0.027*lamb0-0.038;

b0=10/lamb0/exp(12*r/R);

lamb=lamb0*r/R;

x0=[a0;b0];

for i=1:10

lb=[0;0];

ub=[0.5;1];

options=optimset('largescale','off','Algorithm','active-set');

[x,fval]=fmincon('objfun',x0,[],[],[],[],lb,ub,'confun',options);

a=x.(1);

b=x.(2);

xi=x;

if((abs(a-a0)<=0.00001)&&(abs(b-b0)<=0.00001)),break,end

x0=xi;

end

C0=0.21897;

for i=0:10

Re=[100000 500000 1000000 2000000 3000000 4000000 5000000 9000000];

Cl=[1.2424 1.0853 1.001 0.9616 0.9663 0.9693 0.9016 1.0023];

alpha=[8 6 5 4.5 4.5 4.5 4 5];

Re(i)=V1*(1-a)*C0/1.48e-5/sin(phi);

Cl(i)=interpi(Re,Cl,Re(i),'spline');

alpha(i)=interpi(Re,alpha,Re(i),'spline');

S=8*pi*sin(phi)*sin(phi)*(1-a*F)*a*F/cos(phi)/(1-a)^2;

C(i+1)=S*r/B/Cl(i);

C=C(i+1);

if(abs(C-C0)<=0.00001),break,end

C0=C(i+1);

end

theta=phi-alpha(i);

printf(fid,'%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n',a,b);

目标函数

function fval=objfun(x0)

global lamb0;

global B;

global R;

global r;

global phi;

lamb=lamb0*r/R;

phi=atan((1-x0.(1))/(1+x0.(2))/lamb);

f=B*(R-r)/2/r/sin(phi);

F=2*acos(exp(-f))/pi;

fval=-8*x0.(2)*(1-x0.(1))*F*lamb^3/lamb0^2;

条件函数

function [c,ceq]=confun(x0)

global lamb0;

global B;

global R;

global r;

global phi;

lamb=lamb0*r/R;

phi=atan((1-x0.(1))/(1+x0.(2))/lamb);

f=B*(R-r)/2/r/sin(phi);

F=2/pi*acos(exp(-f));

ceq=x0.(1)*(1-x0.(1)*F)-x0.(2)*(1+x0.(2))*lamb^2;

c=[];

程序采用的是Wilson算法，但是在MATLAB中运行的时候总是出问题。我是菜鸟，该程序是边学边做的，可能有出错的地方，望大侠们帮帮吧。。。