## Gibbs energy minimization and the NIST webbook

| categories: optimization | View Comments

Gibbs energy minimization and the NIST webbook

# Gibbs energy minimization and the NIST webbook

John Kitchin

## Contents

In Post 1536 we used the NIST webbook to compute a temperature dependent Gibbs energy of reaction, and then used a reaction extent variable to compute the equilibrium concentrations of each species for the water gas shift reaction.

Today, we look at the direct minimization of the Gibbs free energy of the species, with no assumptions about stoichiometry of reactions. We only apply the constraint of conservation of atoms. We use the NIST Webbook to provide the data for the Gibbs energy of each species.

As a reminder we consider equilibrium between the species , , and , at 1000K, and 10 atm total pressure with an initial equimolar molar flow rate of and .

H2

H2O

CO

CO2

```function nist_wb_constrained_min
```
```T = 1000; %K
R = 8.314e-3; % kJ/mol/K

P = 10; % atm, this is the total pressure in the reactor
Po = 1; % atm, this is the standard state pressure
```

## We define species like this.

We are going to store all the data and calculations in vectors, so we need to assign each position in the vector to a species. Here are the definitions we use in this work.

```1  CO
2  H2O
3  CO2
4  H2```
```species = {'CO' 'H2O' 'CO2' 'H2'};
```

## Heats of formation at 298.15 K

```Hf298 = [
-110.53  % CO
-241.826 % H2O
-393.51  % CO2
0.0]; % H2
```

## Shomate parameters for each species

`   A             B           C          D          E       F          G       H`
```WB = [25.56759  6.096130     4.054656  -2.671301  0.131021 -118.0089 227.3665   -110.5271  % CO
30.09200  6.832514     6.793435  -2.534480  0.082139 -250.8810 223.3967   -241.8264    % H2O
24.99735  55.18696   -33.69137    7.948387 -0.136638 -403.6075 228.2431   -393.5224   % CO2
33.066178 -11.363417  11.432816  -2.772874 -0.158558 -9.980797 172.707974    0.0]; % H2
```

## Shomate equations

```t = T/1000;
T_H = [t; t^2/2; t^3/3; t^4/4; -1/t; 1; 0; -1];
T_S = [log(t); t; t^2/2; t^3/3; -1/(2*t^2); 0; 1; 0];

H = WB*T_H; % (H - H_298.15) kJ/mol
S = WB*T_S/1000; % absolute entropy kJ/mol/K

Gjo = Hf298 + H - T*S; % Gibbs energy of each component at 1000 K
```

## Gibbs energy of mixture

this function assumes the ideal gas law for the activities

```    function G = func(nj)
Enj = sum(nj);
G = sum(nj.*(Gjo'/R/T + log(nj/Enj*P/Po)));
end
```

## Elemental conservation constraints

We impose the constraint that all atoms are conserved from the initial conditions to the equilibrium distribution of species. These constraints are in the form of A_eq x = beq, where x is the vector of mole numbers. CO H2O CO2 H2

```Aeq = [ 1    0    1    0   % C balance
1    1    2    0   % O balance
0    2    0    2]; % H balance

% equimolar feed of 1 mol H2O and 1 mol CO
beq = [1   % mol C fed
2   % mol O fed
2]; % mol H fed
```

## Constraints on bounds

This simply ensures there are no negative mole numbers.

```LB = [0 0 0];
```

## Solve the minimization problem

```x0 = [0.5 0.5 0.5 0.5]; % initial guesses

options = optimset('Algorithm','sqp');
x = fmincon(@func,x0,[],[],Aeq,beq,LB,[],[],options)
```
```Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in
feasible directions, to within the default value of the function tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.

x =

0.4549    0.4549    0.5451    0.5451

```

## Compute mole fractions and partial pressures

The pressures here are in good agreement with the pressures found in Post 1536 . The minor disagreement (in the third or fourth decimal place) is likely due to convergence tolerances in the different algorithms used.

```yj = x/sum(x);
Pj = yj*P;

for i=1:numel(x)
fprintf('%5d%10s%15.4f atm\n',i, species{i}, Pj(i))
end
```
```    1        CO         2.2745 atm
2       H2O         2.2745 atm
3       CO2         2.7255 atm
4        H2         2.7255 atm
```

## Computing equilibrium constants

We can compute the equilibrium constant for the reaction . Compared to the value of K = 1.44 we found at the end of Post 1536 , the agreement is excellent. Note, that to define an equilibrium constant it is necessary to specify a reaction, even though it is not necessary to even consider a reaction to obtain the equilibrium distribution of species!

```nuj = [-1 -1 1 1]; % stoichiometric coefficients of the reaction
K = prod(yj.^nuj)
```
```K =

1.4359

```
```end

% categories: optimization
% tags: Thermodynamics, reaction engineering
```