clcdcalc.m

用MATLAB编写关于飞机设计的matlab程序

      function  [II,KK,CL,CD,RE,MT,SF,CF]=CLCDCALC(PHI,T,Z,ALPHA,I,J,rows,MPH,LDA,C,NU,Acousvel,Radius,SF,CF,Azero,X,RE,MT,MTUC,RPM,REAFD,ACL,BCD)
   
                     %C----T(0:9)相对厚度，Z(0:1,0:25,0:25)翼型参数后改为T(1:10) Z(1:2,1:26,1:26)
      %              Real PHI, CF, SF, T(0:9), Z(0:1,0:25,0:25), ALPHA(0:25,0:25)
      %              Real RD, D1, D2,CL, CD,MPH,C1,RE(0:25,0:25),C(0:9), NU,MT(0:10),MF
      %              Real Acousvel, Azero, P3, X(0:9), MTUC(0:10),ALPHADEG(0:25,0:25)
      %              Integer J, rows, I, II, KK, LDA
      %              Common /logunits/imon,ikey
%
     


      C1 = 360.0;%C1 = 360.0
      P3 = 1000.0;%P3 = 1000.0
      RD = 57.29578;%RD = 57.29578
      MF = 1.46667;%MF = 1.46667
      D1 = 0.0;%D1 = 0.0
      D2 = 0.0;%D2 = 0.0
      II = 2;%II = 1
      KK = 2;%KK = 1
      CF = cos(PHI);%CF = COS(PHI)
      SF = sin(PHI);%SF = SIN(PHI)
      while (T(I)>Z(1,1,KK+1))%1     If (T(I).gt.Z(0,0,KK+1)) Then
             KK = KK+1;       %KK = KK+1
      end                     %Goto 1（找到相对厚度区间下限）
                              %Endif
      if ((ALPHA(I,J)*RD)>Z(1,rows,1))                        %If ((ALPHA(I,J)*RD).gt.Z(0,rows,0)) Then
          II = rows-1; %2     If ((ALPHA(I,J)*RD).gt.Z(0,II+1,0)) Then
                       %II = II +1
                       %Goto 2（找到迎角区间下限）
      else                   %Endif
      %（既迎角大于表中的迎角极限）
           while ((ALPHA(I,J)*RD)>Z(1,II+1,1))
              II = II +1; %II = rows - 1（迎角极限既迎角区间下限）它的rows从0开始
           end 
      end                %Goto 30
                         %Endif
%2     If ((ALPHA(I,J)*RD).gt.Z(0,II+1,0)) Then
%        II = II +1
%        Goto 2
%      Endif
%30    Continue
      D1 = (T(I) - Z(1,1,KK))/(Z(1,1,KK+1)-Z(1,1,KK));%D1 = (T(I) - Z(0,0,KK))/(Z(0,0,KK+1)-Z(0,0,KK))
      D2 = (ALPHA(I,J)*RD - Z(1,II,1))/(Z(1,II+1,1) - Z(1,II,1));%D2 = (ALPHA(I,J)*RD - Z(0,II,0))/(Z(0,II+1,0) - Z(0,II,0))
      CL = Z(1,II,KK) + D1*(Z(1,II,KK+1) - Z(1,II,KK)) +D2*(Z(1,II+1,KK) - Z(1,II,KK)); %CL = Z(0,II,KK) + D1*(Z(0,II,KK+1) - Z(0,II,KK)) + 
                                                                                        %*     D2*(Z(0,II+1,KK) - Z(0,II,KK))
      CD = Z(2,II,KK) + D1*(Z(2,II,KK+1) - Z(2,II,KK)) +D2*(Z(2,II+1,KK) - Z(2,II,KK)); %CD = Z(1,II,KK) + D1*(Z(1,II,KK+1) - Z(1,II,KK)) + 
                                                                                        %*     D2*(Z(1,II+1,KK) - Z(1,II,KK))
      D1 = MPH*MF/SF;%D1 = MPH*MF/SF （得到气流速度）
      if (LDA<0) %If (LDA.lt.0) Then
          D1= ((RPM/60.0)*3.14159*2.0*Radius*X(I))/CF/12.0;%  D1= ((RPM/60.)*3.14159*2.*Radius*X(I))/CF/12.0
      end        %Endif
      RE(I,J) = D1*C(I)/NU/12.0;%RE(I,J) = D1*C(I)/NU/12.0
      if (I==9)%If (I.eq.8) Then（桨尖）
        MT(J) = (D1/X(10))/Acousvel;%MT(J) = (D1/X(9))/Acousvel（得到桨尖马赫数）
      end%Endif
      if (abs(LDA)==1)%If (ABS(LDA).eq.1) Then
          %RETURN
          return    
      end%Endif
      if (abs(LDA)~=3)%If (ABS(LDA).ne.3) Then
          D2 = sqrt(1.0 - D1*D1/(Acousvel*Acousvel));%D2 = SQRT(1.0 - D1*D1/(Acousvel*Acousvel))
          CL = CL/D2;%CL = CL/D2
          CD = CD/D2;%CD = CD/D2
      end%Endif
      if(abs(LDA)==2);%If(ABS(LDA).eq.2) Then
         %RETURN
         return
       end             %End
