# Description of multiple equilibria in solution

Solution chemistry is still widely teaching in the first year post-bac. Typically in pH problems, you are looking for an approaching solutions using the dominant reaction approach. However, you must be aware than today, it's relatively easy to implement numerical algorithm to find the exact solution. 

These numerical approaches are largely used to discuss the validity of the common approximations. Pedagogical software like ChimGéné able someone to simulate E-pH diagrams, i-E curves (<https://chimsoft.com/chimsoft/chimgene>), ...

Your goal is to understand how to describe a chemical equilibrium, without any approximation, using simple numerical recipes.

# A first approach : find the equilibrium of a system

## Weak acid equilibrium in water

$$
\begin{align} HA (aq) &= H^+ (aq) + A^-(aq) & K_A\\ H_2O(l) &= H^+(aq) + HO^- (aq) & K_e\end{align}
$$
**Q. 1** The system is describe by 4 degrees of liberty. Identify them. Write two equilibria and two conservation laws.

**R.1**

**Q.2** Combine the equation in order to get a third order polynomial which describe the problem.

**R.2** 

In order to solve numerical the previous polynomial, called $\mathit{P}$, you might use robust numerical methods. One of the solvers is the dichotomy algorithm.



This algorithm consists in defining an interval $[x_1, x_2]$ where you want to find a solution :

$$
\mathit{P} (x_s) = 0 \quad \text{with } x_s \in [x_1, x_2] = [x_1^{(0)}, x_2^{(0)}]
$$
This n-step solution $x_s^{(n)}$ of the iterative method might converge to the final solution $x_s \pm \epsilon$ . $\epsilon$ is a user-controlled convergence criteria. The iterative transformation is defined as : 

$$
\begin{align}
\forall n > 1, \quad &\begin{cases}
x_1^{(n)} &= \frac{x_1^{(n-1)} + x_2^{(n-1)}}{2} \\
x_2^{(n)} &= x_2^{(n-1)}
\end{cases} \text{ if } P(x_1^{(n-1)})\cdot P(x_1^{(n)})>0 \\
\quad &\begin{cases}
x_2^{(n)} &= \frac{x_1^{(n-1)} + x_2^{(n-1)}}{2} \\
x_1^{(n)} &= x_1^{(n-1)}
\end{cases} \text{ else } 
\end{align}
$$
where $n$ is a given iteration step.

**Q.3** Interpret the transformation. Write a code using the dichotomy algorithm. Benchmark the function on a polynomial of your choice and plot the iterative process.

**R.3** 


```
#R.3
```

*Note :* Implement conditions on the maximum number of iteration and on the value of $x_2$ compared to $x_1$.

Now, let's implement a second iterative method called *Regula falsi* method*.* The process is exactly the same than the one before, but the iterative transformation is different :

$$
\begin{align}

\forall n > 1, \quad &\begin{cases}
x_1^{(n)} &= x_1^{(n-1)} - \frac{x_1^{(n-1)} - x_2^{(n-1)}}{P(x_1^{(n-1)}) - P(x_2^{(n-1)})} P(x_1^{(n-1)})\\
x_2^{(n)} &= x_2^{(n-1)}
\end{cases} \text{ if } P(x_1^{(n-1)})\cdot P(x_1^{(n)})>0 \\
\quad &\begin{cases}
x_2^{(n)} &= x_2^{(n-1)} - \frac{x_2^{(n-1)} - x_1^{(n-1)}}{P(x_2^{(n-1)}) - P(x_1^{(n-1)})} P(x_2^{(n-1)}) \\
x_1^{(n)} &= x_1^{(n-1)}
\end{cases} \text{ else } 
\end{align}
$$
**Q.4** Interpret the transformation. Add this transformation in the iterative solver presented in **Q.3**. Benchmark the function on the polynomial used in **Q.3** and plot the iterative process. Compare the convergence time of the two methods.


```
#R.4
```

**R.4** 

**Q.5** Consider the acetic acid in water $pK_A (CH_3COOH/CH_3COO^-, 25 ^\circ C) = 4.8$ . Find the solution for three concentrations : 

$$
c = 10^{-2}, 10^{-4}, 10^{-6} mol \cdot L^{-1}
$$
**R.5** Compare the results obtained from **Q.3** algorithm with the one obtain with the dominant reaction approximation.

## Another way to solve the problem

The first step is to implement a Newton-Raphson algorithm. This robust method consists in finding the roots of an $N$-variable real function $f$ : 

$$
f(X_s) = 0
$$
In 1D, the transformation is easy to visualize. From the knowledge of the tangent at a given point, the cross-point between $y = 0$ and the tangent, allow you to approach the root of $f$.

 



**Q.6** Using the first order development of f around $x_n$ write a single step transformation $x_n \to x_{n+1}$

**R.6**

As all iterative method, convergence criteria might be defined : 

1. $\lvert f(x_n) \rvert < \epsilon_1$, the residue in between your function and 0;
2. $\lvert x_{n+1} - x_n\rvert < \epsilon_2$, the resolution of your solution, estimating the accuracy of your estimation compared to the exact solution $\lvert x_1^{(\infty)} - x_s\rvert < \epsilon_2$ ;
3. a maximum number of iteration

**Q.7** Implement this transformation for 1D function as in **Q.3**. Benchmark the function on the polynomial used in **Q.3** and plot the iterative process. Compare the convergence time of the previous methods.

**R.7** 

Chemical media are complex system described by : 

* multiple equilibria ; 
* multiple charge/mass conservation law .

In our simple weak acid example, system state belongs to the value of a single vector : 

$$
X = \begin{pmatrix}
[AH] \\
[A^-] \\
[H^+] \\
[HO^-]
\end{pmatrix}
$$
The chemical system are always multicomponents system (at least, one solvent and one solute). So mathematically, any transformation $f$ of a chemical system is multidimensional. You must generalized the Newton-Raphson method in $N$ dimensions. The extension to multiple dimension can be summarize by this block diagram : 



**Q.8** What is the name of the matrix $[\nabla f]$ ?

**R.8** 

**Q.9** Using the system defined in **Q.1**, write the matrix $[\nabla f]$ for a system state described by $X$.

**R.9** 

**Q.10** Develop a code using a N-dimension Newton-Raphson algorithm. Find the equilibrium state of acetic acid **Q.5** in the three cases.


```
#R.10
```

*Note :*  You can write the equilibria relations in log scale.

Before starting with more complex systems, you shall be aware of some limits of the previous approach. Classically, Newton-Raphson might be unstable if your system present too many degrees of freedom. Each result obtained must be interpreted as we say : "avec l'aide de votre $6^e$ sens chimique. Le plus important de tous !" .

# Treat a complex system with different kind of equilibria

## Single equilibrium state for one composition

Using the numerical recipes we study in the previous section, combined with a well-written object-oriented code, we are able to describe most of the single phase chemical equilibria.

Consider a well-known chemical system $Fe(II)/Ce(IV)$ composed of : 

* 1L water solution,
* 0.1 mol of $FeCl_2$, 0.05 mol of $Ce(SO_4 )_2$ and 0.1 mol of $H_2SO_4$,
* restrict the system state to these species :  

$$
H^+, HO^-, Cl^-, Fe^{3+}, Fe^{2+}, Ce^{4+}, Ce^{3+}, SO_4^{2-}, FeSO_4^+
$$
**Q.11** Write the 11 equations which describe the system.

**R.11**

**Q.12** Modify the code **Q.10** for an undermined N number of chemical species described by N (non-)linear equations. Find the equilibrium state of the previous system.


```
#R.12
```

*Data : Thermodynamic information for the Fe/Ce system in water at 25°C*

$$
\begin{align}
pK_e &= 14 \\
\log \beta (FeSO_4^+/Fe^{3+}) &= 2.03\\
E^\circ (Fe^{3+}/Fe^{2+}) &= 0.77 V/ESH \\
E^\circ (Ce^{4+}/Ce^{3+}) &= 1.71 V/ESH
\end{align}
$$
## Multiple equilibrium state for successive composition : titration

Titration reactions are necessarily process controlled by the thermodynamic (the reaction is fast) so they can be describe by the previous formalism. An experimenter prepares this experimental setup.

*Material :*

* $0.1 mol \cdot L^{-1}$ Mohr's salt solution , $0.1 mol \cdot L^{-1}$ titration solution of ceric disulfate,
* Volumetric Pipette, buret, magnetic stirrer, platin electrode, saturated calomel electrode,
* voltmeter.

*Protocol :*

Put $20 cm^{3}$ of ferric solution in a beaker and add a magnetic stirrer. Stire the solution and introduce both electrode on the solution. Connect it to a voltmeter. Fill the burette with the ceric disulfate reaction. Start the burette fall and acquire the potential.



To treat this problem numerically, one might consider it as succession of N equilibria. The complete burette fall is considered in the $[0; 2V_e]$ interval volume.

The volume increment between two equilibria corresponds to equal volume burette fall, defined as : 

$$
\delta V = \frac{2V_{eq}}{N-1}
$$
where $V_{eq}$ is equivalence volume, $N$the number of point in the curve. The iterative method can be simplified by the block diagram below : 



**Q.13** Calling the equilibrium solver developed in **Q.12**, developed an iterative method to obtain the curve : 

$$
E = f(V) \quad [\emptyset]
$$
Save the volume, the electrochemical potential and the concentration of each chemical species in a npz file.


```
#R.13
```

**Q.14** Write a post analysis script which import the outputs. Check the electric potential profile to what you are expecting and benchmark your code.

Implement four methods to determine the volume of end of the titration : 

* tangent method;
* first derivative method;
* second derivative method;
* Gram method.

Compare the value obtained with the one expected. Discuss the interest and precision of each method. 


```
#R.14
```

*Note :* You might have data with different values of $N$ to get an interesting discussion.

**Q.15**  How the pH changing during the titration process ? Add an acid buffer and compare the evolution.


```
#R.15
```

**R.15** 

*Data :* Complementary thermodynamic information in water at 25°C

$$
\begin{align} 
pK_A (H_2SO_4 / HSO_4^-) & < 0 \\
pK_A (HSO_4^- / SO_4^{2-}) & = 2.0
\end{align}
$$
# An original approach : Use semantic tableau

Now, you get some ideas about solution model. Your goal in this section is to implement a more elegant methods. One interesting construction is to create a logic process to get a semantic tableau.

Let's imagine you are studying a fix concentration of iron in a solution with one degree of oxidation : $III$. We choose to describe the solution as composed by :

$$
H^+, Fe^{3+}, HO^-, Fe(OH)_1^{2+}, Fe(OH)_2^{+}, Fe(OH)_4^{-} \quad (\star)
$$
with these equilibria constants in water at 25°C : 

$$
\begin{align}
K_e &= 10^{-14} \\ 
\log \beta_1 &= 11.81 \\  
\log \beta_2 &= 22.33\\
\log \beta_4 &= 34.4
\end{align}
$$
**Q.16** How many concentrations an experimenter should know to determine the concentration of all the species using the mass action law. Why ?

**R.16**

Now, you must simplify the problem !  In order to construct an easiest mathematical problem, the number of variables will be reduce from 6 to 2. We choose the pH (so $[H^+]$) and $[Fe^{3+}]$ as our basis set to generate all the configurations. 

**Q.17** Using the equilibrium constants, rewrite the concentration of the 6 species $(\star)$ in such a way : 

$$
\forall i \in [1, 6], C_i = [H^+]^{n_H} \cdot[Fe^{3+}]^{n_{Fe}} \cdot 10^{-pK}
$$
where $n_H$, $n_{Fe}$ and $pK$ are three parameters to determine.

**R.17** 

**Q.18** Using **Q.17**, complete this tableau : 



What is the meaning of a row ? What is the meaning of the first two columns $H_T$ and $Fe_T$.

**R.18** 

**Q.19** Now we want to convert the problem in a simple matrix problem where the unknown component is **X** : 

$$
X = \left( \log [H^+] \log[Fe^{3+}] \right)
$$
From the tableau in **Q.18**, write :

* a vector (6,) containing the log species concentrations **(C);** 
* a vector (6,) containing the reactions constants **(pK);**
* a matrix (6,2) containing the stoichiometric coefficients **(A)**;
* a vector(2,) containing the conservation law constrain **(T)**. 

Write **C** as a function of **pK**, **A** and **X**.

**R.19** 

To solve the problem, we want to minimize a relation. The mass balance is an ideal quantity to minimize in this problem :

$$
R = {^t A} \cdot 10^C - T
$$
where **T** is the mass contrain, **$^t A$** the transposition of stoichiometric coefficient, **T** the conservation law constrain.

**Q.20** What does mass balance represent ? Using the previously Newton-Raphson method, minimize the residuals in the mass balance. For a total E-5 M of iron, find the reconstruct the distribution of iron species as a function of pH.

 *Note :* To minimize **R** as a function of **X**, you need to find the jacobian and fix an origin.

**R.20** 


```
#R.20 
```

# Open Question

**R.21** Generalize the semantic tableau method for multiphasic system. You can add $Fe(OH)_3(s)$ to extend the previous problem for higher concentration in iron.


```
#R.21
```

*Data :* $pK_s(Fe(OH)_3 (s), 25 ^\circ C) = 15.09$

*Note :* There is as many tableau as there is phases ...

This project is derived from the architecture of *Chimie des solutions : De l'élémentaire aux calculs numériques, Marc Blétry & Marc Presset, Chapter 12 : Modélisation numérique des équilibres en solution idéale.*

To complete, these resources are used : 

*L'Oxydoréduction : Concepts et expériences*, Jean Sarrazin & Michel Verdaguer.

*Chemical Reaction Equilibrium Analysis Theory and Algorithms*,  William R. Smith and Ronald W. Missen

*Solution of simultaneous chemical equilibria in heterogeneous systems : implementation in Matlab*, Scott Smith, 23 September 2019

3986b795111caac3e9210173528e8e20c16d228a
ipynb
Jupyter Notebook
Cours/TP_ProjetsInformatiques_L3/description-of-multiple-equilibria-in-solution.nextjournal.ipynb
jfolivie/Ressources
47ef37c017675798dcb7e1070df2e26106fa9e70
CC0-1.0
Cours/TP_ProjetsInformatiques_L3/description-of-multiple-equilibria-in-solution.nextjournal.ipynb
jfolivie/Ressources
47ef37c017675798dcb7e1070df2e26106fa9e70
CC0-1.0
Cours/TP_ProjetsInformatiques_L3/description-of-multiple-equilibria-in-solution.nextjournal.ipynb
jfolivie/Ressources
47ef37c017675798dcb7e1070df2e26106fa9e70
CC0-1.0
Qwen/Qwen-72B
1. YES
2. YES