Solving function in Matlab

7 Ansichten (letzte 30 Tage)
Haocheng Du
Haocheng Du am 4 Mär. 2020
Beantwortet: David Goodmanson am 5 Mär. 2020
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
  1 Kommentar
darova
darova am 4 Mär. 2020
Please attach original equations

Melden Sie sich an, um zu kommentieren.

Antworten (1)

David Goodmanson
David Goodmanson am 5 Mär. 2020
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

Kategorien

Mehr zu Programming finden Sie in Help Center und File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by