使用Mathematica对化学反应平衡进行模拟
使用Mathematica对化学反应平衡进行模拟
反应平衡动力模拟
对于任何一个溶液中发生的基元可逆反应可以用以下方式表示
graph LR
A[反应物] --k1--> B
B[生成物] --k2--> A\(k_1\)和\(k_2\)是此反应的反应速率常数
反应方程式可表示为
\[ A\leftrightharpoons B\quad K=\frac{k_1}{k_2} \]
或者也可以写成下面这样（我瞎编的）
\[ A\stackrel{k_2}{\underset{k_1}{\leftrightharpoons}} B \]
由此可以写出物质\(A\)和\(B\)的浓度关于时间\(t\)的微分方程
\[ \begin{cases} \frac{dA}{dt}=-k_1A+k_2B\\ \frac{dB}{dt}=+k_1A-k_2B \end{cases} \]
再翻译翻译就是
\[ \begin{cases} A'(t)=-k_1A(t)+k_2B(t) \\ B'(t)=+k_1A(t)-k_2B(t) \end{cases} \]
如果反应的反应物或生成物有多个物质或化学计量数不为\(1\)，则需要带入物质浓度的幂
\[ a_1A_1+a_2A_2+\cdots+a_nA_n\stackrel{k_2}{\underset{k_1}{\leftrightharpoons}}b_1B_1+b_2B_2+\cdots+b_mB_m \]
\[ \begin{cases} \frac{dA_i}{dt}=-k_1\prod_{u=1}^nA_u^{a_u}+k_2\prod_{v=1}^mB_v^{b_v} \\ \frac{dB_j}{dt}=+k_1\prod_{u=1}^nA_u^{a_u}-k_2\prod_{v=1}^mB_v^{b_v} \end{cases}\\ \begin{cases} A_i'(t)=-k_1\prod_{u=1}^nA_u^{a_u}(t)+k_2\prod_{v=1}^mB_v^{b_v}(t) \\ B_j'(t)=+k_1\prod_{u=1}^nA_u^{a_u}(t)-k_2\prod_{v=1}^mB_v^{b_v}(t) \end{cases} \]
Mathematica!
完整程序会放在最后
首先我们需要一个Mathematica
其次需要一些用来测试的化学反应
这里以碳酸的电离平衡为例
\[ \scriptsize \begin{matrix} \mathrm{CO_2}+&\mathrm{H_2O}&\stackrel{k_{12}}{\underset{k_{11}}{\leftrightharpoons}}&\mathrm{H_2CO_3}&&&&&&&\quad K=\frac{k_{11}}{k_{12}}=&\frac{1}{600}&& \\ &&&\mathrm{H_2CO_3}&\stackrel{k_{22}}{\underset{k_{21}}{\leftrightharpoons}}&\mathrm{HCO_3^-}&&+&\mathrm{H^+}&&\quad K=\frac{k_{21}}{k_{22}}=&4.3&\cdot&10^{-7} \\ &&&&&\mathrm{HCO_3^-}&\stackrel{k_{32}}{\underset{k_{31}}{\leftrightharpoons}}\mathrm{CO_3^{2-}}&+&\mathrm{H^+}&&\quad K=\frac{k_{31}}{k_{32}}=&5.6&\cdot&10^{-11} \\ &\mathrm{H_2O}&\stackrel{k_{42}}{\underset{k_{41}}{\leftrightharpoons}}&&&&&&\mathrm{H^+}&+\mathrm{OH^-}&\quad K=\frac{k_{41}}{k_{42}}=&&&10^{-14}&=K_w \end{matrix} \]
初始化
首先，为了避免奇奇怪怪的\(\mathfrak{BUG}\)，先清除全局变量
1 | Clear["Global`*"] |
其次需要声明方程的反应速率常数
1 | b=5; |
T是时间上限，即微分方程自变量的最大值
其中\(k_{ij}\)值与化学方程式中的\(k_{ij}\)有些许不同和错位，亿些细节放到最后再解释
接着就可以开始解微分方程了~
微分方程！
1 | {H2CO3,H,OH,HCO3,CO3,CO2}=NDSolveValue[{ |
NDSolveValue[微分方程列表,所求因变量列表,{自变量,自变量最小值,自变量最大值}]，返回数值解的列表，其元素为所求因变量与自变量的数值关系（大概）
可视化
1 | Plot[{H2CO3[t]},{t,0,T},PlotLegends -> "Expressions"]; |
Plot[因变量列表,{自变量,自变量最小值,自变量最大值},其他参数]，返回绘制的折线图。其他参数中PlotLegends -> "Expressions"可显示图例（Plot[]中只绘制一个因变量时无作用）
以上代码绘制了\(\mathrm{H_2CO_3}\)，\(\mathrm{HCO_3^-}\)，\(\mathrm{CO_3^{2-}}\)，\(\mathrm{H^+}\)，\(\mathrm{OH^-}\)，\(p\mathrm{H}\)，\(p\mathrm{OH}\)以及三个可逆反应的\(Q-K\)关系
\(\mathfrak{BUG}\)
- 必须要加
Clear["Global`*"]，否则会出现无法预料的\(\mathfrak{BUG}\) - \(K=\frac{k_1}{k_2}\)，但对\(k_1\)和\(k_2\)本身的大小并没有要求，即每一组符合条件的\(k_1,k_2\)都可以用于表示同一个方程。对于多个反应，改变其中几个反应的\(k_1,k_2\)（保证\(k_1,k_2\)的比值\(K\)不变），反应达到平衡时各组分浓度不会发生改变，但平衡前各组分浓度变化的过程会发生改变。但\(|k_1|\)，\(|k_2|\)越大，方程越快达到平衡。
- 程序中所有\(k_{ij}\)是经过调试找到的能较快达到平衡且较“好看”的组合。增加了一个参数
b以更好控制\(k\)的数量级。这可能并不符合实际反应原理，但有利于调试 - 加入\(\mathrm{H_2O}\stackrel{k_{42}}{\underset{k_{41}}{\leftrightharpoons}}\mathrm{H^+}+\mathrm{OH^-}\)且画出\(p\mathrm{H},p\mathrm{OH}\)后出现了\(p\mathrm{H}+p\mathrm{OH}\neq 14\)的情况。需要将\(k_{41},k_{42}\)增大使得此反应更快平衡（或增大反应时间
T也能看到效果）这里调整到\(k_{41}=1,k_{42}=10^{14}\)，若再调高则超出了Mathematica的int范围（¿） - 测试中各组分浓度初值尽量不超过\(1(mol/L)\)（我猜是这个单位？），否则可能出现不论如何调整\(k\)和
T也没有反应的情况，原因不明，我猜是因为水解平衡反应不过来了 - 测试中组分浓度越低，效果越好，这大概就是越稀越水解吧
T较大时可视化的图表中的图线会出现抽风的情况，我猜是因为每次\(dt\)采样的间隔较大的原因- 如果认为\(\mathrm{CO_2}\)生成后立即溢出，即视\(c(\mathrm{CO_2})\)在任何时间均为常数，则\(\mathrm{CO_3^{2-}}\)在酸性环境下生成\(\mathrm{CO_2}\uparrow\)的反应可视化效果较好
完整程序
1 | Clear["Global`*"] |
测试输出
