基于速率模型的乙醇胺(MEA)吸收捕集二氧化碳模拟、教学。
刚接触化工模拟的朋友看到MEA吸收CO₂的文献可能会懵——又是传质方程又是反应动力学，密密麻麻的公式看得人头皮发麻。今天咱们换个姿势，用Python把抽象概念变成看得见的代码，手把手拆解这个经典流程。
先来点干货，典型的MEA吸收塔模型核心就两个微分方程：
def absorption_model(t, y): dCO2_gas = -kL_a*(y[0] - H*y[1]) # 气相浓度变化 dCO2_liq = (kL_a*(y[0] - H*y[1]) - rxn_rate)/L # 液相浓度变化 return [dCO2_gas, dCO2_liq]这段代码藏着三个关键参数：kLa是传质系数，H是亨利常数，rxnrate是反应速率。重点看那个负号——气相浓度降低的量正好等于液相增加的量（当然要考虑反应消耗），这就是质量守恒的代码表达。
反应速率计算更有意思，MEA和CO₂的跳舞可以用下面这个函数表示：
def get_rxn_rate(C_CO2, C_MEA): # 二级反应速率公式 k2 = 0.0175 # m³/(mol·s) 30℃数据 return k2 * C_CO2 * C_MEA这里的k2是个需要实测的参数。注意函数参数用的是实时浓度，说明反应速率在吸收过程中是动态变化的——这也是为啥必须用微分方程求解。
接下来上主程序：
from scipy.integrate import solve_ivp import numpy as np # 操作参数 L = 0.5 # 液层高度(m) H = 0.8 # 亨利常数 kL_a = 0.15 # 总传质系数(1/s) C_MEA0 = 300 # 初始MEA浓度(mol/m³) # 求解区间和初值 t_span = (0, 30) # 30秒模拟 y0 = [120.0, 15.0] # 初始气相和液相CO₂浓度 sol = solve_ivp(absorption_model, t_span, y0, method='BDF')这里选用了BDF求解器，适合处理可能出现的刚性问题。跑完程序后画个图看看趋势：
import matplotlib.pyplot as plt plt.figure(figsize=(10,6)) plt.plot(sol.t, sol.y[0], 'r--', label='Gas phase') plt.plot(sol.t, sol.y[1], 'b-', label='Liquid phase') plt.xlabel('Time (s)') plt.ylabel('CO₂ Concentration (mol/m³)') plt.legend() plt.show()不出意外的话会看到气相浓度快速下降，液相浓度先上升后趋缓——这就是传质推动力减弱的典型特征。有意思的是，如果把kL_a调大，曲线会变得更陡，这提示我们在实际操作中可以通过增加气液接触面积来提高吸收效率。
调试时遇到过坑：当MEA浓度过低时，液相浓度可能出现负值。后来加了约束条件：
dCO2_liq = max((kL_a*(y[0] - H*y[1]) - rxn_rate)/L, 0)这种处理虽然粗暴，但能避免计算结果崩掉。不过更严谨的做法应该是建立代数约束，这就得搬出DAE求解器了，下次再展开说。
最后留个思考题：如果把进料浓度设为随时间变化的正弦波，系统响应会呈现什么特征？试着改改初值条件，说不定能发现有趣的动态行为。
