使用scipy与matplotlib绘制勒夏特列原理化学平衡互动图像

存在的错误点 #

• 反应速率方程逻辑完全错误：针对二级基元可逆反应A + B <=> C + D，正反应速率应为k1*Ca*Cb（与两种反应物浓度乘积成正比），逆反应速率应为k2*Cc*Cd（与两种产物浓度乘积成正比），原代码中速率项仅关联了Ca的浓度，不符合反应速率计算规则
• 各物质浓度变化率计算逻辑错误：A、B的浓度变化应为逆反应生成量减去正反应消耗量，C、D的浓度变化应为正反应生成量减去逆反应消耗量，原代码的变化率计算完全不符合反应计量关系
• 求解后浓度赋值错误：原代码错误将产物Cc、Cd的取值指向了反应物Ca、Cb的计算结果，导致输出曲线完全重复
• 存在冗余代码：提前初始化的Ca、Cb、Cc、Cd数组后续被odeint求解结果完全覆盖，属于无效代码

修正后可运行代码 #

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# 对应可逆反应 A + B <=> C + D
k1 = 0.05   # 正反应速率常数
k2 = 0.05   # 逆反应速率常数
tf = 200
dt = 0.2
t = np.arange(0, tf + 0.01, dt)

def dC(C, tm):
    Ca, Cb, Cc, Cd = C
    # 计算正、逆反应速率
    r_forward = k1 * Ca * Cb
    r_reverse = k2 * Cc * Cd
    # 计算各物质浓度变化率
    dCa_dt = -r_forward + r_reverse
    dCb_dt = -r_forward + r_reverse
    dCc_dt = r_forward - r_reverse
    dCd_dt = r_forward - r_reverse
    return [dCa_dt, dCb_dt, dCc_dt, dCd_dt]

# 传入初始浓度：A=1, B=1, C=0, D=0 求解微分方程
C = odeint(dC, [1, 1, 0, 0], t)
Ca = C[:, 0]
Cb = C[:, 1]
Cc = C[:, 2]
Cd = C[:, 3]

# 绘图输出
plt.plot(t, Ca, 'r--', linewidth=2.0, label='反应物A')
plt.plot(t, Cb, 'k--', linewidth=2.0, label='反应物B')
plt.plot(t, Cc, 'b--', linewidth=2.0, label='产物C')
plt.plot(t, Cd, 'm--', linewidth=2.0, label='产物D')
plt.xlabel('反应时间')
plt.ylabel('物质浓度')
plt.legend()
plt.show()

效果说明 #

修正后代码运行后可以得到符合勒夏特列原理的平衡曲线：反应物A、B浓度从初始值1逐渐下降，产物C、D浓度从0逐渐上升，最终四条曲线达到平稳的平衡状态，和示例图的变化趋势完全一致。你可以调整k1/k2的比值、初始浓度参数，观察平衡状态的变化，验证勒夏特列原理。

内容的提问来源于stack exchange，提问作者Fabio Lima