> restart:#"m09_p54"

> read`../therm_chem.m`:with(therm_chem);with(therm_proc):interface(displayprecision=2):

[Ateo, Mf, PCI, PCS, eqEQ, eqMIX, eq_fit, get_hgs_data, hgs_r25, nulist, seqEBE]

Considérese la combustión estequiométrica de hidrógeno con oxígeno puro a presión constante, entrando los gases a 1 MPa y 0 ºC. Se pide:
a) Calor a evacuar para que la salida sea a 2500 K, suponiendo que no hay disociación.
b) Composición de equilibrio a 2500 K si se disocia algo de vapor en sus componentes.
c) Composición de equilibrio a 2500 K si la disociación también genera OH.
d) Temperatura adiabática suponiendo que no hubiera disociación.
e) Temperatura adiabática y composición de equilibrio con disociación.
f) Exceso de hidrógeno que haría que la temperatura de combustión adiabática sin disociación fuese 2500 K.

Datos:

> su1:="H2":su2:="O2":su3:="H2O":dat:=[p1=1e6*Pa_,T1=(0+273.15)*K_,T2=2500*K_];

[p1 = `+`(`*`(0.1e7, `*`(Pa_))), T1 = `+`(`*`(273.15, `*`(K_))), T2 = `+`(`*`(2500, `*`(K_)))]

Eqs. const.:

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

a) Calor a evacuar para que la salida sea a 2500 K, suponiendo que no hay disociación.

El efecto de la presión es nulo con el modelo de gases ideales, y el de la temperatura de entrada es despreciable en el poder calorífico, así que bastará tenerlo en cuenta en el calentamiento de los productos. Tomaremos valores medios usuales para las capacidades térmicas en caliente.

> eq15_6_2;qs:=PCI-c[p,H2O]*(T2-T1);eq0:=H2+(1/2)*O2=H2O;PCImol_:=PCI(eq0);eqcp:=cpComp[6];qs_:=subs(PCI=PCImol_,eqcp,dat,qs);

q[s] = `+`(`*`(a, `*`(PCI)), `-`(`*`(Sum(`*`(x[Com[i]], `*`(c[p, i])), i = 1 .. CP), `*`(`+`(T, `-`(T25))))))
`+`(PCI, `-`(`*`(c[p, H2O], `*`(`+`(T2, `-`(T1))))))
`+`(H2, `*`(`/`(1, 2), `*`(O2))) = H2O
`+`(`/`(`*`(241820.00, `*`(J_)), `*`(mol_)))
c[p, H2O] = `+`(`/`(`*`(47, `*`(J_)), `*`(mol_, `*`(K_))))
`+`(`/`(`*`(137158.05, `*`(J_)), `*`(mol_))) (1)

i.e. si la combustión es completa habría que evacuar 137 kJ/mol (137 kW por cada mol/s de H2 entrante, o de vapor formado).

b) Composición de equilibrio a 2500 K si se disocia algo de vapor en sus componentes.

> eq1:=eqMIX(a*H2+a*(1/2)*O2=[2,6,8]);eqNX;eqBH;eqBO;eqEQ(eq0);eqEQ_:=evalf(subs(p=p1,T=T2,dat,eqEQ(eq0)));sol_:=solve({eqNX,eqBH,eqBO,eqEQ_},{a,x[H2O],x[H2],x[O2]});

`+`(`*`(a, `*`(H2)), `*`(`/`(1, 2), `*`(a, `*`(O2)))) = `+`(`*`(H2, `*`(x[H2])), `*`(H2O, `*`(x[H2O])), `*`(O2, `*`(x[O2])))
1 = `+`(x[O2], x[H2O], x[H2])
0 = `+`(`*`(2, `*`(x[H2O])), `-`(`*`(2, `*`(a))), `*`(2, `*`(x[H2])))
0 = `+`(`-`(a), `*`(2, `*`(x[O2])), x[H2O])
`/`(`*`(x[H2O]), `*`(`^`(x[O2], `/`(1, 2)), `*`(x[H2]))) = `+`(`/`(`*`(0.4809227354e-2, `*`(exp(`+`(`/`(`*`(29085.87925, `*`(K_)), `*`(T)))))), `*`(`^`(`/`(`*`(p0), `*`(p)), `/`(1, 2)))))
`/`(`*`(x[H2O]), `*`(`^`(x[O2], `/`(1, 2)), `*`(x[H2]))) = 1717.157188
{a = .9956452236, x[H2] = 0.8709552895e-2, x[H2O] = .9869356707, x[O2] = 0.4354776447e-2} (2)

i.e. en equilibrio a 2500 K habría en base molar 98,7 %H2O, 0,9 %H2, y 0,4 %O2.

c) Composición de equilibrio a 2500 K si la disociación también genera OH.

