Combustión:

Determinación de aire estequiométrico en la combustión de carbón:

Con la información del análisis último del carbón, determinamos las fórmula molecular del mismo

UA=[84.36 1.89 4.4 0.63 .89];

MW=[12 2 16*2 14*2 32];
n=zeros(1,length(UA));
y=zeros(1,length(UA));
for i=1:length(UA)
n(i)=UA(i)/MW(i);
end
nt=sum(n);
for i=1:length(UA)
y(i)=n(i)/nt;
end

Establecemos el sistema de ecuaciones para la reacción:

1
Ordenando el sistema de ecuaciones

A=[1 0 0 0 0;
0 1 0 0 0;
1 0.5 1 0 -1;
0 0 0 1 -3.76;
0 0 1 0 0]

A = 5×5
1 0 0 0 0
0 1 0 0 0
1 0.5 1 0 -1
0 0 0 1 -3.76
0 0 1 0 0

B=n'

B = 5×1
7.03
0.945
0.1375
0.0225
0.027813

sol=A\B

sol = 5×1
7.03
0.945
0.027813
27.819
7.3928

nt_prod=sol(1)+sol(2)+sol(3)+sol(4)

nt_prod =
35.822

for i=1:4
y_prod(i)=sol(i)/nt_prod;
end

2
 y_prod

y_prod = 1×4
0.19625 0.02638 0.0007764 0.7766

Las moles de oxígeno que ingresan al sistema son:

a_th=sol(5)

a_th =
7.3928

Finalmente, la relación AC (masa de aire/masa combustible) tomando como base para el cálculo 100 kg de
carbón es

Una expresión equivalente es:

AC=a_th*4.76*29/100

AC =
10.205

AC=(a_th/0.21)*29/100

AC =
10.209

Combustión de aire húmedo

3
Asumiremos que el contenido de humedad en el aire de combustión no reacciona, por tanto, analizaremos el
sistema como aire seco y al final, sumaremos tanto a productos como a reactivos la humedad. Usando una
base de cálculo de 100 kmol, la reacción de combustión para el aire seco será:

A=[1 0 0 0;
0 1 0 0;
0 0 1 -3.76;
1 0.5 0 -1]

A = 4×4
1 0 0 0
0 1 0 0
0 0 1 -3.76
1 0.5 0 -1

B=[0.72+0.03 2*0.72+0.09 0.14 0.02+0.03]'

B = 4×1
0.75
1.53
0.14
0.05

sol2=A\B

sol2 = 4×1
0.75
1.53
5.6484
1.465

Ahora a través de CoolProp, estimaremos la humedad específica del aire (kg water/kg dry air)

W=py.CoolProp.CoolProp.HAPropsSI('W','T',273.15+20,'P',101325,'R',0.8)

W =
0.011753

4
También sabemos que hay 1.465*(4.76) kmol de aire seco, con esta información determinaremos las moles de
agua a la entrada.

kg_dry_air=sol2(4)*29/0.21; %29 kg es el MW del aire

kmol_water=W*kg_dry_air/18

kmol_water =
0.13209

En los productos de combustión tenemos de moles de gas producido. Expresando esto en

base seca (kg water/kg dry air) tenemos:

W_prod=((sol2(2)+kmol_water)*18)/((sol2(1)*44)+(sol2(3)*28))

W_prod =
0.15651

Finalmente, la temperatura de condensación será:

T_dew_point=py.CoolProp.CoolProp.HAPropsSI('Tdp','T',273.15+20,'P',101325,'W',W_prod
)-273.15

T_dew_point =
60.331

Aplicación de la ley de Hess y las entalpías de formación

Solución

5
Como se ha discutido, la reacción directa no es posible en la práctica. Adicionalmente, se nos proponen 3
reacciones que podemos utilizar. Una forma de encontrar la solución es a través de matrices.

format shortG
A=[-1 0 0;
0 -1 0;
0 0 -1;
-7/2 -1 -1/2;
2 1 0;
3 0 1];
B=[1 -2 -3 0 0 0]';
sol=mldivide(A,B)';
deltaH_Rx=[-1560; -393.5; -286]

deltaH_Rx = 3×1
-1560
-393.5
-286

deltaH_form=sol*deltaH_Rx

deltaH_form =
-85

Solución
Balance de la reacción de combustión con exceso de aire
En primer lugar vamos a establecer la reacción de combustión completa

6
 clear all
format shortG
A=[1 0 0 0;
0 1 0 0;
1 0.5 0 -1;
0 0 1 -3.76];
B=[3 4 0 0]';
coeff=A\B

coeff = 4×1
3
4
18.8
5

a_th=coeff(4);

