0

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.

Thomas Murdoch
  • 13
  • 2
  • Also, the guesses in the fsolve() function are real results from another sofware (gasEq). Variables e, f, g, h, j k, l should be near those values. – Thomas Murdoch Dec 04 '21 at 22:00
  • You can plug the values return by `fsolve` back into `equations()` to check if the result is approximately zero. I just tried that, and the magnitudes of the components returned by `equations()` range from 0 to 4e-14. So, numerically, it is a solution. – Warren Weckesser Dec 05 '21 at 01:19
  • Is there a way to look for more than 1 solution then? – Thomas Murdoch Dec 05 '21 at 05:32
  • @ThomasMurdoch Do you know what values you are expecting to see? If this is a valid solution (like Warren said) then, if you're expecting something else, it's likely that your equation has multiple roots, and so to find a different root you would need to pass in different initial conditions (as I believe fsolve() is trying to do a linear gradient descent, so will only ever find the closest local minimima/zero) – Richard Dec 05 '21 at 19:37
  • This equation has only one root. The reaction almost does not advance because the constants are very small. The correct values (at least for H2O and CO2) are many orders higher. – azelcer Dec 05 '21 at 19:48

1 Answers1

1

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.

ncica
  • 7,015
  • 1
  • 15
  • 37
azelcer
  • 1,383
  • 1
  • 3
  • 7
  • Good observation! ....@Thomas perhaps check your parameters, as that might be why you're not seeing what you think you should – Richard Dec 05 '21 at 19:38
  • Thanks for the comment. The hard coding of values made it incorrect. Also, some of my equations are not right and I fixed them after posting this. Running into another problem but I think it is more of a chemistry question. Thanks for the help – Thomas Murdoch Dec 05 '21 at 21:46