System Dynamic modelling approach for resolving the thermo-chemistry of wood gasification

For Multiphysics problems that require a thorough understanding of multiple, influential, highly transient process parameters, a System Dynamic model can constitute either an alternative option, or a compact prelude to a more expensive 3-D Finite Element or Finite Volume model. As a rather uncommon example for the application of such a modelling method, this work presents a System Dynamic modelling concept, devised for resolving the thermo-chemistry within a wood gasification reactor. It compares the modelling concept as well as its results to a classic, thermo-chemical solution algorithm based on the minimization of LaGrangian Multipliers for resolving the gasification equilibrium equations. In contrast to the latter, the System Dynamic solver can consider the impact of reaction kinetics as well as molecular mass transfer effects on the gasification equilibrium. Thus the transient production rates of methane, hydrogen, carbon (di-) oxide and water, as well as the residual amounts of pyrolysis gas and oxygen, which occur during the gasification of a wood particle, can be predicted.


INTRODUCTION
To be able to handle the highly unstable wood gasification process on a large industrial scale, either a relatively high amount of maintenance (see e.g.[1]), or a thorough understanding of its thermodynamic basics is required.In order to contribute to achieving the latter, an extensive, single particle, thermo-fluid dynamic gasification model has been developed and presented in [2].
This work focuses on advances in connection with the core part of that previous full particle gasification model: the thermo-chemical wood gas equilibrium solver.
Hereby two alternative more advanced methods to calculate wood gas compositions are investigated: 1.A Gibbs free energy minimization approach based on the method of Lagrangian Multipliers (e.g.: [4], [6], [7], [13]) has been implemented within Matlab [21].The solver has been adapted for the case of wood gas composition calculation and is laid out in chapter 2 of this work.The main advantage of this scheme, as compared to the equilibrium constant solver [2], beside its efficiency, is that only the relevant, chemical species, but not the chemical reactions need to be known in advance [4] in order to resolve the equilibrium.2. A new procedure, based on the principles of System Dynamics [22], the Karlsruhe Physics course ( [5], [18][19][20]) as well as bond graph theory [23], hereby denoted as the "System Dynamics Solver" has been created within Berkeley Madonna [32] and is laid out in chapter 3. The main advantage of this scheme, besides its ability to visualize decisive dependencies, is that kinetic mechanisms as well as thermodynamic aspects are included in the calculations.Thus the impact of reaction-and transport kinetics on chemical equilibrium formation can be investigated much more thoroughly than with comparable solvers.In contrast to the Gibbs free energy minimization approach though, the involved chemical reactions need to be known in advance.

THE WOOD GASIFICATION SYSTEM AND RELATED REACTIONS
The seven main chemical species i, to be considered in the context of thermo-chemical wood gas composition calculations are hydrogen, water, carbon (di)-oxide, methane, oxygen and pyrolysis gas.At temperatures ranging from 300K to 1500K, the wood gas compounds are in intense thermo-chemical interaction.According to [1] and [6], the predominant gasification reactions can be summarized as follows.
Boudouard reaction: Those mechanisms occur in combination with a variety of possible oxidation reactions as seen in Table 1.
In Eqns.1-3 and Table 1, the "(g/s)" notation hints to the fact that those reactions can occur homo-as well as heterogeneously [2].
Undergoing the reactions, cited above and being fed with oxygen from process air as well as pyrolysis gas from cellulose gasification, the seven relevant chemical wood gas species constitute a highly dynamic system of formidable complexity.The wood gas system is sketched in Fig. 1.

THE LAGRANGIAN MULTIPLIER MODEL FOR WOOD GAS COMPOSITION
The method of Lagrangian Multipliers has been used to implement a Gibbs free energy minimization approach according [4], [6], [7] and [13], for wood gas composition equilibrium calculation.This type of solver is based upon two main ideas: the minimization of over all Gibbs free energy within the entire system and the minimization of divergence concerning its atomic species balance.
In the following, the set up behind the Matlab based Lagrangian Multiplier solver for a system of i molecular species and j atomic species, is laid out.In the case of wood 138 System Dynamic modelling approach for resolving the thermo-chemistry of wood gasification The function consists of two components: one expresses the sum of all Gibbs free energies of formation ΔG i of all molecular species i, while the other stands for the j atomic species balances.Thereby λ j and b 0 j are the Lagrangian Multiplier and the input rate of atomic species j respectively, a i, j is the number of atoms j per molecule i and N i is the total number of molecules within the system.