El proceso de combustión se lleva a cabo con un exceso de aire del 50%, por tanto, el

a_th_real=a_th*1.5

a_th_real =
7.5

Asimismo, el análisis de los gases muestra que todo el se convierte en agua pero solo el 90% del C se
convierte en , dejando que el 10% restante se convierta en CO. Teniendo en mente estas relaciones,
procedemos a balancear la ecuación de combustión real.

A=[1 0 1 0 0;
0 1 0 0 0;
1 0.5 0.5 0 1;
0 0 0 1 0;
1/0.9 0 -1/0.1 0 0];
B=[3 4 a_th_real 3.76*a_th_real 0]';
coeff_real=A\B

coeff_real = 5×1
2.7

7
 4
0.3
28.2
2.65

n_comb_real=sum(coeff_real)

n_comb_real =
37.85

for i=1:length(coeff_real)
y_comb(i)=coeff_real(i)/n_comb_real;
end
y_comb

y_comb = 1×5
0.071334 0.10568 0.007926 0.74505 0.070013

Pasamos a calcular la relación AC recordando que:

AC=((a_th_real/0.21)*29)/(3*12+8)

AC =
23.539

El flujo de combustible es de por tanto el flujo de de aire es:

m_dot_air=AC*0.05/60 %Se hace la conversión para obtener el flujo en kg/s

m_dot_air =
0.019616

Determinación de la transferencia de calor considerando que los gases de combustión

abandonan la cámara de combustión a una temperatura de 1500K
El libro de Cengel nos propone determinar

Y nos entrega la siguiente tabla (Que toma los valores de los apéndices A-20 en adelante):

8
Para luego continuar con la determinación del calor de combustión según:

Este método es iterativo y no es tan exacto como la evaluación de la entalpía a partir de relaciones de Cp. A
continuación, explicaremos cómo abordar este problema usando propiedades termodinámicas más generales
en MATLAB. También es cierto que podríamos evaluar las entalpías de las sustancias a determinadas
temperaturas usando, por ejemplo coolprop. Este método podría llegar a ser incluso más exacto pero igual
sería iterativo cuando quisiéramos determinar condiciones como la temperatura de llama adiabática.

Evaluación de la entalpía de una especie en función de la temperatura:

Podemos evaluar la entalpía con la siguiente expresión para la especie i en la reacción:

Así las cosas, para determinar la entalpía de reacción a condiciones estándar tenemos

Y en nuestro caso concreto, considerando la reacción balanceada:

9
Expresando la ecuación para la entalpía de formación estándar de la reacción y considerando que los
coeficientes estequiométricos son negativos para los reactivos y positivos para los productos tenemos:

La ley de Hess nos indica que las entalpías de formación de un elemento en su estado de referencia estable es
cero, esto reduce la expresión a:

Esta expresión resulta fácilmente evaluable en MATLAB a través de un ciclo "for".