> eq2:=eqMIX(a*H2+a*(1/2)*O2=[2,6,8,9]);eqNX;eqBH;eqBO;eq3:=H2+O2=2*OH;eqEQ0_:=evalf(subs(p=p1,T=T2,dat,eqEQ(eq0)));eqEQ1_:=evalf(subs(p=p1,T=T2,dat,eqEQ(eq3)));sol2_:=solve({eqNX,eqBH,eqBO,eqEQ0_,eqEQ1_},{a,x[H2O],x[H2],x[O2],x[OH]})[1];

`+`(`*`(a, `*`(H2)), `*`(`/`(1, 2), `*`(a, `*`(O2)))) = `+`(`*`(H2, `*`(x[H2])), `*`(H2O, `*`(x[H2O])), `*`(O2, `*`(x[O2])), `*`(OH, `*`(x[OH])))
1 = `+`(x[O2], x[H2O], x[H2], x[OH])
0 = `+`(`*`(2, `*`(x[H2O])), `-`(`*`(2, `*`(a))), `*`(2, `*`(x[H2])), x[OH])
0 = `+`(`-`(a), `*`(2, `*`(x[O2])), x[H2O], x[OH])
`+`(H2, O2) = `+`(`*`(2, `*`(OH)))
`/`(`*`(x[H2O]), `*`(`^`(x[O2], `/`(1, 2)), `*`(x[H2]))) = 1717.157188
`/`(`*`(`^`(x[OH], 2)), `*`(x[O2], `*`(x[H2]))) = 1.465857644
{a = .9932784755, x[H2] = 0.9986407405e-2, x[H2O] = .9798354265, x[O2] = 0.3264882919e-2, x[OH] = 0.6913283132e-2} (3)

i.e. se disocia más; ahora, en equilibrio a 2500 K habría en base molar 98,0 %H2O, 1,0 %H2, 0,3 %O2, y 0,7 %OH.

d) Temperatura adiabática suponiendo que no hubiera disociación.

> eq15_7_2;eq0;eqTa:=Ta=T1+PCI/c[p,H2O];eqTa_:=subs(eqcp,PCI=PCImol_,dat,%);

Ta = `+`(T25, `/`(`*`(a, `*`(PCI)), `*`(Sum(`*`(x[Com[i]], `*`(c[p, i])), i = 1 .. CP))))
`+`(H2, `*`(`/`(1, 2), `*`(O2))) = H2O
Ta = `+`(T1, `/`(`*`(PCI), `*`(c[p, H2O])))
Ta = `+`(`*`(5418.256383, `*`(K_))) (4)

i.e. demasiado alta para no haber disociación. Sería algo menor si se tomara un cp-medio más adecuado, algo mayor de los 47 J/(mol·K) usado, pero aun así demasiado alta, i.e. hay disociación.

e) Temperatura adiabática y composición de equilibrio con disociación.

Ahora tenemos 7 incógnitas: a, x[H2O], x[H2], x[O2], x[OH], PCI, y T.

> eq2;eqNX;eqBH;eqBO;eqEQ0_:=evalf(subs(p=p1,dat,eqEQ(eq0)));eqEQ1_:=evalf(subs(p=p1,dat,eqEQ(eq3)));PCI_:=PCI(eq2);eqTa:=T=subs(T25=T1,rhs(eq15_7_2));eqTa_:=subs(Ta=T,T25=T1,cpComp,eq15_7_3);

