Debutant matlab
Résolu
chevrotine22
Messages postés
20
Date d'inscription
Statut
Membre
Dernière intervention
-
rodfaction Messages postés 1 Date d'inscription Statut Membre Dernière intervention -
rodfaction Messages postés 1 Date d'inscription Statut Membre Dernière intervention -
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
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
A voir également:
- Debutant matlab
- Logiciel de programmation pour débutant - Guide
- Logiciel montage vidéo débutant - Guide
- Apprendre le coran pour débutant (+ pdf) - Télécharger - Histoire & Religion
- Platine dj debutant - Forum Enregistrement / Traitement audio
- Comment utiliser un ordinateur pour un débutant - Astuces et Solutions
14 réponses
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));
Error using /
Matrix dimensions must agree.
Error in master (line 26)
P=R * T /(V - b) - ac /(V*(V + b)+b*(V-b));
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 ?
P=R .* T ./(V - b) - a.*c./(V.*(V + b)+b.*(V-b));
ça fonctionne ? Ou ça donne toujours la même erreur ?
Vous n’avez pas trouvé la réponse que vous recherchez ?
Posez votre question
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 ?
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 ?
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
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
Il n'était pas nécessaire de mettre tout ça...
Tu as écris dans le sujet que :
Donc V et T sont des vecteurs. Je demande simplement de quelle longueur ils sont. Donc il faut taper:
Obtiens-tu des nombres différents ?
P.S. : Si tu ne sais pas ce que tu as écrit, qui l'a écrit alors...?
Tu as écris dans le sujet que :
% 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;
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...?
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
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
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 :)
Re-teste ton programme une fois qu'il aura arrangé ça :)
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
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'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;
i=40:10:90;
T(i)=273.15+i
Tr=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 on
axis([0 1600 -40 60])
xlabel('Volume in cm3/mol')
ylabel('pressure in bar')
title('Isotherms for propane')
Matlab te le dit: écrire length(i) ici ça n'a pas de sens.
Tu devrais plutôt écrire:
Mais ça ne résoudra pas ton problème du calcul de l'équation de P car T n'a toujours pas la même longueur que V.
...et je n'ai pas l'impression que tu ais sérieusement vraiment cherché...
Tu devrais plutôt écrire:
i = 40:10:90;
T = 273.15 + i;
Tr = T./Tc;
RT = R.*T;
Mais ça ne résoudra pas ton problème du calcul de l'équation de P car T n'a toujours pas la même longueur que V.
...et je n'ai pas l'impression que tu ais sérieusement vraiment cherché...
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
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
Alors là, soit je ne comprends rien, soit tu m'expliques mal...
Ton but n'est-il pas d'avoir un vecteur pression P en fonction d'un vecteur température T et d'un vecteur volume V ?
Parce que dans ton exemple length(v)=2000 et length(T)=1 ce n'est pas du tout le cas vu que T est une constante.
Et tant donnée cette incompréhension totale, je reviens au base: dans ton équation P, que veux-tu faire varier (les variables) et que souhaites-tu fixer (les constantes) ?
Ton but n'est-il pas d'avoir un vecteur pression P en fonction d'un vecteur température T et d'un vecteur volume V ?
Parce que dans ton exemple length(v)=2000 et length(T)=1 ce n'est pas du tout le cas vu que T est une constante.
Et tant donnée cette incompréhension totale, je reviens au base: dans ton équation P, que veux-tu faire varier (les variables) et que souhaites-tu fixer (les constantes) ?
Mon but est de tracer P en fonction de V donc:
je prend une temperature T et je trace ma courbe
=> P=R*T(i)./(v - b) - a./(v.*(v + b)+b*(v-b));
avec R; b; a sont des constante
et v une variable
ensuite la température T change et je recommence
donc T est une constante
mon but est de tracer plusieurs fois P en fonction de V pour différentes valeurs de T
je prend une temperature T et je trace ma courbe
=> P=R*T(i)./(v - b) - a./(v.*(v + b)+b*(v-b));
avec R; b; a sont des constante
et v une variable
ensuite la température T change et je recommence
donc T est une constante
mon but est de tracer plusieurs fois P en fonction de V pour différentes valeurs de T
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');
end
hold off
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
j'obtiens de simple parabole alors que je dois avoir ca :
https://upload.wikimedia.org/wikipedia/commons/3/3e/Real_Gas_Isotherms.svg
Okay.
Je pense que maintenant ton script est juste. Le problème doit être dans les données maintenant: des températures fausses ou une constante a ou b mal calculée... Il va falloir que tu décortiques tout ça.
Je te conseille donc de prendre une calculette et de faire les calculs étape par étape et de comparer le résultats de ces étapes avec ceux de Matlab. Tu trouveras très vite où est l'erreur comme ça :)
Je pense que maintenant ton script est juste. Le problème doit être dans les données maintenant: des températures fausses ou une constante a ou b mal calculée... Il va falloir que tu décortiques tout ça.
Je te conseille donc de prendre une calculette et de faire les calculs étape par étape et de comparer le résultats de ces étapes avec ceux de Matlab. Tu trouveras très vite où est l'erreur comme ça :)
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
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
Salut Chevrotine,
Saches que ton idée est bonne et oui il y a des boucle while sous MATLAB :
http://www.mathworks.fr/fr/help/symbolic/mupad_ref/while.html
Mais pourquoi faire une condition sur le signe de S ? (Car tu écris while S>0 or S<0 ). Parce que là ça signifie que: Tant que S est positif ou négatif, je reste dans la boucle. Je ne pense pas que tu risques de sortir un jour de la boucle :-P. Donc à mon humble avis, la boucle while ne t'est pas utile.
Par contre une boucle for pourrait peut-être t'aider... ;-)
Et tu as aussi un moyen de calculer les intégrales:
http://www.mathworks.fr/fr/help/matlab/ref/integral.html
Par contre, comment ferais-tu pour récupérer Vmin et Vmax ?
Saches que ton idée est bonne et oui il y a des boucle while sous MATLAB :
http://www.mathworks.fr/fr/help/symbolic/mupad_ref/while.html
Mais pourquoi faire une condition sur le signe de S ? (Car tu écris while S>0 or S<0 ). Parce que là ça signifie que: Tant que S est positif ou négatif, je reste dans la boucle. Je ne pense pas que tu risques de sortir un jour de la boucle :-P. Donc à mon humble avis, la boucle while ne t'est pas utile.
Par contre une boucle for pourrait peut-être t'aider... ;-)
Et tu as aussi un moyen de calculer les intégrales:
http://www.mathworks.fr/fr/help/matlab/ref/integral.html
Par contre, comment ferais-tu pour récupérer Vmin et Vmax ?