Debutant matlab

Résolu/Fermé
chevrotine22 Messages postés 21 Date d'inscription mardi 10 décembre 2013 Statut Membre Dernière intervention 14 septembre 2017 - 14 mai 2014 à 14:18
rodfaction Messages postés 1 Date d'inscription jeudi 28 janvier 2016 Statut Membre Dernière intervention 28 janvier 2016 - 28 janv. 2016 à 12:22
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

14 réponses

JulienJust Messages postés 139 Date d'inscription mardi 25 juin 2013 Statut Membre Dernière intervention 2 septembre 2014 18
14 mai 2014 à 14:40
Bonjour,

Quelle erreur s'affiche pour ce script du coup ?

Cdlt
0
chevrotine22 Messages postés 21 Date d'inscription mardi 10 décembre 2013 Statut Membre Dernière intervention 14 septembre 2017
14 mai 2014 à 14:43
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));
0
JulienJust Messages postés 139 Date d'inscription mardi 25 juin 2013 Statut Membre Dernière intervention 2 septembre 2014 18
Modifié par JulienJust le 14/05/2014 à 15:01
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 ?
0
chevrotine22 Messages postés 21 Date d'inscription mardi 10 décembre 2013 Statut Membre Dernière intervention 14 septembre 2017
14 mai 2014 à 15:03
non y a toujours la meme erreur
0

Vous n’avez pas trouvé la réponse que vous recherchez ?

Posez votre question
JulienJust Messages postés 139 Date d'inscription mardi 25 juin 2013 Statut Membre Dernière intervention 2 septembre 2014 18
14 mai 2014 à 15:20
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 ?
0
chevrotine22 Messages postés 21 Date d'inscription mardi 10 décembre 2013 Statut Membre Dernière intervention 14 septembre 2017
14 mai 2014 à 15:30
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
0
JulienJust Messages postés 139 Date d'inscription mardi 25 juin 2013 Statut Membre Dernière intervention 2 septembre 2014 18
14 mai 2014 à 16:17
Il n'était pas nécessaire de mettre tout ça...

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...?
0
chevrotine22 Messages postés 21 Date d'inscription mardi 10 décembre 2013 Statut Membre Dernière intervention 14 septembre 2017
14 mai 2014 à 18:20
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
0
JulienJust Messages postés 139 Date d'inscription mardi 25 juin 2013 Statut Membre Dernière intervention 2 septembre 2014 18
14 mai 2014 à 20:40
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 :)
0
chevrotine22 Messages postés 21 Date d'inscription mardi 10 décembre 2013 Statut Membre Dernière intervention 14 septembre 2017
17 mai 2014 à 18:09
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'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')
0
JulienJust Messages postés 139 Date d'inscription mardi 25 juin 2013 Statut Membre Dernière intervention 2 septembre 2014 18
18 mai 2014 à 01:03
Matlab te le dit: écrire length(i) ici ça n'a pas de sens.

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é...
0
chevrotine22 Messages postés 21 Date d'inscription mardi 10 décembre 2013 Statut Membre Dernière intervention 14 septembre 2017
18 mai 2014 à 14:48
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
0
JulienJust Messages postés 139 Date d'inscription mardi 25 juin 2013 Statut Membre Dernière intervention 2 septembre 2014 18
19 mai 2014 à 09:26
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) ?
0
chevrotine22 Messages postés 21 Date d'inscription mardi 10 décembre 2013 Statut Membre Dernière intervention 14 septembre 2017
19 mai 2014 à 10:47
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
0
JulienJust Messages postés 139 Date d'inscription mardi 25 juin 2013 Statut Membre Dernière intervention 2 septembre 2014 18
Modifié par JulienJust le 19/05/2014 à 11:42
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
0
chevrotine22 Messages postés 21 Date d'inscription mardi 10 décembre 2013 Statut Membre Dernière intervention 14 septembre 2017
19 mai 2014 à 19:23
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
0
JulienJust Messages postés 139 Date d'inscription mardi 25 juin 2013 Statut Membre Dernière intervention 2 septembre 2014 18
20 mai 2014 à 09:32
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 :)
0
chevrotine22 Messages postés 21 Date d'inscription mardi 10 décembre 2013 Statut Membre Dernière intervention 14 septembre 2017
20 mai 2014 à 22:54
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
0
JulienJust Messages postés 139 Date d'inscription mardi 25 juin 2013 Statut Membre Dernière intervention 2 septembre 2014 18
21 mai 2014 à 09:36
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 ?
0
chevrotine22 Messages postés 21 Date d'inscription mardi 10 décembre 2013 Statut Membre Dernière intervention 14 septembre 2017
21 mai 2014 à 10:40
je sais pas comment récupérer Vmin et Vmax donc si tu a une idée je suis preuneur
0
JulienJust Messages postés 139 Date d'inscription mardi 25 juin 2013 Statut Membre Dernière intervention 2 septembre 2014 18
21 mai 2014 à 10:48
je pense que tu devrais déjà commencer par le début en faisant ce que tu as dit:

"maintenant je dois tracer une droite parallèle a l'axe X, elle va être coupé en 3 points, ce qui fait 2 aires. "
0
rodfaction Messages postés 1 Date d'inscription jeudi 28 janvier 2016 Statut Membre Dernière intervention 28 janvier 2016
Modifié par rodfaction le 28/01/2016 à 13:12
Bonjour Chevrotine22,
Avez vous trouvé une solution ? J'ai exactement le même besoin que vous.
Pouvez vous me partager votre code final svp ?
Merci
0