Wie NaN Problem finden?

WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Partner:

Forum

Option

[Erweitert]

• Diese Seite per Mail weiterempfehlen

Gehe zu:

eh27

Forum-Anfänger


	Beiträge: 12

	Anmeldedatum: 31.03.16

	Wohnort: ---

	Version: ---

Verfasst am: 26.04.2016, 11:13 Titel: Wie NaN Problem finden?

Hallo zusammen.

Ich habe folgendes Problem. Mein Programm sagt mir bei dem Wert pi_druck (= p_Z2 / p_Z1 ) an der Stelle k=68001 NaN obwohl p_Z1(68001) einen normalen Zahlenwert hat und p_Z2(68001) auch einen normalen Zahlenwert hat. Wo entsteht das NaN?

Danke für eure Hilfe.

Code:

ainput = 3;
kolbenwegmax=2*pi*0.057;
b_Z1=0.03;
h_Z1=0.03;
A_kolben=b_Z1*h_Z1;
V_max=kolbenwegmax*A_kolben;
p_Umg=101325; %Umgebungsdruck in pascal
T_Umg=293.75; %Umgebungstemperatur in Kelvin
kappa=1.4;
c_p=1005; %
c_v=718; %
R=287; %

T_w=50; %50 Kelvin ist T_Wand erhöht
VE_Z1_offen=320; %Ventil Zylinder 1->2 auf
VE_Z1_zu=350; %Ventil Zylinder 1 zu

ZS_Z1_offen=350; %Zylindersperre Z1 auf
ZS_Z1_zu=356; %Zylindersperre Z1 zu

Versatz = 330 - VE_Z1_offen; % negativer Versatz von Ventil Z1 offen
Versatz_ges=100*(330 - 10 - Versatz)+1; % Versatz mit 1/100 Genauigkeit

A_VE1=0.00005; % Ventilquerschnitt in m^2
dt=1/(6*2400*100); % dt=1/6n * 100 für 1/100°KW

alpha_BA = 30 + Versatz ; % Brennanfang
alpha_BD = 25; % Brenndauer
m = 2 ; % Formfaktor
m_Luft = 1.33779E-05; % Masse Luft
Q_B_ges = m_Luft * 40.98E06; % Q_max

if isempty(ainput)
ainput=10;
end
alphai=0:0.01:(ainput*360);

%Präallozieren

kw = zeros(size(alphai));
alpha_w_Z1= zeros(size(alphai));
alpha_w_Z2= zeros(size(alphai));
dQw_Z1= zeros(size(alphai)); %Wandwärmeverluste
dQw_Z2= zeros(size(alphai)); %Wandwärmeverluste
dT_Z1= zeros(size(alphai)); %Temperaturänderung
dp_Z1= zeros(size(alphai)); %Druckänderung
dT_Z2= zeros(size(alphai));
dp_Z2= zeros(size(alphai)); %Druckänderung
c_s= zeros(size(alphai)); %Schallgeschwindigkeit
rho_1= zeros(size(alphai)); % Dichte 1
rho_2= zeros(size(alphai)); % Dichte 2
rho= zeros(size(alphai)); %Dichte Strömung 1->2
m1= zeros(size(alphai)); %Masse im Zylinder 1
m2= zeros(size(alphai)); %Masse im Zylinder 2
m_punkt= zeros(size(alphai));
dW_Z1= zeros(size(alphai)); %Volumenänderungsarbeit Z1
dW_Z2= zeros(size(alphai)); %Volumenänderungsarbeit Z2
dm= zeros(size(alphai)); %Masseveränderung 1->2
dH= zeros(size(alphai)); %Enthalpieänderung
VE_Z1_o= zeros(size(alphai));
VE_Z1_z= zeros(size(alphai));
ZS_Z1_o= zeros(size(alphai));
ZS_Z1_z= zeros(size(alphai));
s_k= zeros(size(alphai));
s_k2= zeros(size(alphai));
pi_druck= zeros(size(alphai));
p_Z1= zeros(size(alphai));
T_Z1= zeros(size(alphai));
p_Z2= zeros(size(alphai));
T_Z2= zeros(size(alphai));
V_Z1= zeros(size(alphai));
V_Z2= zeros(size(alphai));
p_Z1_hK= zeros(size(alphai));
T_Z1_hK= zeros(size(alphai));

pi_druck_Z1=zeros(size(alphai));

dT_Z1_vK = zeros(size(alphai));
dp_Z1_vK = zeros(size(alphai));
dT_Z1_hK = zeros(size(alphai));
dp_Z1_hK = zeros(size(alphai));
p_0 = zeros(size(alphai));
w = zeros(size(alphai));
qexp = zeros(size(alphai));
q = zeros(size(alphai));
dq = zeros(size(alphai));
dQ_B = zeros(size(alphai));