nu=[-1 -a_th_real -3.76*a_th_real coeff_real'];

Hf_i=[-118910 0 0 -393509 -241818 -110525 0 0];
for i=1:length(nu)
nu_Hf_i(i)=nu(i)*Hf_i(i);
end
DeltaHf_Rx=sum(nu_Hf_i)

DeltaHf_Rx =
-1.944e+06

Esta entalpía de formación que acabamos de determinar es equivalente al LHV del combustible (propano). En
un escenario en donde los productos de combustión salen a las mismas condiciones de temperatura que los
reactivos, liberando toda la energía a los alrededores. En nuestro caso, parte de esa energía no es liberada y
termina aumentando la temperatura de los gases de combustión en un proceso más cercano a la realidad.

La estimación de la entalpía de las especies individuales así como la entalpía de reacción a una determinada
temperatura se puede generalizar de la siguiente manera:

Simplificando esta expresión obtenemos:

En donde, se calcula como:

Para que esta notación tenga sentido, debemos entender que los coeficientes estequiométricos son negativos
para los reactivos y positivos para los productos. Por esta razón cuando definimos la matriz "nu" asignamos

10
el signo negativo a los reactivos. Esta convención de signos también nos indicará que si obtenemos un
negativo significa que la reacción exotérmica y por tanto estará liberando energía al ambiente.

Para este cálculo seguiremos el siguiente ejemplo:

Creación de matriz de coeficientes para el .

CP_matrix=[31.74e-9 -15.72e-5 30.48e-2 -4.04;

1.312e-9 -0.7155e-5 1.520e-2 25.48;
-2.873e-9 0.8081e-5 -0.1571e-2 28.9;
7.469e-9 -3.501e-5 5.981e-2 22.26;
-3.595e-9 1.055e-5 0.1923e-2 32.24;
-2.222e-9 0.5372e-5 0.1675e-2 28.16;
-2.873e-9 0.8081e-5 -0.1571e-2 28.9;
1.312e-9 -0.7155e-5 1.520e-2 25.48];
[m,n]=size(CP_matrix);
for i=1:m
for j=1:n
DeltaCP(i,j)=nu(i)*CP_matrix(i,j);
end
end
DeltaCP

DeltaCP = 8×4
-3.174e-08 0.0001572 -0.3048 4.04
-9.84e-09 5.3663e-05 -0.114 -191.1
8.1019e-08 -0.00022788 0.044302 -814.98
2.0166e-08 -9.4527e-05 0.16149 60.102
-1.438e-08 4.22e-05 0.007692 128.96
-6.666e-10 1.6116e-06 0.0005025 8.448
-8.1019e-08 0.00022788 -0.044302 814.98
3.4768e-09 -1.8961e-05 0.04028 67.522

DeltaCP_reac_i=DeltaCP(1:3,:)

DeltaCP_reac_i = 3×4
-3.174e-08 0.0001572 -0.3048 4.04
-9.84e-09 5.3663e-05 -0.114 -191.1
8.1019e-08 -0.00022788 0.044302 -814.98

DeltaCP_prod_i=DeltaCP(4:8,:)

DeltaCP_prod_i = 5×4
2.0166e-08 -9.4527e-05 0.16149 60.102
-1.438e-08 4.22e-05 0.007692 128.96
-6.666e-10 1.6116e-06 0.0005025 8.448
-8.1019e-08 0.00022788 -0.044302 814.98
3.4768e-09 -1.8961e-05 0.04028 67.522

for i=1:n
DeltaCP_reac(i)=sum(DeltaCP_reac_i(:,i));
end
for i=1:n
DeltaCP_prod(i)=sum(DeltaCP_prod_i(:,i));
end

11
 DeltaCP_reac

DeltaCP_reac = 1×4
3.9439e-08 -1.7022e-05 -0.3745 -1002

DeltaCP_prod

DeltaCP_prod = 1×4
-7.2422e-08 0.00015821 0.16566 1080

Tref=25+273.15;
Treac=7+273.15;
Tprod=1500;
del_H_reac=polyint(DeltaCP_reac)

del_H_reac = 1×5
9.8596e-09 -5.6739e-06 -0.18725 -1002 0

del_H_prod=polyint(DeltaCP_prod)

del_H_prod = 1×5
-1.8106e-08 5.2736e-05 0.08283 1080 0

H_reac=polyval(del_H_reac,Treac)-polyval(del_H_reac,Tref)

H_reac =
19994

H_reac=diff(polyval(del_H_reac,[Tref Treac]))

H_reac =
19994

H_prod=polyval(del_H_prod,Tprod)-polyval(del_H_prod,Tref)

H_prod =
1.5621e+06

H_prod=diff(polyval(del_H_prod,[Tref Tprod]))

H_prod =
1.5621e+06

Finalmente, la entalpía de reacción, considerando la temperatura de entrada de productos y de reactivos es:

Del_H_Rx=DeltaHf_Rx+H_prod+H_reac

Del_H_Rx =
-3.6191e+05

DeltaHf_Rx

DeltaHf_Rx =
-1.944e+06

Como comentario final, es importante notar que es mayor en magnitud que ya que parte de la
energía que debería ser liberada a los alrededores se consume elevando la temperatura de los gases de
combustión.

12
Temperatura de llama adiabática
Como se mencionó en la clase. Si la combustión es adiabática no hay liberación de calor hacia los alrededores
y por tanto . Esta condición hará que los gases de combustión alcancen la máxima temperatura
posible a las condiciones de la reacción

Del ejercicio anterior hemos determinado todas las magnitudes pero será entonces la temperatura de
llama adiabática que desamos encontrar para que se satisfaga que:

Expresemos el valor del lado izquierdo de la ecuación:

DeltaHf_Rx+H_reac

ans =
-1.924e+06

fun=(DeltaHf_Rx+H_reac)+polyval(del_H_prod,1500)

fun =
-31290

pol_prod=[del_H_prod(:,1:4) DeltaHf_Rx+H_reac]

pol_prod = 1×5
-1.8106e-08 5.2736e-05 0.08283 1080 -1.924e+06

Tsol=roots(pol_prod)

Tsol = 4×1 complex

5227.1 + 0i
-1918 + 3111.4i
-1918 - 3111.4i
1521.7 + 0i

fun=(DeltaHf_Rx+H_reac)+polyval(del_H_prod,Tsol(4,1))

fun =
6.9849e-10

13