%
%************** The Reynolds Nun\mber Coefficient Correction  *************** 
% 
       R0 = Z(1,2,KK);%R0 = Z(0,1,KK)KK列第一个升力系数
       R1 = Z(1,2,1);%R1 = Z(0,1,0)最小角度
       R2 = R0;%R2 = R0
       R3 = R1;%R3 = R1
       R4 = Z(2,2,KK);%R4 = Z(1,1,KK)KK列第一个阻力系数
       for JJ = 2:rows%Do 5 JJ = 1, rows
            if (R0>Z(1,JJ,KK))%If (R0.gt.Z(0,JJ,KK)) Then
                R0 = Z(1,JJ,KK);%R0 = Z(0,JJ,KK)
                R1 = Z(1,JJ,1);%R1 = Z(0,JJ,0)
            end%Endif找到KK列的最小升力系数和此时的角度
            if (R2<Z(1,JJ,KK))%If (R2.lt.Z(0,JJ,KK)) Then
                R2 = Z(1,JJ,KK);%R2 = Z(0,JJ,KK)
                R3 = Z(1,JJ,1);%R3 = Z(0,JJ,0)
            end%Endif找到KK列的最大升力系数和此时的角度
            if (R4>Z(2,JJ,KK))%If (R4.gt.Z(1,JJ,KK)) Then
                R4 = Z(2,JJ,KK);%R4 = Z(1,JJ,KK)
            end%Endif找到KK列的最最小阻力系数
       end%5     Continue
       Azero = (R2 - R0)/(R3 - R1); %Azero = (R2 - R0)/(R3 - R1)     Azero即为升力线斜率
%C
%C	NOTE: THIS PORTION OF THE CODE IS NOT COMPATABLE
%C		WITH OLDER DATA FILES
%C
      R5 = REAFD;%R5 = REAFD
      R6 = ACL;%R6 = ACL
      R7 = BCD;%R7 = BCD	
      R8 = R2*((RE(I,J)/R5)^R6 - 1.0);%R8 = R2*((RE(I,J)/R5)**R6 - 1.0)
      R9 = R4 * ((R5/RE(I,J))^R7 - 1.0);%R9 = R4 * ((R5/RE(I,J))**R7 - 1.0)
      CD = CD + R9;%CD = CD + R9
      ALPHADEG(I,J)=ALPHA(I,J)*RD;%ALPHADEG(I,J)=ALPHA(I,J)*RD
%C
%C	NEW CL CORRECTION ROUTINE
%C
      if (ALPHADEG(I,J)>R3)%If (ALPHADEG(I,J).GT.R3) Then
          CL=CL+R8;%CL=CL+R8
          return%RETURN
      end%Endif
      if (ALPHADEG(I,J)>=(R3+(R8/Azero)))%If (ALPHADEG(I,J).GE.(R3+(R8/AZERO))) Then
          CL=CL+R8*((ALPHADEG(I,J)-(R3+(R8/Azero)))/(R3-(R3+(R8/Azero))));%CL=CL+R8*((ALPHADEG(I,J)-(R3+(R8/AZERO)))/
                                                                          %&   (R3-(R3+(R8/AZERO))))
          return%RETURN
      end %Endif
      return%RETURN
      %End
%C