Solving function in Matlab

3 次查看（过去 30 天）
Haocheng Du
Haocheng Du 2020-3-4
回答： David Goodmanson 2020-3-5
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 个）

David Goodmanson
David Goodmanson 2020-3-5
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

类别

在 Help Center 和 File Exchange 中查找有关 Fluid Dynamics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by