MATLAB CODE FOR TURBOFAN ANALYSIS
clc clear all % Input parameters Ta=220; %ambient temperature (K) M=0.85;% M=0.85; % Flight Mach number To4= 1700; % Turbine Inlet temperature (K) Q_R= 42480*1000; % Lower heating value (KJ/kg) Po8_o2= 1.5; % fan pressure ratio r=1.4; c_p=1005; u_i=M*sqrt(r*287*Ta); %===================================================================== % Efficiencies of all the components eta_d =0.9; eta_f =0.9; eta_c =0.9; eta_cc =0.9; eta_t =0.9; eta_n =0.9; %===================================================================== % Diffuser To2 = Ta*(1+((r-1)/2)*M^2); Po2_a= (1+(eta_d*((To2/Ta)-1)))^(r/(r-1)); % Fan To8=To2*(1+((1/eta_f)*((Po8_o2^((r-1)/r)-1)))); % fan outlet temperature % Bypass duct nozzle Po8_a=Po8_o2*Po2_a; % Po8/Pa Po8_c=(1-((1/eta_n)*((r-1)/(r+1))))^(r/(1-r)); %Po8/Pc crictical pr.ratio if Po8_a if Po8_a <= Po8_c %Condition for choking u_ef=sqrt(2*eta_n*(r/(r-1))*287*To8*((1-(1/Po8_a))^((r-1)/r))); else Tcf=To8*(2/(r+1)); u_ef=sqrt(r*287*Tcf); F=(287*Tcf*Po8_c/u_ef)*((1/Po8_c)-(1/Po8_a)); end %===================================================================== Po4_o3=eta_cc; %combustion chamber pr.ratio for B=5:10 for B=5:10 for rp=2:100 for rp=2:100 % Compressor outlet temperature To3(length(B),rp)= To8*((rp^((r-1)/r)-1)/eta_c+1); % Combustion Chamber fuel-air ratio f(length(B),rp)=(To4-To3(length(B),rp))/((Q_R/c_p)-To4); % Turbine outlet temperature To5(B,rp)=To4-(To3(length(B),rp)-To8)-(B*(To8-To2)); % Pressure ratio across turbine Po5_o4(B,rp)=(1-((1-(To5(B,rp)/To4))*(1/eta_t)))^(r/(r-1)); % Nozzle Po6_a(B,rp)=Po2_a*rp*Po4_o3*Po5_o4(B,rp)*Po8_o2; %Po6/Pa Po6_c=(1-((1/eta_n)*((r-1)/(r+1))))^(r/(1-r)); % Po6/Pc %============================================================================ %Unchoked nozzle if(Po6_a(B,rp) if (Po6_a(B,rp) <= Po6_c) %Condition for choking
u_e(B,rp)=sqrt(2*eta_n*c_p*To5(B,rp)*(1-(1/Po6_a(B,rp)^((r1)/r)))); F_s(B,rp)=(1+f(length(B),rp))*u_e(B,rp)+B*u_ef-(1+B)*u_i+F; TSFC(B,rp)=f(length(B),rp)/F_s(B,rp); %============================================================================ % Choked nozzle else Tc(B,rp)=To5(B,rp)*(2/(r+1)); u_e(B,rp)=sqrt(r*287*Tc(B,rp)); F_sp(B,rp)= (287*Tc(B,rp)*Po6_c/u_e(B,rp))*((1/Po6_c)(1/Po6_a(B,rp)))+F; F_s(B,rp)=(1+f(length(B),rp))*u_e(B,rp)+B*u_ef(1+B)*u_i+F_sp(B,rp); TSFC(B,rp)=f(length(B),rp)/F_s(B,rp); end %============================================================================ vo(B,rp)=((1+f(length(B),rp))*u_e(B,rp)^2)-u_i^2+(B*(u_ef^2-u_i^2)); eta(B,rp)= 2*F_s(B,rp)*u_i/vo(B,rp); %Propulsive efficiency etat(B,rp)=vo(B,rp)/(2*f(length(B),rp)*Q_R); %Thermal efficiency etao(B,rp)=eta(B,rp)*etat(B,rp); %Overall efficiency end end for B=5:10 for B=5:10 plot ((1:rp),TSFC(B,:)); hold on plot ((1:rp),TSFC(B,:)); hold on plot ((30:45),eta(B,30:45)); hold on plot ((30:45),etat(B,30:45)) hold on plot ((30:45),etao(B,30:45)) end
