Solution Procedure

The preceding sections utilized equilibrium constants and the first law of thermodynamics to construct a set of five simultaneous non-linear equations to model combustion of hydrogen and oxygen. The full set of equations are reproduced here. For convenience, these functions are labeled $f_1$ through $f_5$:^[37],[39]

This system can be solved using a Newton-Raphson method analogous to the scheme described in Section 2.2.^[15] To initiate the algorithm, an initial guess must be provided for the value of each of the unknown variables. These initial guesses will be stored in the vector $A_i$:

$$ A_i = \left[\begin{matrix} a_i \\ b_i\\c_i\\d_i\\T_i\end{matrix}\right]$$

Next, Equations $f_1$ through $f_5$ must be differentiated with respect to the five unknown variables of interest: a,b,c,d, and T. Therefore there are $5 \times 5= 25 $ partial derivatives, $F^\prime$ to calculate:

$$ F^\prime = \left[\begin{matrix} \frac{\partial f_1}{\partial a} & \frac{\partial f_1}{\partial b} & \frac{\partial f_1}{\partial c} & \frac{\partial f_1}{\partial d} & \frac{\partial f_1}{\partial T} \\ \frac{\partial f_2}{\partial a} & \frac{\partial f_2}{\partial b} & \frac{\partial f_2}{\partial c} & \frac{\partial f_2}{\partial d} & \frac{\partial f_2}{\partial T} \\\frac{\partial f_3}{\partial a} & \frac{\partial f_3}{\partial b} & \frac{\partial f_3}{\partial c} & \frac{\partial f_3}{\partial d} & \frac{\partial f_3}{\partial T} \\\frac{\partial f_4}{\partial a} & \frac{\partial f_4}{\partial b} & \frac{\partial f_4}{\partial c} & \frac{\partial f_4}{\partial d} & \frac{\partial f_4}{\partial T} \\\frac{\partial f_5}{\partial a} & \frac{\partial f_5}{\partial b} & \frac{\partial f_5}{\partial c} & \frac{\partial f_5}{\partial d} & \frac{\partial f_5}{\partial T}\end{matrix}\right]$$

These derivatives are evaluated numerically using a centered finite difference scheme and the values of $A_i$ guessed above. Here is an example:

Any convenient value can be used for $\Delta a$, provided it is sufficiently small. The smaller $\Delta a$ becomes, the more accurate the derivative will be. Theoretically, as $\Delta a$ approaches zero, the derivative will contain no error.

Next, each of the five equations is evaluated using the guess $A_i$ as an input. The results are stored in the vector $F$:

$$ F = -1\times\left[\begin{matrix} f_1\left(A_i\right) \\ f_2\left(A_i\right)\\f_3\left(A_i\right)\\f_4\left(A_i\right)\\f_5\left(A_i\right)\end{matrix}\right]$$

This reduces the original system of non-linear equations to a linear system:

The matrix $B$ in the above equation is used to store intermediate values that will help reach the final answer. Equation 3.3.1 can be solved for $B$ using any standard technique for systems of linear equations like Gauss elimination with partial pivoting or LU decomposition.^[15] Once $B$ is determined, it can be used to refine the initial guess for the solution of the system:

$$ A_{i+1}=A_{i}+B = \left[\begin{matrix} a_{i+1} \\ b_{i+1}\\c_{i+1}\\d_{i+1}\\T_{i+1}\end{matrix}\right]$$

After solving for the refined solution vector, the entire process is repeated using the new values from $A_{i+1}$ as the initial guess. This procedure is iterated until the solution vector $A$ converges to within a desired tolerance. One possible convergence criteria could be $\text{abs}\left(T_{i+1}-T_{i}\right) < 1 $. This says that when the change in flame temperature between each iteration is less than 1 K, the algorithm has found a solution.

Once the solution converges, the values of $a$, $b$, $c$, and $d$ are plugged in to Equations 3.1.7 - 3.1.12 to solve for the mole fractions of each species in the exhaust mixture. We now have the composition of the combustion gas, along with its pressure and temperature. These results can plugged in the equations found in Sections 2.4 and 2.6 to solve for the relevant thermodynamic and transport properties.

Benchmark Case

In order to test the accuracy of the model, its output is compared against results from NASA’s Chemical Equilibrium with Applications (CEA) program.^{[43] [85]} A benchmark case is established with the following input conditions:

Table 3.3.1

As shown in the following plots, the flame temperature and equilibrium mole fractions converge to their final result within 4 iterations.

Figure 3.3.1 Converging Flame Temperature

Figure 3.3.2 Converging Mole Fractions

The output of this program is in excellent agreement with the results from NASA CEA. The thermodynamic properties have a relative error of less than 2%. A significant amount of error does exist on some of the predicted mole fractions (up to 97% error for $O$). This error is not particularly concerning though. When the chemical reaction reaches equilibrium at 3625 K, the population of $O$ molecules has been almost completely consumed. For all practical purposes, $Z_{O}=0.$ CEA predicts a mole fraction of 0.00205, while this model predicts 0.00407. The large amount of relative error between these near-zero mole fractions has little impact on the thermodynamic and transport properties. Relative error for species with significant mole fractions like $H_2O$ and $H_2$ is much better, 2.62% and 1.71% respectively. The following table contains the complete comparison.

Table 3.3.2: Computational Results

Previous Next

Top