>	restart;#"m07_p25"

Por un conducto circula a presión ambiente una mezcla equimolar de propano/butano. Se pide:

a) Intervalo de temperatura en el que podría haber mezcla bifásica si no se conocieran las concentraciones.

b) Intervalo de temperatura en el que habrá mezcla bifásica.

c) Entalpía y capacidad térmica y densidad en función de la temperatura.

Datos:

>	read"../therm_eq.m":read"../therm_proc.m":with(therm_proc):

>	su1:="C3H8":su2:="C4H10":dat1:=[x[B]=0.5];

Esquema:

>

>

Eqs. const.:

>	Pdat:=get_gas_data(su1),get_liq_data(su1):Bdat:=get_gas_data(su2),get_liq_data(su2):dat0:=Const,SI2,SI1:pv1:=proc(T) global su1;get_pv_data(su1);RETURN(pv(T)):end:pv2:=proc(T) global su1;get_pv_data(su2);RETURN(pv(T)):end:

Genérico:

>	eqNV:=xv1+xv2=1;eqNL:=xl1+xl2=1;eqC1:=xv1*xv0+xl1*xl0=x01;eqC2:=xv2*xv0+xl2*xl0=x02;eqE1:=xv1/xl1=pv1/p;eqE2:=xv2/xl2=pv2/p;sol1:=solve({eqNV,eqNL,eqE1,eqE2},{xv1,xv2,xl1,xl2});sol2:=solve({eqC1,eqC2},{xv0,xl0});

a) Intervalo de temperatura en el que podría haber mezcla bifásica si no se conocieran las concentraciones.

Para mezclas ideales, es el comprendido entre las temperaturas de ebullición de los compuestos puros.

Diagrama T-x a p0=100 kPa.

>

n:=19:T1_:=evalf(subs(dat0,solve(p0=pv1(T),T)));T2_:=evalf(subs(dat0,solve(p0=pv2(T),T)));for i from 0 to n do T||i:=T1_+(T2_-T1_)*(i/n);xv1||i:=subs(sol1,p=p0,pv1=pv1(T||i),pv2=pv2(T||i),dat0,SI0,xv1);xl1||i:=subs(sol1,p=p0,pv1=pv1(T||i),pv2=pv2(T||i),dat0,SI0,xl1);od:plot({[seq([xv1||i,T||i/K_],i=0..n)],[seq([xl1||i,T||i/K_],i=0..n)]},x[B]=0..1,'T'=200..300,color=black,axes=boxed);#for i from 1 to n do print(T||i,xv1||i,xl1||i):od:

b) Intervalo de temperatura en el que habrá mezcla bifásica.

>	x1_:=subs(dat1,x[B]);Mm:='x[P]*M[P]+x[B]*M[B]';Mm_:=subs(x[P]=1-x[B],M[P]=M,M[B]=M_,Pdat,M_=M,Bdat,dat1,Mm);T[dew_]:=fsolve(x1_=subs(sol1,p=p0,pv1=pv1(T),pv2=pv2(T),dat0,SI0,xv1),T=100..500)*K_;T[boil_]:=fsolve(x1_=subs(sol1,p=p0,pv1=pv1(T),pv2=pv2(T),dat0,SI0,xl1),T=100..500)*K_;

c) Entalpía, capacidad térmica y densidad en función de la temperatura.

Cálculo de hmix, cpmix y rhomix para fase liq. (T<Tboil), bifásico, y fase gas (T>Tdew):

Elegimos como referencia h=0 a p0=100 kPa y T25=25 ºC (ambos gaseosos), aunque saldrán casi todos valores negativos.

Mezcla gaseosa:

>

cpmix_gas:='x[P]*cpv1+x[B]*cpv2';cpv1:=subs(Pdat,c[p]*M);cpv2:=subs(Bdat,c[p]*M);cpmix_gas_:=subs(x[P]=1-x[B],dat1,dat0,cpmix_gas);hmix_gas:='cpmix_gas*(T-T25)';hmix_gas_:=subs(x[P]=1-x[B],dat1,dat0,hmix_gas):hmix_gas_dew:=subs(T=T[dew_],dat0,%);vmix_gas:='x[P]*vv1+x[B]*vv2';vv:=R[u]*T/p0;vv1:=subs(dat0,Pdat,R[u]*T/p0):vv2:=subs(dat0,Bdat,R[u]*T/p0):vmix_gas_:=subs(x[P]=1-x[B],dat1,dat0,vmix_gas):rhomix_gas_:=Mm_/vmix_gas_:rhomix_gas_dew_:=subs(T=T[dew_],rhomix_gas_);

Mezcla bifásica:

>

