2014-03-05 21 views
2

我正在尝试使用Matlab对气体进行一些模拟,并且希望绘制从Matrix M中获取的预测数据。此数据应该显示为单行,但看起来它正在被绘制并填充在两个值之间。我哪里错了? 预先感谢为什么我的matlab图填写?

%Basic One Dimensional Simulation 

%Set Parameters 
Y = 1.32; %spec heat ratio 
po = 103200; %pressure at outlet 
pi = 300000; %pressure at inlet 
dx = 0.3; %length of control volume 
cbar = 340.29; %average speed of sound, m/s 
Fr = 4; %Frequency of wave 
Cd = 0.7; %Discharge coefficient 
Tu = 1800; %Rough temperature of LNG combustion 
Am = 1; %amplitude of wave 
A = 1; %Cross sectional area of exhaust chamber 
Ao = 0.4; %Cross sectional area of exhaust exit 
Ai = 0.3; %Cross sectional area of exhaust inlet 

%Define simulation setting time constants 
Ts = 0.0001; %Sampling interval, sec 
T = 1;  %simulation time, sec 
mt = T/Ts; %Number of iterations 


%Initial values 
%Mass flow rates 
Mg = 0.0110352; %Mass flow rate of gas 
Ma = 0.189189; %Mass flow rate of air 
Mi = Mg + Ma; %Total Mass flow rate at inlet 
M4 = Mi; 
M3 = M4; 
M2 = M3; 
M1 = M2; 

p1 = po; %Pressure, Pa 
p2 = p1; 
p3 = p2; 
p4 = p3; 

%Simulation  
M = zeros(mt); 
N = zeros(mt); 
for j = 0:1:mt 
    t=j*Ts; 

Mi = Mg + Ma; %Mass flow rate based on input flow rates 
a = (Ma/0.172); 
Rexhaust = (84373 + 70841*a)/(290.67 + 210.67*a); 

%Predicting Mass Flow Rate at exit 
Cm1 = (((2*Y)/(Rexhaust*(Y-1)))*(((po/p1)^(2/Y))*((po/p1)^((Y+1)/Y))))^(1/2); 
Cm2 = ((2/(Y+1))^(1/(Y-1)))*(2*Y/Rexhaust*(Y+1))^(1/2); 

Pratio = (2/(Y+1))^(Y/(Y-1)); 

if po/p1 > Pratio; 
Mo = Cd*Cm1*Ao*p1/(Tu)^(1/2); 
else 
Mo = Cd*Cm2*Ao*p1/(Tu)^(1/2); 
end 

%Position 1 
dpdt1 = ((Rexhaust*Tu)/(dx*A))*(Mo-M1); 
drhodt1 = dpdt1/(Rexhaust*Tu); 
p1 = dpdt1*Ts; 
M1 = (dpdt1*A*dx)/(Rexhaust*Tu); 


%Position 2 
dpdt2 = ((Rexhaust*Tu)/(dx*A))*(M1-M2); 
drhodt2 = dpdt2/(Rexhaust*Tu); 
p2 = dpdt2*Ts; 
M2 = (dpdt2*A*dx)/(Rexhaust*Tu); 


%Position 3 
dpdt3 = ((Rexhaust*Tu)/(dx*A))*(M2-M3); 
drhodt3 = dpdt3/(Rexhaust*Tu); 
p3 = dpdt3*Ts; 
M3 = (dpdt3*A*dx)/(Rexhaust*Tu); 


%Position 4 
%Predicting Mass Flow Rate at inlet 
Cm1 = (((2*Y)/(Rexhaust*(Y-1)))*(((p4/pi)^(2/Y))*((p4/pi)^((Y+1)/Y))))^(1/2); 
Cm2 = ((2/(Y+1))^(1/(Y-1)))*(2*Y/Rexhaust*(Y+1))^(1/2); 

Pratio = (2/(Y+1))^(Y/(Y-1)); 

if p4/pi > Pratio; 
M4 = Cd*Cm1*Ai*pi/(Tu)^(1/2); 
else 
M4 = Cd*Cm2*Ai*pi/(Tu)^(1/2); 
end 
dpdt4 = ((Rexhaust*Tu)/(dx*A))*(M3-M4); 
drhodt4 = dpdt4/(Rexhaust*Tu); 
p4 = dpdt4*Ts; 
M4 = (dpdt4*A*dx)/(Rexhaust*Tu);  


%Save current values to Matrix M 
d = [t Mo M1 M2 M3 M4 Mi dpdt1 dpdt2 dpdt3 dpdt4 drhodt1 drhodt2 drhodt3 drhodt4]; 
for i = 1:1:15 
    M(j+1,i) = d(i); 
end 

%Save current p to Matrix N 
e = [t po p1 p2 p3 p4];    
for k=1:1:6 
    N(j+1,k) = e(k); 
end 


end 


%FIGURES 
F1 = figure(1); 
subplot(3,1,1); 
plot(M(:, 1), M(:, 2), 'y'); %outlet 
hold on 
plot(M(:, 1), M(:, 3), 'g'); %position 1 
plot(M(:, 1), M(:, 4), 'r'); %position 2 
plot(M(:, 1), M(:, 5), 'b'); %position 3 
plot(M(:, 1), M(:, 6), 'm'); %position 4 
plot(M(:, 1), M(:, 7), 'c'); %inlet 
grid on 
hold off 
legend('outlet','M1','M2','M3','M4','Mi'); 
title('MASS FLOW RATES'); 

subplot(3,1,2); 
plot(M(:, 1), M(:, 8), 'g'); %position 1 
hold on 
plot(M(:, 1), M(:, 9), 'r'); %position 2 
plot(M(:, 1), M(:, 10), 'b'); %position 3 
plot(M(:, 1), M(:, 11), 'm'); %position 4 
grid on 
hold off 
legend('dpdt1','dpdt2','dpdt3','dpdt4'); 
title('RATE OF CHANGE OF PRESSURES'); 

subplot(3,1,3); 
plot(M(:, 1), M(:, 12), 'g'); %position 1 
hold on 
plot(M(:, 1), M(:, 13), 'r'); %position 2 
plot(M(:, 1), M(:, 14), 'b'); %position 3 
plot(M(:, 1), M(:, 15), 'm'); %position 4 
grid on 
hold off 
legend('drhodt1','drhodt2','drhodt3','drhodt4'); 
title('RATE OF CHANGE OF DENSITIES'); 


F2 = figure(2); 
plot(N(:, 1), N(:, 2), 'g'); 
hold on 
plot(N(:, 1), N(:, 2), 'r'); 
plot(N(:, 1), N(:, 2), 'b'); 
plot(N(:, 1), N(:, 2), 'm'); 
hold off`enter code here` 
legend('p1','p2','p3','p4'); 
title('PRESSURES'); 

回答

0

填充区域是由于数据的性质(即高频分量)。下图显示了相同的数据,plot(M(:, 1), M(:, 6), 'm');,但以不同的比例(右侧放大)。

enter image description here

如果你把M(1:2:end, 1)M(2:2:end, 1),你将有当地的极值。

+0

谢谢!更好地降低频率! – PogFrog

相关问题