Problem with inf value

19 views (last 30 days)
Plamen Bonev on 20 May 2021
Commented: Plamen Bonev on 20 May 2021
I am having inf values, while I am supposed to be having 1.0396e+01 value... (which is apparently not at all an infinite value)...
I have 4 loops in my code : and I am trying to get the ratios condGaz(i,j,k) / condGaz(i,j,g), where k=1:3 ; g=1:3.
Here is a part of my code :
----------- Before this, I have 3 loops i (for Temperature), j (for Pressure), k=1:3 and this loop is my final loop ------------
for g=1:3
if k~=g
Phigk(i,j,k,g) = condGaz(i,j,k) / condGaz (i,j,g);
else
Phigk(i,j,k,g) =0 ;
end
end
end
end
end
Here are the condGaz(i,j,k) values.
And here are the values that I obtain for my ratio Phigk (i,j,k,g) = condGaz(i,j,k) / condGaz (i,j,g);
I get the impression that for g=1, I get all values right...
As soon as it gets to g=2 I get 1 value right and the other is INF
And for g=3, I get all the values INF....
I am not supposed to be having INF values because I manually did the ratio of for example Phigk(1,1,2,3) = condGaz(1,1,2) / condGaz(1,1,3) = 1.23e+00. (which is INF value if you see the results with the loops for Phigk(:,:,2,3)...
I tried format shortE and stuff... And it did not work... So I am confused about the sense of this INF value.... It is supposed to give INF value when the result is too big... But here we are talking values of 1.23e+00.... Can you please help me?
Plamen Bonev on 20 May 2021
Here is the full version of my code :
format shortE
%Constante des gaz parfaits en J/mol.K
Rg=8.314;
%---------------------- Températures en K -------------------------
%T=[200 : 100 : 1000];
%T=[1000 : 100 : 2000];
T=[600 : 200 : 2000];
% -------------------- Pressions en Pa -----------------------------
%------------ P = [0.01 Mpa ; 0.1 Mpa] ---------------------
% P=[0.1e5:0.1e5:1e5]
% ----------- P = [0.1 MPa (1 bar) ; 10 MPa (10 bar)] ---------
% P=[1e5:2e5:100e5];
% ----------- P = [40MPa, 45MPa] -----------
P=[400e5:50e5:450e5];
% ----------- P = [0.1MPa, 1MPa] -----------
% P=[1e5:1e5:10e5];
% P=[0.1e5 : 1e5 : 200e5];
% ---------- Echelle log pour la Pression (de 1e4 Pa (0.01 MPa) à 1e7 Pa (10MPa) avec 1000 points entre les
% deux) ---------
% P=logspace(4, 7, 1000)
%He indice 1
%Kr indice 2
%Xe indice 3
% ------------------- Nombre de gaz considéré -----------------------
NbGaz=[1, 2, 3];
%--------------Paramètres critiques pour chaque gaz (pour le modèle avancé) ----------
%Masses molaires en kg/mol
M=[4.003e-3, 83.798e-3, 131.293e-3];
%Pressions critiques Pc en Pa
Pc=[0.2275e6, 5.51e6, 5.84e6];
%Températures critiques Tc en K
Tc=[5.2, 209.4, 289.7];
%Densités critiques rhoc en kg/m^3
rhoc=[69.64, 908.4, 1110];
%Coefficients directeurs pour viscosité dynamique
Anu=[3.063e-7, 6.963e-7, 7.568e-7];
%Coefficients de température pour viscosité dynamique en K
Tnu=[-21.33, 71.07, 112.31];
%Exposants pour viscosité dynamique
n=[0.724, 0.667, 0.655];
%--------------- Constantes a et b pour chaque gas (pour le modèle
%simplifié) ---------------
%He indice 1
%Kr indice 2
%Xe indice 3
aURGAP=[176.32e-5, 9.8e-5, 4.6e-5];
bURGAP=[0.77227, 0.8007, 0.8425];
aROCHE=[284e-5, 8.11e-5, 5.14e-5];
bROCHE=[0.7, 0.83, 0.83];
aBISON=[263.9e-5, 8.247e-5, 4.351e-5];
bBISON=[0.7085, 0.8363, 0.8616];
%------------ Initialisation des vecteurs/matrices -----------------
visc0=zeros(length(T),1);
cond0=zeros(length(T),1);
Vetoile=zeros(length(NbGaz),1);
condPseudoCr=zeros(length(NbGaz),1);
rho=zeros(length(T), length(P), length(NbGaz));
B2=zeros(length(T), length(NbGaz));
B3=zeros(length(T), length(NbGaz));
theta=zeros(length(T), length(NbGaz));
coef1=zeros(length(T), length(P));
coef2=zeros(length(T), length(P), length(NbGaz));
coef3=zeros(length(T), length(P), length(NbGaz));
racines=cell(length(T), length(P), length(NbGaz));
vraiRacines=cell(length(T), length(P), length(NbGaz));
B=cell(length(T), length(P), length(NbGaz));
imagi=cell(length(T), length(P), length(NbGaz));
rhor=zeros(length(T), length(P), length(NbGaz));
psi=zeros(length(T), length(P), length(NbGaz));
condGaz=zeros(length(T), length(P), length(NbGaz));
condGazURGAP=zeros(length(T), length(NbGaz));
condGazROCHE=zeros(length(T), length(NbGaz));
condGazBISON=zeros(length(T), length(NbGaz));
rapportURGAP=zeros(length(T), length(P), length(NbGaz));
pourcentageAugm=zeros(length(T), length(P), length(NbGaz));
diffURGAP=zeros(length(T), length(P), length(NbGaz));
diffBISON=zeros(length(T), length(P), length(NbGaz));
diffROCHE=zeros(length(T), length(P), length(NbGaz));
Psigk=zeros(length(T), length(P), length(NbGaz), length(NbGaz));
Phigk=zeros(length(T), length(P), length(NbGaz), length(NbGaz));
for k=1:numel(NbGaz)
%Volumes molaires caractéristiques
Vetoile(k)=Rg*Tc(k)/Pc(k);
condPseudoCr(k)=(0.201e-4*((Tc(k))^(0.277))*((M(k))^(-0.465)))/((0.291*Vetoile(k))^(0.415));
for i=1:numel(T)
%----------Conductivité de chaque gaz (modèle simplifié)---------
%Données a et b URGAP
condGazURGAP(i,k)=aURGAP(k)*(T(i))^(bURGAP(k));
%Données a e b ROCHE
condGazROCHE(i,k)=aROCHE(k)*(T(i))^(bROCHE(k));
%Données a et b BISON
condGazBISON(i,k)=aBISON(k)*(T(i))^(bBISON(k));
%----------Modèle avancé---------
theta(i,k)=(T(i))/(Tc(k));
visc0(i,k)=Anu(k)*(T(i)-Tnu(k))^(n(k));
cond0(i,k)=visc0(i,k)*((15*Rg)/(4*M(k)));
if k==2 | k==3
B2(i,k)=(-102.6+(102.732-0.001*theta(i,k)-0.44*((theta(i,k))^(-1.22)))*tanh(4.5*sqrt(theta(i,k))))*Vetoile(k);
B3(i,k)=(0.0757+(-0.0862+(-3.6e-5)*theta(i,k)+0.0237*((theta(i,k))^(-0.059)))*tanh(0.84*theta(i,k)))*((Vetoile(k))^2);
else
B2(i,k)=(8.4-0.0018*T(i)+(115/(sqrt(T(i))))-835/(T(i)))*(10^(-6));
B3(i,k)=0;
end
for j=1:numel(P)
%coef1, 2 et 3 corréspondent aux coefficients devant rho
coef1(i,j,k)=(Rg*T(i))/(P(j));
coef2(i,j,k)=(B2(i,k))*((Rg*T(i))/(P(j)));
coef3(i,j,k)=(B3(i,k))*((Rg*T(i))/(P(j)));
% On cherche les racines de l'équation cubique d'état
racines{i,j,k}=((roots([coef3(i,j,k), coef2(i,j,k), coef1(i,j,k), -1])));
% On récupère uniquement les racines réelles
vraiRacines{i,j,k}=racines{i,j,k}(imag(racines{i,j,k})==0);
% La densité corréspond à la valeur maximale de chaque racine trouvée
% pour chaque jeu de (P, T, gaz)
rho(i,j,k)=(max(vraiRacines{i,j,k}))*M(k);
% Densité réduite
rhor(i,j,k)=(rho(i,j,k))/rhoc(k);
% Psi
% psi(i,j,k)=0.645*rhor(i,j,k)+0.331*((rhor(i,j,k))^2)+0.0368*((rhor(i,j,k))^3)-0.0128*((rhor(i,j,k))^4);
psi(i,j,k)=0.645+0.331*((rhor(i,j,k))^1)+0.0368*((rhor(i,j,k))^2)-0.0128*((rhor(i,j,k))^3);
% % Calcul final de la conductivité d'un gaz
condGaz(i,j,k)=cond0(i,k)+(condPseudoCr(k)*psi(i,j,k)*(rhor(i,j,k)));
%RapportsURGAP
rapportURGAP(i,j,k)=(condGaz(i,j,k))/(condGazURGAP(i,k));
%Calcul d'erreur entre le modèle avancé/modèle simple
% Erreur URGAP / avancé (en %)
diffURGAP(i,j,k)=(condGazURGAP(i,k)-condGaz(i,j,k))*10^2;
% Erreur BISON/ avancé (en %)
diffBISON(i,j,k)=(condGazBISON(i,k)-condGaz(i,j,k))*10^2;
% Erreur ROCHE/ avancé (en %)
diffROCHE(i,j,k)=(condGazROCHE(i,k)-condGaz(i,j,k))*10^2;
% POUR GUILHERME : LE PROBLEME COMMENCE ICI.
for g=1:numel(NbGaz)
disp(['i,j,k,g = ' num2str([i j k g]) ' ; condGaz (i,j,g) = ' num2str(condGaz(i,j,g))])
Phigk(i,j,k,g)=(condGaz(i,j,k))/(condGaz(i,j,g));
end
end
end
end

Walter Roberson on 20 May 2021
condGaz(i,j,k)=cond0(i,k)+(condPseudoCr(k)*psi(i,j,k)*(rhor(i,j,k)));
That sets condGaz(I,j,:) according to the largest k value that has been seen so far. For example at k = 1, condGaz(i,j,1) will be set, and condGaz(i,j,2) will not have been set yet.
for g=1:numel(NbGaz)
disp(['i,j,k,g = ' num2str([i j k g]) ' ; condGaz (i,j,g) = ' num2str(condGaz(i,j,g))])
but numel(NbGaz) is greater than k = 1, so in this inner loop within for k, as soon as g is greater than the current k, you index condGaz(i,j,:) elements that are still 0, and that gets you a division by 0.
Perhaps that loop needs to be postponed until after the other arrays have been fully built.
Plamen Bonev on 20 May 2021
Thank you, Walter, you solved it all for me.
Cheers.