Generate Random Numbers to satisfy a variable on the left and right hand side

Greetings and Hi !

Supposedly if i would like to generate a random number for M1 and M2 which would then be inputted into the coming equations.

R = 287; %Gas constant
T01= 339; %Stagnation Temperature or total inlet temperature
Patm = 101325;
P01 = [124638.717 126432.080 128225.442 130018.805 131812.168 135398.894 137192.256637168 142572.345132743 142572.345132743 147952.433628319 155125.884955752 165886.061946903 167679.424778761 183819.690265487 201753.318584071 225067.035398230 234033.849557522 239413.938053097 248380.752212389 253760.840707965 263445 273577.500000000 283710 293842.500000000 303975 309041.250000000 310054.500000000 310561.125000000 311067.750000000 311574.375000000 314107.500000000];
A1 = 0.0025686; %Area inlet volute
kc = 1.4; %specific heat ratio
kc0 = 0.2; %specific heat ratio at total state
kc1 = 2.498; %specific heat ratio at inlet volute
kc2 = 3.498; %specific heat ratio at inlet rotor

%Mach Number 1 (Guess)
M1 = rand(1,31);% enter the vector M1
rho01 = P01./(R*T01);
rho1 = rho01./(1 + kc0 * M1.^2).^kc1;

for iter_1 = 1:length(M1)
    rho01 = P01./(R*T01); %Stagnation density
    rho1 = rho01./(1 + kc0 * M1.^2).^kc1; %flow density inlet volute
    T1 = T01./(1 + kc0 * M1.^2); %static temperature inlet volute
    P1 = P01./(1 + kc0 *M1.^2).^kc2; %static pressure inlet volute
    C1 = M1.*(kc*R*T1).^0.5; %absolute velocity inlet
    mdot1 = rho1.*A1.*C1; %mass flow rate inlet volute
end

%% Station 2
M2 = rand(1,31);% enter the value M2
%M2 = [0.320 0.333 0.346 0.359 0.371 0.395 0.407 0.441 0.441 0.473 0.514 0.571 0.580 0.660 0.744 0.848 0.888 0.912 0.953 0.977 1.023 1.074 1.130 1.192 1.271 1.324 1.336 1.344 1.351 1.359 1.411]

T02 = 339; %Stagnation temperature or total temperature at rotor inlet
N = 47503; %Rotor speed
omega = (2*pi/60)*N; %Rotor 100% speed angular velocity
ktp_volute= 0.9; %Coeeficient Total Pressure Loss Volute
SW = 0.95; % Swirl Coefficient
B2 = 1; %Inlet rotor blockage factor
r1 = 0.0739; %Inlet volute radius
r2max = 0.047576; %Rotor shroud tip radius at leading edge
r2min = 0.036006; %Rotor hub radius at leading edge
rrms2 = (r2max + r2min)/2; %Average between r2max and r2min
U2 = rrms2 * omega; %Blade speed rotor tip

for iter_2 = 1:length(M2)
    
Ctheta2_exitvolute = C1.*SW.*(r1 / rrms2 ); %Tangential absolute velocity inlet rotor
P2 = P01./(((1 + kc0 * M2.^2).^kc2).*(1+ktp_volute)-ktp_volute); %Static pressure inlet rotor
T2 = T02./(1 + kc0 * M2.^2); %static temperature inlet rotor
P02 = P2.*(1 + kc0 * M2.^2).^kc2; %Total pressure inlet rotor
rho2 = P2./(R*T2); %Stagnation density

% using continuity equation, mass flow rate inlet rotor = mass flow rate inlet volute
% mdot1 = mdot2
mdot2 = mdot1; %mass flow inlet rotor
leh = 0.018; %rotor leading height
A2 = pi * (r2max +r2min) * leh;%Area Inlet Rotor

%Absolute velocity inlet rotor
alpha2 = (atand((rho2./rho1).*((A2/rrms2)./( A1/r1)).*SW)); 

C2 = M2.*(kc*R*T2).^0.5; %absolute velocity inlet volute
Rad = deg2rad(alpha2);
Ctheta2_inletrotor = C2.*sin(Rad); %Tangential absolute velocity inlet rotor
Ctheta2_exitvolute = Ctheta2_inletrotor ; %DeltaCtheta2

end

