>	restart:#"m04_p05".

Se va a hacer uso de la ecuación de van der Waals con las constantes determinadas a partir de las condiciones en el punto crítico, en lugar de por ajuste en una zona más amplia. Se pide:
a) Calcular a, b y R en función de Tcr, pcr y vcr.
b) Calcular a, b y vcr en función de Tcr, pcr y R. Esto es mejor porque así la R corresponde a la de los gases ideales y porque la vcr tiene gran incertidumbre experimental.
c) Expresión del factor de compresibilidad en función de las magnitudes reducidas.
d) Expresar la presión reducida en función de las otras magnitudes reducidas.
e) Puntos en los que las isotermas son horizontales en el diagrama Z-p y temperatura de Boyle, y comparación con la del diagrama generalizado de compresibilidad.
f) Curva de inversión y temperatura de inversión, y comparación con la del diagrama generalizado de compresibilidad.
g) Repetir el problema para otras ecuaciones de estado:Dietericci, Redlich-Kwong, Berthelot).
Datos:

>	#eq1:=Z=v/(v-b)-a/(v*R*T);#van der Waals

>	#eq1:=Z=(v/(v-b))*exp(-a/(v*R*T));#Dieterici

>	eq1:=Z=v/(v-b)-a/((v+b)*R*T^(3/2));#Redlich-Kwong

>	#eq1:=Z=v/(v-b)-a/(v*R*T^2);#Berthelot

>	eq2:=Z=p*v/(R*T);

a) Calcular a, b en función de Tcr y pcr.

>	assign(eq1):p:=collect(simplify(subs(v_b=v-b,collect(subs(v=v_b+b,solve(eq2,p)),a))),{T,a});

>	ppri:=diff(p,v):psec:=diff(ppri,v):

>	### WARNING: allvalues now returns a list of symbolic values instead of a sequence of lists of numeric values sol1:=solve(subs({T=Tcr,v=vcr},{p=pcr,ppri=0,psec=0}),{a,b,R}):sol1_:=allvalues(%)[1];sol1:=evalf(sol1_,3);

b) Calcular a, b y vcr en función de Tcr, pcr y R. Esto es mejor porque así la R corresponde a la de los gases ideales y porque la vcr tiene gran incertidumbre experimental.

>	sol1:=solve(subs({T=Tcr,v=vcr},{p=pcr,ppri=0,psec=0}),{a,b,vcr}):sol1_:=allvalues(%)[1];sol1:=evalf(sol1_,3);

c) Expresión del factor de compresibilidad en función de las magnitudes reducidas.

>	eqDefTR:=TR=T/Tcr;eqDefpR:=pR='p'/pcr;eqDefvR:=vR=v/vcr;Zcr:='pcr*vcr/(R*Tcr)';Zcr_:=subs(sol1_,Zcr);Zcr__:=evalf(Zcr_,2);

d) Expresar la presión reducida en función de las otras magnitudes reducidas.

>	Z_:=collect(subs(T=TR*Tcr,p=pR*pcr,v=vR*vcr,sol1,Z),vcr);eqpR:=pR='(Z/Zcr)*TR/vR';pR:=expand((Z_/Zcr_)*TR/vR);#Z__:=convert(map(simplify,{op(Z_)}),whattype(Z_));

Sólo para la de van der Waals.

>	#eqvR:=vR='(Z__/Zcr_)*TR/pR__';eqvR;eqZ3:=simplify(Z__*(vR-1/3)*vR=Z_*(vR-1/3)*vR);eqZ3:=expand(subs(eqvR,0=lhs(eqZ3)-rhs(eqZ3)));eqZ3_:=expand(eqZ3/coeff(rhs(eqZ3),Z__,3));

>	plot({seq(subs(TR=TR10/10,pR),TR10=6..15)},vR=0.34..3,'pR'=0..2,color=black);plot({seq(subs(TR=TR10/10,[pR,Z_,vR=0.34..10]),TR10=6..15)},'pR'=0..2,'Z'=0..1.1,color=black);

e)Puntos en los que las isotermas son horizontales en el diagrama Z-p y temperatura de Boyle, y comparación con la del diagrama generalizado de compresibilidad.

>	eqDef_TB:=Diff('Z','pR')[TR]=0;Diff('Z','pR')[TR]=Diff('Z',vR)[TR]/Diff('pR',vR)[TR];eqDef_TB:=Diff('Z',vR)[TR]=0;Diff('Z',vR)=diff(Z_,vR);TB_vR:=subs(Tcr=1,expand(solve(diff(Z_,vR)=0,TR)));TB:=evalf(convert(asympt(TB_vR,vR),polynom));'TB'=evalf(subs(vR=1e99,TB),2);

La del diagrama generalizado de compresibilidad es aprox TR=2,5.

El lugar geométrico de los mínimos es:

>	Z_Tmin:=simplify(subs(TR=TB_vR,Z_)):pR_Tmin:=simplify(subs(TR=TB_vR,pR)):

f) Curva de inversión y temperatura de inversión, y comparación con la del diagrama generalizado de compresibilidad.

>

eqJT:=mu=Diff(T,'p')[h];eqdh:=dh=c['p']*dT+(alpha*T-1)*v*dp;eqmu0:=T*Diff(v,T)['p']-v=0;eqmu0:=Diff('Z',TR)['pR']=0;eqmu0:=Diff('Z',TR)['vR']+Diff('Z',vR)['TR']*Diff(vR,TR)['pR']=0;dZ_dTR_vR:=diff(Z_,TR):dZ_dvR_TR:=diff(Z_,vR):dpR_dTR_vR:=diff(pR,TR):dpR_dvR_TR:=diff(pR,vR):dvR_dTR_pR:='-dpR_dTR_vR/dpR_dvR_TR';eqmu0:='dZ_dTR_vR+dZ_dvR_TR*dvR_dTR_pR=0';TI_vR:=(solve(subs(Tcr=1,pcr=1,vcr=1,eqmu0),TR))[1]:TI_:=subs(vR=1e9,TI_vR);

En el diagrama de correcciones de h se atisba TI=5.

La curva de inversión, en paramétricas, será:

>	pR_TI:=subs(TR=TI_vR,pR):

>	plot(subs(Tcr=1,[TI_vR,pR_TI,vR=.334..100]),Temp=0..TI_,pres=0..20,numpoints=100);

Respecto a las constantes a y b,

>	eqmu:=mu[JT]='(T*(-Diff(p,T)/Diff(p,v))-v)/cp';eqmu:=value(eqmu):eqmu0_lin:=0=expand(convert(series(simplify(expand(subs(v=1/rho,rhs(eqmu)))),rho,1),polynom));

>