>	restart:#"m16_p09"

>	read`../therm_eq.m`:read`../therm_chem.m`:with(therm_chem);with(therm_proc):

Determinar la composición de equilibrio y la temperatura alcanzada en la combustión estequiométrica, adiabática e isobárica de propano en aire y en oxígeno puro, en los siguientes supuestos:
a) Suponiendo que sólo aparece CO2 y H2O.
b) Suponiendo que además aparece CO, H2 y O2.
Datos:

>	su1:="Aire":su2:="H2O":fuel:="C3H8":dat:=[Mf_=0.044*kg_/mol_]:

Eqs. const.:

>	Adat:=get_gas_data(su1):Wdat:=get_gas_data(su2),get_liq_data(su2):get_pv_data(su2):dat:=op(dat),op(subs(g=g0,[Const])),Adat,SI2,SI1:

a) Suponiendo que sólo aparece CO2, H2O

Con O2 puro

>	eqM:=eqMIX(a*C3H8+b*O2=[4,6]);sol1_O2:=evalf(solve(subs(dat,{eqNX,eqBC,eqBH,eqBO}),{a,b,x[Comp[4]],x[Comp[6]]}));PCI_:=PCI(subs(sol1_O2,eqM/a));Ta_:=subs(cpComp_,sol1_O2,dat,T25+a*PCI_/sum(delta[i]*x[Comp[i]]*c[p,Comp[i]],i=1..C_));subs(sol1_O2,eqM/a);

i.e., con O2 puro y sin disociaciones saldría a 6000 K que es incompatible con no haber disociaciones.

a) Suponiendo que sólo aparece CO2 y H2O

Con aire

>	eqM:=eqMIX(a*C3H8+b*(c21*O2+c79*N2)=[3,4,6]);sol1_A:=solve(subs(dat,{eqNX,eqBC,eqBH,eqBO,eqBN}),{a,b,x[Comp[3]],x[Comp[4]],x[Comp[6]]});eqA0:=Ateo(fuel);PCI_:=PCI(subs(sol1_A,eqM/a));Ta_:=subs(cpComp_,sol1_A,dat,T25+a*PCI_/sum(delta[i]*x[Comp[i]]*c[p,Comp[i]],i=1..C_));

b) Suponiendo que además aparece CO, H2 y O2.

Con O2 puro. Pruebo T hasta que Qs=0.

>	eq1:=subs(p=p0,eqEQ(CO+H2O=CO2+H2));eq2:=subs(p=p0,eqEQ(CO+(1/2)*O2=CO2));

Con O2 puro.

>

eqM:=eqMIX(a*C3H8+5*a*O2=[2,4,6,7,8]);n:=7:arr:=array(1..n,1..7):for i from 1 to n do T:=evalf( 1500+2500*(i/n))*K_;sol1_:=fsolve(subs(dat,{eqNX,eqBC,eqBH,eqBO,eq1,eq2}),{a,x[Comp[2]],x[Comp[4]],x[Comp[6]],x[Comp[7]],x[Comp[8]]},{a=0..1,x[Comp[2]]=0..1,x[Comp[4]]=0..1,x[Comp[6]]=0..1,x[Comp[7]]=0..1,x[Comp[8]]=0..1});PCI_:=PCI(subs(sol1_,eqM/a));eqBE_:=q[s]/a=subs(cpComp_,sol1_,dat,PCI_-sum(delta[ii]*x[Comp[ii]]*c[p,Comp[ii]],ii=1..C_)*(T-T25)/a);Qs_:=subs(sol1_,rhs(eqBE_)/a);arr[i,1]:=T/K_;arr[i,2]:=subs(SI0,Qs_)/1e6;arr[i,3]:=subs(sol1_,x[Comp[2]]);arr[i,4]:=subs(sol1_,x[Comp[4]]);arr[i,5]:=subs(sol1_,x[Comp[6]]);arr[i,6]:=subs(sol1_,x[Comp[7]]);arr[i,7]:=subs(sol1_,x[Comp[8]]);od: pl2:=pla(arr,2,1):pl3:=pla(arr,3,1):pl4:=pla(arr,4,1):pl5:=pla(arr,5,1):pl6:=pla(arr,6,1):pl7:=pla(arr,7,1):plot([pl2,pl3,pl4,pl5,pl6,pl7],Temperature=0..4000,x=0..1,color=black);

