基于轴向扩散模型拟合RTD数据估算Peclet数的Python实现求助

架构老林

2026-7-15

Python实现闭-闭边界轴向扩散模型拟合RTD数据求Peclet数

核心思路

闭-闭边界条件下的轴向扩散模型，其停留时间分布密度函数$E(t)$的解析解为无限级数形式，我们可以通过截断级数（取前N项即可满足收敛精度）来近似，再用非线性拟合工具将模型与实验数据匹配，求解Peclet数（$Pe = \frac{uL}{D}$，其中$u$为流速，$L$为反应器长度，$D$为轴向扩散系数）。

步骤1：明确模型表达式

闭-闭边界的$E(t)$公式为：
$$
E(t) = \frac{1}{\bar{t}} \sum_{n=1}^{\infty} \frac{4Pe}{(4n^2\pi2 + Pe^2)} \exp\left( -\frac{Pe(1 + \frac{4n^2\pi2}{Pe^2})}{4} \cdot \left( \frac{t}{\bar{t}} - 1 \right)^2 \right)
$$
其中$\bar{t}$为平均停留时间，可通过实验数据积分计算（$\bar{t} = \frac{\int t \cdot C(t) dt}{\int C(t) dt}$），也可由反应器体积和流量直接推导（$\bar{t} = V/Q$）。

步骤2：Python实现方案

导入依赖库

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

定义轴向扩散模型函数

取前10项级数足够满足工程计算精度：

def axial_diffusion_closed_closed(t, Pe, t_mean):
    """
    闭-闭边界轴向扩散模型的E(t)计算
    参数:
        t: 时间数组
        Pe: 待拟合的Peclet数
        t_mean: 平均停留时间
    返回:
        E(t): 停留时间分布密度数组
    """
    E = np.zeros_like(t)
    theta = t / t_mean  # 无量纲时间
    for n in range(1, 11):  # 取前10项级数
        term = (4 * Pe) / (4 * n**2 * np.pi**2 + Pe**2)
        exponent = -Pe * (1 + (4 * n**2 * np.pi**2) / Pe**2) * (theta - 1)**2 / 4
        E += term * np.exp(exponent)
    E /= t_mean  # 归一化到时间维度
    return E

处理实验数据

假设你有实验测得的时间数组t_exp和对应浓度数组C_exp，先将浓度转换为$E(t)$（归一化处理）：

# 替换成你的实际实验数据
t_exp = np.array([10, 20, 30, 40, 50, 60, 70, 80, 90, 100])
C_exp = np.array([0.1, 0.5, 1.2, 1.8, 1.5, 1.0, 0.6, 0.3, 0.1, 0.05])

# 计算E(t)：浓度除以积分面积（总物料量）实现归一化
area = np.trapz(C_exp, t_exp)
E_exp = C_exp / area

# 计算平均停留时间t_mean（也可直接用反应器参数V/Q计算）
t_mean = np.trapz(t_exp * C_exp, t_exp) / area

拟合模型求解Peclet数

使用curve_fit进行非线性拟合，初始值建议用你Excel得到的结果，或经验值（层流Pe≈2，湍流Pe>100）：

# 初始猜测值：Pe设为10，t_mean用计算值
initial_guess = [10, t_mean]

# 执行拟合
params, covariance = curve_fit(axial_diffusion_closed_closed, t_exp, E_exp, p0=initial_guess)

# 提取结果
Pe_fit = params[0]
t_mean_fit = params[1]

print(f"拟合得到的Peclet数: {Pe_fit:.2f}")
print(f"拟合得到的平均停留时间: {t_mean_fit:.2f}")

可视化拟合效果

# 生成拟合曲线的时间点
t_fit = np.linspace(min(t_exp), max(t_exp), 100)
E_fit = axial_diffusion_closed_closed(t_fit, Pe_fit, t_mean_fit)

# 绘图对比
plt.scatter(t_exp, E_exp, label='实验数据', color='red')
plt.plot(t_fit, E_fit, label='拟合曲线', color='blue')
plt.xlabel('时间 t')
plt.ylabel('E(t)')
plt.title('RTD数据拟合轴向扩散模型（闭-闭边界）')
plt.legend()
plt.grid(True)
plt.show()

注意事项

若拟合不收敛，可调整Pe的初始值，或增加级数项数（比如把循环改成range(1,21)）。
确保$E(t)$的积分接近1，否则说明数据归一化有误。
若已知准确的$\bar{t}$（由反应器几何参数计算），可在模型中固定该值，仅拟合Pe，提升拟合稳定性。

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