for k = 1:numel(alphai)
if alphai(k)>360
kw_ungerundet=alphai(k)/360;
kw(k)=alphai(k)-(floor(kw_ungerundet)*360);
VE_Z1_o(k) = 360+VE_Z1_offen*(floor(kw_ungerundet));
VE_Z1_z(k) = 360+VE_Z1_zu*(floor(kw_ungerundet));
ZS_Z1_o(k) = 360+ZS_Z1_offen*(floor(kw_ungerundet));
ZS_Z1_z(k) = 360+ZS_Z1_zu*(floor(kw_ungerundet));
test(k) = 36001+36001*(floor(kw_ungerundet));
else
kw(k)=alphai(k);
VE_Z1_o(k)=VE_Z1_offen;
VE_Z1_z(k)=VE_Z1_zu;
ZS_Z1_o(k)=ZS_Z1_offen;
ZS_Z1_z(k)=ZS_Z1_zu;
test(k)=36001;
end
s_k(k)=kolbenwegmax*(kw(k)/360); %Kolbenweg

if k > Versatz_ges
s_k2(k)=(kolbenwegmax*((kw(k-Versatz_ges))/360)); %Kolbenweg 2
kw2(k)=(kw(k-Versatz_ges));
else
s_k2(k)=0.5;
kw2(k)=0;
end

V_Z1(k)=A_kolben*(s_k(k)-kolbenwegmax)*(-1);
V_Z2(k)=A_kolben*(s_k2(k));

dV_Z1=-8.953539062722324e-09; %diff(V_Z1);
dV_Z2=8.953539062722324e-09;
if alphai(k)< 1
p_Z1(k)=p_Umg/((V_Z1(k)/V_max)^kappa);
T_Z1(k)=T_Umg/((p_Umg/p_Z1(k))^((kappa-1)/kappa));
p_Z1_hK(k)= 101325; %p_Z1 hinter Kolben
T_Z1_hK(k)= 293.15; %T_Z1 hinter Kolben
p_Z2_hK(k)= 101325; %p_Z2 hinter Kolben
T_Z2_hK(k)= 293.15; %T_Z2 hinter Kolben

elseif alphai(k)>=1 && alphai(k)<VE_Z1_offen
p_Z1(k)=p_Z1(k-1)/((V_Z1(k)/V_Z1(k-1))^kappa);
T_Z1(k)=T_Z1(k-1)/((p_Z1(k-1)/p_Z1(k))^((kappa-1)/kappa));
dH(k)=0;
dQw_Z1(k)=0;
p_Z2(k)=101325;
T_Z2(k)=293.15;
p_Z1_hK(k)= 101325; %p_Z1 hinter Kolben
T_Z1_hK(k)= 293.15; %T_Z1 hinter Kolben

elseif alphai(k)>=VE_Z1_o(k) && alphai(k)<=VE_Z1_z(k)
p_Z1(k)=p_Z1(k-1)+dp_Z1(k-1);
T_Z1(k)=T_Z1(k-1)+dT_Z1(k-1);
p_Z2(k)=p_Z2(k-1)+dp_Z2(k-1);
T_Z2(k)=T_Z2(k-1)+dT_Z2(k-1);
p_Z1_hK(k)= 101325; %p_Z1 hinter Kolben
T_Z1_hK(k)= 293.15; %T_Z1 hinter Kolben

elseif alphai(k)>=ZS_Z1_o(k) && alphai(k)<=ZS_Z1_z(k)
p_Z1(k)=p_Z1(k-1)+dp_Z1_vK(k-1);
T_Z1(k)=T_Z1(k-1)+dT_Z1_vK(k-1);
p_Z1_hK(k)= p_Z1_hK(k-1)+dp_Z1_hK(k-1); %p_Z1 hinter Kolben
T_Z1_hK(k)= T_Z1_hK(k-1)+dT_Z1_hK(k-1); %T_Z1 hinter Kolben

elseif alphai(k)>ZS_Z1_z(k) && alphai(k)<=test(k)
p_Z1(k)=p_Z1(k-1);
T_Z1(k)=T_Z1(k-1);
p_Z1_hK(k)= p_Z1_hK(k-1)+dp_Z1_hK(k-1); %p_Z1 hinter Kolben
T_Z1_hK(k)= T_Z1_hK(k-1)+dT_Z1_hK(k-1); %T_Z1 hinter Kolben

