氢气燃烧速率模拟 · Heishan Press

在氢气燃烧过程中，有三种成分的浓度值得关注：

氧气（反应物的代表）
氢原子（中间物的代表，决定了反应速率）
水（产物）

根据对各步基元反应的分析，可以得出：反应的整体速率，由最慢的一步反应 $H+O_2\to OH+O$ 所决定，因此氢原子浓度至关重要。考虑到氢原子的增殖与终止，可以列出一个一阶线性常微分方程表征氢原子浓度的变化： $\frac{\mathrm{d}c_H}{\mathrm{d}\tau}=w_H+fc_H-gc_H$ 其中 $w_H$ 表征了氢原子的初始生成速率，$f$ 表示了增殖效应，而 $g$ 表征了销毁效应。

氢原子浓度分析

本问题提法如下：已知 $\omega_H=0.01$，$f=3$。

分别求 $g=2.9, 3, 3.1$ 的时候，$c_H$ 随时间变化的曲线。
$f=3$，而 $t<30$ 时，$g=2.9$；$t\geq30$ 时，$g=0.1(50-c_H)$。求 $c_H$ 随时间变化的曲线。

此处将两部分合而为一，将总计 $4$ 条曲线绘制在同一幅图中。

解法

课上建议采用最为简便的 Euler 法数值求解这一常微分方程问题。此处打算对这一做法进行最简便的一种优化：将 Euler 法与后退 Euler 法结合，组成一个预估-校正公式。对一个一般的一阶线性方程边值问题 $\begin{cases} y'(x)=f(x,y)\\y(x_0)=y_0 \end{cases}$ 列出如下所示的迭代格式： $\begin{cases} p_i=y_i+f(x_i,y_i)\cdot\Delta x\\ y_{i+1}=y_i+f(x_{i+1},p_i)\cdot\Delta x \end{cases}$ 其中 $p_i$ 相当于是用 Euler 法对 $y_{i+1}$ 做的事先预估，再将其代入到后退 Euler 法中以避免隐式方程的求解。

以下便遵循这个步骤来求解问题。

程序源码与结果

import numpy as np
import matplotlib.pyplot as plt

以上引入了必要的 Numpy 库与 Matplotlib 库。

def init():
    '初始化各基本参数'
    global t, dt, c, w
    t = [0]
    dt = 0.05
    c = [0]
    w = .01

以上的 init 函数用于在各次模拟之前初始化计算所用的变量与参数。

f = 3  # 增殖因子
g_list = [2.9,3,3.1]  # 销毁因子

for num in range(4):
    init()
    for i in range(750):
        t = t + [t[i] + dt]
        if num == 3:
            if t[i] < 30: g = 2.9
            else: g = 0.1 * (50 - c[i])
        else: g = g_list[num]
        p = c[i] + (0.01 + (f - g) * c[i]) * dt
        c = c + [c[i] + (0.01 + (f - g) * p) * dt]
    if num == 3: plt.plot(t, c)
    else: plt.plot(t, c, linewidth=1, linestyle='--')


plt.title('$c_H$ versus $t$')
plt.xlabel('$t$')
plt.ylabel('$c_H$')
plt.legend(['$g=2.9$', '$g=3.0$', '$g=3.1$', '$g=0.1(50-c_H)$'])
plt.show()

从结果可以看出，当 $g=2.9$ 时生成速率较高，氢原子浓度呈指数上升；当 $g=3.0$ 时增殖与销毁持平，氢原子仅按初始生成速率线性生成；当 $g=3.1$ 时销毁速率高于增殖速率，氢原子浓度几乎没有变化。

特别地，若服从问题 (2) 中提出的销毁因子变化规律，氢原子的浓度将在 $g$ 的突变点处发生转折，迅速降至接近于 $0$ 的水平。推敲此种规律的提出背景，可能来自于反应物的输入突然停止、反应温度骤降等诸多可能。作为体现总包反应速率的中间产物，氢原子浓度的骤降直接意味着反应的终结。

对规律的修正

