debutant matlab Résolu/Fermé

Question

Bonjour, 
je suis sous windows,
et j'essaye de crée un programme matlab pour tracer les isotherme d'une équation d'État

voici mon code

% isotherm(nomgaz) d'Andrews

nomgaz = 'ch4';
eos='PengRobinson'
[u,w,Omega,m,R,Tc,Pc]=proprigaz(nomgaz,'PengRobinson');
nomgaz
ac=0.15724*(R*Tc)^2/Pc
b = 0.07780*R*Tc/Pc
Pc
Tc

% Nombre d'isotherme nisot
nisot = 6;
nech=2000;

% Volume molaire V
pas= 15;
V=0.001:pas:pas*nech
P=zeros(1,nech);
TK=273.15;
dT=20;
T=Tc+dT:TK-dT*(nisot-1);
Tr=T./Tc;
RT=R.*T;

P=R * T /(V - b) - ac /(V*(V + b)+b*(V-b));

% plotting isotherms for T
figure(1)
h=plot(V,P);
set(h,'color',rand(1,3),'linewidth',1);
hold on
axis([0 1600 -40 60])
xlabel('Volume in cm3/mol')
ylabel('pressure in bar')
title('Isotherms for propane')


je n'obtiens rien alors que je devrais obtenir une graphe de ce type
https://upload.wikimedia.org/wikipedia/commons/3/3e/Real_Gas_Isotherms.svg

merci d'avance
Chevrotine22

JulienJust · Answer

Bonjour,

Quelle erreur s'affiche pour ce script du coup ?

Cdlt

chevrotine22 · Answer

il affiche:

Error using  / 
Matrix dimensions must agree.

Error in master (line 26)
P=R * T /(V - b) - ac /(V*(V + b)+b*(V-b));

JulienJust · Answer

Peut-être devrais-tu taper:

P=R .* T ./(V - b) - a.*c./(V.*(V + b)+b.*(V-b));

ça fonctionne ? Ou ça donne toujours la même erreur ?

chevrotine22 · Answer

non y a toujours la meme erreur

JulienJust · Answer

Bon apparemment c'est clairement un problème de dimension donc ma première solution est futile, désolé.

Je te conseille en fait de voir si tout tes vecteurs sont de la même taille, c'est-à-dire si :

T et V sont de même taille ?


P.S: ac, b et R sont bien constants ?

chevrotine22 · Answer

alors je met ci joint la fonction d'ou vient mes variables

ac, b et R sont constant
T et une temperature
V un volume
je ne vois pas de quel taille tu parle
(je debute)

function [u,w,Omega,m,R,Tc,Pc]=proprigaz(nomgaz,eos) 
% Propriétées de corps purs
%  P en MPa
%  T en K
R = 83.14413;    %  J/(mol K)  expression en kg litre^2 /s^3 /mol K

% choix du gaz
if(strcmp(nomgaz,'ar'))
% argon a=1.355;   b=0.0320;
Teb=87.30;
Pc=150.87;
Tc=48.98;
elseif(strcmp(nomgaz,'ch4'))
% méthane a=2.303;   b=0.0431;
Teb=111.67;
Pc=45.99;
Tc=190.56;
Omega=0.011;
elseif(strcmp(nomgaz,'c2h6'))
% éthane a=2.303;   b=0.0431;
Tc = 305.32;
Pc = 48.72;
Omega= 0.099;
elseif(strcmp(nomgaz,'c3h8'))
% propane a=2.303;   b=0.0431;
Teb=231.02;
Tc = 369.83;
Pc = 42.48;
Omega= 0.152;
elseif(strcmp(nomgaz,'c4h10'))
% butane a=;   b=;
Tc = 425.12;
Pc = 37.96;
Omega= 0.200;
elseif(strcmp(nomgaz,'h2'))
% a=0.2452;   b=0.0265;
Teb=20.28;
Tc=32.97;
Pc=12.93;
elseif(strcmp(nomgaz,'h2o'))
% a=5.537;   b=0.0305;
Teb=373.2;
Tc=647.14;
Pc=220.6;
elseif(strcmp(nomgaz,'o2'))
% a=1.382;   b=0.0319;
Teb=90.20;
Pc=154.59;
Tc=50.43;
elseif(strcmp(nomgaz,'co'))
% a=1.472;    b=0.0395;
Teb=81.7
Tc=132.91;
Pc=34.99;
elseif(strcmp(nomgaz,'co2'))
% a=3.658;    b=0.0429;
Teb=194.6
Tc=304.13;
Pc=73.75;
end

% choix de l'eos
if(strcmp(eos,'VanderWals'))
    u=0; w=0;
elseif(strcmp(eos,'RedllichKwong'))
    u=1; w=0;
elseif(strcmp(eos,'PengRobinson'))
    u=2; w=-1;
    if(Omega<0.49)
    m = 0.37464 + Omega*(1.54226 - 0.26992*Omega)
    else
    m = 0.379642+ Omega*(1.48503 - Omega*(0.164423+0.0166666*Omega));
    end
elseif(strcmp(eos,'Harmens'))
    u=3; w=-2;
end

end

JulienJust · Answer

Il n'était pas nécessaire de mettre tout ça...

Tu as écris dans le sujet que :

% Volume molaire Vpas= 15;V=0.001:pas:pas*nechP=zeros(1,nech);TK=273.15;dT=20;T=Tc+dT:TK-dT*(nisot-1);Tr=T./Tc;RT=R.*T; 