hmix:='xv0*(xv1*hv1+xv2*hv2)+xl0*(xl1*hl1+xl2*hl2)';hv1:='cpv1*(T-T25)';hv1_:=subs(dat0,hv1):hv2:='cpv2*(T-T25)';hv2_:=subs(dat0,hv2):hl1:='-hlv251+cl1*(T-T25)';cl1:=subs(Pdat,c*M);hlv251:=subs(Pdat,T=T25,dat0,hlv(T)*M);hl1_:=subs(Pdat,dat0,hl1):hl2:='-hlv252+cl2*(T-T25)';cl2:=subs(Bdat,c*M);hlv252:=subs(Bdat,T=T25,dat0,hlv(T)*M);hl2_:=subs(Bdat,dat0,hl2):cpmix:=Diff('hmix',T);cpmix:=Dhmix/DT;vmix:='xv0*(xv1*vv1+xv2*vv2)+xl0*(xl1*vl1+xl2*vl2)';vl1:=subs(Pdat,M/rho);vl2:=subs(Bdat,M/rho);

hay que calcular las xi para cada T, y con ellas la hamix y la vmix:

>

for i from 0 to n do T||i:=T[dew_]+(T[boil_]-T[dew_])*(i/n);xv1||i:=subs(sol1,p=p0,pv1=pv1(T||i),pv2=pv2(T||i),dat0,SI0,xv1);xl1||i:=subs(sol1,p=p0,pv1=pv1(T||i),pv2=pv2(T||i),dat0,SI0,xl1);xv0||i:=subs(sol2,sol1,x01=x1_,x02=1-x1_,p=p0,pv1=pv1(T||i),pv2=pv2(T||i),dat0,SI0,xv0);hmix||i:=subs(T=T||i,dat0,xv0||i*(xv1||i*hv1+(1-xv1||i)*hv2)+(1-xv0||i)*(xl1||i*hl1+(1-xl1||i)*hl2));vmix||i:=subs(T=T||i,dat0,xv0||i*(xv1||i*vv1+(1-xv1||i)*vv2)+(1-xv0||i)*(xl1||i*vl1+(1-xl1||i)*vl2));rhomix||i:=Mm_/vmix||i;od:

y ahora el cp=Dh/DT (rellenamos el último para que tengan el mismo rango para dibujar):

>	for i from 0 to n-1 do cpmix||i:=(hmix||(i+1)-hmix||i)/(T||(i+1)-T||i);od:cpmix||n:=cpmix||(n-1):

Mezcla líquida:

>

hmix_liq:='(x[P]*cl1+x[B]*cl2)*(T-T25)-x[P]*hlv251-x[B]*hlv252';hlv251:=subs(Pdat,T=T25,dat0,hlv(T)*M):hlv252:=subs(Bdat,T=T25,dat0,hlv(T)*M):hmix_liq_:=subs(x[P]=1-x[B],dat1,dat0,hmix_liq):hmix_liq_boil:=subs(T=T[boil_],dat0,%);cpmix_liq:='x[P]*cl1+x[B]*cl2';cpmix_liq_:=subs(x[P]=1-x[B],dat1,dat0,cpmix_liq);vmix_liq:='x[P]*vl1+x[B]*vl2';vmix_liq_:=subs(x[P]=1-x[B],dat1,dat0,vmix_liq);rhomix_liq_:=Mm_/vmix_liq_:rhomix_liq_boil_:=subs(T=T[boil_],rhomix_liq_);

>

n_:=1:hgas:=[seq([T[dew_]+(T2_-T[dew_])*(i/n_),subs(T=T[dew_]+(T2_-T[dew_])*(i/n_),dat0,hmix_gas_)],i=0..n_)]:hliq:=[seq([T1_+(T[boil_]-T1_)*(i/n_),subs(T=T1_+(T[boil_]-T1_)*(i/n_),dat0,hmix_liq_)],i=0..n_)]:rhogas:=[seq([T[dew_]+(T2_-T[dew_])*(i/n_),subs(T=T[dew_]+(T2_-T[dew_])*(i/n_),dat0,rhomix_gas_)],i=0..n_)]:rholiq:=[seq([T1_+(T[boil_]-T1_)*(i/n_),subs(T=T1_+(T[boil_]-T1_)*(i/n_),dat0,rhomix_liq_)],i=0..n_)]:cpgas:=[seq([T[dew_]+(T2_-T[dew_])*(i/n_),subs(T=T[dew_]+(T2_-T[dew_])*(i/n_),dat0,cpmix_gas_)],i=0..n_)]:cpliq:=[seq([T1_+(T[boil_]-T1_)*(i/n_),subs(T=T1_+(T[boil_]-T1_)*(i/n_),dat0,cpmix_liq_)],i=0..n_)]:

>	hbif:=[seq([T||i,hmix||i],i=0..n)]:cpbif:=[seq([T||i,cpmix||i],i=0..n)]:rhobif:=[seq([T||i,rhomix||i],i=0..n)]:

>	plot(subs(SI0,{hliq,hbif,hgas}),T_K=T1_/K_..T2_/K_,h_J_mol=-30000..0,color=black);

>	plot(subs(SI0,{cpliq,cpbif,cpgas}),T_K=T1_/K_..T2_/K_,cp_J_molK=0..2500,color=black);

>	plot(subs(SI0,{rholiq,rhobif,rhogas}),T_K=T1_/K_..T2_/K_,rho_kg_m3=0..600,color=black);

Nótese que la densidad cae muy bruscamente en cuanto empieza la ebullición.

>