I am trying to compute the final composition and temperature of a mixture of syn-gas and air. They enter in at 300 K and 600 K respectively. The syn-gas is a mixture of CO and H2 the proportions of which I am varying from 3:1 to 1:3. Later on, this ratio is fixed and additional nitrogen gas is introduced. Lastly, heat loss is accounted for and its effect on temperature/composition is calculated. Focusing on the first part, I am having a hard time balancing the system of non-linear equations. The general equation for the chemical reaction is as follows:
aCO + bH2 + c*(O2 + 79/21 * N2) + dN2 = eCO + fH2 + gO2 + hN2 + jCO2 + kH2O + lNO
From conservation of species:
Carbon: a = e + j
Oxygen: a + 2c = e + 2g + 2*j + k + l
Hydrogen: 2b = 2f + 2*k
Nitrogen: 2*(79/21)c + 2d = 2*h + l
Since there are three compounds, there are three partial pressure equilibrium values called K_p. K_p is a function of temperature and from empirical data is a known constant.
K_p_NO = (X_NO * P) / (sqrt(X_N2*P)*sqrt(X_O2 * P))
K_p_H2O = (X_H2O * P) / (sqrt(P*X_O2)X_H2P)
K_p_CO2 = (X_CO2 * P) / (sqrt(P*X_O2)X_COP)
Where X is the mole fraction. E.G. X_H2O = k/(e+f+g+h+j+k+l)
The variables e,f,g,h,j,k,l are 7 unknowns and there are 7 equations. Variables a, b, c, and d are varied manually and are treated as known values. Using scipy, I implemented fsolve() as shown below:
from scipy.optimize import fsolve # required library
import sympy as sp
import scipy
# known values hard coded for testing
a = 0.25
b = 0.75
c = a + b
d = 0
kp_NO = 0.00051621
kp_H2O = 0.0000000127
kp_CO2 = 0.00000001733
p = 5 # pressure
dec = 10 # decimal point precision
# Solving the system of equations
def equations(vars):
(e, f, g, h, j, k, l) = vars
f1 = e + j - a
f2 = e + 2*g + 2*j + k + l - a - 2*c
f3 = f + k - b
f4 = 2*h + l - 2*d - (2*79/21)*c
f5 = kp_NO - (l/sp.sqrt(c*(79/21)*c))
f6 = kp_H2O - (k/(b*sp.sqrt((c*p)/(e + f + g + h + j + k + l))))
f7 = kp_CO2 - (j/(a*sp.sqrt((c*p)/(e + f + g + h + j + k + l))))
return[f1, f2, f3, f4, f5, f6, f7]
e, f, g, h, j, k, l = scipy.optimize.fsolve(equations, (0.00004, 0.00004, 0.49, 3.76, 0.25, 0.75, 0.01))
# CO, H2, O2, N2, CO2, H2O, NO
print(e, f, g, h, j, k, l)
The results are
0.2499999959640893 0.7499999911270915 0.999499382628763 3.761404150987935 4.0359107126181326e-09 8.872908576472292e-09 0.001001221833654118
Notice that e = 0.24999995 but a = 0.25. It seems that the chemical reaction is progressing little if at all. Why am I getting my inputs back as results? I know something is wrong because in some cases, the chemical coefficients are negative.
Things I've tried: Triple checked my math/definitions. Used nsolve() from sympy, nonlinsolve() from scipy, other misc. solvers.
I do not see anything wrong in your code, albeit I would have solved it differently. This question should be moved to Chemistry stackexchange: the Kp values you are using are wrong. The value for water formation at 500K is around 7.65E22 (when pressure is expresed in bar), I am pretty sure that the constant for CO to CO2 is also much higher.
I would have written this as a comment but I do not have enough reputation.