However what i need from M1 and M2 is to generate numbers in which when inputted will give the same value between Ctheta2_exitvolute and Ctheta2_inletrotor

[Torsten]

Then prescribe M1 and solve for M2 (or vice versa) using the condition

Ctheta2_exitvolute = Ctheta2_inletrotor

[Danish Iqhwan Rohaizad]

Hi there Torsten ! Yes I have assigned both of my Ctheta's to each other. However, what Im seeking help is the coding for variable M1 and M2. How can I code these two variables to generate a random number that would allow both of the Ctheta's to be of the same value ?

[Torsten]

That's what I tried to explain.

The relation

Ctheta2_exitvolute = Ctheta2_inletrotor   (1)

implicitly defines a relation between M1 and M2.

Thus given M1 (or M2), M2 (or M1) is not arbitrary random, but fixed by relation (1).

[Torsten]

@Danish Iqhwan Rohaizad

Should the relation between M1 and M2 to satisfy

Ctheta2_exitvolute = Ctheta2_inletrotor

depend on P01 ?

[Danish Iqhwan Rohaizad]

Actually if you look through my coding, you can see that certain variables are used in another. So in the case of these two variables (M1 and M2), they are used in the equation C1 and C2 respectively. The values obtain from the random generated in M1 and M2 would lead to DeltaP = Ctheta_exitvolute - Ctheta_inletrotor having a value of 0. That is what Im trying to seek

Answer 1

[Torsten]

Here is a code with which you get M2 versus M1 curves for the vector P01 of pressures for which the relation

Ctheta_exitvolute - Ctheta_inletrotor = 0

is satisfied.

Interestingly, the curves don't depend on P01. That's why I asked whether the M1-M2 relation should depend on P01.

Be a little patient - the code takes quite a while.

P01 = [124638.717 126432.080 128225.442 130018.805 131812.168 135398.894 137192.256637168 142572.345132743 142572.345132743 147952.433628319 155125.884955752 165886.061946903 167679.424778761 183819.690265487 201753.318584071 225067.035398230 234033.849557522 239413.938053097 248380.752212389 253760.840707965 263445 273577.500000000 283710 293842.500000000 303975 309041.250000000 310054.500000000 310561.125000000 311067.750000000 311574.375000000 314107.500000000];
m1 = 0:0.001:0.6;
m2 = 0:0.001:5;
for k=1:numel(P01)
    p01 = P01(k);
    for i=1:numel(m1)
        M1 = m1(i);
        resl = fun(m2(1),M1,p01);
        flag = 0;
        for j=1:numel(m2)-1
            resr = fun(m2(j+1),M1,p01);
            if resl*resr <=0
                sol_m2(i,k) = (m2(j)+m2(j+1))/2;
                flag=1;
                break
            end
            resl = resr;
        end  
        if flag==0 
            sol_m2(i,k)=NaN;
        end
        %if mod(i,100)==0
        % i
        %end
    end
    k
end
plot(m1,sol_m2)
end


function res = fun(x,M1,P01)

R = 287; %Gas constant
T01= 339; %Stagnation Temperature or total inlet temperature
Patm = 101325;
%P01 = [124638.717 126432.080 128225.442 130018.805 131812.168 135398.894 137192.256637168 142572.345132743 142572.345132743 147952.433628319 155125.884955752 165886.061946903 167679.424778761 183819.690265487 201753.318584071 225067.035398230 234033.849557522 239413.938053097 248380.752212389 253760.840707965 263445 273577.500000000 283710 293842.500000000 303975 309041.250000000 310054.500000000 310561.125000000 311067.750000000 311574.375000000 314107.500000000];
A1 = 0.0025686; %Area inlet volute
kc = 1.4; %specific heat ratio
kc0 = 0.2; %specific heat ratio at total state
kc1 = 2.498; %specific heat ratio at inlet volute
kc2 = 3.498; %specific heat ratio at inlet rotor

rho01 = P01./(R*T01);
rho1 = rho01./(1 + kc0 * M1.^2).^kc1;

T1 = T01./(1 + kc0 * M1.^2); %static temperature inlet volute
P1 = P01./(1 + kc0 *M1.^2).^kc2; %static pressure inlet volute
C1 = M1.*(kc*R*T1).^0.5; %absolute velocity inlet
mdot1 = rho1.*A1.*C1; %mass flow rate inlet volute

