1% Tutorial 4:   Chemical Equilibrium
2%
3%   Topics:
4%     - the equilibrate method
5%     - specifying fixed TP, HP, UV, SV, or SP
6%     - checking reaction rates of progress
7%
8help tut4
9
10% To set a gas mixture to a state of chemical equilibrium, use the
11% 'equilibrate' method.
12%
13g = GRI30('None');
14set(g,'T',1200.0,'P',oneatm,'X','CH4:0.95,O2:2,N2:7.52')
15equilibrate(g,'TP')
16
17% The above statement sets the state of object 'g' to the state of
18% chemical equilibrium holding temperature and pressure
19% fixed. Alternatively, the specific enthalpy and pressure can be held
20% fixed:
21
22disp('fixed H and P:');
23set(g,'T',1200.0,'P',oneatm,'X','CH4:0.95,O2:2.0,N2:7.52');
24equilibrate(g,'HP')
25
26
27% Other options are
28%     'UV'   fixed specific internal energy and specific volume
29%     'SV'   fixed specific entropy and specific volume
30%     'SP'   fixed specific entropy and pressure
31
32disp('fixed U and V:');
33set(g,'T',1200.0,'P',oneatm,'X','CH4:0.95,O2:2,N2:7.52');
34equilibrate(g,'UV')
35
36disp('fixed S and V:');
37set(g,'T',1200.0,'P',oneatm,'X','CH4:0.95,O2:2,N2:7.52');
38equilibrate(g,'SV')
39
40disp('fixed S and P:');
41set(g,'T',1200.0,'P',oneatm,'X','CH4:0.95,O2:2,N2:7.52');
42equilibrate(g,'SP')
43
44% How can you tell if 'equilibrate' has correctly found the
45% chemical equilibrium state? One way is verify that the net rates of
46% progress of all reversible reactions are zero.
47
48% Here is the code to do this:
49set(g,'T',2000.0,'P',oneatm,'X','CH4:0.95,O2:2,N2:7.52');
50equilibrate(g,'TP')
51rf = rop_f(g);
52rr = rop_r(g);
53format short e;
54for i = 1:nReactions(g)
55   if isReversible(g,i)
56      disp([i, rf(i), rr(i), (rf(i) - rr(i))/rf(i)]);
57   end
58end
59
60
61% You might be wondering how 'equilibrate' works. (Then again, you might
62% not, in which case you can go on to the next tutorial now.) Method
63% 'equilibrate' invokes Cantera's chemical equilibrium solver, which
64% uses an element potential method. The element potential method is
65% one of a class of equivalent 'nonstoichiometric' methods that all
66% have the characteristic that the problem reduces to solving a set of
67% M nonlinear algebraic equations, where M is the number of elements
68% (not species). The so-called 'stoichiometric' methods, on the other
69% hand, (including Gibbs minimization), require solving K nonlinear
70% equations, where K is the number of species (usually K >> M). See
71% Smith and Missen, "Chemical Reaction Equilibrium Analysis" for more
72% information on the various algorithms and their characteristics.
73%
74% Cantera uses a damped Newton method to solve these equations, and
75% does a few other things to generate a good starting guess and to
76% produce a reasonably robust algorithm. If you want to know more
77% about the details, look at the on-line documented source code of
78% Cantera C++ class 'ChemEquil' at https://cantera.org.
79
80%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
81clear all
82cleanup
83%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
84%   end of tutorial 4
85%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
86