else
p_Z1(k)=p_Z1(k-1)/((V_Z1(k)/V_Z1(k-1))^kappa);
T_Z1(k)=T_Z1(k-1)/((p_Z1(k-1)/p_Z1(k))^((kappa-1)/kappa));
dH(k)=1;
dQw_Z1(k)=0;
p_Z1_hK(k)= 101325; %p_Z1 hinter Kolben
T_Z1_hK(k)= 293.15; %T_Z1 hinter Kolben
end
p_Z1(36002)= p_Z1(35990);
V_Z1(36001)=V_Z1(36000);
pi_druck(k)=p_Z2(k)/p_Z1(k);
if pi_druck(k)<0.52
pi_druck(k)=0.52;
end
pi_druck_Z1(k)=p_Z1_hK(k)/p_Z1(k);
if pi_druck_Z1(k)<0.52
pi_druck_Z1(k)=0.52;
end
if pi_druck_Z1 >= 1

end
c_s(k) = sqrt(2 * c_p * T_Z1(k) * (1 - (pi_druck(k)^((kappa-1) / kappa))));
c_s_Z1(k) = sqrt(2 * c_p * T_Z1(k) * (1 - (pi_druck_Z1(k)^((kappa-1)/kappa))));%Schallgeschwindigkeit
rho_1(k) = p_Z1(k) / (T_Z1(k) * R); % Dichte 1
rho_1_hK(k) = p_Z1_hK(k) / (T_Z1_hK(k) * R); % Dichte 1 hinter Kolben
rho_2(k) = p_Z2(k) / (T_Z2(k) * R); % Dichte 2
rho(k)= rho_1(k) * (pi_druck(k)^(1/kappa)); %Dichte Strömung 1->2
rho_Z1(k) = rho_1(k) * (pi_druck_Z1(k)^(1/kappa)); %Dichte Strömung 1->1 hinter KOlben
m1(k) = rho_1(k) * V_Z1(k); %Masse im Zylinder 1
m1_hK(k) = rho_1_hK(k) * (V_max-V_Z1(k)); %Masse im Zylinder 1
m2(k)= rho_2(k) * V_Z2(k); %Masse im Zylinder 2
m_punkt(k) = c_s(k) * rho(k) * A_VE1;
m_punkt_Z1(k) = c_s_Z1(k) * rho_Z1(k) * A_kolben; % Strom Z1 - Z1 hK
dW_Z1(k) = dV_Z1 * p_Z1(k); %Volumenänderungsarbeit Z1
dW_Z1_hK(k) = -dV_Z1 * p_Z1_hK(k); %Volumenänderungsarbeit Z1 hK
dW_Z2(k) = dV_Z2 * p_Z2(k); %Volumenänderungsarbeit Z2
dm(k) = m_punkt(k) * dt; %Masseveränderung 1->2
dm_Z1(k) = m_punkt_Z1(k) * dt; %Masseveränderung 1->1 hK
dH(k) = dm(k) * T_Z1(k) * c_p; %Enthalpieänderung
dH_Z1(k) = dm_Z1(k) * T_Z1(k) * c_p; %Enthalpieänderung 1-> 1 hK