T02 = 339; %Stagnation temperature or total temperature at rotor inlet
N = 47503; %Rotor speed
omega = (2*pi/60)*N; %Rotor 100% speed angular velocity
ktp_volute= 0.9; %Coeeficient Total Pressure Loss Volute
SW = 0.95; % Swirl Coefficient
B2 = 1; %Inlet rotor blockage factor
r1 = 0.0739; %Inlet volute radius
r2max = 0.047576; %Rotor shroud tip radius at leading edge
r2min = 0.036006; %Rotor hub radius at leading edge
rrms2 = (r2max + r2min)/2; %Average between r2max and r2min
U2 = rrms2 * omega; %Blade speed rotor tip
Ctheta2_exitvolute = C1.*SW.*(r1 / rrms2 ); %Tangential absolute velocity inlet rotor
P2 = P01./(((1 + kc0 * x.^2).^kc2).*(1+ktp_volute)-ktp_volute); %Static pressure inlet rotor
T2 = T02./(1 + kc0 * x.^2); %static temperature inlet rotor


P02 = P2.*(1 + kc0 * x.^2).^kc2; %Total pressure inlet rotor
rho2 = P2./(R*T2); %Stagnation density
% using continuity equation, mass flow rate inlet rotor = mass flow rate inlet volute
% mdot1 = mdot2
mdot2 = mdot1; %mass flow inlet rotor
leh = 0.018; %rotor leading height
A2 = pi * (r2max +r2min) * leh;%Area Inlet Rotor
%Absolute velocity inlet rotor
alpha2 = (atand((rho2./rho1).*((A2/rrms2)./( A1/r1)).*SW)); 
C2 = x.*(kc*R*T2).^0.5; %absolute velocity inlet volute
Rad = deg2rad(alpha2);
Ctheta2_inletrotor = C2.*sin(Rad); %Tangential absolute velocity inlet rotor
res = Ctheta2_exitvolute - Ctheta2_inletrotor;
end

[Danish Iqhwan Rohaizad]

Hi there Torsten, though i appreciate for your coding however im not sure whether yours are relevant to mine since as for my codes, I need every variable to show its values as it will be used for the calculations of the next stage (Not included in the above codes). And as for the updated codes, I decided to use the convergence concept in order to achieve the same value for Ctheta2_exitvolute and Ctheta2_inletrotor and focused at 1 point only.

%% Calculation

%Station 1

R = 287; %Gas constant

T01= 339; %Stagnation Temperature or total inlet temperature

Patm = 101325;

P01 = 124638.717;

A1 = 0.0025686; %Area inlet volute

kc = 1.4; %specific heat ratio

kc0 = 0.2; %specific heat ratio at total state

kc1 = 2.498; %specific heat ratio at inlet volute

kc2 = 3.498; %specific heat ratio at inlet rotor

%% Station 2

M2 = rand(1)

T02 = 339; %Stagnation temperature or total temperature at rotor inlet

N = 47503; %Rotor speed

omega = (2*pi/60)*N; %Rotor 100% speed angular velocity

ktp_volute= 0.9; %Coeeficient Total Pressure Loss Volute

SW = 0.95; % Swirl Coefficient

B2 = 1; %Inlet rotor blockage factor

r1 = 0.0739; %Inlet volute radius

r2max = 0.047576; %Rotor shroud tip radius at leading edge

r2min = 0.036006; %Rotor hub radius at leading edge

rrms2 = (r2max + r2min)/2; %Average between r2max and r2min

U2 = rrms2 * omega; %Blade speed rotor tip

for iter_2 = 1:length(M2)

P2 = P01./(((1 + kc0 * M2.^2).^kc2).*(1+ktp_volute)-ktp_volute); %Static pressure inlet rotor

T2 = T02./(1 + kc0 * M2.^2); %static temperature inlet rotor

P02 = P2.*(1 + kc0 * M2.^2).^kc2; %Total pressure inlet rotor

rho2 = P2./(R*T2); %Stagnation density

leh = 0.018; %rotor leading height

A2 = pi * (r2max +r2min) * leh;%Area Inlet Rotor

m1 = rand() %Changing M1 to get value DeltaCtheta2

for iter_1 = 1:length(m1)

rho01 = P01./(R*T01); %Stagnation density

rho1 = rho01./(1 + kc0 * m1.^2).^kc1; %flow density inlet volute