Step 2 -Minimization of the LaGrange function
A minimization of the total LaGrange function yields minimum Gibbs free energy of the system, as well as the best possible solution for a closed atomic balance.Thus L is derived with respect to the two types of variables (N i an λ j ) to be sought as seen in Eqn. 5. (5)

Step 3 -Setting up the minimization criteria equation system
In combination with the introduction of the molecular species balance and the total number of molecules N tot , Eqn. 5 can be used to derive three main types of minimization criteria equations, Eqns.6 to 8. Eqns.6 to 8 constitute a system of j + i + 1 equations for i unknown numbers of molecules, j unknown Lagrangian Multipliers and the unknown total number of molecules.
Inserting for the case of wood gas equilibrium, a system of 11 (7 + 3 + 1) equations emerges.By stating that π j =-λ j /RT and by slightly rewriting the type I criterion from Eqn. 6, the full system of minimization criteria equations for wood gasification can be written as seen in Eqn. 9.

Step 4 -Solving the minimization criteria equation system by Newton Raphson
The minimization criteria equation system, shown in Eqn.6 to 9 is non-linear in nature and cannot be solved algebraically.Thus the solution has to be achieved numerically.Following the recommendation of e.g.[4], the multi-dimensional Newton Raphson [24] scheme is applied to linearize the system by a first order approximation.The basic idea behind the well known Newton Raphson scheme is shown in Eqn.10. (10) In Eqn. 10 x k is the solution vector of system f and iteration k, also shown in Eqn.11, s k is the deviation vector between x k and x k + 1 , also seen in Eqn. 12 and ∇ xk is the Jacobian of system f with respect to x k , also seen in Eqn. 13. (11) Combining Eqn. 6 to 8 with Eqn. 10 to 13, provides the full system of j + i + 1, linearized equations, as seen in Eqn.14. Hereby the first, second and third row represent a series of i, j and 1 equations respectively.

Step 5 -Application to wood gasification
At this point, all that remains to be done is to insert into Eqn.14 for the special case of wood gasification equilibrium calculation.This then yields the full, linearized equation system of 11 equations and 11 unknowns, of which Eqn. 15 shows the left-and Eqn.16 the right hand side in matrix notation.
142 System Dynamic modelling approach for resolving the thermo-chemistry of wood gasification In the Matlab implementation, a reasonable starting point x 0 is chosen.Then a conjugate gradient solver [25] is used to solve the system seen in Eqn. 15 and Eqn.16, such that s 0 is yielded, according to Eqn. 12.The Newton Raphson scheme is then iterated using x k + 1 = x k + s k until a sufficiently converged solution for x k and thus for all N i , λ j and N tot is found.

Properties of the Lagrangian Multiplier solver
In contrast to the iterative, equilibrium constant based solver presented in [2] and the System Dynamic solver, laid out in chapter 3, the Lagrangian Multiplier solver does not require any apriori knowledge of chemical reactions.It is relatively fast, efficient and well suited for resolving the equilibrium state of any thermochemical system with relatively large numbers of molecular species Y i .Its memory requirements are proportional to Y i 2 and its speed in relation to Y i is bound by the properties of the hereby chosen conjugate gradient solver [25].
The yielded result corresponds to the thermodynamic equilibrium composition of the gas mixture, at a global minimum of total Gibbs free energy within the system.Kinetic aspects or a temporal resolution of reaction mechanisms cannot be considered though.This is a major disadvantage concerning the thorough research of gasification processes.However, the solver is very well suited to efficiently investigate thermodynamic wood gas equilibrium compositions.Fig 2 shows some exemplary results (wood gas composition x i ) produced by the Lagrangian Multiplier solver, as temperature (left) and the ratio of oxygen to carbon atoms R O/C (right) are varied along the x-axis.
In the following, the Lagrangian Multiplier solver is mainly used to cross-reference and verify the System Dynamic solver, presented in chapter 3.

THE SYSTEM DYNAMIC MODEL FOR WOOD GAS COMPOSITION
The System Dynamic model for wood gas composition determination is based upon ideas and concepts formulated e.g. by Callen [30], Burkhardt [19], Herrmann [5], Maurer [20], Fuchs [26], Job [27] and Falk [28] (in arbitrary order) as well as Gibbs fundamental equation [29] seen in Eqn. 17.It thus relates to the concepts described within The Karlsruhe Physics course [18] and System Physics [19,20] and was created in full awareness of related criticism [31].

