Solving function in Matlab

Hi, I'm currently trying solve the shock tube problem with a matlab code I wrote. With several assumed initial conditions, my codes will calculate W3 and W3_P based on the shock Mach number Ms that I arbitrarily defined. The objective is to let W3 = W3_P, so I tried while loop and function handle but none worked. Any good suggestions?

function shock_3

Ms = 3; %This is a starting value

%State 1 conditions
P1 = 101.325;    % Assuming P1 = 1atm = 101.325kPa at SSL;
Ra = 287;        % R for air is 287;
T1 = 300;        % Assuming the room temperature is 27C
C1 = sqrt(1.4*Ra*T1);
Rho1 = P1/(Ra*T1);
ushock = Ms*C1;
u1 = ushock;

%State 5 conditions
P5 = 500*101.325;% storage pressure.
Rh = 4124;
T5 = 300;
C5 = sqrt(1.4*Rh*T5);
Rho5 = P5/(Rh*T5);
k = P5/(Rho5^(1.4));

P2 = P1*(2.8*(Ms^2-1)/2.4)+P1;
u2 = u1 - C1*2*(Ms^2-1)/(Ms*2.4);
W3 = ushock - u2;
P3 = P2;

Rho3 = (P3/k)^(1/1.4);
C3 = sqrt(1.4*P3/Rho3);
W3_P = 2*(C5-C3)/0.4;

%So obviously W3 is not the same with W3_P. I tried fsolve and a while loop but none worked. Any suggestions would help. Thanks in advance.

end

[darova]

Please attach original equations

Answer 1

[David Goodmanson]

Hi HD,

Never a bad idea to do a plot. I modified the code a bit, so it does vector in with Ms, vector out.

Could you explain the physical significance of Ra and Rh?

Ms = 3:.01:8;
P1 = 101.325;    % Assuming P1 = 1atm = 101.325kPa at SSL;
Ra = 287;        % R for air is 287;
T1 = 300;        % Assuming the room temperature is 27C
C1 = sqrt(1.4*Ra*T1);
Rho1 = P1/(Ra*T1);
ushock = Ms*C1;
u1 = ushock;
%State 5 conditions
P5 = 500*101.325;% storage pressure.
Rh = 4124;
T5 = 300;
C5 = sqrt(1.4*Rh*T5);
Rho5 = P5/(Rh*T5);
k = P5/(Rho5^(1.4));
P2 = P1*(2.8*(Ms.^2-1)/2.4)+P1;
u2 = u1 - C1*2*(Ms.^2-1)./(Ms*2.4);
W3 = ushock - u2;
P3 = P2;
Rho3 = (P3/k).^(1/1.4);
C3 = sqrt(1.4*P3./Rho3);
W3_P = 2*(C5-C3)/0.4;
plot(Ms,W3,Ms,W3_P)
grid on

You can read the value off the plot, or for a more accurate number:

Ms0 = fzero(@fun,6)
function y = fun(Ms)
P1 = 101.325;    % Assuming P1 = 1atm = 101.325kPa at SSL;
Ra = 287;        % R for air is 287;
T1 = 300;        % Assuming the room temperature is 27C
C1 = sqrt(1.4*Ra*T1);
Rho1 = P1/(Ra*T1);
ushock = Ms*C1;
u1 = ushock;
%State 5 conditions
P5 = 500*101.325;% storage pressure.
Rh = 4124;
T5 = 300;
C5 = sqrt(1.4*Rh*T5);
Rho5 = P5/(Rh*T5);
k = P5/(Rho5^(1.4));
P2 = P1*(2.8*(Ms.^2-1)/2.4)+P1;
u2 = u1 - C1*2*(Ms.^2-1)./(Ms*2.4);
W3 = ushock - u2;
P3 = P2;
Rho3 = (P3/k).^(1/1.4);
C3 = sqrt(1.4*P3./Rho3);
W3_P = 2*(C5-C3)/0.4;

y = W3 - W3_P
end

Ms0 = 6.5405