T1 = T01./(1 + kc0 * m1.^2); %static temperature inlet volute

P1 = P01./(1 + kc0 *m1.^2).^kc2; %static pressure inlet volute

C1 = m1.*(kc*R*T1).^0.5; %absolute velocity inlet

mdot1 = rho1.*A1.*C1; %mass flow rate inlet volute

Ctheta2_exitvolute = C1.*SW.*(r1 / rrms2 ); %Tangential absolute velocity inlet rotor

% using continuity equation, mass flow rate inlet rotor = mass flow rate inlet volute

% mdot1 = mdot2

mdot2 = mdot1; %mass flow inlet rotor

%Absolute velocity inlet rotor

alpha2 = (atand((rho2./rho1).*((A2/rrms2)./( A1/r1)).*SW));

Rad = deg2rad(alpha2);

%Absolute velocity inlet rotor

C2 = M2(iter_2).*(kc*R*T2).^0.5;

Ctheta2_inletrotor = C2.*sin(Rad);

%DeltaCtheta

%syms Ctheta2_exitvolute Ctheta2_inletrotor

%DeltaCtheta2 = solve([Ctheta2_exitvolute - Ctheta2_inletrotor == 0],Ctheta2_exitvolute)

tolerance = Ctheta2_exitvolute - Ctheta2_inletrotor

while isempty(tolerance) || tolerance == 0 || tolerance > 1

C1 = m1.*(kc*R*T1).^0.5 ;

Ctheta2_exitvolute = (m1.*(kc*R*T1).^0.5).*SW.*(r1 / rrms2 ) ;

accuracy = m1; %abs(Ctheta2_exitvolute - Ctheta2_inletrotor)

while tolerance ~= 0

%m1 = m1 + 0.01; %want to continue iterate until achieve a value that corresponds with DeltaCtheta2

C1 = m1.*(kc*R*T1).^0.5;

accuracy = abs(m1);

if (Ctheta2_exitvolute < Ctheta2_inletrotor)

m1 = (m1 + (m1+0.01))/2

Ctheta2_exitvolute = m1.*(kc*R*T1).^0.5

fprintf('Criteria NOT satisfied with value Mach Number 1 of %f\n', m1);

elseif Ctheta2_exitvolute > Ctheta2_inletrotor

m1 = (m1 + (m1 - 0.01))/2

Ctheta2_exitvolute = m1.*(kc*R*T1).^0.5;

fprintf('Criteria NOT satisfied with value Mach Number 1 of %f\n', m1);

else abs(Ctheta2_exitvolute == Ctheta2_inletrotor) %<= eps(Cheta2_inletrotor)

%Continue

fprintf('Criteria satisfied with value Mach Number 1 of %f\n', m1);

break

end

C1

m1

tolerance

end

if Ctheta2_exitvolute == Ctheta2_inletrotor

fprintf('Delta_Ctheta2 satisfied with value of %f\n', Ctheta2_exitvolute, Ctheta2_inletrotor);

end

end

end

end

However, I am still unable to achieve the desired value 0 and was hoping the help of this community to mitigate on my problem

Much Thanks,

Danish

[Torsten]

The task remains the same:

Given M2, you search for m1 to make Ctheta2_exitvolute == Ctheta2_inletrotor.

Do we agree so far ?

[Danish Iqhwan Rohaizad]

yes

[Walter Roberson]

However, I am still unable to achieve the desired value 0

It is possible that your goal cannot be reached. You are using floating point expressions, and calculating values in two different ways, but you are expecting the two values to be bit for bit identical . That is what it means to compare A == B -- that A and B are bit for bit identical . Exact equality.

For example, adjacent representable numbers

format long g
A = hex2num('3fe63c093c802fe9')

A = 
         0.694828622975817

B = hex2num('3fe63c093c802fe8')

B = 
         0.694828622975817

fprintf('A = %.15g\n', A)

A = 0.694828622975817

fprintf('B = %.15g\n', B)

B = 0.694828622975817

They look the same

A == B

ans = logical
   0

But they compare different

A-B

ans = 
      1.11022302462516e-16

There is an actual difference

fprintf('%.99g\n', A)

0.69482862297581704513760314512182958424091339111328125

fprintf('%.99g\n', B)

0.69482862297581693411530068260617554187774658203125

Because of this, you should never be using == on floating point numbers that are calcuated in different ways. Instead you should be calculating

if abs(A - B) < SomeTolerance

Answer 2

[Torsten]

ok.

This code should return M1, given M2 and a starting guess for M1, such that the pair (M1,M2) satisfies

Ctheta2_exitvolute == Ctheta2_inletrotor

M2 = 2.0;
M10 = 0.5;
M1 = fsolve(@(x)fun(x,M2),M10)
fun(M1,M2)

function res = fun(m1,M2)
%% Calculation
%Station 1
R = 287; %Gas constant
T01= 339; %Stagnation Temperature or total inlet temperature
Patm = 101325;
P01 = 124638.717;
A1 = 0.0025686; %Area inlet volute
kc = 1.4; %specific heat ratio
kc0 = 0.2; %specific heat ratio at total state
kc1 = 2.498; %specific heat ratio at inlet volute
kc2 = 3.498; %specific heat ratio at inlet rotor
%% Station 2
%M2 = rand(1) 
T02 = 339; %Stagnation temperature or total temperature at rotor inlet
N = 47503; %Rotor speed
omega = (2*pi/60)*N; %Rotor 100% speed angular velocity
ktp_volute= 0.9; %Coeeficient Total Pressure Loss Volute
SW = 0.95; % Swirl Coefficient
B2 = 1; %Inlet rotor blockage factor
r1 = 0.0739; %Inlet volute radius
r2max = 0.047576; %Rotor shroud tip radius at leading edge
r2min = 0.036006; %Rotor hub radius at leading edge
rrms2 = (r2max + r2min)/2; %Average between r2max and r2min
U2 = rrms2 * omega; %Blade speed rotor tip
%for iter_2 = 1:length(M2)

P2 = P01./(((1 + kc0 * M2.^2).^kc2).*(1+ktp_volute)-ktp_volute); %Static pressure inlet rotor
T2 = T02./(1 + kc0 * M2.^2); %static temperature inlet rotor
P02 = P2.*(1 + kc0 * M2.^2).^kc2; %Total pressure inlet rotor
rho2 = P2./(R*T2); %Stagnation density
leh = 0.018; %rotor leading height
A2 = pi * (r2max +r2min) * leh;%Area Inlet Rotor
%m1 = rand() %Changing M1 to get value DeltaCtheta2
%for iter_1 = 1:length(m1)

rho01 = P01./(R*T01); %Stagnation density
rho1 = rho01./(1 + kc0 * m1.^2).^kc1; %flow density inlet volute
T1 = T01./(1 + kc0 * m1.^2); %static temperature inlet volute
P1 = P01./(1 + kc0 *m1.^2).^kc2; %static pressure inlet volute
C1 = m1.*(kc*R*T1).^0.5; %absolute velocity inlet
mdot1 = rho1.*A1.*C1; %mass flow rate inlet volute
Ctheta2_exitvolute = C1.*SW.*(r1 / rrms2 ); %Tangential absolute velocity inlet rotor

% using continuity equation, mass flow rate inlet rotor = mass flow rate inlet volute
% mdot1 = mdot2
mdot2 = mdot1; %mass flow inlet rotor
%Absolute velocity inlet rotor
alpha2 = (atand((rho2./rho1).*((A2/rrms2)./( A1/r1)).*SW)); 
Rad = deg2rad(alpha2);

%Absolute velocity inlet rotor
%C2 = M2(iter_2).*(kc*R*T2).^0.5; 
C2 = M2.*(kc*R*T2).^0.5; 
Ctheta2_inletrotor = C2.*sin(Rad);

res = Ctheta2_exitvolute - Ctheta2_inletrotor;
%DeltaCtheta
%syms Ctheta2_exitvolute Ctheta2_inletrotor 
%DeltaCtheta2 = solve([Ctheta2_exitvolute - Ctheta2_inletrotor == 0],Ctheta2_exitvolute)

%tolerance = Ctheta2_exitvolute - Ctheta2_inletrotor 

%while isempty(tolerance) || tolerance == 0 || tolerance > 1

%C1 = m1.*(kc*R*T1).^0.5 ;

%Ctheta2_exitvolute = (m1.*(kc*R*T1).^0.5).*SW.*(r1 / rrms2 ) ;

%accuracy = m1; %abs(Ctheta2_exitvolute - Ctheta2_inletrotor) 

%while tolerance ~= 0

%m1 = m1 + 0.01; %want to continue iterate until achieve a value that corresponds with DeltaCtheta2 
%C1 = m1.*(kc*R*T1).^0.5; 
%accuracy = abs(m1);

%if   (Ctheta2_exitvolute < Ctheta2_inletrotor) 

%m1 = (m1 + (m1+0.01))/2
%Ctheta2_exitvolute = m1.*(kc*R*T1).^0.5
%fprintf('Criteria NOT satisfied with value Mach Number 1 of %f\n',   m1);

%elseif Ctheta2_exitvolute > Ctheta2_inletrotor

%m1 = (m1 + (m1 - 0.01))/2
%Ctheta2_exitvolute = m1.*(kc*R*T1).^0.5;
%fprintf('Criteria NOT satisfied with value Mach Number 1 of %f\n', m1);

%else abs(Ctheta2_exitvolute == Ctheta2_inletrotor) %<= eps(Cheta2_inletrotor)
%Continue
%fprintf('Criteria satisfied with value Mach Number 1 of %f\n', m1);

%break

%end
%C1
%m1
%tolerance
%end
%if Ctheta2_exitvolute == Ctheta2_inletrotor  
%fprintf('Delta_Ctheta2 satisfied with value of %f\n', Ctheta2_exitvolute, Ctheta2_inletrotor); 
% end
%end
%end
%end
end

[Danish Iqhwan Rohaizad]

I cannot express how thankful I am to you Torsten ! Mind if i ask on how to proceed if lets say there are other M variables (M6_rel & M7) which are at different stages ?

Should I create another function in the same script ? Any help would be very much appreciated

[Torsten]

Mind if i ask on how to proceed if lets say there are other M variables (M6_rel & M7) which are at different stages ?

Should I create another function in the same script ?

In the above case, you have 1 equation (Ctheta2_exitvolute - Ctheta2_inletrotor) to determine 1 unknown (M1).

If you have another equation for another unknown, you can either integrate this equation into the same script if the new unknown somehow depends on M1 and M2, or - if it is independent of M1 and M2 - use a second script.

If you want to determine both unknowns simultaneously, the call to fsolve would then comprise an initial vector with two elements:

M = fsolve(@(x)fun(x,M2,M7),[M10,M6_rel0])

and "fun" would return a vector "res" also with two elements:

function res = fun(m,M2,M7)
  M1 = m(1);
  M6_rel = m(2);
  ...
  res(1) = ...;
  res(2) = ...;
end

[Danish Iqhwan Rohaizad]

Hi there Torsten, thanks again for your help. For my codes, I decied to put two functions in a single script. However due to this, MATLAB only considers on the first function and not on the second. Attached below is the updated codes.Any advise is much appreciated,1

P01 = 124638.717;

A1 = 0.0025686; %Area inlet volute

kc = 1.4; %specific heat ratio

kc0 = 0.2; %specific heat ratio at total state

kc1 = 2.498; %specific heat ratio at inlet volute

kc2 = 3.498; %specific heat ratio at inlet rotor

%Station2

T02 = 339; %Stagnation temperature or total temperature at rotor inlet

N = 47503; %Rotor speed

omega = (2*pi/60)*N; %Rotor 100% speed angular velocity

ktp_volute= 0.9; %Coeeficient Total Pressure Loss Volute

SW = 0.95; % Swirl Coefficient

B2 = 1; %Inlet rotor blockage factor

r1 = 0.0739; %Inlet volute radius

P2 = P01./(((1 + kc0 * M2.^2).^kc2).*(1+ktp_volute)-ktp_volute); %Static pressure inlet rotor

T2 = T02./(1 + kc0 * M2.^2); %static temperature inlet rotor

P02 = P2.*(1 + kc0 * M2.^2).^kc2; %Total pressure inlet rotor

rho2 = P2./(R*T2); %Stagnation density

leh = 0.018; %rotor leading height

r2max = 0.047576; %Rotor shroud tip radius at leading edge

r2min = 0.036006; %Rotor hub radius at leading edge

rrms2 = (r2max + r2min)/2; %Average between r2max and r2min

r6max = 0.039325; %rotor shroud tip radius at trailing edge

r6min = 0.015535; %rotor hub radius at trailing edge

rrms6 = (r6min + r6max )/2; %average between r6max and r6min

U2 = rrms2 * omega; %Blade speed rotor tip

U6 = rrms6 * omega; %Blade speed exit rotor

B6 = 1.00; %Exit rotor blockage factor

Zr = 12; %Number of Blades

beta_b2 = 20 ; %Blade angle at inlet rotor (in degrees)

beta_6 = -40 * (pi/180) %Relative Flow Angle at exit rotor (degree to rad)

A6_C = pi * (r6max^2 - r6min^2); %Area Exit Rotor

cr = 0.2; %coefficient loss for Mach number < 0.75

A2 = pi * (r2max +r2min) * leh;%Area Inlet Rotor

rho01 = P01./(R*T01); %Stagnation density

rho1 = rho01./(1 + kc0 * m1.^2).^kc1; %flow density inlet volute

T1 = T01./(1 + kc0 * m1.^2); %static temperature inlet volute

P1 = P01./(1 + kc0 *m1.^2).^kc2; %static pressure inlet volute

C1 = m1.*(kc*R*T1).^0.5; %absolute velocity inlet

mdot1 = rho1.*A1.*C1; %mass flow rate inlet volute

Ctheta2_exitvolute = C1.*SW.*(r1 / rrms2 ) %Tangential absolute velocity inlet rotor

% using continuity equation, mass flow rate inlet rotor = mass flow rate inlet volute

mdot2 = mdot1 %mass flow inlet rotor

%Absolute velocity inlet rotor

alpha2 = (atand((rho2./rho1).*((A2/rrms2)./( A1/r1)).*SW));

Rad = deg2rad(alpha2);

%Absolute velocity inlet rotor

C2 = M2.*(kc*R*T2).^0.5;

Ctheta2_inletrotor = C2.*sin(Rad)

res = Ctheta2_exitvolute - Ctheta2_inletrotor %DeltaCtheta2

Cm2 = C2.* cos(Rad);

mdot2_calculation = rho2.* (Ctheta2_inletrotor./tan(Rad)).* A2.* B2 %Mass flow rate at inlet rotor in terms of calculation

Delta_mdot = mdot2 - mdot2_calculation %Difference of mass flow rate at rotor

Theta2 = atand(Cm2./Ctheta2_exitvolute); %Angle opposite to alpha2

%Theta2 + alpha2

Theta2_alpha2 = Theta2 + alpha2;

%Relative velocity inlet rotor

x2 = U2.^2 + C2.^2 - 2.*U2.*C2.*cosd(Theta2);

W2 = sqrt(x2); %Relative Velocity Inlet Rotor

Wtheta2 = U2 - abs(Ctheta2_exitvolute); %Tangential relative velocity inlet rotor

% Recheck and Reuse

W2_rechecking = sqrt(Wtheta2.^2 + Cm2.^2); %Validating the value of W2 using hypotenuse

% Relative Flow Angle @ Inlet Rotor, beta2

beta2 = (Wtheta2<0) .* (abs(acosd(Cm2./W2_rechecking)) + (Wtheta2 > 0) .* (-acosd(Cm2./W2_rechecking))); % One-Line Statement

beta2_fcn = @(Wtheta2) (Wtheta2<0) .*(abs(acosd(Cm2./W2_rechecking)) + (Wtheta2 > 0) .* (-acosd(Cm2./W2_rechecking))); % Anonymous Function

M2_rel = W2_rechecking./ sqrt(kc * R * T2 ); %Relative Mach Number Inlet Rotor

T02_rel = T2.*(1 + kc0 * M2_rel.^2); %Relative Total Temperature Inlet Rotor

P02_rel = P2.*(1 + kc0 * M2_rel.^2).^kc2;

end

[Danish Iqhwan Rohaizad]

sorry this is the initial part of my code

M2 = 0.320;

M10 = rand(1);

M1 = fsolve(@(x)fun_1(x,M2),M10)

fun_1(M1,M2)

function main()

fun_1(m1,M2)

fun_2(m2,m6W)

end

function res = fun_1(m1,M2)

%% Calculation

%Station 1

R = 287; %Gas constant

T01= 339; %Stagnation Temperature or total inlet temperature

Patm = 101325;

[Torsten]

I also see only one function (fun_1). Where is the second (fun_2) ?

[Walter Roberson]

I decied to put two functions in a single script. However due to this, MATLAB only considers on the first function and not on the second.

