Invalid expression. Check for missing multiplication operator, missing or unbalanced delimiters, or other syntax error.

M_0= Mach00
P_0= Ambien pressure(P0 in kPa))101 %kPascals
T_0= Ambient temperature (T0 in K)298 %Kelvin
G=   Gamma1.4
cp=  Specific heat1004 %J/kg*K
m_0= Mass flow rate100 %kg/s
A_1= Alpha 10
OPR= Overall pressure ratio20
M_1= Design Mach0.5
U_t= Tip velocity (U in m/s)450
R_d= Degree of reaction (R)0.5
s=   Pitchline Solidity1
S=   Stator Solidity1.5
i=   Angle of incidence3.7
t_c= Thickness to chord ratio0.065
rho= Density at sea level1.225
dv=  Dynamic viscosity at sea level0.00001789

rht=Hub-to-tip ratio0.50.30.50.5
dh= De haller criterion0.720.720.750.72
e_c=Polytropic efficiency0.90.850.920.9


Tt_1 = T_0;
Pt_1 = P_0;

T_1 = (Tt_1/(1+((G-1)/2)*(M_1^2)))
a_1 = sqrt((G-1)*cp*T_1)
Cz_1= M_1 * a_1
Cz_2 = Cz_1
P_1 = (Pt_1/(1+((G-1)/2)*(M_1^2)^(G/G-1)))
R = ((cp*(G-1))/G)
Rm= .5
pa = (P_1/(R*T_1))
A1= (m_0/(pa*Cz_1))
rt_1 = sqrt(A1/ 3.14*(1-rht^2))
rh_1 = (rt_1*.5)
Pt_2 = 20*Pt_1
tau = (OPR^((G-1)/(e_c*G)))
Tte = Tt_1*tau




Te = (Tte-((Cz_1^2)/(2*cp)))
a_e = sqrt((G-1)*cp*Te)
M_2 = (Cz_1/a_e)
P_2 = (Pt_2/(1+((G-1)/2)*(M_2)^2)^(G/(G-1))) %Exit pressure
p2 = (P_2/(R*Te))
A_2 = (m_0/(p2*Cz_2))

Blade height design

 % Constant tip radius
rh_2 = sqrt(rt_1^2-(A_2/3.14))
rm_a = ((rh_2 + rt_1)/2) 
B_h1 = (A_2/(2*3.14*rm_a))

% Constant pitchline radius
rm_b = ((rh_1 + rt_1)/2)
B_h2 = (A_2/(2*3.14*rm_b))

%Constant hub radius
rt_2 = sqrt(((rh_1^2)+(A_2/3.14)))
rm_c = ((rh_1 + rt_2)/2) 
B_h3 = (A_2/ (2*3.14*rm_c))

Pitchline Calculations

rm_1 = ((rh_1 + rt_1)/2)
w = (U_t/rt_1) %rad/s
V_m = U_t*(rm_1/rt_1)
B1 = atand(V_m/Cz_2)
V_r= sqrt(V_m^2 + Cz_2^2)
M_m = (V_r/a_1)
C01= Cz_1* tan(A_1)
C02= ((2*V_m*Rm)-(1+C01))



%To prove design choice is feasible
W_1m=((V_m^2)+(Cz_2^2))^.5
W_2m = Cz_2
DCt= (C01+C02)
Dr_m = 1- (W_2m/W_1m) + ((abs(DCt))/(2*s*W_1m))

if Dr_m > .6
    x = uifigure
uialert (x, 'Enter a new degree of reaction', 'Unacceptable')
end

Validate with De Haller criterion

W2m = dh*W_1m
B2 = acosd(W_2m/W2m) %Beta2
W02 = (Cz_2*tand(B2))  %Relative swirl at rotor exit
C02m = (V_m-W02)  %Absolute swirl at rotor exit
Ym = (C02m/V_m)   %Stage loading factor
Alpha_2 = atand(C02m/Cz_2) % Absolute flow angle at rotor exit 
C2m = sqrt((Cz_2^2)+(C02m^2))   %Absolute velocity at rotor exit
cr = Cz_2/C2m

if cr < .72; disp 'Error'
end

Dr = 1- (W2m/W_1m) + ((abs(C02m-C01))/(2*s*W_1m)) %D factor rotor
Ds = 1- (Cz_2/C2m) + ((abs(C02m-C01))/(2*S*C2m)) %D factor stator
Rm = 1- ((C01+C02m)/(2*V_m)) %Pitchline degree of reaction

%Rotor blade design at pitchline
%To calculate chord lenght, calculate density and kinematic viscosity at
%16 Km altitude highest

rho1 = (.136)*(pa) %Density at 16km
v = (dv/rho) %Kinematic viscosity
v1 = v*5.849 %Kinematic viscocity at 16 Km
Re = 300000
cl = ((Re*v1)/W_1m) %meters


d = ((abs(B2-B1))/(4*(sqrt(s)))+2) %Deviation angle
k1m = B1 - i %Blade leading edge angle
k2m = B2 - d %Blade trailing edge angle
psi = k1m - k2m %Camber angle
Y = (B1 - (2/psi) -i) %Stagger angle
N = ceil(((2*3.14*rm_1)/cl))
% Hay que encontrar como convertir cl a cm para que te de 33

%Stator blade design at pitchline
Tt2m = (Tt_1 + ((V_m*(C02m-C01))/cp))  %Total temperature 2
T2m = Tt2m - ((C2m^2)/(2*cp)) % Static temperature
a_2 = sqrt(G*R*T2m); %Speed of sound 
M_2 = C2m/a_2 %Mach at strator entrance

Choose a Naca series 65 cascade profile

NACA 65 Series cascade properties

Stall Margin

DT_t = Tt2m-Tt_1   %Total stagnation temperature rise for compressor
DT = Tte - Tt_1 %Temperature rise
Ns = ceil(DT/DT_t)

[Voss]

Where does the error happen? The code runs fine for me.

[Jan]

I cannot find a problem also.

Answer 1

[Nipun]

Hi Jamyl,

I understand that you are facing an error while running the shared code snippet in MATLAB R2021b.

I tried running the shared code in MATLAB R2021b and it ran to completion without any errors. It will be helpful if you can point to the line numbers or a function where the error occurs.

Regards,

Nipun