Donc V et T sont des vecteurs. Je demande simplement de quelle longueur ils sont. Donc il faut taper:

length(T);  % Pour la taille de T (= nombre de coordonnées dans le vecteur T)
length(V);  % Pour la taille de V (= nombre de coordonnées dans le vecteur V)

Obtiens-tu des nombres différents ? 


P.S. : Si tu ne sais pas ce que tu as écrit, qui l'a écrit alors...?

chevrotine22 · Answer

il y a effectivement une difference

Length(T)=0
Length(V)=2000

Je n'ai pas ecris tout ce code, c'est mon "prof" qui nous l'a donné et a demandé a ce qu'on le finissent, j'ai juste ajouté la fonction finale P et le code pour les graphes

JulienJust · Answer

Okay dans ce cas le prof a du oublié quelque chose car T est nulle :-/ voilà où est le problème. 

Re-teste ton programme une fois qu'il aura arrangé ça :)

chevrotine22 · Answer

bonjour j'ai modifié mon code mais rien y fait. 
maintenant j'ai ce message:

??? Undefined function or method 'lenght' for input arguments of type 'double'.

Error in ==> master at 26
lenght(i);
 
voici le code
% isotherm(nomgaz) d'Andrewsnomgaz = 'ch4';eos='PengRobinson'[u,w,Omega,m,R,Tc,Pc]=proprigaz(nomgaz,'PengRobinson');nomgazac=0.15724*(R*Tc)^2/Pcb = 0.07780*R*Tc/PcPcTc% Nombre d'isotherme nisotnisot = 6;nech=2000;% Volume molaire Vpas= 15;V=0.001:pas:pas*nechP=zeros(1,nech);TK=273.15;dT=20;i=40:10:90;T(i)=273.15+iTr=T(i)./Tc;RT=R.*T(i);lenght(i);P=R*T(i)./(V - b) - ac./(V.*(V + b)+b*(V-b));figure(1)h=plot(v,P);set(h,'color',rand(1,3),'linewidth',1);hold onaxis([0 1600 -40 60])xlabel('Volume in cm3/mol')ylabel('pressure in bar')title('Isotherms for propane')

chevrotine22 · Answer

y a pas besoin que V et i est la meme taille, dans un exemple que j'ai 
length(v)=2000
et length(T)=1

comme ca c'est bon mais la courbe a pas la bonne forme

% isotherm(nomgaz) d'Andrews

nomgaz = 'ch4';
eos='PengRobinson'
[u,w,Omega,m,R,Tc,Pc]=proprigaz(nomgaz,'PengRobinson');
nomgaz
ac=0.15724*(R*Tc)^2/Pc
b = 0.07780*R*Tc/Pc
Pc
Tc

% Nombre d'isotherme nisot
nisot = 6;
nech=2000;

for i=40:10:90

v=0.001:pas:pas*nech;

T(i)=273.15+i;

Tre = T(i)/Tc;

a = 0.45724*(R*Tc)^2/Pc*(1 + m*(1 - sqrt(Tre)))^2;

P=R*T(i)./(v - b) - a./(v.*(v + b)+b*(v-b));

figure(1)
h=plot(v,P);
set(h,'color',rand(1,3),'linewidth',1);
hold on
axis([0 1600 -40 60])
xlabel('Volume en cm3/mol')
ylabel('pression en bar')
title('Isothermes pour propane')
end

JulienJust · Answer

Okay je comprends mieux maintenant. Fais un test avec cette modification et dis moi ce que ça donne s'il te plait:

nisot = 6;nech = 2000;v = 0.001:pas:pas*nech;n = 40:10:90;T = zeros(1,length(n));figure(1);for ii = 1:length(T);  T(ii) = 273.15 + ii;  Tre = T(ii)/Tc;  a = 0.45724*(R*Tc)^2/Pc*(1 + m*(1 - sqrt(Tre)))^2;  P = R*T(ii)./(v - b) - a./(v.*(v + b)+b*(v-b));  h = plot(v,P);  set(h,'color',rand(1,3),'linewidth',2); grid on; hold on;  %axis([0 1600 -40 60]);  xlabel('Volume en cm3/mol');  ylabel('pression en bar');  title('Isothermes pour propane');endhold off

chevrotine22 · Answer

ca fonctionne mais les courbes n'on pas la bonne forme:

j'obtiens de simple parabole alors que je dois avoir ca :
https://upload.wikimedia.org/wikipedia/commons/3/3e/Real_Gas_Isotherms.svg

chevrotine22 · Answer

c'est bon ca marche,
maintenant je dois tracer une droite parallèle a l'axe X, elle va etre coupé en 3 points, ce qui fait 2 aires.
Je dois calculer la surface de ces 2 aires pour différentes valeurs de P (axe Y) jusqu'à que les 2 aires soit égaux.
La formule c'est:
int(v1,v2) Pdv =S=[R*T*ln(V-b)-(a:(2*b*sqrt(2)))*ln((V+b*(1-sqrt(2)))/(V+b*(1+sqrt(2))))]

je sais pas si on peux faire des boucles while avec matlab mais si oui je ferais un truc du genre

while S>0 or S<0

P=0:0,0001:-40 %c'est un exemple

% pour chaque valeur de P, V1 = min de v et V2 = max de v
%calculer l'intégrale entre les points V1 et V2

Debutant matlab

14 réponses

Discussions similaires