Error with vpa solve

29 次查看（过去 30 天）
Alessio Falcone
Alessio Falcone 2021-11-9
编辑： Alessio Falcone 2021-11-9
Hello everyone.
The program computes for the first Tc but then gives me an error. Can you help me?
Tb=112; %K
Z=0.27;
R=0.0821;
DH=8200; %J/mol
for Tc=120:0.01:200
%clapeyron Pc
Pc=exp((DH/8.314)*((1/Tb)-(1/Tc)));
syms a b c Vc;
X=R*Tc/(Vc-b)-(a/Tc)/((Vc+c)^2)-Pc;
eqn1=R*Tc/(Vc-b)-(a/Tc)/((Vc+c)^2)-Pc==0;
eqn2=diff(X,Vc)==0;
eqn3=diff(X,Vc,2)==0;
eqn4=(Pc*Vc)/(R*Tc)-Z==0;
sol=vpasolve([eqn1,eqn2,eqn3,eqn4],[a,b,c,Vc]);
Vc=sol.Vc;
a=sol.a
b=sol.b
c=sol.c
syms V
P_star=1;
Pb=R*Tb/(V-b)-(a/Tb)/((V+c)^2);
Z=R*Tb/(V-b)-(a/Tb)/((V+c)^2)-P_star==0;
V_soluzione=real(vpasolve(Z,V));
Vliq=min(V_soluzione);
Vvs=max(V_soluzione);
Y=Pb-P_star;
integralearee=real(vpaintegral(Y,Vliq,Vvs))
if integralearee<=0
break
end
end
vol_liq_sat=Vliq
vol_vap_sat=Vvs
  1 个评论
Steven Lord
Steven Lord 2021-11-9
Please show us the the full and exact text of the error and/or warning messages (all the text displayed in red and/or orange in the Command Window.)

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

David Hill
David Hill 2021-11-9
Tc=fzero(@I,126);
[~,Vliq,Vvs]=I(Tc);
function [integralearee,Vliq,Vvs]=I(Tc)
Tb=112; %K
Z=0.27;
R=0.0821;
DH=8200; %J/mol
syms a b c Vc V;
Pc=exp((DH/8.314)*((1/Tb)-(1/Tc)));
X=R*Tc/(Vc-b)-(a/Tc)/((Vc+c)^2)-Pc;
eqn1=R*Tc/(Vc-b)-(a/Tc)/((Vc+c)^2)-Pc==0;
eqn2=diff(X,Vc)==0;
eqn3=diff(X,Vc,2)==0;
eqn4=(Pc*Vc)/(R*Tc)-Z==0;
sol=vpasolve([eqn1,eqn2,eqn3,eqn4],[a,b,c,Vc]);
a=sol.a;
b=sol.b;
c=sol.c;
P_star=1;
Pb=R*Tb/(V-b)-(a/Tb)/((V+c)^2);
Z=R*Tb/(V-b)-(a/Tb)/((V+c)^2)-P_star==0;
V_soluzione=real(vpasolve(Z,V));
Vliq=min(V_soluzione);
Vvs=max(V_soluzione);
Y=Pb-P_star;
integralearee=real(vpaintegral(Y,Vliq,Vvs));
end
  1 个评论
Alessio Falcone
Alessio Falcone 2021-11-9
编辑：Alessio Falcone 2021-11-9
Hi thank you very much for the reply.
How do I set that the program must stop when 'integralearee' is equal to 0 or very close to it? furthermore 'Vliq' must be approximately equal to 0.036. Then I have to show the Tc found on the display, thank you very much

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Assumptions 的更多信息

标签

产品


版本

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by