SYSTEM DYNAMIC METHODS IN PHYSICS
The very basis of applying System Dynamic methodology within physics is Gibbs fundamental equation seen in Eqn. 17. (17) Eqn. 17 shows that the global change of energy dE within any system, stems from a change of its basic, conservative quantities, such as not exclusively entropy S, volume V, molecules N, mass m, charge Q, direction i momentum P i or direction i normal angular momentum L i as well as from the extent of their denoted, respective potentials temperature T, pressure p, chemical potential μ, gravitational potential ψ, electrical potential ϕ, direction i velocity v i or direction i normal angular velocity ω i .As a consequence, the relevant variables of any physical system are distinguished to be either general conservative quantities Ψ on the one hand or general potentials ϕ on the other.In addition to this, the general system capacity κ can be introduced in order to relate the filling content Ψ to the respective filling height ϕ as seen in Eqn.18. (18) The relation between this point of view and the general differential balance equation of quantity Ψ (Eqn.19), for any physical system, can be graphically interpreted in System Dynamic terms as seen in Fig. 3. (19) dE TdS pdV dN dQ v dP dL dm ... In Eqn.19 and Fig. 3, u is a flow velocity in and out of the system, Γ is a diffusion coefficient and s Ψ is a source term of quantity Ψ.
Table 2 sums up some examples of the cited concept, relating it to some, selected fields of physics.

SYSTEM DYNAMIC METHODS TO RESOLVE CHEMICAL EQUILIBRIA
Based on the concepts, described in chapter 3.1, a System Dynamic interprtation of reacting chemical systems can be achieved, as is shown in the following.
Inserting into Eqn.19 for the special case of a chemically reacting system, the general balance equation becomes Eqn.20. (20) Therby the conservative quantity is the volume specific number of molecules of species i aka.concentration c i .The system volume is expressed by V sys and the reaction rate, related to reaction j is I R, j .
If fluxes across system boundaries are disregarded, the molecular fluxes I Ni of species i within a chemically reactive system are caused solely by the chemical reactions.Thus Eqn.20 simplifies to Eqn. 21.
144 System Dynamic modelling approach for resolving the thermo-chemistry of wood gasification  With those simplifications, System Dynamic terminology can describe any thermochemical system on the basis of the principles stated within Table 3.
According to Gibbs fundamental equation (Eqn.17) the chemical potential μ can be identified to be the decisive potential variable within chemistry, besides the molar composition of a mixture x i .Note that, starting with Table 3, the chemical potential of species i is in the following denoted as its molar Gibbs free energy ΔG i .
The System Dynamic scheme, needed to resolve the dynamic behavior and equilibrium composition x i, End of any chemical system of i species and j reactions can now be described on the basis of a graphical template, depicted in Fig. 4.
In Fig. 4 the green rectangles symbolize the system content of the conservative quantities: Molecules N i and related Gibbs free energy ΔG i .The full arrows denote fluxes of the conservative quantities, while dashed arrows stand for information flow within the model.The potentials are x i and ΔG i and the capacities are the total number of molecules N tot as well as the number of molecules of each species N i .Table 4 sums up all necessary relations to correctly interpret Fig. 4.
In Table 4, Θ stands for the standard state of pure compound i, standard pressure p 0 and temperature T, while γ j, i is the stoichiometric coefficient of species i, within reaction j, ΔH i is the molar enthalpy of formation, ΔS i is the molar entropy of formation and K R, j is the reaction specific, kinetic relation.
Table 3: Comparison of conservative quantities Ψ, potentials ϕ and capacities κ, including units.Adaption to the field of thermo-chemistry, in analogy to Table 2 Figure 4: Graphical template for the system dynamic scheme to resolve dynamics and equilibrium of any chemically reacting system.
Combining the symbolism of Fig. 4. with the relations from Table 4, yields the ordinary differential equations to solve for the number of molecules i and total Gibbs free energy of species i against reaction time t (see Eqn. 22 and Eqn.23 respectively).(22) (23)

