Matlab Code into Simulink Help

Question

ameen 2015-5-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/215757-matlab-code-into-simulink-help

编辑： Walter Roberson 2015-5-8

I want to put the following code into Simulink. I tried MATLAB Function block but it gives me many errors. May be a S-Function block will help. I want to have a block with 4 inputs (J, P_D, BARatio, NosBlades) and to have 3 outputs (TT, TQ, TE). Here is the original code

J=0.4;
P_D=0.875; 
BARatio=0.55;  
NosBlades=4;    
%%%%%%%%%%
C2(2)=0.208;
C2(3)=0.241;
C2(4)=0.263;
C2(5)=0.276;
C2(6)=0.279;
C2(7)=0.269;
C2(8)=0.241;
C2(9)=0.184;
%%%%%%%%%
PA(2)=0.888;
PA(3)=1.008;
PA(4)=1.055;
PA(5)=1.060;
PA(6)=1.039;
PA(7)=1.000;
PA(8)=0.948;
PA(9)=0.888;
%%%%%%%%
VV=0;
LA=0.97;
RT=0.045;
LM=0.8*11.3;
LI=0.8*12.67;
%%%%%%%%%
for i=2:9
      C1(i)=C2(i)*BARatio*4.0/(NosBlades*0.5);
      CP=NosBlades*C1(i);                         %Nc/D
      X(i)=i/10;                                  %x/R
      AI(i)=0;                                    %Assume initial Angle of Attack=0
      PR(i)=P_D*PA(i);                            %Pitch/diameter ratio at blade element    
    TD=RT-(RT*0.935)*X(i);
    TC(i)=TD/C1(i);    
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      counter=1;
      SubConverged=0;
      while and(counter <2000, SubConverged==0)
          counter=counter+1;
          HA=atan(PR(i)/(pi*X(i)))/pi*180-AI(i);  
          H1=tan(HA/180*pi);                      %tan phi
          AA=J/(pi*X(i));                         %tan psi (undisturbed flow angle?)
          EI=AA/H1;                               %Ideal Efficiency
          EA(i)=0.9*EI;                           %Local Efficiency Estimate
          AE=0.0;                                 %tan(gamma) =atan(Cd(i)/Cl(i)) assume =0.0
          KF=X(i)*H1;                             %Lambda parameter for Goldstein correction
          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
          SF=NosBlades/(2.0*KF)-0.5;
          F1=cosh(SF);
          F2=cosh(SF*X(i));
          F3=F2/F1;
          F4=acos(F3);
          K=2.0*F4/pi;
          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
          GK(i)=K;                                
          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
          counter2=0;
          SubSubConverged=0;
          while and(counter2<2000, SubSubConverged==0)
              counter2=counter2+1;
              FA(i)=(1-EI)/(EI+AA*AA/EA(i));                      %axial flow factor
              KT(i)=pi*J^2*X(i)*K*FA(i)*(1+FA(i));                %local dKT/dx
              AE=tan(AE/180*pi);                                  %tan(gamma)
              FT(i)=1-EI*(1+FA(i));                               %a' tangential flow factor
              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
              CL(i)=KT(i)*cos(HA/180*pi);
              CL(i)=CL(i)*4/(pi^2)/(X(i)^2);
              CL(i)=CL(i)/((1-FT(i))^2*(1-H1*AE)*CP);
              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
              F6=0.0107+(AI(i)+1.0)*(-0.0015+AI(i)*(0.0015+0.000965*(AI(i)-1.0)));
              F7=-0.03833+(AI(i)+1.0)*(0.0133+AI(i)*(-0.015-0.01166*(AI(i)-1.0)));
              F8=0.8193+(AI(i)+1.0)*(-0.0138+AI(i)*(0.0903+0.079*(AI(i)-1.0)));
              F9=-3.076+(AI(i)+1.0)*(-0.0728+AI(i)*(-0.3162-0.2437*(AI(i)-1.0)));
              CD(i)=F6+TC(i)*(F7+(TC(i)-0.06)*(F8+F9*(TC(i)-0.12)));
              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
              AG=CD(i)/CL(i);                                     
              AE=atan(AG)/pi*180;                                 %gamma
              EZ=AA/tan((HA+AE)/180*pi);                          %
              if (abs(EZ-EA(i)) < 0.001)                          
                  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                  K2(2)=1.0+2.857*(BARatio-0.4)*(BARatio-0.6)*(BARatio-0.9);
                  K2(3)=1.19+(BARatio-0.4)*(BARatio-0.6)*(0.267+(BARatio-0.9)*0.1665);
                K2(4)=1.3+(BARatio-0.4)*(0.1+(BARatio-0.6)*(1+(BARatio-0.9)*0.43));
                K2(5)=1.54+(BARatio-0.4)*(0.15+(BARatio-0.6)*(1.1+0.286*(BARatio-0.9)));
                K2(6)=1.67+(BARatio-0.4)*(0.55+(BARatio-0.6)*(0.1667+(BARatio-0.9)*2.9465));
                K2(7)=1.8+(BARatio-0.4)*(0.75+(BARatio-0.6)*(0.3+(BARatio-0.9)*2.835));
                K2(8)=1.8+(BARatio-0.4)*(1.0+(BARatio-0.6)*(1.333+(BARatio-0.9)*1.905));
                K2(9)=1.75+(BARatio-0.4)*(1.25+(BARatio-0.6)*(1.5+(BARatio-0.9)*3.55));
              U1=-0.65*KF*KF+1.1*KF+0.664;
              U2=(0.85+(KF-0.3)*(-4.0+(KF-0.4)*(15.42-47.95*(KF-0.5))));
              U2=-0.09+(KF-0.2)*U2;
              U3=(1.375+(KF-0.3)*(-3.75+(KF-0.4)*(20.85-75.7875*(KF-0.5))));
              U3=-0.2+(KF-0.2)*U3;
              K1=U1+(BARatio-0.4)*(U2+U3*(BARatio-0.8));
              CC(i)=K1*K2(i);
                  if VV==1
                      MC(i)=MT(i)*TC(i);
                      AC=(CC(i)*CL(i)/LI-MC(i))*LM/0.1097+(MC(i)*LI/LA);
                  elseif VV==5
                      MC(i)=0.5*TC(i);
                      AC=(CC(i)*CL(i)/LI-MC(i))*LM/0.1097+(CL(i)/LA);
                      MT(i)=MC(i)/TC(i);
                  else
                      AC=CL(i)/LA;
                      MC(i)=CL(i)/LI;
                      MC(i)=MC(i)*CC(i);
                      if (MC(i)<0.5*TC(i))
                          MT(i)=MC(i)/TC(i);
                      else
                          VV=5;
                          MC(i)=0.5*TC(i);
                          AC=(CC(i)*CL(i)/LI-MC(i))*LM/0.1097+(CL(i)/LA);
                          MT(i)=MC(i)/TC(i);
                      end              
                  end
                 SubSubConverged=1;                                %Set subsub loop converged flag to 1
              else
                  EA(i)=EZ;                                        %update local efficiency estimate
              end          
          end %subsubloop
          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
          if (abs((AC-AI(i))/(AC))<0.1) 
                             %subloop converged difference between angle of attack assumed AI and calculated angle of attack AC <0.1*AC
              KQ(i)=4.935*J*X(i)*X(i)*X(i)*K*FT(i)*(1+FA(i));      %calculate dkq/dx
              SubConverged=1;                                      %set subloop converged flag to 1
          else
              AI(i)=(AC+AI(i))/2;                                  %New Guess of Angle of Attack AI(i)
          end
      end %subloop
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  end     %End Main Loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X(10)=1.0;
KT(10)=0.0;
KQ(10)=0.0;
FA(10)=0.0;
FT(10)=0.0;
GK(10)=0.0;
MT(10)=MT(9);
PR(10)=PR(9);
AI(10)=AI(9);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Numerical Integration of Total Values
TT=(KT(2)+2.0*(KT(4)+KT(6)+KT(8))+4.0*(KT(3)+KT(5)+KT(7)+KT(9)));   
TT=0.1*TT/3.0;                                                      
TQ=(KQ(2)+2.0*(KQ(4)+KQ(6)+KQ(8))+4.0*(KQ(3)+KQ(5)+KQ(7)+KQ(9)));   
TQ=0.1*TQ/3.0;                                                      
TE=J*TT/(2.0*3.1416*TQ);