📄 mainprogram.m
字号:
end
%---41尾 %41 Continue
%***********************************拟合桨叶厚度*****************************
% *** Curve Fit Blade Thickness ***
%
for I = 1:(N+1) %Do 42 I = 0,N
F(2,I) = T(I);%F(1,I) = T(I)
end %42 Continue
%**************************************************************************
SPLINE;%Call SPLINE(N, F)
%**************************************************************************
for I = 1:10 %Do 43 I = 0,9
FX = X(I); %FX = X(I)
[FY]=FNEVAL(FX, P3, F);%Call FNEVAL(FX, FY, P3, F)
T(I) = FY; %T(I) = FY
end %43 Continue
%拟合桨叶角
% *** Curve Fit Blade Angles ***
%
for I = 1:(N+1) %Do 44 I = 0,N
F(2,I) = BETA(I);%F(1,I) = Beta(I)
end %44 Continue
%**************************************************************************
SPLINE;%Call SPLINE(N, F)
%**************************************************************************
for I = 1:10 %Do 45 I = 0,9
FX = X(I); %FX = X(I)
[FY]=FNEVAL(FX, P3, F);%Call FNEVAL(FX, FY, P3, F)
BETA(I) = FY; %Beta(I) = FY
end %45 Continue
%
fprintf('**********拟合完毕*********\n')
%
% 调入翼型数据
% *** Load airfoil data from file ***
%
fprintf('**********Reading Cl & Cd data for the--- %s---airfoil.\n',AFds(AF,:));%Write(imon,46) AFds(AF)
%46 FORMAT(1X,'Reading Cl & Cd data for the ',A20,' airfoil.')
if AF==1
flatplat;
end
if AF==2 %检查AF调入相应的翼型文件
ClarkY;
end
if AF==3
RAF6;
end
if AF==4
NACA44Lo;
end
if AF==5
Symmetri;
end %OPEN (UNIT=43, FILE = AFds(AF), STATUS='OLD')
AFNAME; %Read(43, 446) AFNAME
%446 format(1x,A20)
%c调入(I = 0,ROWS Z(0,I,J), J=0,COLS)
ROWS=rows; %￥￥￥￥￥￥￥￥￥在此就把ROWS和COLS增加1，则后面CHECKDATA调用时就无须再次增加1
COLS=cols; %Read(43,*) ROWS,COLS
DUMMYNAME1;%Read(43,446) DUMMYNAME
for I = 1:ROWS %Do 50 I = 0,ROWS
for J=1:COLS%Read(43,*) (Z(0,I,J), J=0,COLS)
Z(1,I,J)=Lift(I,J); %升力系数矩阵
end
end %50 Continue
DUMMYNAME2;%Read(43,446) DUMMYNAME
for I = 1:ROWS%Do 51 I = 0,ROWS
for J=1:COLS%Read(43,*) (Z(1,I,J), J = 0,COLS)
Z(2,I,J)=Drag(I,J); %阻力系数矩阵
end
end %51 Continue
% 继续录入
REAFD;
ACL;
BCD; %READ(43,*) REAFD,ACL,BCD
%CLOSE(unit=43)关闭翼型文件
if (RA==3) % 采用该方式-----Analysis including induced velocity and tip losses
%**************************** Load Goldstein Constants ********************
%
fprintf('Reading tip loss data for- %4i -blades',NB);
if (NB==2)
%cWrite(imon,*) 'Reading tip loss data for', NB, ' blades. '
GoldK2;%OPEN(UNIT=55, FILE = 'GoldK2.d',STATUS='OLD')
end
if (NB==3)
%cWrite(imon,*) 'Reading tip loss data for', NB, ' blades. '
GoldK3;%OPEN(UNIT=55, FILE = 'GoldK3.d',STATUS='OLD')
end
if (NB==4)
%cWrite(imon,*) 'Reading tip loss data for', NB, ' blades. '
GoldK4;%OPEN(UNIT=55, FILE = 'GoldK4.d',STATUS='OLD')
end
Goldstr;%Read(55,447) Goldstr
%447 format(1x,A3)
for I = 1:10 %Do 53 I = 0,9
for J=1:12 %Read(55,*) (GK(I,J), J=0,11)
GK(I:J);%(读入桨叶损失文件)
end
end%53 Continue
%CLOSE(unit=55)---关闭桨叶损失文件
end
%
N=9;% N=9;
% C
[T,Fzero]=CHECKDATA(N,Z,BETA,X,T,ROWS,LDA,MPH,RPM,A,Radius,S28,COLS,Jmin,Jmax,H);%CALL CHECKDATA(N,Z,Beta,X,T,ROWS,LDA,MPH,RPM,A,Radius,
%& S28,COLS,Jmin,Jmax,F1,F2,F3,Fzero,H)
%
% *** Calculate Solidity and Pitch ***
%
Kgold = 1.0;% Kgold = 1.0
fprintf('\n');% Write(imon,*)
for I = 1:10% Do 60 I = 0,9
S(I) = NB*C(I)/(3.14159*Radius);% S(I) = Float(NB)*C(I)/(3.14159*Radius)
GP(I) = 2.0*3.14159*Radius*X(I)*tan(BETA(I)*DR);% GP(I) = 2.0*3.14159*Radius*X(I)*TAN(BETA(I)*DR)
end % Continue
%
% *** Start J loop, Calculate PHI开始计算 J PHI***
%
fprintf('********Running calculations for ADVANCE RATIO, J, and RADIAL POSITION, X:*****\n');% Write(imon,*) 'Running calculations for ADVANCE RATIO, J, and RADIAL POSITION, X:'
N=9;
% Real(J)强制转换为实型 ，创建进距比矩阵Jar(J)，11个元素
for J = 1:11%Do 61 J = 0,10
Jar(J) = Jmin + (J-1)*(Jmax-Jmin)/10.0;% Jar(J) = Jmin + Real(J)*(Jmax-Jmin)/10.0
%fprintf('%4i,%6f \n',J,Jar(J)); %c Write(imon,555) J,Jar(J)
%c555 FORMAT(1x,'J (', i4, ') = ', f5.2)
for I=1:9
M=0;
PHI=atan(Jar(J)/(X(I)*3.14159));
ALPHA(I,J)=BETA(I)*DR-PHI;%此处攻角单位是弧度
Theta=0.0;
if ((RA==2)|(RA==3))
FLAG=-1;
while (FLAG<=0)
SF=sin(PHI);
CF=cos(PHI);
[II,KK,CL,CD,RE,MT,SF,CF]=CLCDCALC(PHI,T,Z,ALPHA,I,J,ROWS,MPH,LDA,C,NU,A,Radius,SF,CF,Azero,X,RE,MT,MTUC,RPM,REAFD,ACL,BCD);
Azero=(Z(1,II+1,KK) - Z(1,II,KK))/((Z(1,II+1,1)-Z(1,II,1))*DR);%攻角区间上的升力线斜率 1/rad
if (RA==2)
Kgold=1.0;%不计桨叶尖损失，只考虑诱导速度修正
end
if (RA==3)
kk=2;
while (SF>=GK(1,KK+1))
KK=KK+1;
end%考虑桨叶尖损失，拟合出损失系数
Kgold=GK(I+1,KK) + (SF-GK(1,KK))*(GK(I+1,KK+1)-GK(I+1,KK))/(GK(1,KK+1)-GK(1,KK));
end
Thetaold=Theta;
Thetanew = S(I)*CL/(8.0*Kgold*X(I)*SF+S(I)*Azero);
Deltheta=Thetanew-Thetaold;
if (abs(Deltheta)>0.7)
Theta=Thetaold+(0.5-0.45*(M/30.0))*Deltheta/(abs(Deltheta))*0.7;
else
Theta=Thetaold+ (0.5-0.45*(M/30.0))*Deltheta;
end
PHI=atan(Jar(J)/(X(I)*3.14159)) + Theta;
ALPHA(I,J) = BETA(I) * DR - PHI;
M = M +1;
FLAG=1;
if (RA>1)&(abs(sin(PHI) - SF)>C2)&(M<30)
FLAG=-1;
continue
end
end
end
SF=sin(PHI);
CF=cos(PHI);
[II,KK,CL,CD,RE,MT,SF,CF]=CLCDCALC(PHI,T,Z,ALPHA,I,J,ROWS,MPH,LDA,C,NU,A,Radius,SF,CF,Azero,X,RE,MT,MTUC,RPM,REAFD,ACL,BCD);
if (M==30)
F5=1;
end
%
D1 = cos(Theta)/cos(PHI - Theta);%D1 = COS(Theta)/COS(PHI - Theta)
D1 = D1*D1;% D1 = D1*D1
D2 = (CL*SF + CD*CF)*D1;% D2 = (CL*SF + CD*CF)*D1
D1 = (CL*CF - CD*SF)*D1;% D1 = (CL*CF - CD*SF)*D1
DT(I,J) = Czero*X(I)*X(I)*S(I);% DT(I,J) = Czero*X(I)*X(I)*S(I)
DQ(I,J) = DT(I,J)*X(I)*D2/2.0;% DQ(I,J) = DT(I,J)*X(I)*D2/2.0
DT(I,J) = DT(I,J) * D1;% DT(I,J) = DT(I,J) * D1
end%62 Continue 62循环结束，还有61循环
%C
%
%
% *** Integrate to get Cp(j) and Ct(j) ***
%
for I = 1:(N+1)%Do 70 I = 0,N
F(1,I) = X(I);% F(0,I) = X(I)
F(2,I) = DQ(I,J);% F(1,I) = DQ(I,J)
end %70 Continue
%***********************************************************
%调用子例行程序INTEGR(N, F, FY)
[FY]=INTEGR(N,F);
%**********************************************************
CP(J) = FY*2.0*3.14159; %CP(J) = FY*2.0*3.14159
for I = 1:N+1 %￥71 Do 71 I = 0,N (内部小循环)
F(2,I) = DT(I,J);%F(1,I) = DT(I,J)
end %71 Continue
%************************************************
[FY]=INTEGR(N, F);
CT(J) = FY;%CT(J) = FY
%**********************************************************
if ((AF~=1)&(CP(J)~=0.0))%If ((AF.ne.1).and.(CP(J).ne.0.0)) Then
ETA(J) = CT(J)*Jar(J)/CP(J);%ETA(J) = CT(J)*Jar(J)/CP(J)
end %Endif
if (ETA(J)<0.0)% IF (ETA(J).LT.0.0) THEN
fprintf('EFFICIENCY NEGATIVE\n');%c WRITE(IMON,*)'EFFICIENCY NEGATIVE'
end%ENDIF
end%61 Continue 大循环61
fprintf(' THRUST, POWER, EFFICIENCY, AND VELOCITIES \n')%c Write(L,*) ' THRUST, POWER, EFFICIENCY, AND VELOCITIES '
Jar; %c Write(L,211) ( Jar(J), J = 0,10 )
CT;%c Write(L,211) ( Ct(J), J = 0,10 )
CP;%c Write(L,211) ( Cp(J), J = 0,10 )
ETA;%c Write(L,211) ( eta(J), J = 0,10 )
if (LDA>0)%c if (LDA.gt.0) then
528.0*MPH./(Jar*Radius);%c write(L,211) (528.0*MPH/(Jar(J)*Radius), J=0,10)
else %c else if (LDA.lt.0) then
Jar*Radius*RPM/528.0;%c write(L,211) (Jar(J)*Radius*RPM/528.0, J=0,10)
end%c end if
if (Fzero==1)%c If (Fzero.eq.1) Then
fprintf('*Thickness values limited by available Cl and Cd data for the selected airfoil section.\n');%c Write(L,*) 'Thickness values limited by available Cl and Cd data'
%c Write(L,*) ' for the selected airfoil section.'
end%c Endif
if ((F1==1)|(F2==1))%c If ((F1.eq.1).or.(F2.eq.1)) Then
fprintf('* Advance ratio values limited by available Cl and Cd data for the selected airfoil section.\n');%c Write(L,*) '* Advance ratio values limited by available Cl and Cd'
%c Write(L,*) ' data for the selected airfoil section. '
end%c Endif
if (F5==1)%c If (F5.eq.1) Then
fprintf('* Warning: Not all values converged within program tolerance in loop to find angle corrections\n');%c Write(L,*) '* Warning: Not all values converged within program'
%c Write(L,*) ' tolerance in loop to find angle corrections.'
end%c End If
%99 Continue
%c CLOSE(UNIT=3)
%200 FORMAT(11F7.3)
%201 FORMAT(F12.2, 9F7.3)
%202 FORMAT(11F7.3)
%203 FORMAT(I3, 4f10.3)
%204 FORMAT(I3,f7.3)
%205 FORMAT(F5.2,3F7.3)
%206 FORMAT(I3, 2f7.3)
%207 FORMAT(10f7.3)
%208 FORMAT(I3, 6f9.3)
%209 FORMAT(A19,9f7.2)
%210 FORMAT(A19,9f7.3)
%c211 FORMAT(A4,11F7.3)
%211 format(11f7.3)
%212 FORMAT(A4,11f7.0)
%213 FORMAT(A8, f5.2, 10f7.2)
%214 FORMAT(f4.2, f9.2, 10f7.2)
%215 FORMAT(f4.2, f9.3, 10f7.3)
%216 FORMAT(f4.2, e8.0, 10e8.0)
%C
% return
%9999 WRITE(IMON,*)'CURVE FIT FAILED - CHECK YOUR DATA!'
% STOP
% END⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -