How to assign multiple variables

Greetings everyone, Im currently working on this coding but seems to have a bit problem. So, basically im trying to code so that the values of my P01 variable are assigned to each M1 and M2. For different P01 value, the M1 and M2 would yield different values as per calculated using their respective equations However, i dont know whats wrong with my codes but it seems matlab only calculated on the first values of P01 only. Any help would be very much appreciated.

The value of P01 are as follows:

[124638.717 126432.080 128225.442 130018.805 131812.168 135398.894 137192.257 142572.345 142572.345 147952.434 155125.885 165886.062 167679.425 183819.690 201753.319 225067.035 234033.850 239413.938 248380.752 253760.841 263445.000 273577.500 283710.000 293842.500 303975.000 309041.250 310054.500 310561.125 311067.750 311574.375 314107.500]

%% Calculation
%Station 1
R = 287; %Gas constant
T01= 339; %Stagnation Temperature or total inlet temperature
Patm = 101325;
P01 ='VarName10';  %Total Inlet Pressure
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
M2 = 0.3 %guess value

%Mach Number 1 (Guess)
k=0;
for k= 1:length(P01)
    M1 = 0.1 %guess value of M1
    M2 = M2+0.001;
    
    Delta_Ctheta2=0.1;
while Delta_Ctheta2<-0.001 | Delta_Ctheta2>0.001
    if Delta_Ctheta2<-0.001; 
        M1=M1-0.0001;
    elseif Delta_Ctheta2>0.001;
        M1=M1+0.0001;

        end
    
    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

%% Station 2
%BOUNDARY CONDITION @ STATION 2
T02 = T01; %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

%CALCULATION @ STATION 2---------------------------------------------------

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)); 
Rad = alpha2 * pi/180;
C2 = M2.*(kc*R*T2).^0.5; %absolute velocity inlet volute
Ctheta2_inletrotor = C2.*sin(Rad);
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  
while delta_mdot <-0.001 | delta_mdot > 0.001
    if delta_mdot <-0.001; 
        M2=M2-0.0001;
        
    elseif delta_mdot >0.001
        M2=M2+0.0001;
        
        if mdot2_calculation == mdot2
           break
        end
    end   
    %Absolute velocity inlet rotor
alpha2 = (atand((rho2./rho1).*((A2/rrms2)./( A1/r1)).*SW)); 
Rad = alpha2 * pi/180;
C2 = M2.*(kc*R*T2).^0.5; %absolute velocity inlet volute
Ctheta2_inletrotor = C2.*sin(Rad);
Cm2 = C2.* cos(Rad);  
mdot2=mdot1;
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  

end

C1 = M1.*(kc*R*T1).^0.5;
Ctheta2_exitvolute = C1.*SW.*(r1 / rrms2 );
C2 = M2.*(kc*R*T2).^0.5;
Ctheta2_inletrotor = C2.*sin(Rad);
Delta_Ctheta2= Ctheta2_exitvolute - Ctheta2_inletrotor;

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

%STORED PARAMETERS---------------------------------------------------------
Ctheta2_exitvolute100(k,1)=Ctheta2_exitvolute
Ctheta2_inletrotor100(k,1)=Ctheta2_inletrotor
k=k+1;
break
end
end

[Torsten]

Do you try to determine two unknows (M1, M2) with one equation (Ctheta2_exitvolute==Ctheta2_inletrotor) ?

This will not work. You will have to prescribe one of them (either M1 or M2) and determine the other for a fixed inlet pressure P01.

Why don't you use the code where "fsolve" is invoked ?

[Danish Iqhwan Rohaizad]

The two unknowns are given an initial guess, however, what im trying to achieve is that the codes will keep on iterating using the value of the initial guess and will keep on iterating until achieving a value that does not just satisfy Delta_Ctheta2 = 0 but other Delta's as well(which is not provided on the above codes). I was thinking something like a while loop in loop

[Torsten]

For your own benefit, write down the variables you want to solve for and the equations that have to be fulfilled.

The number of both should be equal.

Then call "fsolve" prescribing the equations which are to be used to solve for the variables.

Your changes to M1 and M2 in the code above lead to nowhere because you don't update your M1 and M2 variables purposefully and directed towards a root like it's done with Newton's method, secant method or something similar.

Answer 1

[Ishu]

Hi Danish,

I ran your code at my end and tried to reprocude the error. In the code you are calculating "Ctheta2_exitvolute" and "Ctheta2_inletrotor" and storing these calculated values in "Cthete2_exitvolute100(k,1)" and "Ctheta2_inletrotor100(k,1)" respectively.

As you can see your "Ctheta2_intetrotor" is a 1x9 vector and "Ctheta2_inletrotor100(k,1)" can store only one variable at each iteration. So that is why it is giving you an error.

I am assuming that you want to store one calculated variable each iteration, hence you can make these changes in your code:

rho01 = P01(k)./(R*T01); %Stagnation density   % change P01  -> P01(k)

P1 = P01(k)./(1 + kc0 *M1.^2).^kc2; %static pressure inlet volute  % change P01  -> P01(k)

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

Hope it helps.