MATLAB can find a named function under the following circumstances:

• the function name matches the name of a variable in the current workspace or the parameter list to the current function, and that variable is a function handle
• the current function is a nested function and the function name matches the name of a variable in a nesting function or the parameter list to a nesting function
• the function name matches the name of a function inside the current .m file (and that function is not in a deeper nesting level)
• the function name matches the name of a .m file somewhere on the MATLAB path (there are rules about which path entry takes priority)
• the "dominant" parameter is an object that has a method in its heirarchy with the name of the function
• a class has been loaded which has the name of the function as a static class method

Notice that nothing in this list includes the possibility of having a single external .m file that contains multiple functions that you intend to address directly by name. If you have main.m and fun_1.m then putting function fun_2 inside fun_1.m does not give MATLAB a way to know that it should look inside fun_1.m to find fun_2.m

You can put your fun_1 and fun_2 inside your main.m or you can put your fun_1 in fun_1.m and put your fun_2 in fun_2.m

[Danish Iqhwan Rohaizad]

Sorry there Torsten, this is the rest of my code

%% Station 3

function res2 = fun_2(m2,m6W)

M6_rel = rand(1);

M6rel0 = rand(1);

m6W = fsolve(@(x)fun_2(x,M6_rel),M6rel0);

fun_2(M6_rel,m6W)

%optimum incidence angle

i_opt = atand((-1.98 * tan(Rad))./(Zr - 1.98));

%Deviation from optimum incidence

diff = beta2 - i_opt - beta_b2;

Rad_diff = deg2rad(diff); %Diff in radians

%Absolute diff

abs_diff = abs(Rad_diff);

%coefficient incidence loss W2 if diff > pi/4

Kinc1 = 0.7.*(0.5 + (Rad_diff) - (pi/4));

%coefficient incidence loss W2 if diff < pi/4

Kinc2 = 0.200.*sin(Rad_diff).^2;

Rad_Kpass2 = deg2rad(beta2 - i_opt); %Relative Flow Angle @ Inlet Rotor - Optimum incidence angle (in radians)

Kpass2 = cr.*(cos(abs(Rad_Kpass2)).^2)./2; %coefficient internal passage loss W6

Kpass6 = cr/2; %Coefficient Internal Passage

T06_rel = T02_rel - kc0./(kc*R) * (U2.^2 - U6.^2); %Relative Total Temperature Exit Rotor

T6 = T06_rel./ (1+(kc0*m2.^2));

W6 = m2.*sqrt(kc * R * T6);

Wcr2 = sqrt((2 * kc * R * T02_rel)./(kc+1));

Wcr6 = Wcr2.* sqrt(T06_rel/T02_rel);

P06_idrel = P02_rel.* (T06_rel/T02_rel).^kc2;

% Ideal Relative Velocity @ Exit Rotor,W6_ideal

if diff > pi/4

W6_ideal = sqrt(((W6.^2).* (1+Kpass6)) + ((Kpass2 + Kinc1).* (W2_rechecking.^2)));

else diff < pi/4

W6_ideal = sqrt(((W6.^2).* (1+Kpass6)) + ((Kpass2 + Kinc2).* (W2_rechecking.^2)));

end

%Static pressure at exit rotor, P6

if abs_diff > pi/4

P6 = (1 - ((W6_ideal./Wcr6).^2).* ((kc-1)./(kc+1))).^(1/kc2).* P06_idrel;

else abs_diff < pi/4

P6 = (1 - ((W6_ideal./Wcr6).^2).* ((kc-1)./(kc+1))).^(1/kc2).* P06_idrel;

end

rho6 = P6./(R * T6); %flow density exit rotor

m6W = rho6.*A6_C.*(W6 * cos(beta_6)); %mass flow rate exit rotor - perpendicular to W

res2= m6W - mdot2_calculation %Delta mass flow rate

Delta_mdot_26W = res2

% using continuity equation, mass flow rate exit rotor (perpendicular to W)= mass flow rate exit rotor (perpendicular to Cm6)

% m6W = m6C

m6C = m6W;

Cm6 = m6C./(rho6 * A6_C);

Wm6 = Cm6; %meriodional relative velocity exit rotor

end

[Danish Iqhwan Rohaizad]

Hi there Walter, if lets say some of the variables in fun_2 is related to fun_1. How should I progress from there ?