EXAMPLE AND VALIDATION: HETEROGENEOUS WATER GAS SHIFT REACTION EQUILIBRIUM
As an example of the presented modeling methodology, the System Dynamic scheme shall be applied to the case of the heterogeneous Water Gas Shift reaction (WGS), as seen in Eqn.24.Furthermore this example shall be used to verify the functionality of the modeling procedure.
(24) 3.3.1.System Dynamic model of the Water Gas Shift reaction system Fig. 5 shows the graphical interpretation of the complete, System Dynamic model of the Water Gas Shift reaction system.
In Fig. 5 the rectangles are containers or integrators of the relevant conservative quantities, which are the numbers of molecules N CO , N H2O and N H2 (red), as well as the respective Gibbs free energies G i (green).Underlying relations (blue spots) correspond to those presented in Table 4, with reaction j=WGS and species i = CO, H 2 O and H 2 .Note that the extent of species N i corresponds to the capacity of the denoted Gibbs free energy 146 System Dynamic modelling approach for resolving the thermo-chemistry of wood gasification By adding a relatively simple sub-model, shown in Fig. 6, and by using I G, WGS = ΣI G, i , the total reduction of Gibbs free energy within the entire system ΔG WGS can be monitored from starting conditions to equilibrium.

Results and Validation
A validation of the Water Gas Shift model, described in chapter 3.3.1,was achieved as follows.
First, the model was realized within the System Dynamic software Berkeley Madonna [32].Arbitrary starting conditions were assumed and the molar composition as well as ΔG WGS of the Water Gas shift system was then plotted over time, as seen in Fig. 7. Thereby the constant result towards the right of the time axis, at minimum ΔG WGS , corresponds to the equilibrium state of the system.Note that the total Gibbs free Energy reduction of -8441.6J depends on the system parameters p and T, as well as on the starting conditions, the initial numbers of molecules N i, 0 .Validation of the model could then be achieved by retrieving the exact same results for the equilibrium molar composition x i, End and ΔG WGS , through analytical calculation.The analytical calculation scheme to retrieve those results is shown in Table 5.

APPLICATION TO WOOD GASIFICATION
System Dynamic methodology and its application to resolving the thermodynamics and kinetics of a chemical system, has been introduced in chapters 3.1 and 3.2.Its capacity has then been demonstrated on the validated example of the Water Gas Shift reaction within chapter 3.3.On this basis, the concept is now applied to a more complex, thermo-chemical system, featuring i = 7 species and j = 9, parallel types of reactions: the example of wood gasification.
A full System Dynamic model for resolving the thermo-chemistry of wood gasification has been created within the software Berkeley Madonna.It features all i = 7 relevant chemical species, as well as all j = 9 chemical reactions, cited in chapter 1.2.In the following, only some exemplary outtakes of the model are presented, since the general, underlying methodology has already been discussed.

Species Balance
One very basic aspect of the model is setting up all i = 7 species balances, which can be achieved rather efficiently by using a graphical programming tool such as Berkeley Madonna.The sub-model for considering the species containers, as well as their respective, reactive molar fluxes is shown in Fig. 8.Note that within Fig. 8 all molar fluxes, associated with the various oxidation reactions, shown in Table 1, are summarized within one, single oxidation flux.The graphical interpretations, depicted within Fig. 8 can be summarized and solved as an equation system, given by Eqn.22.

CALCULATION OF REACTION KINETICS
Another model aspect to be pointed out in this context is the calculation of reaction kinetics.
Fig. 9 presents the exemplary sub-model for calculating the reaction kinetics of the Boudouard reaction.It highlights the thermodynamic impact on chemical flux direction, rather than the implementation of any specific rate law.Thus Fig. 9 depicts the steps from thermodynamic standard values, calculated on the basis of NASA technical memo 4513 [8], to Gibbs free energies of formation of species i within the mixture, towards the Gibbs free energy of Boudouard reaction, considering relations from Table 4.