i.e. Tad=3630 K (donde Qs=0).

>

T:=3630*K_;sol1_:=fsolve(subs(dat,{eqNX,eqBC,eqBH,eqBO,eq1,eq2}),{a,x[Comp[2]],x[Comp[4]],x[Comp[6]],x[Comp[7]],x[Comp[8]]},{a=0..1,x[Comp[2]]=0..1,x[Comp[4]]=0..1,x[Comp[6]]=0..1,x[Comp[7]]=0..1,x[Comp[8]]=0..1});PCI_:=PCI(subs(sol1_,eqM/a));eqBE_:=q[s]/a=subs(cpComp_,sol1_,dat,PCI_-sum(delta[ii]*x[Comp[ii]]*c[p,Comp[ii]],ii=1..C_)*(T-T25)/a);Qs_:=subs(sol1_,rhs(eqBE_)/a);

Con aire.

>

Digits:=15:eqM:=eqMIX(a*C3H8+a*rhs(eqA0)*(c21*O2+c79*N2)=[2,3,4,6,7,8]);n:=5:arr:=array(1..n,1..8):for i from 1 to n do T:=evalf( 1500+2000*(i/n))*K_;sol1_:=fsolve(subs(dat,{eqNX,eqBC,eqBH,eqBO,eqBN,eq1,eq2}),{a,x[Comp[2]],x[Comp[3]],x[Comp[4]],x[Comp[6]],x[Comp[7]],x[Comp[8]]},{a=0..1,x[Comp[2]]=0..1,x[Comp[3]]=0..1,x[Comp[4]]=0..1,x[Comp[6]]=0..1,x[Comp[7]]=0..1,x[Comp[8]]=0..1});PCI_:=PCI(subs(sol1_,eqM/a));eqBE_:=q[s]/a=subs(cpComp_,sol1_,dat,PCI_-sum(delta[ii]*x[Comp[ii]]*c[p,Comp[ii]],ii=1..C_)*(T-T25)/a);Qs_:=subs(sol1_,rhs(eqBE_)/a);arr[i,1]:=T/K_;arr[i,2]:=subs(SI0,Qs_)/1e6;arr[i,3]:=subs(sol1_,x[Comp[2]]);arr[i,4]:=subs(sol1_,x[Comp[4]]);arr[i,5]:=subs(sol1_,x[Comp[6]]);arr[i,6]:=subs(sol1_,x[Comp[7]]);arr[i,7]:=subs(sol1_,x[Comp[8]]);arr[i,8]:=subs(sol1_,x[Comp[3]]);od: pl2:=pla(arr,2,1):pl3:=pla(arr,3,1):pl4:=pla(arr,4,1):pl5:=pla(arr,5,1):pl6:=pla(arr,6,1):pl7:=pla(arr,7,1):pl8:=pla(arr,8,1):plot([pl2,pl3,pl4,pl5,pl6,pl7,pl8],Temperature=0..4000,x=0..1,color=black);

i.e. la Tad es 2265 K (donde Qs=0):

>

Digits:=4:T:=2265*K_;sol1_:=fsolve(subs(dat,{eqNX,eqBC,eqBH,eqBO,eqBN,eq1,eq2}),{a,x[Comp[2]],x[Comp[3]],x[Comp[4]],x[Comp[6]],x[Comp[7]],x[Comp[8]]},{a=0..1,x[Comp[2]]=0..1,x[Comp[3]]=0..1,x[Comp[4]]=0..1,x[Comp[6]]=0..1,x[Comp[7]]=0..1,x[Comp[8]]=0..1});PCI_:=PCI(subs(sol1_,eqM/a));eqBE_:=q[s]/a=subs(cpComp_,sol1_,dat,PCI_-sum(delta[ii]*x[Comp[ii]]*c[p,Comp[ii]],ii=1..C_)*(T-T25)/a);Qs_:=subs(sol1_,rhs(eqBE_)/a);

>