crc_NH3_gibbs := [[298.15, -16.407 * 10^3], [300, -16.223 * 10^3], [400, -5.980 * 10^3], [500, 4.764 * 10^3], [600, 15.846 * 10^3], [700, 27.161 * 10^3], [800, 38.639 * 10^3], [900, 50.231 * 10^3], [1000, 61.903 * 10^3], [1100, 73.629 * 10^3], [1200, 85.392 * 10^3], [1300, 97.177 * 10^3], [1400, 108.975 * 10^3], [1500, 120.779 * 10^3]]:

plots:-display( plot( crc_NH3_gibbs, style = point ), plot( Property( Gmolar, "NH3(g)", "temperature" = T ), T = 298.15 .. 1500 ) )

g_NO  := T -> Property( "Gmolar", "NO(g)", "temperature" = T ):
g_O2  := T -> Property( "Gmolar", "O2(g)", "temperature" = T ):
g_NO2 := T -> Property( "Gmolar", "NO2(g)", "temperature" = T ):

Gibbs_rn := T ->  2 * g_NO2(T) - 2 * g_NO(T) - g_O2(T):
plot( Gibbs_rn, 300 * Unit(K) .. 1000 * Unit(K), labels = [ "Temperature (K)", "Gibbs Energy (J/mol)" ], labeldirections = [ horizontal, vertical ] )

fsolve( Gibbs_rn( T ) = 0, T = 1000 * Unit( K ) )

$765.77⟦K⟧$

The equilibrium composition is found by minimizing the Gibbs energy of the combustion products (formulated as a series of equations constructed via the method of Lagrange multipliers).

The adiabatic flame temperature is found by equating the enthalpy of the reactants and the enthalpy of the products.

h_CO  := Property("Hmolar", "CO(g)",  "temperature" = T ):
h_CO2 := Property("Hmolar", "CO2(g)", "temperature" = T ):
h_O   := Property("Hmolar", "O(g)",   "temperature" = T ):
h_O2  := Property("Hmolar", "O2(g)",  "temperature" = T ):

g_CO  := Property("Gmolar", "CO(g)",  "temperature" = T ):
g_O2  := 0:
g_CO2 := Property("Gmolar", "CO2(g)", "temperature" = T ):
g_O   := Property("Gmolar", "O(g)",   "temperature" = T ):

h_r_CO2 := Property( Hmolar, "CO2(g)", temperature = 298.15 );
h_r_CO  := Property( Hmolar, "CO(g)",  temperature = 298.15 );
h_r_O   := Property( Hmolar, "O(g)",   temperature = 298.15 );
h_r_O2  := Property( Hmolar, "O2(g)",  temperature = 298.15 );

$\mathrm{h\_r\_CO2}≔\mathrm{-393510.0001}$

$\mathrm{h\_r\_CO}≔\mathrm{-110535.1957}$

$\mathrm{h\_r\_O}≔249175.0027$

$\mathrm{h\_r\_O2}≔0.$

h_f_CO  := Property( "HeatOfFormation", "CO(g)"  );
h_f_O2  := Property( "HeatOfFormation", "O2(g)"  );
h_f_CO2 := Property( "HeatOfFormation", "CO2(g)" );
h_f_O   := Property( "HeatOfFormation", "O(g)"  )

$\mathrm{h\_f\_CO}≔\mathrm{-110535.196}$

$\mathrm{h\_f\_O2}≔0.$

$\mathrm{h\_f\_CO2}≔\mathrm{-393510.000}$

$\mathrm{h\_f\_O}≔249175.003$

R := 8.314:

con_1 := 2 * n_1 + n_2 + 2 * n_3 + n_4 = 4:

con_2 := n_1 + n_2 = 2:

n_t := n_1 + n_2 + n_3 + n_4:

gibbs := n_1 * ( g_CO2 + R * T * ln( n_1 / n_t) ) + n_2 * ( g_CO + R * T * ln( n_2 / n_t ) ) + n_3 * ( g_O2 + R * T * ln( n_3 / n_t ) ) + n_4 * ( g_O + R * T * ln( n_4 / n_t ) ):

eqComposition_1 := L_1 * diff( lhs( con_1 ), n_1 ) + L_2 * diff( lhs( con_2 ), n_1 ) = diff( gibbs, n_1 ):
eqComposition_2 := L_1 * diff( lhs( con_1 ), n_2 ) + L_2 * diff( lhs( con_2 ), n_2 ) = diff( gibbs, n_2 ):
eqComposition_3 := L_1 * diff( lhs( con_1 ), n_3 ) + L_2 * diff( lhs( con_2 ), n_3 ) = diff( gibbs, n_3 ):
eqComposition_4 := L_1 * diff( lhs( con_1 ), n_4 ) + L_2 * diff( lhs( con_2 ), n_4 ) = diff( gibbs, n_4 ):

flameTemp := 2 * h_f_CO + h_f_O2 = n_1 * h_CO2 + n_2 * h_CO + n_3 * h_O2  + n_4 * h_O :

sol := fsolve( { eqComposition_1, eqComposition_2, eqComposition_3, eqComposition_4, flameTemp, con_1, con_2 }, { L_1 = -10000, L_2 = -10000, T = 3000, n_1 = 0.1, n_2 = 0.1, n_3 = 0.1, n_4 = 0.1 } )

$\mathrm{sol}≔\left\{\mathrm{L\_1}=\mathrm{-23096.57164}\,\mathrm{L\_2}=\mathrm{-368734.5094}\,T=2975.362560\,\mathrm{n\_1}=1.130266439\,\mathrm{n\_2}=0.8697335610\,\mathrm{n\_3}=0.3841072729\,\mathrm{n\_4}=0.1015190152\right\}$

eval(T, sol)

$2975.362560$

eval( seq( n_||i / ( n_1 + n_2 + n_3 + n_4 ), i = 1 .. 4 ), sol )

$0.4547209870,0.3499051990,0.1545313850,0.04084242900$