仔细推敲氢原子从生成、增殖到全部销毁（尽管从数学上与实在上来说，其都不可能全部销毁）的过程，可以看出：在反应物没有持续供应的情况下，导致氢原子最终散尽的不是 $g$ 值的上升，而应该是 $f$ 值的下降——更进一步的说，是反应物的耗尽，导致 $f$ 从 $3$ 很快的掉到 $1$ 左右——与消耗持平，变为不分支的链反应。

为了模拟 $f$ 值的下降，考虑以下两种规律：

$f$ 值在 $10$ 秒内从 $3$ 线性的降低到 $0$；
$f$ 值以 $\tau=\frac{20}{\mathrm{e}}$ 的时间常数指数下降，趋向于 $0$ 。（该常数由尝试得出，可参考以下的结果）

以下对这两种过程分别进行模拟。首先是线性下降：

g = 1 # 销毁因子固定为 1

# 线性下降规律
init()
for i in range(1000):
    t = t + [t[i] + dt]
    if t[i] <= 10: f = 3 - .3 * t[i]
    else: f = 0
    p = c[i] + (w + (f - g) * c[i]) * dt
    c = c + [c[i] + (w + (f - g) * p) * dt]
plt.plot(t, c)

# 指数下降规律
init()
tau = 20 / np.e
for i in range(1000):
    t = t + [t[i] + dt]
    f = 3 * np.exp(- t[i] / tau)
    p = c[i] + (w + (f-g) * c[i]) * dt
    c = c + [c[i] + (w + (f - g) * p) * dt]
plt.plot(t, c)

plt.title('$c_H$ versus $$t$$ (with $$f$$ decreasing)')
plt.xlabel('$t$')
plt.ylabel('$c_H$')
plt.legend(['$f=0.3(10-t)$', '$f=3\exp(-t\cdot\mathrm{e}/20)$'])
plt.show()

从绘图结果可以看出，经过精心取定的两个规律大体上产生了相近的结果：在 $t=10$ 以内的时间之中，氢原子浓度很快地攀升，又很快的由于反应物的耗尽而回落到接近于 $0$ 的水平。

产物浓度分析

以上仅分析了中间物——氢原子的浓度变化规律。除此以外，产物——水的浓度变化，也是非常重要的一个指标，值得研究。

根据对各基元反应情况的分析，水的生成速率可按 $\frac{\mathrm{d}c_{H_2O}}{\mathrm{d}t}=2kc_Hc_{O_2}$ 计算，其中反应常数 $k$ 在等温条件下维持不变；而 $c_{O_2}$ 为反应物之一的消耗情况，若我们认为反应时氧气过量，且反应物均按指数规律衰减——与之前氢原子浓度分析时的最后一种情景吻合，则其将很快以指数规律衰减到一个正数（将氢气耗尽后的剩余部分）$c_{O_2,\text{res}}$。以下不妨对反应做这样的假设：

氢原子浓度按上一节中最后 $f$ 按 $\tau=20/\mathrm{e}$ 指数衰减的结果来考虑，不计入中间物与反应物之间的相互影响（用相同的衰变规律表示值这一影响足矣）；
氧气浓度按同样的时间常数衰减，设初值为 $10$，末值为 $5$；
常数 $2k=K$ 设为 $0.002$。

接下来据以上假设进行模拟——这种模拟只能推算出产物浓度的总体变化规律，具体的数值并无参考价值。

# 部分参数直接用上面已有的结果

c_w = [0]
for i in range(1000):
    c_o = 5 + 5 * np.exp(-t[i] / tau)
    c_w = c_w + [c_w[i] + .002 * c[i] * c_o]

plt.title('Concentration of water over time')
plt.xlabel('$t$')
plt.ylabel('$c_{H_2O}$')
plt.plot(t, c_w)
plt.show()

可以看到，在这种假设下，水的浓度呈现「由慢转快——稳定生成——由快转慢——停止生成」的增长速率变化规律，这与我们的预期相符合。以上也可以说明，之前对氢原子浓度变化的假设基本合理，应该与某种情形之下的反应条件相吻合——只是具体的细节需要修正与改进。

这是一次燃烧学课程作业的内容。