4 Stage Axial Flow Compressor Design using MATLAB
The main objective of this project is to design a 4-stage- Axial compressor with inlet guide vane, having constant tip diameter and using free vortex distribution i.e. swirling flow with constant angular momentum. Also, for the same volume flow rate and inlet conditions, results for the two different fluids- Air and Hydrogen has to be compared.
Compressor has been designed using MATLAB code, which uses just a handful of parameters such as ambient conditions, and based on some assumptions calculates flow conditions at various nodes from hub to tip and inlet to outlet of the wheel throughout the compressor.
Based on the thermodynamic laws, calculation is carried out along each streamline.Based on the required Pressure ratio and volume flow rate set by user, calculations are first carried out iteratively along mean-line. Then using Free vortex law, flow variables are calculated along the blade from hub to tip. Followed by circumferential velocity profile assumption, using which flow properties are calculated from inlet to outlet part of the wheel.
Based on the blade angles calculated at each point and basic profile assumption, blade shape along pressure side is calculated. Then using assumed Thickness distribution, suction surfaces are plotted from inlet to outlet. Based on the blade profile, boundary layer thickness is calculated at pressure and suction side of the blades.
Various Limitation parameters are part of the code, which warns the user, about the Mach numbers which have exceeded a safe value. If such limiting parameters are not checked then, in practical scenario, there could be shocks or flow separation that might drastically reduce the efficiency of compressor.All the Mach numbers are limited to be less than 1, as supersonic flows can lead to shocks in the compressor. Higher the Mach number, higher will be the probability of shock formation. Also to avoid diffusion, De-Haller number is limited to 0.72.
Project Plan:
Design Approach: Mean Line Design Approach
Results:
Meridional Sectional Veiw designed based on matalb code after final calculations for both fluids:
Mach number contour plot for both fluids, to check for Mach number limitations avoiding supersonic conditions, as blades are designed based on subsonic design criteria:
Velocity plots on high and low pressure sides of the blade:
3-D representation of the rotor blade using matlab code:
MATLAB CODE for the same is given below:
clear all
close all
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Given parameters%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Dimensional
R=287; %specific gas constant used watervapour--%467 air---%287 H2--4124
kvisc=1.48e-5;%the viscosity of air is 1.48×10?5 m2/s at given temp.
D0=0.125;%(m)
rt=D0/2;%radius of tip-constant for all calculations -as COD design method is use
%Thermal parameters
cp=1.005*10^3; % specific heat at constant pressure taken constant for all calculations %water vapour--1.93e3;;;air---1.005*10^3.....H2-14.31e3
gamma=1.41; %gamma-specific heat ratios %water vapour --1.33; air---1.4; H2-1.405
gammas=2*gamma/(gamma+1);
T01=283.15; % static temperature before igv inlet in kelvin
P01=1.2281*10^3; %inlet static pressure before igv in pa
H01=cp*T01; %static enthalpy at igv inlet
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%chosen parameters%%%%%%%%%%%%%%%%%%%%%%%%%%
N=[48000 53000 56000 59000];%(rpm) %%optimized answers for air--[48000 53000 56000 59000];;; 2nd case H2--[96000 106000 112000 118000]; H2_best-[84000 84000 84000 84000]
Rht=[0.65 0.6565]; %ratio of hub to tip ratio at inlet
set_VFlow=1.5; %m3/s set valued--1.5
Rm=1.02; %rms slope line
alpha5=13;%optimized answers for air--13;11
alpha9=11;%optimized answers for air--11;11
alpha13=11;%optimized answers for air--11;11
alpha17=0;%optimized answers for air--0;0
AR_stator=1.55; %aspect ratio for stator 2.25
AR_rotor=1.25; %aspect ratio for rotor 2
AR_IGV=1.55;%aspect ratio for rotor 2.25
alpha1=18;% initial guess for degreees turning angle given by stator to the flow %optimized answers for air--18
phi=0.7; %flow coefficient %0.5 for increasing pressure ratio, phi is increased %optimized answers for air--0.7; 2nd case H2-0.35; --0.4
PSI=[0.5 0.505 0.5 0.505 0.5 0.505 0.5 0.505] ; % initial guess forLoading coefficient
Cm1=100; %initial guess for meridional velocity
PR_comp=[1.1 1.095 1.09 1.09]; %optimized answers for air -[1.49 1.45 1.4 1.39]; H2--[1.025 1.02 1.02 1.02];;; 2nd H2--[1.035 1.04 1.045 1.055];;;;;; best H2--[1.1 1.095 1.09 1.09]
%%%%%%%%%%%%%%%%%%%%%%%%%%%Mean-Line-Design method%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%% Pressure loop Iteration %%%%%%%%%%%%%%%%%%%%%%%%
n_PR=0; %setting no. of steps
PR_rel_error_level=1;
conv_PR=0;
while conv_PR ==0
if abs(PR_rel_error_level)<1e-4
conv_PR=1;
end
for m =1:2
%%%%%%%%%%%%Calculation for inlet profile/Inlet Guide vane %%%%%%%%%%%%%%%%%%%%%%%%%%%
%Iteration for Cm0
Cm1_old=Cm1;
rt1=rt;
Cm1_error=1;
damp=0.1;
itr_cm=0;
%%%%%% Iteration to correct Volume flow rate to set volume flow rate%%%%%%
%%%%%%%%Hub-tip ratio is varied to correct to the set volume flow rate%%%%
conv_VF=0;
VF_rel_error_level=1;
while conv_VF ==0
if abs(VF_rel_error_level)<1e-4
conv_VF=1;
end
for j=1:2
rh1=rt1*Rht(j);
while(Cm1_error>1e-2) % correcting Axial velocity to set flow coefficient
Cm1_old=Cm1;
C1=Cm1_old/cosd(alpha1);
rmean1=((rh1^2+rt^2)/2)^0.5;
rmean2=rmean1*Rm; %rms radius at Rotor
umean2=2*pi*rmean2*N(1)/60;%rotation velocity at Rotor
Cm1=umean2*phi;
Cm1_error=abs((Cm1/Cm1_old)-1);
Cm1=Cm1+(Cm1_old-Cm1)*damp;
C_circum1=C1*sind(alpha1);
itr_cm=itr_cm+1;
end
%Calculated variables for IGV outlet
H1=H01-(C1^2)/2;
T1=H1/cp;
a1=(gamma*R*T1)^0.5;
a1s=(gammas*R*T01)^0.5;% a1*
M1=C1/a1;
Mc1s=C1/a1s;% Mc1* for IGV outlet
P1=P01/((1+((gamma-1)/2)*M1^2)^(gamma/(gamma-1)));
Dens0=P01/(R*T01);
Dens1=Dens0/((1+((gamma-1)/2)*M1^2)^(1/(gamma-1)));
Area1=pi*(rt1^2-rh1^2);
massflow=Dens1*Area1*Cm1;
Volume_flow(j)=massflow/Dens1; %Volume flow rate--m3/s
end
if conv_VF==1
break;
end
if conv_VF==0
VFlow_slope=(Volume_flow(1)-Volume_flow(2))/(Rht(1)-Rht(2));
VF_rel_error_level=(set_VFlow/Volume_flow(1))-1;
VF_error_level=set_VFlow-Volume_flow(1);
Rht(1)=Rht(1)+(VF_error_level/VFlow_slope)*0.01;
Rht(2)=Rht(1)*1.01;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Stage1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%Calculation for Rotor Inlet at mean%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rt2=rt1;
Cm2=Cm1; %meridional velocity kept constant from stator to Rotor
C_circum2=C_circum1/Rm;
alpha2=atand(C_circum2/Cm2);
C2=Cm2/cosd(alpha2); %absolute velocity calcualted
W_circum2=(umean2-C_circum2); %relative circum velocity
Beta2=atand(W_circum2/Cm2); %blade angle aka relative flow angle
W2=Cm2/cosd(Beta2); %relative velocity at blade inlet
%Stagnation properties at inlet remains constant from stator out to rotor
%in
H02=H01;
P02=P01;
T02=T01;
Dens02=Dens0;
%Static properties at inlet
H2=H02-C2^2/2; %Static heat enthalpy
T2=H2/cp; %Temperature at inlet
a2=(gamma*R*T2)^0.5; % sonic velocity at inlet
a2s=(gammas*R*T02)^0.5; %critical sonic velocity
T02r=T02+(W2^2-C2^2)/(2*cp); %calculated for relative critical sonic velocity
a2sr=(gammas*R*T02r)^0.5;% relative critical sonic velocity
M2=C2/a2; %Mach number calculated at inlet
Mc2s=C2/a2s; %critical mach no. at inlet
Mw2s=W2/a2sr; %critical relative mach no. at inlet
P2=P02/((1+((gamma-1)/2)*M2^2)^(gamma/(gamma-1))); % static pressure at inlet
Dens2=Dens02/((1+((gamma-1)/2)*M2^2)^(1/(gamma-1)));
I2=H2+W2^2/2-umean2^2/2; % inlet rothalpy
%Dimensions for rotor inlet
area2=massflow/(Dens2*Cm2); %Rotor inlet area
rh2=((pi*rt^2-area2)/pi)^0.5; % Hub unlet radius
%%%%%%%%%%%%%%%%%%%%%Rotor Oulet Calculations%%%%%%%%%%%%%%%%%%%%%%%%
rt3=rt2; %tip radius is constant inlet to outlet
Cm3=Cm2; %meridional velocity is constant inlet to outlet
umean3=umean2; % rotational velocity constant inlet to outlet
C_circum3=umean3*PSI(m)+C_circum2; %circumferential velocity calculated using PSI at outlet
C3=(Cm3^2+C_circum3^2)^0.5; % absolute velocity calculated at outlet
alpha3=atand(C_circum3/Cm3); %absolute angle calculated at outlet
W_circum3=umean3-C_circum3; % relative circumferential velocity calc at oulet
Beta3=atand(W_circum3/Cm3); %blade angle calculated at oulet
W3=Cm3/cosd(Beta3); %relative velocity calc. at outlet
%static & total thermal properties
I3=I2; %Rothalpy remains constant
H3=I3-W3^2/2+umean3^2/2; %static enthalpy calculated using Rothalpy
T3=H3/cp; % Temperature at rotor outlet
H03=H3+C3^2/2; % stagnation enthalpy at rotor outlet
a3=(gamma*R*T3)^0.5; % sonic velocity at rotor outlet
M3=C3/a3; % Mach no. at rotor outlet
T03=T3*(1+((gamma-1)*M3^2/2)); % total temp at rotor outlet
a3s=(gammas*R*T03)^0.5;
Mc3s=C3/a3s;
Mu3s=umean3/a3s;
Iso_effi=0.8; % isoentropic efficiency chosen
P03=P02*(1+Iso_effi*(T03-T02)/T02)^(gamma/(gamma-1)); % total pressure calc. using assumed efficiency
Dens03=P03/(R*T03);
P3=P03/(1+((gamma-1)*M3^2/2))^(gamma/(gamma-1)); % Static pressure calc. at rotor outlet
Dens3=Dens03/((1+((gamma-1)/2)*M3^2)^(1/(gamma-1)));
%Rotor Oulet dimensions
area3=massflow/(Dens3*Cm3);
rh3=((pi*rt3^2-area3)/pi)^0.5;
dH=W3/W2;
rmean3=rmean2;
%%%%%%%%%%%%%%%%%Stator Inlet Calculations%%%%%%%%%%%%%%%%%%%%%%%%%%
rt4=rt;
Cm4=Cm3;
C_circum4=C_circum3/Rm;
alpha4=atand(C_circum4/Cm4);
C4=Cm4/cosd(alpha4);
%Stagnation enthalpies
H04=H03;
P04=P03;
T04=T03;
Dens04=Dens03;
%Static enthalpies
H4=H04-C4^2/2;
T4=H4/cp;
a4=(gamma*R*T4)^0.5;
a4s=(gammas*R*T04)^0.5;
M4=C4/a4;
Mc4s=C4/a4s;
P4=P04/(1+((gamma-1)*M4^2/2))^(gamma/(gamma-1));
Dens4=Dens04/((1+((gamma-1)/2)*M4^2)^(1/(gamma-1)));
%Stator dmensions
area4=massflow/(Dens4*Cm4);
rh4=((pi*rt4^2-area4)/pi)^0.5;
rmean4=rmean3*Rm;
%%%%%%%%%%%%%%%%Stator Outlet calculations%%%%%%%%%%%%%%%%%%%%%%%%
%alpha5=alpha2; %for reaction--0.5
Cm5=Cm4;
C5=Cm5/cosd(alpha5);
C_circum5=C5*sind(alpha5);
T05=T04;
H05=H04;
P05=P04;
Dens05=Dens04;
H5=H05-C5^2/2;
T5=H5/cp;
a5=(gamma*R*T5)^0.5;
a5s=(gammas*R*T05)^0.5;
M5=C5/a5;
Mc5s=C5/a5s;
P5=P05/(1+((gamma-1)*M5^2/2))^(gamma/(gamma-1));
Dens5=Dens05/((1+((gamma-1)/2)*M5^2)^(1/(gamma-1)));
%Stator dimensions-outlet
rt5=rt4;
area5=massflow/(Dens5*Cm5);
rh5=((pi*rt5^2-area5)/pi)^0.5;
rmean5=rmean4;
React1=(H3-H2)/(H05-H02); % Stage 1 reaction
PR_stage1=P5/P1;
PR(m) = PR_stage1; %total pressure ratio
if conv_PR==1
break
end
end %end of 2 calculation for pressure loop
if conv_PR==0
PR_slope=(PR(1)-PR(2))/(PSI(1)-PSI(2));
PR_rel_error_level=(PR_comp(1)/PR(1))-1;
PR_error_level=PR_comp(1)-PR(1);
PSI(1)=PSI(1)+(PR_error_level/PR_slope)*0.01;
PSI(2)=PSI(1)*1.01;
n_PR=n_PR+1;
end
end
n_PR=0; %setting no. of steps
PR_rel_error_level2=1;
conv_PR2=0;
while conv_PR2 ==0
if abs(PR_rel_error_level2)<1e-4
conv_PR2=1;
end
for m =3:4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%Stage 2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rotor2 inlet
%%%%%%%%%%%%%%%%%%Calculation for Rotor Inlet2 at mean%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rmean6=rmean5*Rm;
umean6=2*pi*rmean6*N(2)/60;%rotation velocity at Rotor
rt6=rt5;
Cm6=Cm5; %meridional velocity kept constant from stator to Rotor
C_circum6=C_circum5/Rm;
alpha6=atand(C_circum6/Cm6);
C6=Cm6/cosd(alpha6); %absolute velocity calcualted
W_circum6=(umean6-C_circum6); %relative circum velocity
Beta6=atand(W_circum6/Cm6); %blade angle aka relative flow angle
W6=Cm6/cosd(Beta6); %relative velocity at blade inlet
%Stagnation properties at inlet remains constant from stator out to rotor
%in
H06=H05;
P06=P05;
T06=T05;
Dens06=Dens05;
%Static properties at inlet
H6=H06-C6^2/2; %Static heat enthalpy
T6=H6/cp; %Temperature at inlet
a6=(gamma*R*T6)^0.5; % sonic velocity at inlet
a6s=(gammas*R*T06)^0.5; %critical sonic velocity
T06r=T06+(W6^2-C6^2)/(2*cp); %calculated for relative critical sonic velocity
a6sr=(gammas*R*T06r)^0.5;% relative critical sonic velocity
M6=C6/a6; %Mach number calculated at inlet
Mc6s=C6/a6s; %critical mach no. at inlet
Mw6s=W6/a6sr; %critical relative mach no. at inlet
P6=P06/((1+((gamma-1)/2)*M6^2)^(gamma/(gamma-1))); % static pressure at inlet
Dens6=Dens06/((1+((gamma-1)/2)*M6^2)^(1/(gamma-1)));
I6=H6+W6^2/2-umean6^2/2; % inlet rothalpy
%Dimensions for rotor inlet
area6=massflow/(Dens6*Cm6); %Rotor inlet area
rh6=((pi*rt6^2-area6)/pi)^0.5; % Hub unlet radius
%%%%%%%%%%%%%%%%%%%%%Rotor Oulet Calculations%%%%%%%%%%%%%%%%%%%%%%%%
rmean7=rmean6;
%Velocity triangles
rt7=rt6; %tip radius is constant inlet to outlet
Cm7=Cm6; %meridional velocity is constant inlet to outlet
umean7=umean6; % rotational velocity constant inlet to outlet
C_circum7=umean7*PSI(m)+C_circum6; %circumferential velocity calculated using PSI at outlet
C7=(Cm7^2+C_circum7^2)^0.5; % absolute velocity calculated at outlet
alpha7=atand(C_circum7/Cm7); %absolute angle calculated at outlet
W_circum7=umean7-C_circum7; % relative circumferential velocity calc at oulet
Beta7=atand(W_circum7/Cm7); %blade angle calculated at oulet
W7=Cm7/cosd(Beta7); %relative velocity calc. at outlet
%static & total thermal properties
I7=I6; %Rothalpy remains constant
H7=I7-W7^2/2+umean7^2/2; %static enthalpy calculated using Rothalpy
T7=H7/cp; % Temperature at rotor outlet
H07=H7+C7^2/2; % stagnation enthalpy at rotor outlet
a7=(gamma*R*T7)^0.5; % sonic velocity at rotor outlet
M7=C7/a7; % Mach no. at rotor outlet
T07=T7*(1+((gamma-1)*M7^2/2)); % total temp at rotor outlet
a7s=(gammas*R*T07)^0.5;
Mc7s=C7/a7s;
Mu7s=umean7/a7s;
Iso_effi=0.8; % isoentropic efficiency chosen
P07=P06*(1+Iso_effi*(T07-T06)/T06)^(gamma/(gamma-1)); % total pressure calc. using assumed efficiency
Dens07=P07/(R*T07);
P7=P07/(1+((gamma-1)*M7^2/2))^(gamma/(gamma-1)); % Static pressure calc. at rotor outlet
Dens7=Dens07/((1+((gamma-1)/2)*M7^2)^(1/(gamma-1)));
%Rotor Oulet dimensions
area7=massflow/(Dens7*Cm7);
rh7=((pi*rt7^2-area7)/pi)^0.5;
dH2=W7/W6;
%%%%%%%%%%%%%%%%%Stator2 Inlet Calculations%%%%%%%%%%%%%%%%%%%%%%%%%%
rt8=rt7;
Cm8=Cm7;
C_circum8=C_circum7/Rm;
alpha8=atand(C_circum8/Cm8);
C8=Cm8/cosd(alpha8);
%Stagnation enthalpies
H08=H07;
P08=P07;
T08=T07;
Dens08=Dens07;
%Static enthalpies
H8=H08-C8^2/2;
T8=H8/cp;
a8=(gamma*R*T8)^0.5;
a8s=(gammas*R*T08)^0.5;
M8=C8/a8;
Mc8s=C8/a8s;
P8=P08/(1+((gamma-1)*M8^2/2))^(gamma/(gamma-1));
Dens8=Dens08/((1+((gamma-1)/2)*M8^2)^(1/(gamma-1)));
%Stator dmensions
area8=massflow/(Dens8*Cm4);
rh8=((pi*rt8^2-area8)/pi)^0.5;
rmean8=rmean7*Rm;
%%%%%%%%%%%%%%%%Stator2 Outlet calculations%%%%%%%%%%%%%%%%%%%%%%%%
%alpha9=alpha6;
Cm9=Cm8;
C9=Cm9/cosd(alpha9);
C_circum9=C9*sind(alpha9);
T09=T08;
H09=H08;
P09=P08;
Dens09=Dens08;
H9=H09-C9^2/2;
T9=H9/cp;
a9=(gamma*R*T9)^0.5;
a9s=(gammas*R*T09)^0.5;
M9=C9/a9;
Mc9s=C9/a9s;
P9=P09/(1+((gamma-1)*M9^2/2))^(gamma/(gamma-1));
Dens9=Dens09/((1+((gamma-1)/2)*M9^2)^(1/(gamma-1)));
%Stator dimensions-outlet
rt9=rt8;
area9=massflow/(Dens9*Cm9);
rh9=((pi*rt9^2-area9)/pi)^0.5;
rmean9=rmean8;
React2=(H7-H6)/(H09-H06);% Stage 2 reaction
PR_stage2=P9/P6;
PR(m) = PR_stage2; %total pressure ratio
if conv_PR2==1
break
end
end %end of 2 calculation for pressure loop
if conv_PR2==0
PR_slope2=(PR(3)-PR(4))/(PSI(3)-PSI(4));
PR_rel_error_level2=(PR_comp(2)/PR(3))-1;
PR_error_level2=PR_comp(2)-PR(3);
PSI(3)=PSI(3)+(PR_error_level2/PR_slope2)*0.01;
PSI(4)=PSI(3)*1.01;
n_PR=n_PR+1;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Stage 3%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n_PR=0; %setting no. of steps
PR_rel_error_level3=1;
conv_PR3=0;
while conv_PR3 ==0
if abs(PR_rel_error_level3)<1e-4
conv_PR3=1;
end
for m =5:6
% Rotor3 inlet
%%%%%%%%%%%%%%%%%%Calculation for Rotor Inlet3 at mean%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rmean10=rmean9*Rm;
umean10=2*pi*rmean10*N(3)/60;%rotation velocity at Rotor
rt10=rt9;
Cm10=Cm9; %meridional velocity kept constant from stator to Rotor
C_circum10=C_circum9/Rm;
alpha10=atand(C_circum10/Cm10);
C10=Cm10/cosd(alpha10); %absolute velocity calcualted
W_circum10=(umean10-C_circum10); %relative circum velocity
Beta10=atand(W_circum10/Cm10); %blade angle aka relative flow angle
W10=Cm10/cosd(Beta10); %relative velocity at blade inlet
%Stagnation properties at inlet remains constant from stator out to rotor
%in
H010=H09;
P010=P09;
T010=T09;
Dens010=Dens09;
%Static properties at inlet
H10=H010-C10^2/2; %Static heat enthalpy
T10=H10/cp; %Temperature at inlet
a10=(gamma*R*T10)^0.5; % sonic velocity at inlet
a10s=(gammas*R*T010)^0.5; %critical sonic velocity
T010r=T010+(W10^2-C10^2)/(2*cp); %calculated for relative critical sonic velocity
a10sr=(gammas*R*T010r)^0.5;% relative critical sonic velocity
M10=C10/a10; %Mach number calculated at inlet
Mc10s=C10/a10s; %critical mach no. at inlet
Mw10s=W10/a10sr; %critical relative mach no. at inlet
P10=P010/((1+((gamma-1)/2)*M10^2)^(gamma/(gamma-1))); % static pressure at inlet
Dens10=Dens010/((1+((gamma-1)/2)*M10^2)^(1/(gamma-1)));
I10=H10+W10^2/2-umean10^2/2; % inlet rothalpy
%Dimensions for rotor inlet
area10=massflow/(Dens10*Cm10); %Rotor inlet area
rh10=((pi*rt10^2-area10)/pi)^0.5; % Hub unlet radius
%%%%%%%%%%%%%%%%%%%%%Rotor Oulet Calculations%%%%%%%%%%%%%%%%%%%%%%%%
rmean11=rmean10;
%Velocity triangles
rt11=rt10; %tip radius is constant inlet to outlet
Cm11=Cm10; %meridional velocity is constant inlet to outlet
umean11=umean10; % rotational velocity constant inlet to outlet
C_circum11=umean10*PSI(m)+C_circum10; %circumferential velocity calculated using PSI at outlet
C11=(Cm11^2+C_circum11^2)^0.5; % absolute velocity calculated at outlet
alpha11=atand(C_circum11/Cm11); %absolute angle calculated at outlet
W_circum11=umean11-C_circum11; % relative circumferential velocity calc at oulet
Beta11=atand(W_circum11/Cm11); %blade angle calculated at oulet
W11=Cm11/cosd(Beta11); %relative velocity calc. at outlet
%static & total thermal properties
I11=I10; %Rothalpy remains constant
H11=I11-W11^2/2+umean11^2/2; %static enthalpy calculated using Rothalpy
T11=H11/cp; % Temperature at rotor outlet
H011=H11+C11^2/2; % stagnation enthalpy at rotor outlet
a11=(gamma*R*T11)^0.5; % sonic velocity at rotor outlet
M11=C11/a11; % Mach no. at rotor outlet
T011=T11*(1+((gamma-1)*M11^2/2)); % total temp at rotor outlet
a11s=(gammas*R*T011)^0.5;
Mc11s=C11/a11s;
Mu11s=umean11/a11s;
Iso_effi=0.8; % isoentropic efficiency chosen
P011=P010*(1+Iso_effi*(T011-T010)/T010)^(gamma/(gamma-1)); % total pressure calc. using assumed efficiency
Dens011=P011/(R*T011);
P11=P011/(1+((gamma-1)*M11^2/2))^(gamma/(gamma-1)); % Static pressure calc. at rotor outlet
Dens11=Dens011/((1+((gamma-1)/2)*M11^2)^(1/(gamma-1)));
%Rotor Oulet dimensions
area11=massflow/(Dens11*Cm11);
rh11=((pi*rt11^2-area11)/pi)^0.5;
dH3=W11/W10;
%%%%%%%%%%%%%%%%%Stator3 Inlet Calculations%%%%%%%%%%%%%%%%%%%%%%%%%%
rt12=rt11;
Cm12=Cm11;
C_circum12=C_circum11/Rm;
alpha12=atand(C_circum12/Cm12);
C12=Cm12/cosd(alpha12);
%Stagnation enthalpies
H012=H011;
P012=P011;
T012=T011;
Dens012=Dens011;
%Static enthalpies
H12=H012-C12^2/2;
T12=H12/cp;
a12=(gamma*R*T12)^0.5;
a12s=(gammas*R*T012)^0.5;
M12=C12/a12;
Mc12s=C12/a12s;
P12=P012/(1+((gamma-1)*M12^2/2))^(gamma/(gamma-1));
Dens12=Dens012/((1+((gamma-1)/2)*M12^2)^(1/(gamma-1)));
%Stator dmensions
area12=massflow/(Dens12*Cm12);
rh12=((pi*rt12^2-area12)/pi)^0.5;
rmean12=rmean11*Rm;
%%%%%%%%%%%%%%%%Stator Outlet calculations%%%%%%%%%%%%%%%%%%%%%%%%
%alpha13=alpha10;
Cm13=Cm12;
C13=Cm13/cosd(alpha13);
C_circum13=C13*sind(alpha13);
T013=T012;
H013=H012;
P013=P012;
Dens013=Dens012;
H13=H013-C13^2/2;
T13=H13/cp;
a13=(gamma*R*T13)^0.5;
a13s=(gammas*R*T013)^0.5;
M13=C13/a13;
Mc13s=C13/a13s;
P13=P013/(1+((gamma-1)*M13^2/2))^(gamma/(gamma-1));
Dens13=Dens013/((1+((gamma-1)/2)*M13^2)^(1/(gamma-1)));
%Stator dimensions-outlet
rt13=rt12;
area13=massflow/(Dens13*Cm13);
rh13=((pi*rt13^2-area13)/pi)^0.5;
rmean13=rmean12;
React3=(H11-H10)/(H013-H010);% Stage 3 reaction
PR_stage3=P13/P10;
PR(m) = PR_stage3; %total pressure ratio
if conv_PR3==1
break
end
end %end of 2 calculation for pressure loop
if conv_PR3==0
PR_slope3=(PR(5)-PR(6))/(PSI(5)-PSI(6));
PR_rel_error_level3=(PR_comp(3)/PR(5))-1;
PR_error_level3=PR_comp(3)-PR(5);
PSI(5)=PSI(5)+(PR_error_level3/PR_slope3)*0.01;
PSI(6)=PSI(5)*1.01;
n_PR=n_PR+1;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Stage 4%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% n_PR=0; %setting no. of steps
PR_rel_error_level4=1;
conv_PR4=0;
while conv_PR4 ==0
if abs(PR_rel_error_level4)<1e-4
conv_PR4=1;
end
for m =7:8
% Rotor4 inlet
%%%%%%%%%%%%%%%%%%Calculation for Rotor Inlet4 at mean%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rmean14=rmean13*Rm;
umean14=2*pi*rmean14*N(4)/60;%rotation velocity at Rotor
rt14=rt13;
Cm14=Cm13; %meridional velocity kept constant from stator to Rotor
C_circum14=C_circum13/Rm;
alpha14=atand(C_circum14/Cm14);
C14=Cm14/cosd(alpha14); %absolute velocity calcualted
W_circum14=(umean14-C_circum14); %relative circum velocity
Beta14=atand(W_circum14/Cm14); %blade angle aka relative flow angle
W14=Cm14/cosd(Beta14); %relative velocity at blade inlet
%Stagnation properties at inlet remains constant from stator out to rotor
%in
H014=H013;
P014=P013;
T014=T013;
Dens014=Dens013;
%Static properties at inlet
H14=H014-C14^2/2; %Static heat enthalpy
T14=H14/cp; %Temperature at inlet
a14=(gamma*R*T14)^0.5; % sonic velocity at inlet
a14s=(gammas*R*T014)^0.5; %critical sonic velocity
T014r=T014+(W14^2-C14^2)/(2*cp); %calculated for relative critical sonic velocity
a14sr=(gammas*R*T014r)^0.5;% relative critical sonic velocity
M14=C14/a14; %Mach number calculated at inlet
Mc14s=C14/a14s; %critical mach no. at inlet
Mw14s=W14/a14sr; %critical relative mach no. at inlet
P14=P014/((1+((gamma-1)/2)*M14^2)^(gamma/(gamma-1))); % static pressure at inlet
Dens14=Dens014/((1+((gamma-1)/2)*M14^2)^(1/(gamma-1)));
I14=H14+W14^2/2-umean14^2/2; % inlet rothalpy
%Dimensions for rotor inlet
area14=massflow/(Dens14*Cm14); %Rotor inlet area
rh14=((pi*rt14^2-area14)/pi)^0.5; % Hub unlet radius
%%%%%%%%%%%%%%%%%%%%%Rotor Oulet Calculations%%%%%%%%%%%%%%%%%%%%%%%%
rmean15=rmean14;
%Velocity triangles
rt15=rt14; %tip radius is constant inlet to outlet
Cm15=Cm14; %meridional velocity is constant inlet to outlet
umean15=umean14; % rotational velocity constant inlet to outlet
C_circum15=umean15*PSI(m)+C_circum14; %circumferential velocity calculated using PSI at outlet
C15=(Cm15^2+C_circum15^2)^0.5; % absolute velocity calculated at outlet
alpha15=atand(C_circum15/Cm15); %absolute angle calculated at outlet
W_circum15=umean15-C_circum15; % relative circumferential velocity calc at oulet
Beta15=atand(W_circum15/Cm15); %blade angle calculated at oulet
W15=Cm15/cosd(Beta15); %relative velocity calc. at outlet
%static & total thermal properties
I15=I14; %Rothalpy remains constant
H15=I15-W15^2/2+umean15^2/2; %static enthalpy calculated using Rothalpy
T15=H15/cp; % Temperature at rotor outlet
H015=H15+C15^2/2; % stagnation enthalpy at rotor outlet
a15=(gamma*R*T15)^0.5; % sonic velocity at rotor outlet
M15=C15/a15; % Mach no. at rotor outlet
T015=T15*(1+((gamma-1)*M15^2/2)); % total temp at rotor outlet
a15s=(gammas*R*T015)^0.5;
Mc15s=C15/a15s;
Mu15s=umean15/a15s;
Iso_effi=0.8; % isoentropic efficiency chosen
P015=P014*(1+Iso_effi*(T015-T014)/T014)^(gamma/(gamma-1)); % total pressure calc. using assumed efficiency
Dens015=P015/(R*T015);
P15=P015/(1+((gamma-1)*M15^2/2))^(gamma/(gamma-1)); % Static pressure calc. at rotor outlet
Dens15=Dens015/((1+((gamma-1)/2)*M15^2)^(1/(gamma-1)));
%Rotor Oulet dimensions
area15=massflow/(Dens15*Cm15);
rh15=((pi*rt15^2-area15)/pi)^0.5;
dH4=W15/W14;
%%%%%%%%%%%%%%%%%Stator3 Inlet Calculations%%%%%%%%%%%%%%%%%%%%%%%%%%
rt16=rt15;
Cm16=Cm15;
C_circum16=C_circum15/Rm;
alpha16=atand(C_circum16/Cm16);
C16=Cm16/cosd(alpha16);
%Stagnation enthalpies
H016=H015;
P016=P015;
T016=T015;
Dens016=Dens015;
%Static enthalpies
H16=H016-C16^2/2;
T16=H16/cp;
a16=(gamma*R*T16)^0.5;
a16s=(gammas*R*T016)^0.5;
M16=C16/a16;
Mc16s=C16/a16s;
P16=P016/(1+((gamma-1)*M16^2/2))^(gamma/(gamma-1));
Dens16=Dens016/((1+((gamma-1)/2)*M16^2)^(1/(gamma-1)));
%Stator dmensions
area16=massflow/(Dens16*Cm16);
rh16=((pi*rt16^2-area16)/pi)^0.5;
rmean16=rmean15*Rm;
%%%%%%%%%%%%%%%%Stator Outlet calculations%%%%%%%%%%%%%%%%%%%%%%%%
alpha17=0;
Cm17=Cm16;
C17=Cm17/cosd(alpha17);
C_circum17=C17*sind(alpha17);
T017=T016;
H017=H016;
P017=P016;
Dens017=Dens016;
H17=H017-C17^2/2;
T17=H17/cp;
a17=(gamma*R*T17)^0.5;
a17s=(gammas*R*T017)^0.5;
M17=C17/a17;
Mc17s=C17/a17s;
P17=P017/(1+((gamma-1)*M17^2/2))^(gamma/(gamma-1));
Dens17=Dens017/((1+((gamma-1)/2)*M17^2)^(1/(gamma-1)));
%Stator dimensions-outlet
rt17=rt16;
area17=massflow/(Dens17*Cm17);
rh17=((pi*rt17^2-area17)/pi)^0.5;
rmean17=rmean16;
React4=(H15-H14)/(H017-H014);% Stage 4 reaction
PR_stage4=P17/P14;
PR(m) = PR_stage4; %total pressure ratio
if conv_PR4==1
break
end
end %end of 2 calculation for pressure loop
if conv_PR4==0
PR_slope4=(PR(7)-PR(8))/(PSI(7)-PSI(8));
PR_rel_error_level4=(PR_comp(4)/PR(7))-1;
PR_error_level4=PR_comp(4)-PR(7);
PSI(7)=PSI(7)+(PR_error_level4/PR_slope4)*0.01;
PSI(8)=PSI(7)*1.01;
n_PR=n_PR+1;
end
end
Pressure_ratio_static=P17/P1;
% %Height and chord length calculation
% IGV_avg_height=rt-(rh1+rh2)*0.5;
% IGV_chordlength=IGV_avg_height/AR_IGV;
%
% Rotor_avg_height=rt-(rh2+rh3)*0.5;
% R_chordlength=Rotor_avg_height/AR_rotor;
%
% stator_avg_height=rt-(rh3+rh4)*0.5;
% S_chordlength=stator_avg_height/AR_stator;
%
m=20;
%%%%Calculations for Vortex laws i.e. calculation for other streamlines%%%
%%%%%%%%%%%%rt_global=rt; % tip radius constant everywhere.%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Cm_global=Cm1; %axial/meridional velocity constant everywhere
r_station1=linspace(rh1,rt1,m);
r_station2=linspace(rh2,rt2,m);
r_station3=linspace(rh3,rt3,m);
r_station4=linspace(rh4,rt4,m);
r_station5=linspace(rh5,rt5,m);
r_station6=linspace(rh6,rt6,m);
r_station7=linspace(rh7,rt7,m);
r_station8=linspace(rh8,rt8,m);
r_station9=linspace(rh9,rt9,m);
r_station10=linspace(rh10,rt10,m);
r_station11=linspace(rh11,rt11,m);
r_station12=linspace(rh12,rt12,m);
r_station13=linspace(rh13,rt13,m);
r_station14=linspace(rh14,rt14,m);
r_station15=linspace(rh15,rt15,m);
r_station16=linspace(rh16,rt16,m);
r_station17=linspace(rh17,rt17,m);
%%%%%%%loop for calculating variables from hub to tip%%%%%%%%
%%%%%%%%%%%%%%%Variables are calculated from hub to tip%%%%%%%%%%%%%%%%%
H01_station1=H01*ones(1,m);
for i=1:m
%%%%%%%%%%%%%Inlet guide vane/Inlet profile calculation%%%%%%%%%%%%%%%
C_circum1_station1(i)=(C_circum1*rmean1)/(r_station1(i));
alpha1_station1(i)=atand(C_circum1_station1(i)/Cm_global);
C1_station1(i)=Cm_global/cosd(alpha1_station1(i));
%Thermal variables and Mach number
H1_station1(i)=H01_station1(i)-(C1_station1(i)^2)/2;
T1_station1(i)=H1_station1(i)/cp;
a1_station1(i)=(gamma*R*T1_station1(i))^0.5;
a1s_station1(i)=(gammas*R*T01)^0.5;
M1_station1(i)=C1_station1(i)/a1_station1(i);
Mc1s_station1(i)=C1_station1(i)/a1s_station1(i);
P1_station1(i)=P01/((1+((gamma-1)/2)*M1_station1(i)^2)^(gamma/(gamma-1)));
Dens1_station1(i)=P1_station1(i)/(R*T1_station1(i));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Stage 1%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%Rotor Inlet calculation for all radial coordinates%%%%%%%%%%%%%%%
C_circum2_station2(i)=(C_circum2*rmean2)/(r_station2(i));
alpha2_station2(i)=atand(C_circum2_station2(i)/Cm_global);
C2_station2(i)=Cm_global/cosd(alpha2_station2(i));
umean2_station2(i)=(umean2)*r_station2(i)/rmean2;
phi_station2(i)=Cm_global/umean2_station2(i);
W_circum2_station2(i)=(umean2_station2(i)-C_circum2_station2(i));
Beta2_station2(i)=atand(W_circum2_station2(i)/Cm_global);
W2_station2(i)=Cm_global/cosd(Beta2_station2(i));
%Static properties at inlet
H02_station2(i)=H01_station1(i);
H2_station2(i)=H02_station2(i)-C2_station2(i)^2/2;
T2_station2(i)=H2_station2(i)/cp;
a2_station2(i)=(gamma*R*T2_station2(i))^0.5;
a2s_station2(i)=(gammas*R*T02)^0.5;
M2_station2(i)=C2_station2(i)/a2_station2(i);
Mc2s_station2(i)=C2_station2(i)/a2s_station2(i);
T02r_station2(i)=T02+(W2_station2(i)^2-C2_station2(i)^2)/(2*cp);
a2sr_station2(i)=(gammas*R*T02r_station2(i))^0.5;
Mw2s_station2(i)=W2_station2(i)/a2sr_station2(i);
P2_station2(i)=P02/((1+((gamma-1)/2)*M2_station2(i)^2)^(gamma/(gamma-1)));
Dens2_station2(i)=P2_station2(i)/(R*T2_station2(i));
I2_station2(i)=H2_station2(i)+W2_station2(i)^2/2-umean2_station2(i)^2/2;
%%%%%%%%%%%%%%%%%%%%%%%Rotor Oulet Calculations for all radial points%%%%%%%%%%%%%%%%%%%%%%
umean3_station3(i)=(umean3)*r_station3(i)/rmean3;
C_circum3_station3(i)=(C_circum3*rmean3)/(r_station3(i));
C3_station3(i)=(Cm_global^2+ C_circum3_station3(i)^2)^0.5;
alpha3_station3(i)=atand(C_circum3_station3(i)/Cm_global);
W_circum3_station3(i)=umean3_station3(i)-C_circum3_station3(i);
Beta3_station3(i)=atand(W_circum3_station3(i)/Cm_global);
W3_station3(i)=Cm_global/cosd(Beta3_station3(i));
PSI_station(i)=((C_circum3_station3(i)-C_circum2_station2(i)))/(umean3_station3(i));
%%%%static & total thermal properties
I3_station3(i)=I2_station2(i);
H3_station3(i)=I3_station3(i)-W3_station3(i)^2/2+umean3_station3(i)^2/2;
T3_station3(i)=H3_station3(i)/cp;
H03_station3(i)=H3_station3(i)+C3_station3(i)^2/2;
a3_station3(i)=(gamma*R*T3_station3(i))^0.5;
M3_station3(i)=C3_station3(i)/a3_station3(i);
T03_station3(i)=T3_station3(i)*(1+((gamma-1)*M3_station3(i)^2/2));
a3s_station3(i)=(gammas*R*T03_station3(i))^0.5;
Mc3s_station3(i)=C3_station3(i)/a3s_station3(i);
Mu3s_station3(i)=umean3_station3(i)/a3s_station3(i);
T03r_station3(i)=T03+(W3_station3(i)^2-C3_station3(i)^2)/(2*cp);
a3sr_station3(i)=(gammas*R*T03r_station3(i))^0.5;
Mw3s_station3(i)=W3_station3(i)/a3sr_station3(i);
P03_station3(i)=P02*(1+Iso_effi*(T03_station3(i)-T02)/T02)^(gamma/(gamma-1));
P3_station3(i)=P03_station3(i)/(1+((gamma-1)*M3_station3(i)^2/2))^(gamma/(gamma-1));
dH_station(i)=W3_station3(i)/W2_station2(i);
Dens3_station3(i)=P3_station3(i)/(R*T3_station3(i));
%%%%%%%%%%%%%%%%%Stator Inlet Calculations%%%%%%%%%%%%%%%%%%%%%%%%%%
C_circum4_station4(i)=(C_circum4*rmean4)/(r_station4(i));
alpha4_station4(i)=atand(C_circum4_station4(i)/Cm_global);
C4_station4(i)=Cm_global/cosd(alpha4_station4(i));
%stagnation properties
H04_station4(i)=H03_station3(i);
P04_station4(i)=P03_station3(i);
T04_station4(i)=T03_station3(i);
%Static enthalpies
H4_station4(i)=H04_station4(i)-C4_station4(i)^2/2;
T4_station4(i)=H4_station4(i)/cp;
a4_station4(i)=(gamma*R*T4_station4(i))^0.5;
a4s_station4(i)=(gammas*R*T04_station4(i))^0.5;
M4_station4(i)=C4_station4(i)/a4_station4(i);
Mc4s_station4(i)=C4_station4(i)/a4s_station4(i);
P4_station4(i)=P04_station4(i)/(1+((gamma-1)*M4_station4(i)^2/2))^(gamma/(gamma-1));
Dens4_station4(i)=P4_station4(i)/(R*T4_station4(i));
%%%%%%%%%%%%%%%%Stator Outlet calculations%%%%%%%%%%%%%%%%%%%%%%%%
C_circum5_station5(i)=(C_circum5*rmean5)/(r_station5(i));
alpha5_station5(i)=atand(C_circum5_station5(i)/Cm_global);
T05_station5(i)=T04_station4(i);
C5_station5(i)=Cm_global/cosd(alpha5_station5(i));
H05_station5(i)=H04_station4(i);
P05_station5(i)=P04_station4(i);
H5_station5(i)=H05_station5(i)-C5_station5(i)^2/2;
T5_station5(i)=H5_station5(i)/cp;
a5_station5(i)=(gamma*R*T5_station5(i))^0.5;
a5s_station5(i)=(gammas*R*T05_station5(i))^0.5;
M5_station5(i)=C5_station5(i)/a5_station5(i);
Mc5s_station5(i)=C5_station5(i)/a5s_station5(i);
P5_station5(i)=P05_station5(i)/(1+((gamma-1)*M5_station5(i)^2/2))^(gamma/(gamma-1));
Dens5_station5(i)=P5_station5(i)/(R*T5_station5(i));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Stage2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%Rotor Inlet calculation for all radial coordinates%%%%%%%%%%%%%%%
C_circum6_station6(i)=(C_circum6*rmean6)/(r_station6(i));
alpha6_station6(i)=atand(C_circum6_station6(i)/Cm_global);
C6_station6(i)=Cm_global/cosd(alpha6_station6(i));
umean6_station6(i)=(umean6)*r_station6(i)/rmean6;
phi_station6(i)=Cm_global/umean6_station6(i);
W_circum6_station6(i)=(umean6_station6(i)-C_circum6_station6(i));
Beta6_station6(i)=atand(W_circum6_station6(i)/Cm_global);
W6_station6(i)=Cm_global/cosd(Beta6_station6(i));
%Static properties at inlet
H06_station6(i)=H05_station5(i);
T06_station6(i)=T05_station5(i);
H6_station6(i)=H06_station6(i)-C6_station6(i)^2/2;
T6_station6(i)=H6_station6(i)/cp;
a6_station6(i)=(gamma*R*T6_station6(i))^0.5;
a6s_station6(i)=(gammas*R*T06)^0.5;
M6_station6(i)=C6_station6(i)/a6_station6(i);
Mc6s_station6(i)=C6_station6(i)/a6s_station6(i);
P06_station6(i)=P05_station5(i);
T06r_station6(i)=T06_station6(i)+(W6_station6(i)^2-C6_station6(i)^2)/(2*cp);
a6sr_station6(i)=(gammas*R*T06r_station6(i))^0.5;
Mw6s_station6(i)=W6_station6(i)/a6sr_station6(i);
P6_station6(i)=P06_station6(i)/((1+((gamma-1)/2)*M6_station6(i)^2)^(gamma/(gamma-1)));
Dens6_station6(i)=P6_station6(i)/(R*T6_station6(i));
I6_station6(i)=H6_station6(i)+W6_station6(i)^2/2-umean6_station6(i)^2/2;
%%%%%%%%%%%%%%%%%%%%%%%Rotor Oulet Calculations for all radial points%%%%%%%%%%%%%%%%%%%%%%
umean7_station7(i)=(umean7)*r_station7(i)/rmean7;
C_circum7_station7(i)=(C_circum7*rmean7)/(r_station7(i));
C7_station7(i)=(Cm_global^2+ C_circum7_station7(i)^2)^0.5;
alpha7_station7(i)=atand(C_circum7_station7(i)/Cm_global);
W_circum7_station7(i)=umean7_station7(i)-C_circum7_station7(i);
Beta7_station7(i)=atand(W_circum7_station7(i)/Cm_global);
W7_station7(i)=Cm_global/cosd(Beta7_station7(i));
PSI_station7(i)=((C_circum7_station7(i)-C_circum6_station6(i)))/(umean7_station7(i));
%% %static & total thermal properties
I7_station7(i)=I6_station6(i);
H7_station7(i)=I7_station7(i)-W7_station7(i)^2/2+umean7_station7(i)^2/2;
T7_station7(i)=H7_station7(i)/cp;
H07_station7(i)=H7_station7(i)+C7_station7(i)^2/2;
a7_station7(i)=(gamma*R*T7_station7(i))^0.5;
M7_station7(i)=C7_station7(i)/a7_station7(i);
T07_station7(i)=T7_station7(i)*(1+((gamma-1)*M7_station7(i)^2/2));
a7s_station7(i)=(gammas*R*T07_station7(i))^0.5;
Mc7s_station7(i)=C7_station7(i)/a7s_station7(i);
Mu7s_station7(i)=umean7_station7(i)/a7s_station7(i);
T07r_station7(i)=T07_station7(i)+(W7_station7(i)^2-C7_station7(i)^2)/(2*cp);
a7sr_station7(i)=(gammas*R*T07r_station7(i))^0.5;
Mw7s_station7(i)=W7_station7(i)/a7sr_station7(i);
P07_station7(i)=P06_station6(i)*(1+Iso_effi*(T07_station7(i)-T06_station6(i))/T06_station6(i))^(gamma/(gamma-1));%%%check
P7_station7(i)=P07_station7(i)/(1+((gamma-1)*M7_station7(i)^2/2))^(gamma/(gamma-1));
dH2_station(i)=W7_station7(i)/W6_station6(i);
Dens7_station7(i)=P7_station7(i)/(R*T7_station7(i));
%%%%%%%%%%%%%%%%%Stator Inlet Calculations%%%%%%%%%%%%%%%%%%%%%%%%%%
C_circum8_station8(i)=(C_circum8*rmean8)/(r_station8(i));
alpha8_station8(i)=atand(C_circum8_station8(i)/Cm_global);
C8_station8(i)=Cm_global/cosd(alpha8_station8(i));
%stagnation properties
H08_station8(i)=H07_station7(i);
P08_station8(i)=P07_station7(i);
T08_station8(i)=T07_station7(i);
%Static enthalpies
H8_station8(i)=H08_station8(i)-C8_station8(i)^2/2;
T8_station8(i)=H8_station8(i)/cp;
a8_station8(i)=(gamma*R*T8_station8(i))^0.5;
a8s_station8(i)=(gammas*R*T08_station8(i))^0.5;
M8_station8(i)=C8_station8(i)/a8_station8(i);
Mc8s_station8(i)=C8_station8(i)/a8s_station8(i);
P8_station8(i)=P08_station8(i)/(1+((gamma-1)*M8_station8(i)^2/2))^(gamma/(gamma-1));
Dens8_station8(i)=P8_station8(i)/(R*T8_station8(i));
%%%%%%%%%%%%%%%%Stator Outlet calculations%%%%%%%%%%%%%%%%%%%%%%%%
C_circum9_station9(i)=(C_circum9*rmean9)/(r_station9(i));
alpha9_station9(i)=atand(C_circum9_station9(i)/Cm_global);
T09_station9(i)=T08_station8(i);
C9_station9(i)=Cm_global/cosd(alpha9_station9(i));
H09_station9(i)=H08_station8(i);
P09_station9(i)=P08_station8(i);
H9_station9(i)=H09_station9(i)-C9_station9(i)^2/2;
T9_station9(i)=H9_station9(i)/cp;
a9_station9(i)=(gamma*R*T9_station9(i))^0.5;
a9s_station9(i)=(gammas*R*T09_station9(i))^0.5;
M9_station9(i)=C9_station9(i)/a9_station9(i);
Mc9s_station9(i)=C9_station9(i)/a9s_station9(i);
P9_station9(i)=P09_station9(i)/(1+((gamma-1)*M9_station9(i)^2/2))^(gamma/(gamma-1));
Dens9_station9(i)=P9_station9(i)/(R*T9_station9(i));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% stage3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%Rotor Inlet calculation for all radial coordinates%%%%%%%%%%%%%%%
H010_station10(i)=H09_station9(i);
C_circum10_station10(i)=(C_circum10*rmean10)/(r_station10(i));
alpha10_station10(i)=atand(C_circum10_station10(i)/Cm_global);
C10_station10(i)=Cm_global/cosd(alpha10_station10(i));
umean10_station10(i)=(umean10)*r_station10(i)/rmean10;
phi_station10(i)=Cm_global/umean10_station10(i);
W_circum10_station10(i)=(umean10_station10(i)-C_circum10_station10(i));
Beta10_station10(i)=atand(W_circum10_station10(i)/Cm_global);
W10_station10(i)=Cm_global/cosd(Beta10_station10(i));
%Static properties at inlet
T010_station10(i)=T09_station9(i);
H10_station10(i)=H010-C10_station10(i)^2/2;
T10_station10(i)=H10_station10(i)/cp;
a10_station10(i)=(gamma*R*T10_station10(i))^0.5;
a10s_station10(i)=(gammas*R*T010)^0.5;
M10_station10(i)=C10_station10(i)/a10_station10(i);
Mc10s_station10(i)=C10_station10(i)/a10s_station10(i);
P010_station10(i)=P09_station9(i);
T010r_station10(i)=T010_station10(i)+(W10_station10(i)^2-C10_station10(i)^2)/(2*cp);
a10sr_station10(i)=(gammas*R*T010r_station10(i))^0.5;
Mw10s_station10(i)=W10_station10(i)/a10sr_station10(i);
P10_station10(i)=P010_station10(i)/((1+((gamma-1)/2)*M10_station10(i)^2)^(gamma/(gamma-1)));
Dens10_station10(i)=P10_station10(i)/(R*T10_station10(i));
I10_station10(i)=H10_station10(i)+W10_station10(i)^2/2-umean10_station10(i)^2/2;
%%%%%%%%%%%%%%%%%%%%%%%Rotor Oulet Calculations for all radial points%%%%%%%%%%%%%%%%%%%%%%
umean11_station11(i)=(umean11)*r_station11(i)/rmean11;
C_circum11_station11(i)=(C_circum11*rmean11)/(r_station11(i));
C11_station11(i)=(Cm_global^2+ C_circum11_station11(i)^2)^0.5;
alpha11_station11(i)=atand(C_circum11_station11(i)/Cm_global);
W_circum11_station11(i)=umean11_station11(i)-C_circum11_station11(i);
Beta11_station11(i)=atand(W_circum11_station11(i)/Cm_global);
W11_station11(i)=Cm_global/cosd(Beta11_station11(i));
PSI_station11(i)=((C_circum11_station11(i)-C_circum10_station10(i)))/(umean11_station11(i));
%% %static & total thermal properties
I11_station11(i)=I10_station10(i);
H11_station11(i)=I11_station11(i)-W11_station11(i)^2/2+umean11_station11(i)^2/2;
T11_station11(i)=H11_station11(i)/cp;
H011_station11(i)=H11_station11(i)+C11_station11(i)^2/2;
a11_station11(i)=(gamma*R*T11_station11(i))^0.5;
M11_station11(i)=C11_station11(i)/a11_station11(i);
T011_station11(i)=T11_station11(i)*(1+((gamma-1)*M11_station11(i)^2/2));
a11s_station11(i)=(gammas*R*T011_station11(i))^0.5;
Mc11s_station11(i)=C11_station11(i)/a11s_station11(i);
Mu11s_station11(i)=umean11_station11(i)/a11s_station11(i);
T011r_station11(i)=T011_station11(i)+(W11_station11(i)^2-C11_station11(i)^2)/(2*cp);
a11sr_station11(i)=(gammas*R*T011r_station11(i))^0.5;
Mw11s_station11(i)=W11_station11(i)/a11sr_station11(i);
P011_station11(i)=P010_station10(i)*(1+Iso_effi*(T011_station11(i)-T010_station10(i))/T010_station10(i))^(gamma/(gamma-1));%%%check
P11_station11(i)=P011_station11(i)/(1+((gamma-1)*M11_station11(i)^2/2))^(gamma/(gamma-1));
dH3_station(i)=W11_station11(i)/W10_station10(i);
Dens11_station11(i)=P11_station11(i)/(R*T11_station11(i));
%%%%%%%%%%%%%%%%%Stator Inlet Calculations%%%%%%%%%%%%%%%%%%%%%%%%%%
C_circum12_station12(i)=(C_circum12*rmean12)/(r_station12(i));
alpha12_station12(i)=atand(C_circum12_station12(i)/Cm_global);
C12_station12(i)=Cm_global/cosd(alpha12_station12(i));
%stagnation properties
H012_station12(i)=H011_station11(i);
P012_station12(i)=P011_station11(i);
T012_station12(i)=T011_station11(i);
%Static enthalpies
H12_station12(i)=H012_station12(i)-C12_station12(i)^2/2;
T12_station12(i)=H12_station12(i)/cp;
a12_station12(i)=(gamma*R*T12_station12(i))^0.5;
a12s_station12(i)=(gammas*R*T012_station12(i))^0.5;
M12_station12(i)=C12_station12(i)/a12_station12(i);
Mc12s_station12(i)=C12_station12(i)/a12s_station12(i);
P12_station12(i)=P012_station12(i)/(1+((gamma-1)*M12_station12(i)^2/2))^(gamma/(gamma-1));
Dens12_station12(i)=P12_station12(i)/(R*T12_station12(i));
%%%%%%%%%%%%%%%%Stator Outlet calculations%%%%%%%%%%%%%%%%%%%%%%%%
C_circum13_station13(i)=(C_circum13*rmean13)/(r_station13(i));
alpha13_station13(i)=atand(C_circum13_station13(i)/Cm_global);
T013_station13(i)=T012_station12(i);
C13_station13(i)=Cm_global/cosd(alpha13_station13(i));
H013_station13(i)=H012_station12(i);
P013_station13(i)=P012_station12(i);
H13_station13(i)=H013_station13(i)-C13_station13(i)^2/2;
T13_station13(i)=H13_station13(i)/cp;
a13_station13(i)=(gamma*R*T13_station13(i))^0.5;
a13s_station13(i)=(gammas*R*T013_station13(i))^0.5;
M13_station13(i)=C13_station13(i)/a13_station13(i);
Mc13s_station13(i)=C13_station13(i)/a13s_station13(i);
P13_station13(i)=P013_station13(i)/(1+((gamma-1)*M13_station13(i)^2/2))^(gamma/(gamma-1));
Dens13_station13(i)=P13_station13(i)/(R*T13_station13(i));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Stage 4%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%Rotor Inlet calculation for all radial coordinates%%%%%%%%%%%%%%%
C_circum14_station14(i)=(C_circum14*rmean14)/(r_station14(i));
alpha14_station14(i)=atand(C_circum14_station14(i)/Cm_global);
C14_station14(i)=Cm_global/cosd(alpha14_station14(i));
umean14_station14(i)=(umean14)*r_station14(i)/rmean14;
phi_station14(i)=Cm_global/umean14_station14(i);
W_circum14_station14(i)=(umean14_station14(i)-C_circum14_station14(i));
Beta14_station14(i)=atand(W_circum14_station14(i)/Cm_global);
W14_station14(i)=Cm_global/cosd(Beta14_station14(i));
%Static properties at inlet%
H014_station14(i)=H013_station13(i);
T014_station14(i)=T013_station13(i);
H14_station14(i)=H014_station14(i)-C14_station14(i)^2/2;
T14_station14(i)=H14_station14(i)/cp;
a14_station14(i)=(gamma*R*T14_station14(i))^0.5;
a14s_station14(i)=(gammas*R*T014)^0.5;
M14_station14(i)=C14_station14(i)/a14_station14(i);
Mc14s_station14(i)=C14_station14(i)/a14s_station14(i);
P014_station14(i)=P013_station13(i);
T014r_station14(i)=T014_station14(i)+(W14_station14(i)^2-C14_station14(i)^2)/(2*cp);
a14sr_station14(i)=(gammas*R*T014r_station14(i))^0.5;
Mw14s_station14(i)=W14_station14(i)/a14sr_station14(i);
P14_station14(i)=P014_station14(i)/((1+((gamma-1)/2)*M14_station14(i)^2)^(gamma/(gamma-1)));
Dens14_station14(i)=P14_station14(i)/(R*T14_station14(i));
I14_station14(i)=H14_station14(i)+W14_station14(i)^2/2-umean14_station14(i)^2/2;
%%%%%%%%%%%%%%%%%%%%%%%Rotor Oulet Calculations for all radial points%%%%%%%%%%%%%%%%%%%%%%
umean15_station15(i)=(umean15)*r_station15(i)/rmean15;
C_circum15_station15(i)=(C_circum15*rmean15)/(r_station15(i));
C15_station15(i)=(Cm_global^2+ C_circum15_station15(i)^2)^0.5;
alpha15_station15(i)=atand(C_circum15_station15(i)/Cm_global);
W_circum15_station15(i)=umean15_station15(i)-C_circum15_station15(i);
Beta15_station15(i)=atand(W_circum15_station15(i)/Cm_global);
W15_station15(i)=Cm_global/cosd(Beta15_station15(i));
PSI_station15(i)=((C_circum15_station15(i)-C_circum14_station14(i)))/(umean15_station15(i));
%% %static & total thermal properties
I15_station15(i)=I14_station14(i);
H15_station15(i)=I15_station15(i)-W15_station15(i)^2/2+umean15_station15(i)^2/2;
T15_station15(i)=H15_station15(i)/cp;
H015_station15(i)=H15_station15(i)+C15_station15(i)^2/2;
a15_station15(i)=(gamma*R*T15_station15(i))^0.5;
M15_station15(i)=C15_station15(i)/a15_station15(i);
T015_station15(i)=T15_station15(i)*(1+((gamma-1)*M15_station15(i)^2/2));
a15s_station15(i)=(gammas*R*T015_station15(i))^0.5;
Mc15s_station15(i)=C15_station15(i)/a15s_station15(i);
Mu15s_station15(i)=umean15_station15(i)/a15s_station15(i);
T015r_station15(i)=T015_station15(i)+(W15_station15(i)^2-C15_station15(i)^2)/(2*cp);
a15sr_station15(i)=(gammas*R*T015r_station15(i))^0.5;
Mw15s_station15(i)=W15_station15(i)/a15sr_station15(i);
P015_station15(i)=P014_station14(i)*(1+Iso_effi*(T015_station15(i)-T014_station14(i))/T014_station14(i))^(gamma/(gamma-1));%%%check
P15_station15(i)=P015_station15(i)/(1+((gamma-1)*M15_station15(i)^2/2))^(gamma/(gamma-1));
dH4_station(i)=W15_station15(i)/W14_station14(i);
Dens15_station15(i)=P15_station15(i)/(R*T15_station15(i));
%%%%%%%%%%%%%%%%%Stator Inlet Calculations%%%%%%%%%%%%%%%%%%%%%%%%%%
C_circum16_station16(i)=(C_circum16*rmean16)/(r_station16(i));
alpha16_station16(i)=atand(C_circum16_station16(i)/Cm_global);
C16_station16(i)=Cm_global/cosd(alpha16_station16(i));
%stagnation properties
H016_station16(i)=H015_station15(i);
P016_station16(i)=P015_station15(i);
T016_station16(i)=T015_station15(i);
%Static enthalpies
H16_station16(i)=H016_station16(i)-C16_station16(i)^2/2;
T16_station16(i)=H16_station16(i)/cp;
a16_station16(i)=(gamma*R*T16_station16(i))^0.5;
a16s_station16(i)=(gammas*R*T016_station16(i))^0.5;
M16_station16(i)=C16_station16(i)/a16_station16(i);
Mc16s_station16(i)=C16_station16(i)/a16s_station16(i);
P16_station16(i)=P016_station16(i)/(1+((gamma-1)*M16_station16(i)^2/2))^(gamma/(gamma-1));
Dens16_station16(i)=P16_station16(i)/(R*T16_station16(i));
%%%%%%%%%%%%%%%%Stator Outlet calculations%%%%%%%%%%%%%%%%%%%%%%%%
C_circum17_station17(i)=(C_circum17*rmean17)/(r_station17(i));
alpha17_station17(i)=atand(C_circum17_station17(i)/Cm_global);
T017_station17(i)=T016_station16(i);
C17_station17(i)=Cm_global/cosd(alpha17_station17(i));
H017_station17(i)=H016_station16(i);
P017_station17(i)=P016_station16(i);
H17_station17(i)=H017_station17(i)-C17_station17(i)^2/2;
T17_station17(i)=H17_station17(i)/cp;
a17_station17(i)=(gamma*R*T17_station17(i))^0.5;
a17s_station17(i)=(gammas*R*T017_station17(i))^0.5;
M17_station17(i)=C17_station17(i)/a17_station17(i);
Mc17s_station17(i)=C17_station17(i)/a17s_station17(i);
P17_station17(i)=P017_station17(i)/(1+((gamma-1)*M17_station17(i)^2/2))^(gamma/(gamma-1));
Dens17_station17(i)=P17_station17(i)/(R*T17_station17(i));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%IGV inlet to IGV Outlet calculations from Hub to Tip%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%m=20; %number of nodes from inlet to outlet at mean position
t=linspace(0,1,m); %nodal points
cont_x_igv=0.6; %control point on x -->IGV
t_station=meshgrid(t);
Cm_station0_1=Cm_global*ones(m,m);
C_station_circum0_1=zeros(m,m);
alpha_station_0_1=zeros(m,m);
C_station_0_1=zeros(m,m);
x_station_circum0_1=zeros(m,m);
contx_station0_1=cont_x_igv;
C_circum0_station0=zeros(1,m);
C_circum1_station1=C_circum2_station2;
cont_y_station0_1=(C_circum0_station0+C_circum1_station1)/2;
for i=1:m
for j=1:m
C_station_circum0_1(i,j)=((1-t_station(i,j))^2)*(C_circum0_station0(i))+2*(1-t_station(i,j))*t_station(i,j)*(cont_y_station0_1(i))+t_station(i,j)^2*(C_circum1_station1(i));
x_station_circum0_1(i,j)=((1-t_station(i,j))^2)*(0)+2*(1-t_station(i,j))*t_station(i,j)*(contx_station0_1)+t_station(i,j)^2*(t_station(i,j));
C_station_0_1(i,j)=(Cm_global^2+C_station_circum0_1(i,j)^2)^0.5;
alpha_station_0_1(i,j)=-atand(C_station_circum0_1(i,j)/Cm_global);
alpha_station_0_1_m(i,j)=tand(alpha_station_0_1(i,j));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%to check plot%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%plot(x_station_circum0_1(1,:),C_station_circum0_1(1,:));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%plotting pressure side of the blade at mean%%%%%%%%%%%%%%%
x_station_circum0_1_mid=zeros(m,m);
y_station_circum0_1_mid=zeros(m,m);
y_station_circum0_1=zeros(m,m);
for i=1:m
for j=2:m
x_station_circum0_1_mid(i,j)=((x_station_circum0_1(i,j)+x_station_circum0_1(i,j-1))/2);
y_station_circum0_1_mid(i,j)= alpha_station_0_1_m(i,j-1)*(x_station_circum0_1_mid(i,j)-x_station_circum0_1_mid(i,j-1))+y_station_circum0_1_mid(i,j-1);
y_station_circum0_1(i,j)=alpha_station_0_1_m(i,j)*(-x_station_circum0_1_mid(i,j)+x_station_circum0_1(i,j))+y_station_circum0_1_mid(i,j);
end
IGV_Stagger(i)=atand((y_station_circum0_1(i,m)-y_station_circum0_1(i,1))/(x_station_circum0_1(i,m)-x_station_circum0_1(i,1)));
end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%to check plot%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot(x_station_circum0_1(5,:),y_station_circum0_1(5,:)) %pressure side at mean
% hold on
% plot([x_station_circum0_1(5,1) x_station_circum0_1(5,m)], [y_station_circum0_1(5,1) y_station_circum0_1(5,m)],'color','r')
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% %%%%%%%%%%%%%%%%%%%%%%Plotting suction side at mean%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%assuming thickness distribution using bezier curve%%%%%%%%%%%%%
Thickness0=0;
Thick_y_igv=0.04;
Thickness1=0;
Thick_x_igv=cont_x_igv;
for i=1:m
for j=1:m
Thickness_station_0_1(i,j)=((1-t_station(i,j))^2)*(Thickness0)+2*(1-t_station(i,j))*t_station(i,j)*(Thick_y_igv)+t_station(i,j)^2*(Thickness1);
x_station_circum0_1_s(i,j)=((1-t_station(i,j))^2)*(0)+2*(1-t_station(i,j))*t_station(i,j)*(Thick_x_igv)+t_station(i,j)^2*(t_station(i,j));
y_station_circum0_1_s(i,j)=y_station_circum0_1(i,j)+ Thickness_station_0_1(i,j);
end
end
% %plotting suction side with blade geometry
% plot(x_station_circum0_1_s(5,:), y_station_circum0_1_s(5,:),'color','b') %suction side at mean
% hold on
% plot(x_station_circum0_1(5,:),y_station_circum0_1(5,:),'color','r') %pressure side at mean
% hold on
% plot([x_station_circum0_1(5,1) x_station_circum0_1(5,m)], [y_station_circum0_1(5,1) y_station_circum0_1(5,m)],'--')
for i=1:m
for j=2:m %calculating mid points for suction side
x_station_circum0_1_mid_s(i,j)=((x_station_circum0_1_s(i,j)+ x_station_circum0_1_s(i,j-1))/2);
y_station_circum0_1_mid_s(i,j)=y_station_circum0_1_mid(i,j)+Thickness_station_0_1(i,j);
end
end
%
% %%%%%%%%%%%%%%%%%finding blade angles and velocities at suction side%%%%%%%
C_station_circum0_1_s=zeros(m,m);
for i=1:m
for k=2:m
alpha_station0_1_m_s(i,k-1)=(y_station_circum0_1_mid_s(i,k)-y_station_circum0_1_mid_s(i,k-1))/(x_station_circum0_1_mid_s(i,k)-x_station_circum0_1_mid_s(i,k-1));
alpha_station_0_1_s(i,k-1)=atand(alpha_station0_1_m_s(i,k-1));
C_station_circum0_1_s(i,k-1)=Cm_global*alpha_station0_1_m_s(i,k-1);
C_station_0_1_s(i,k-1)=(Cm_global^2+C_station_circum0_1_s(i,k-1)^2)^0.5;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%% Loop for stages 1-4 for Rotor%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Umean_stage(1,:)=umean2_station2;
Umean_stage(2,:)=umean6_station6;
Umean_stage(3,:)=umean10_station10;
Umean_stage(4,:)=umean14_station14;
Controly_stage_r(1,:)=(C_circum2_station2+C_circum3_station3)/4;
Controly_stage_r(2,:)=(C_circum6_station6+C_circum7_station7)/4;
Controly_stage_r(3,:)=(C_circum10_station10+C_circum11_station11)/4;
Controly_stage_r(4,:)=(C_circum14_station14+C_circum15_station15)/4;
I_stage_ir(1,:)=I2_station2;
I_stage_ir(2,:)=I6_station6;
I_stage_ir(3,:)=I10_station10;
I_stage_ir(4,:)=I14_station14;
%
% Dens_avg_r(1,:)=(Dens2_station2+Dens3_station3)/2;
% Dens_avg_r(2,:)=(Dens6_station6+Dens7_station7)/2;
% Dens_avg_r(3,:)=(Dens10_station10+Dens11_station11)/2;
% Dens_avg_r(4,:)=(Dens14_station14+Dens15_station15)/2;
%
C_circum_stage_ir(1,:)=C_circum2_station2; %inlet circum vel for all rotors
C_circum_stage_ir(2,:)=C_circum6_station6;
C_circum_stage_ir(3,:)=C_circum10_station10;
C_circum_stage_ir(4,:)=C_circum14_station14;
C_circum_stage_or(1,:)=C_circum3_station3; %outlet circum vel for all rotors
C_circum_stage_or(2,:)=C_circum7_station7;
C_circum_stage_or(3,:)=C_circum11_station11;
C_circum_stage_or(4,:)=C_circum15_station15;
for s= 1:4
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%Rotor inlet to Rotor Outlet calculations from Hub to Tip%%%%%%%%%
% %%%%%%%%%Velocity triangle calculation from Rotor inlet to Oulet%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
t_station=meshgrid(t);
umean_station_2_3(:,:,s)=meshgrid(Umean_stage(s,:))'; %Umean_stage(s,:)
Cm_station_2_3(:,:,s)=Cm_global*ones(m,m);
C_circum_station_2_3(:,:,s)=zeros(m,m);
alpha_station_2_3(:,:,s)=zeros(m,m);
Beta_station_2_3(:,:,s)=zeros(m,m);
W_station_2_3(:,:,s)=zeros(m,m);
C_station_2_3(:,:,s)=zeros(m,m);
W_station_circum2_3(:,:,s)=zeros(m,m);
x_station_circum2_3(:,:,s)=zeros(m,m);
cont_x=[0.4 0.4 0.4 0.4]; %control point in x for C_circum
cont_x_station2_3=cont_x(s);
cont_y_station2_3(s,:)=Controly_stage_r(s,:);%Controly_stage_r(s,:)
for i=1:m
for j=1:m
C_circum_station_2_3(i,j,s)=((1-t_station(i,j))^2)*(C_circum_stage_ir(s,i))+2*(1-t_station(i,j))*t_station(i,j)*(cont_y_station2_3(s,i))+t_station(i,j)^2*(C_circum_stage_or(s,i));
x_station_circum2_3(i,j,s)=((1-t_station(i,j))^2)*(0)+2*(1-t_station(i,j))*t_station(i,j)*(cont_x_station2_3)+(t_station(i,j)^2)*(t_station(i,j));
C_station_2_3(i,j,s)=(Cm_station_2_3(j)^2+C_circum_station_2_3(i,j,s)^2)^0.5;
alpha_station_2_3(i,j,s)=atand(C_circum_station_2_3(i,j,s)/Cm_station_2_3(j));
W_station_circum2_3(i,j,s)=umean_station_2_3(i,j,s)-C_circum_station_2_3(i,j,s);
W_station_2_3(i,j,s)=(Cm_station_2_3(j)^2+W_station_circum2_3(i,j,s)^2)^0.5;
Beta_station_2_3(i,j,s)=atand(W_station_circum2_3(i,j,s)/Cm_station_2_3(j));
Beta_station_2_3_m(i,j,s)=tand(Beta_station_2_3(i,j,s));
Reynold_2_3(i,j,s)=W_station_2_3(i,j,s)*x_station_circum2_3(i,j,s)/kvisc;
Boundarylayer_2_3(i,j,s)=0.37*x_station_circum2_3(i,j,s)/(Reynold_2_3(i,j,s))^0.2;
H_stage_2_3(i,j,s)=I_stage_ir(s,i)-W_station_2_3(i,j,s)^2/2+Umean_stage(s,i)^2/2;
T_stage_2_3(i,j,s)=H_stage_2_3(i,j,s)/cp;
H0_stage_2_3(i,j,s)=H_stage_2_3(i,j,s)+C_station_2_3(i,j,s)^2/2;
T0_stage_2_3(i,j,s)=H0_stage_2_3(i,j,s)/cp;
T0r_stage_2_3(i,j,s)=T0_stage_2_3(i,j,s)+(W_station_2_3(i,j,s)^2-C_station_2_3(i,j,s)^2)/(2*cp);
asr_stage_2_3(i,j,s)=(gammas*R*T0r_stage_2_3(i,j,s))^0.5;
Mws_stage_2_3(i,j,s)=W_station_2_3(i,j,s)/asr_stage_2_3(i,j,s);
as_stage_2_3(i,j,s)=(gammas*R*T0_stage_2_3(i,j,s))^0.5;
Mcs_stage_2_3(i,j,s)=C_station_2_3(i,j,s)/as_stage_2_3(i,j,s);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%Just to check plot%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% figure(1)
% plot(x_station_circum2_3(1,:,2),C_circum_station_2_3(1,:,2),'color','g'); %hub
% hold on
% plot(x_station_circum2_3(m/2,:,2),C_circum_station_2_3(m/2,:,2),'color','b');%mid
% hold on
% plot(x_station_circum2_3(m,:,2),C_circum_station_2_3(m,:,2),'color','k');%tip
% xlabel("Along Streamline-Inlet to Outlet")
% ylabel("Cu*r")
% legend("Hub side", "Mean",'Tip side')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%plotting pressure side of the blade at mean%%%%%%%%%%%%%%%
x_station_circum2_3_mid(:,:,s)=zeros(m,m);
y_station_circum2_3_mid(:,:,s)=zeros(m,m);
y_station_circum2_3(:,:,s)=zeros(m,m);
for i=1:m
for j=2:m
x_station_circum2_3_mid(i,j,s)=((x_station_circum2_3(i,j,s)+x_station_circum2_3(i,j-1,s))/2);%+x_station_circum2_3(i,j-1);
y_station_circum2_3_mid(i,j,s)=Beta_station_2_3_m(i,j-1,s)*(x_station_circum2_3_mid(i,j,s)-x_station_circum2_3_mid(i,j-1,s))+y_station_circum2_3_mid(i,j-1,s);
y_station_circum2_3(i,j,s)=Beta_station_2_3_m(i,j,s)*(-x_station_circum2_3_mid(i,j,s)+x_station_circum2_3(i,j,s))+y_station_circum2_3_mid(i,j,s);
end
Rotor_Stagger(i,s)=atand((y_station_circum2_3(i,m,s)-y_station_circum2_3(i,1,s))/(x_station_circum2_3(i,m,s)-x_station_circum2_3(i,1,s)));
end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Plot checking%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot(x_station_circum2_3(1,:,1),y_station_circum2_3(1,:,1)) %pressure side at mean
% hold on
% plot([x_station_circum2_3(1,1,1) x_station_circum2_3(1,m,1)], [y_station_circum2_3(1,1,1) y_station_circum2_3(1,m,1)],'--','color','b')
% hold on
%contourf(Mws_stagebystage_2_3x,Mws_stagebystage_2_3y,Mws_stagebystage_2_3_rotor)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%Plotting suction side -rotor hub to tip%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%assuming thickness distribution using bezier curve%%%%%%%%%
%
% x_station_circum2_3_s(:,:,s)=zeros(m,m);
% Thickness_station2_3(:,:,s)=zeros(m,m);
% y_station_circum2_3_s(:,:,s)=zeros(m,m);
Thickness2=0;
Thick_y=0.08*ones(m);
Thickness3=0;
Thick_x=0.5;
for i=1:m
for j=1:m
x_station_circum2_3_s(i,j,s)=((1-t_station(i,j))^2)*(0)+2*(1-t_station(i,j))*t_station(i,j)*(Thick_x)+t_station(i,j)^2*(t_station(i,j));
Thickness_station2_3(i,j,s)=((1-t_station(i,j))^2)*(Thickness2)+2*(1-t_station(i,j))*t_station(i,j)*(Thick_y(i))+t_station(i,j)^2*(Thickness3);
y_station_circum2_3_s(i,j,s)=y_station_circum2_3(i,j,s)+Thickness_station2_3(i,j,s);
end
end
% %%plotting suction side with blade geometry
% plot(x_station_circum2_3_s(1,:,1), y_station_circum2_3_s(1,:,1),'--','color','g') %suction side at mean
% hold on
% plot(x_station_circum2_3(1,:,1),y_station_circum2_3(1,:,1),'color','r') %pressure side at mean
% hold on
% plot([x_station_circum2_3(1,1,1) x_station_circum2_3(1,m,1)], [y_station_circum2_3(1,1,1) y_station_circum2_3(1,m,1)],'--')
%
for i=1:m
for j=2:m %calculating mid points for suction side
x_station_circum2_3_mid_s(i,j,s)=((x_station_circum2_3_s(i,j,s)+ x_station_circum2_3_s(i,j-1,s))/2);%+x_circum2_3(j-1); %calculating midpoints
y_station_circum2_3_mid_s(i,j,s)=y_station_circum2_3_mid(i,j,s)+Thickness_station2_3(i,j,s); % change thickness(j) to Thickness_station2_3(i,j)
end
end
%
%
% %%%%%%%%%%%%%%finding blade angles and velocities at suction side%%%%%%%%%%
for i=1:m
for k=2:m
Beta_station_2_3_m_s(i,k-1,s)=(y_station_circum2_3_mid_s(i,k,s)-y_station_circum2_3_mid_s(i,k-1,s))/(x_station_circum2_3_mid_s(i,k,s)-x_station_circum2_3_mid_s(i,k-1,s));
Beta_station_2_3_s(i,k-1,s)=atand(Beta_station_2_3_m_s(i,k-1,s));
W_station_circum2_3_s(i,k-1,s)=Cm_global*Beta_station_2_3_m_s(i,k-1,s);
W_station_2_3_s(i,k-1,s)=(Cm_global^2+W_station_circum2_3_s(i,k-1,s)^2)^0.5;
C_station_circum2_3_s(i,k-1,s)=umean_station_2_3(i,k-1,s)-W_station_circum2_3_s(i,k-1,s);
alpha_station_2_3_s(i,k-1,s)=atand(C_station_circum2_3_s(i,k-1,s)/Cm_global);
C_station_2_3_s(i,k-1,s)=(Cm_global^2+C_station_circum2_3_s(i,k-1,s)^2)^0.5;
Reynold_2_3_s(i,k-1,s)=W_station_2_3_s(i,k-1,s)*x_station_circum2_3_s(i,k-1,s)/kvisc;
Boundarylayer_2_3_s(i,k-1,s)=0.37*x_station_circum2_3_s(i,k-1,s)/(Reynold_2_3_s(i,k-1,s))^0.2;
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%% Loop for stages 1-4 for Stator%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Controly_stage_s(1,:)=(C_circum4_station4+C_circum5_station5)/4;
Controly_stage_s(2,:)=(C_circum8_station8+C_circum9_station9)/4;
Controly_stage_s(3,:)=(C_circum12_station12+C_circum13_station13)/4;
Controly_stage_s(4,:)=(C_circum16_station16+C_circum17_station17)/1.5;
C_circum_stage_is(1,:)=C_circum4_station4; %inlet circum vel for all stators
C_circum_stage_is(2,:)=C_circum8_station8;
C_circum_stage_is(3,:)=C_circum12_station12;
C_circum_stage_is(4,:)=C_circum16_station16;
C_circum_stage_os(1,:)=C_circum5_station5; %outlet circum vel for all stators
C_circum_stage_os(2,:)=C_circum9_station9;
C_circum_stage_os(3,:)=C_circum13_station13;
C_circum_stage_os(4,:)=C_circum17_station17;
for s =1:4
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % %%%%%%%%%Stator inlet to Stator Outlet calculations from Hub to Tip%%%%%%%%%
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Cm_station4_5(:,:,s)=Cm_global*ones(m,m);
C_station_circum4_5(:,:,s)=zeros(m,m);
alpha_station_4_5(:,:,s)=zeros(m,m);
C_station_4_5(:,:,s)=zeros(m,m);
x_station_circum4_5(:,:,s)=zeros(m,m);
cont_x_s(s)=0.4;
contx_station4_5(s)=cont_x_s(s);
cont_y_station4_5(s,:)=Controly_stage_s(s,:);%Controly_stage_s(s,:)
for i=1:m
for j=1:m
C_station_circum4_5(i,j,s)=((1-t_station(i,j))^2)*(C_circum_stage_is(s,i))+2*(1-t_station(i,j))*t_station(i,j)*(cont_y_station4_5(s,i))+t_station(i,j)^2*(C_circum_stage_os(s,i));
x_station_circum4_5(i,j,s)=((1-t_station(i,j))^2)*(0)+2*(1-t_station(i,j))*t_station(i,j)*(contx_station4_5(s))+t_station(i,j)^2*(t_station(i,j));
C_station_4_5(i,j,s)=(Cm_global^2+C_station_circum4_5(i,j,s)^2)^0.5;
alpha_station_4_5(i,j,s)=-atand(C_station_circum4_5(i,j,s)/Cm_global);
alpha_station_4_5_m(i,j,s)=tand(alpha_station_4_5(i,j,s));
end
end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%to check plot%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot(x_station_circum4_5(1,:,s),C_station_circum4_5(1,:,s));
% hold on
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% %%%%%%%%%%%%%%%%%plotting pressure side of the blade at mean%%%%%%%%%%%%%%%
x_station_circum4_5_mid(:,:,s)=zeros(m,m);
y_station_circum4_5_mid(:,:,s)=zeros(m,m);
y_station_circum4_5(:,:,s)=zeros(m,m);
for i=1:m
for j=2:m
x_station_circum4_5_mid(i,j,s)=((x_station_circum4_5(i,j,s)+x_station_circum4_5(i,j-1,s))/2);%+x_station_circum4_5(i,j-1);
y_station_circum4_5_mid(i,j,s)= alpha_station_4_5_m(i,j-1,s)*(x_station_circum4_5_mid(i,j,s)-x_station_circum4_5_mid(i,j-1,s))+y_station_circum4_5_mid(i,j-1,s);
y_station_circum4_5(i,j,s)=alpha_station_4_5_m(i,j,s)*(-x_station_circum4_5_mid(i,j,s)+x_station_cir- Mehrdad Jafarpourhello, how can i get completly this code ? [email protected]
- Mehrdad Jafarpourpls help me i need this code , pls gmail me [email protected]
- Cfan0804Nice work! Could you send me the complete code to me? Gmail: [email protected]