`+`(`*`(a, `*`(H2)), `*`(`/`(1, 2), `*`(a, `*`(O2)))) = `+`(`*`(H2, `*`(x[H2])), `*`(H2O, `*`(x[H2O])), `*`(O2, `*`(x[O2])), `*`(OH, `*`(x[OH])))
1 = `+`(x[O2], x[H2O], x[H2], x[OH])
0 = `+`(`*`(2, `*`(x[H2O])), `-`(`*`(2, `*`(a))), `*`(2, `*`(x[H2])), x[OH])
0 = `+`(`-`(a), `*`(2, `*`(x[O2])), x[H2O], x[OH])
`/`(`*`(x[H2O]), `*`(`^`(x[O2], `/`(1, 2)), `*`(x[H2]))) = `+`(`*`(0.1520811222e-1, `*`(exp(`+`(`/`(`*`(29085.87925, `*`(K_)), `*`(T)))))))
`/`(`*`(`^`(x[OH], 2)), `*`(x[O2], `*`(x[H2]))) = `+`(`*`(65.32726365, `*`(exp(`+`(`-`(`/`(`*`(9492.422423, `*`(K_)), `*`(T))))))))
`+`(`/`(`*`(241820.00, `*`(x[H2O], `*`(J_))), `*`(mol_)), `-`(`/`(`*`(39460.00, `*`(x[OH], `*`(J_))), `*`(mol_))))
T = `+`(T1, `/`(`*`(a, `*`(PCI)), `*`(Sum(`*`(x[Com[i]], `*`(c[p, i])), i = 1 .. CP))))
T = `+`(T1, `/`(`*`(a, `*`(`+`(`/`(`*`(241820.00, `*`(x[H2O], `*`(J_))), `*`(mol_)), `-`(`/`(`*`(39460.00, `*`(x[OH], `*`(J_))), `*`(mol_)))))), `*`(`+`(`/`(`*`(34, `*`(x[H2], `*`(J_))), `*`(mol_, `*`... (5)

Es más fácil resolver en función de la T las 6 primeras, y comprobar cuando se cumple la séptima

> N:=5:for ii from 1 to N do Ti:=evalf(3000*K_+1000*K_*(ii/N));sol_:=solve(evalf(subs(T=Ti,{eqNX,eqBH,eqBO,eqEQ0_,eqEQ1_})),{a,x[H2O],x[H2],x[O2],x[OH]})[1];subs(sol_,cpComp,T25=T1,dat,eq15_7_3);print('Ti'=Ti,' Ti-Ta'=Ti-rhs(%));od:Ti:=3765*K_;sol_:=solve(evalf(subs(T=Ti,{eqNX,eqBH,eqBO,eqEQ0_,eqEQ1_})),{a,x[H2O],x[H2],x[O2],x[OH]})[1];PCI__=subs(sol_,PCI_);subs(sol_,cpComp,T25=T1,dat,eq15_7_3);

Ti = `+`(`*`(3200., `*`(K_))), `+`(Ti, `-`(Ta)) = `+`(`-`(`*`(1525.357663, `*`(K_))))
Ti = `+`(`*`(3400., `*`(K_))), `+`(Ti, `-`(Ta)) = `+`(`-`(`*`(1030.218136, `*`(K_))))
Ti = `+`(`*`(3600., `*`(K_))), `+`(Ti, `-`(Ta)) = `+`(`-`(`*`(483.196206, `*`(K_))))
Ti = `+`(`*`(3800., `*`(K_))), `+`(Ti, `-`(Ta)) = `+`(`*`(102.121803, `*`(K_)))
Ti = `+`(`*`(4000., `*`(K_))), `+`(Ti, `-`(Ta)) = `+`(`*`(708.055687, `*`(K_)))
`+`(`*`(3765, `*`(K_)))
{a = .9044535730, x[H2] = .1236134881, x[H2O] = .7133607190, x[O2] = 0.2806706114e-1, x[OH] = .1349587317}
PCI__ = `+`(`/`(`*`(167179.4175, `*`(J_)), `*`(mol_)))
Ta = `+`(`*`(3767.329115, `*`(K_))) (6)

i.e.  la solución es Ta=3765 K, x[H2O]=71,3 %, x[H2]=12,4 %, x[O2]=2,8 %, x[OH]=13,5 %, y PCI=167 kJ/mol (en vez de los 242 kJ/mol sin disociación). Experimentalmente resulta Ta=3500 K. La discrepancia con nuerstro modelo es que hemos usado valores medios de cp típicos de combustión de 300 K a 2000 K, que son más bajos que los de 300 K a 3250 K.

f) Exceso de hidrógeno que haría que la temperatura de combustión adiabática sin disociación fuese 2500 K.

Sea s el exceso de H2.

> eq4:=s*H2+H2+(1/2)*O2=H2O+s*H2;PCI_:=PCImol_;eqTa:=T2=T1+PCI/(1*c[p,H2O]+s*c[p,H2]);eqTa_:=subs(cpComp,PCI=PCImol_,dat,SI0,%);s_:=solve(%,s);r[OF]:=mO/mF=((1/2)/(1+s))*subs(Odat,M)/subs(Fdat,M);subs(s=s_,%);

`+`(`*`(s, `*`(H2)), H2, `*`(`/`(1, 2), `*`(O2))) = `+`(`*`(H2, `*`(s)), H2O)
`+`(`/`(`*`(241820.00, `*`(J_)), `*`(mol_)))
T2 = `+`(T1, `/`(`*`(PCI), `*`(`+`(`*`(s, `*`(c[p, H2])), c[p, H2O]))))
2500 = `+`(273.15, `/`(`*`(241820.00), `*`(`+`(`*`(34, `*`(s)), 47))))
1.811554570
`/`(`*`(mO), `*`(mF)) = `+`(`/`(`*`(8.000000000), `*`(`+`(1, s))))
`/`(`*`(mO), `*`(mF)) = 2.845400934 (7)

i.e. para que la Ta sea 2500 K han de entrar s+1=2,8 moles de H2 en vez de 1 mol de H2; la relación másica oxidante/combustible será por tanto rOF=2,8 en vez de rOF=8 para estequiometría. En los cohetes LH2/LOX se suele operar con rOF=6 porque mejora el impulso específico en vacío y rebaja la Ta en la cámara de combustión desde los 3500 K estequiométricos a unos 3250 K.

>