if kw2(k) >= alpha_BA && kw2(k) <= 55
qexp(k) = -6.908 * ((kw2(k)-alpha_BA)/alpha_BD)^(m+1); % Zeile im Exp
q(k) = 1-exp(qexp(k));
dq(k) = q(k) - q(k-1);
dQ_B(k) = dq(k) * Q_B_ges;
m2_Start = m2(35000);
else
dq(k) = 0;
dQ_B(k) = 0;
end
if kw2(k) >= alpha_BA && kw2(k) <= 350
T_Z2(k) = -(((p_Z2(k-1)*dV_Z2)-(m2_Start*c_p*T_Z2(k-1))-dQ_B(k-1)-dQw_Z2(k-1))/(m2_Start*c_p));
p_Z2(k) = (m2_Start * R * T_Z2(k-1))/V_Z2(k);
end
if k > 50000 && kw2(k) >0
T_Z2(k) = 330;
p_Z2(k) = 150000;
end
if kw2(k) > 350 && kw2(k) <= 360
T_Z2(k) = 330;
p_Z2(k) = 150000;
end
w_Z2(k) = (2.28+0.308)*16 + (6.22E-3 * (V_max*T_Umg)/(V_max*p_Umg)*(p_Z2(k)-p_0(k)));
if alphai(k) >= 1
p_0(k)=p_Z2(k)*(V_Z2(k-1)/V_Z2(k))^kappa;
else
p_0(k)=0;
end
alpha_w_Z1(k) = 127.93 * (0.03^(-0.2)) * ((p_Z1(k) * 1e-5)^0.8) * (T_Z1(k)^(-0.53)) * (41.408^.8);
alpha_w_Z1_hK(k) = 127.93 * (0.03^(-0.2)) * ((p_Z1_hK(k) * 1e-5)^0.8) * (T_Z1_hK(k)^(-0.53)) * (41.408^.8);
alpha_w_Z2(k) = 127.93 * ( 0.03^(-0.2)) * ((p_Z2(k) * 1e-5)^0.8) * (T_Z2(k)^(-0.53)) * (w_Z2(k)^.8);
dQw_Z1(k) = (2 * b_Z1 * h_Z1 + 4 * b_Z1 * ( kolbenwegmax - s_k(k))) * alpha_w_Z1(k) * ( 273.15 + T_w - T_Z1(k)) * dt; %Wandwärmeverluste
dQw_Z1_hK(k) = (2 * b_Z1 * h_Z1 + 4 * b_Z1 * (s_k(k))) * alpha_w_Z1_hK(k) * (273.15 + T_w - T_Z1_hK(k)) * dt; %Wandwärmeverluste
dQw_Z2(k) = (2 * b_Z1 * h_Z1 + 4 * b_Z1 * (s_k2(k))) * alpha_w_Z2(k) * ( 273.15 + T_w - T_Z2(k)) * dt; %Wandwärmeverluste
dT_Z1(k) = - (dH(k) + (c_v * T_Z1(k) * (-dm(k))) + dW_Z1(k) - dQw_Z1(k)) / (m1(k) * c_v); %Temperaturänderung
dp_Z1(k) = ((m1(k) * R * dT_Z1(k)) + ( R * T_Z1(k) * (-dm(k))) - dW_Z1(k)) / V_Z1(k); %Druckänderung
dT_Z2(k) = (dH(k) - (c_v * T_Z2(k) * dm(k)) - dW_Z2(k) + dQw_Z2(k)) / (m2(k)*c_v);
dp_Z2(k) = (( m2(k) * R * dT_Z2(k)) + ( R * T_Z2(k) * dm(k)) - dW_Z2(k)) / V_Z2(k); %Druckänderung
dT_Z1_vK(k) = -(dH_Z1(k)+(c_v*T_Z1(k)*(-dm_Z1(k)))+dW_Z1(k)-dQw_Z1(k))/(m1(k)*c_v); %Temperaturänderung 1->1hK
dp_Z1_vK(k) = ((m1(k) * R * dT_Z1_vK(k)) + (R * T_Z1(k) * (-dm_Z1(k))) - dW_Z1(k)) / V_Z1(k); %Druckänderung 1-1 hK
dT_Z1_hK(k) = (dH_Z1(k) - ( c_v * T_Z1_hK(k) * dm_Z1(k)) - dW_Z1_hK(k) + dQw_Z1_hK(k)) / ( m1_hK(k) * c_v);
dp_Z1_hK(k) = ((m1_hK(k) * R * dT_Z1_hK(k)) + ( R * T_Z1_hK(k) * dm_Z1(k)) - dW_Z1_hK(k)) / (V_max-V_Z1(k)); %Druckänderung

end

Funktion ohne Link?

Winkow

Moderator


	Beiträge: 3.842

	Anmeldedatum: 04.11.11

	Wohnort: Dresden

	Version: R2014a 2015a

Verfasst am: 26.04.2016, 11:26 Titel:

Zitat:

obwohl p_Z1(68001) einen normalen Zahlenwert hat und p_Z2(68001) auch einen normalen Zahlenwert hat.

das stimmt nicht. p_Z2 ist bei 68001 nan. jedenfalls dann wenn es zur berechneung von pi_druck kommt. du änderst p_Z2 ja danch nochmal. das hat aber keinen einfluss auf pi_druck. das lässt sich mit dem debugger recht schnell rausfinden.
grüße
_________________

richtig Fragen

eh27

Themenstarter

Forum-Anfänger


	Beiträge: 12

	Anmeldedatum: 31.03.16

	Wohnort: ---

	Version: ---

Verfasst am: 26.04.2016, 14:34 Titel:

Vielen Dank.

Kannst du mir kurz erklären wie man das mit dem debugger rausfinden kann?

Dankeschön!

Winkow

Moderator


	Beiträge: 3.842

	Anmeldedatum: 04.11.11

	Wohnort: Dresden

	Version: R2014a 2015a

Verfasst am: 26.04.2016, 14:42 Titel:

einfach set conditional breakpoint an die stelle setzen und die bedingung auf k==68001 setzen.
_________________

richtig Fragen


Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:	Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen

Impressum | Nutzungsbedingungen | Datenschutz | FAQ |

RSS

Hosted by:

Copyright © 2007 - 2025 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.