VALIDATION, COMPARISON AND RESULTS
In this chapter the validity of the two models, presented within chapters 2 and 3 is proven, by vice versa comparison of the results.Furthermore the applicability and the potential of the System Dynamic model, in terms of researching wood gasification processes, is demonstrated.
Measurement of local, dynamic gas compositions within a thermal reactor is nearly impossible.Thus a reasonable way to verify gas composition calculations is to compare the results of completely different methods among each other.Fig. 10 shows such a comparison and thus validates the two solvers, presented within this paper.The gas composition x i and the total Gibbs free energy G of a wood gas system are thereby plotted against the kinetic ratio R OxCH4 .The kinetic ratio R OxCH4 is here defined as the rate constant of methane oxidation over the sum of rate constants of all occurring oxidation reactions.The output of the Lagrangian Multiplier -Gibbs free energy minimization solver remains unaffected by kinetic aspects since it just returns the thermodynamic equilibrium composition of least 150 System Dynamic modelling approach for resolving the thermo-chemistry of wood gasification Figure 9: Sub-model for calculating Boudouard reaction kinetics on basis of [8] and Table 4 Gibbs free energy (dashed lines).The System Dynamic solver however, can consider these effects and yields varying gas output concentrations (full lines), as methane oxidation kinetics change.In addition to gas output compositions, the System Dynamic solver also returns a result for the total Gibbs free energy of the system (orange).At the very spot, where minimum total Gibbs free energy is found (black vertical line), the results of the System Dynamic solver and the Lagrangian Multiplier solver intersect.A vice versa validation of the solvers is thus achieved.Any discussion regarding specific kinetic rate laws of wood gasification reactions, such as described e.g. in [3], [9] and [17], has been omitted within this paper.The reason for this is that thorough investigation of a complex thermo-chemical process, such as wood gasification cannot be conducted by nailing any one single type of rate law.It much rather requires a broader view of a wide range of possible kinetic constellations.As shown in Fig. 10, the System Dynamic solver can provide just that point of view.The results demonstrate the relatively significant impact of kinetic issues on output gas composition.Thus the fact that a static equilibrium calculation does not suffice to study the full depth of the wood gasification mechanism, is underlined.Even more so, because measured wood gas compositions do in praxis actually vary over relatively wide ranges of composition [1].
In addition to that and in contrast to Gibbs free energy minimization solvers (e.g.: Lagrangian Multiplier), the System Dynamic solver can be used to study dynamically changing gas compositions.One such example is shown in Fig. 11.The result plotted there, shows the composition of a reacting wood gas system, with fixed atomic balance ratios, and fixed temperature, from starting conditions to equilibrium formation.The average molar Gibbs free energy ΔG, is plotted out as well.Its minimum obviously corresponds with compositional equilibrium formation.

CONCLUSION AND OUTLOOK
Within this work, two wood gas composition solvers have been presented and compared: The Lagrangian Multiplier solver and the System Dynamic solver.While the Lagrangian Multiplier solver has been programmed on the basis of templates from literature [6], the System Dynamic wood gasification model, is a novelty.Table 6 sums up a simplified comparison between the two types of solvers.
While the System Dynamic solver cannot compete in terms of efficiency, speed nor memory requirement, it does allow studying wood gasification processes on a deeper, more diverse level, than ordinary Gibbs free energy minimization tools.
In the course of ongoing work, the presented System Dynamic solver will be fully incorporated into a full wood particle gasification model.This new program will be set up, using [2] as role model, but will be wholly based upon System Dynamic methodology.

Figure 1 :
Figure 1: Sketch of wood gas species (composition x i ) in thermo-chemical interaction.System is fed with oxygen flux I O2 from process air and pyrolysis gas flux I COnHm from gasified cellulose

Figure 2 :
Figure 2: Wood gas equilibrium composition against temperature (left) and ratio of oxygen to carbon atoms RO/C (right).

Figure 3 :
Figure 3: System dynamic interpretation of general balances: container (=Integrator) of with in and out fluxes, filling height ϕ (=potential) and cross section κ (=capacity)

Figure 5 :Figure 6 :
Figure 5: Graphical interpretation of system dynamic model of the water gas shift reaction system

Figure 8 :
Figure 8: Sub model for considering species balances and respective molar fluxes from chemical reactions

Figure 10 :
Figure 10: Wood gas compositions x i and total gibbs free energy G against the ratio of kinetic constant of methane oxidation over sum of all kinetic oxidation constants, R OxCH4 .Comparison of results from lagrangian multiplier solver (dashed lines, *) and system dynamic solver (full lines); T = 810K, R O/C = 1.5, R H/C = 2.88

Figure 11 :
Figure 11: Example of dynamically changing wood gas composition and average molar gibbs free energy against time.T = 810 K, R O/C = 1.5, R H/C = 2.88; results are produced by the system dynamic solver

Table 2 :
Comparison of some conservative quantities potentials ϕ, and capacities κ, including units

Table 4 :
Relations in addition to Fig.3

Table 6 :
Comparison between the Lagrangian Multiplier solver and the system dynamic solver