氢气燃烧速率模拟

• Report

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

  • 氧气(反应物的代表)
  • 氢原子(中间物的代表,决定了反应速率)
  • 水(产物)

根据对各步基元反应的分析,可以得出:反应的整体速率,由最慢的一步反应 所决定,因此氢原子浓度至关重要。考虑到氢原子的增殖与终止,可以列出一个一阶线性常微分方程表征氢原子浓度的变化: 其中 表征了氢原子的初始生成速率,$f$ 表示了增殖效应,而 表征了销毁效应。

氢原子浓度分析

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

  1. 分别求 的时候,$c_H$ 随时间变化的曲线。
  2. $f=3$,而 时,$g=2.9$;$t\geq30$ 时,$g=0.1(50-c_H)$。求 随时间变化的曲线。

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

解法

课上建议采用最为简便的 Euler 法数值求解这一常微分方程问题。此处打算对这一做法进行最简便的一种优化:将 Euler 法与后退 Euler 法结合,组成一个预估-校正公式。对一个一般的一阶线性方程边值问题 列出如下所示的迭代格式: 其中 相当于是用 Euler 法对 做的事先预估,再将其代入到后退 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()

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

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

对规律的修正

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

为了模拟 值的下降,考虑以下两种规律:

  1. 值在 秒内从 线性的降低到 $0$;
  2. 值以 的时间常数指数下降,趋向于 。(该常数由尝试得出,可参考以下的结果)

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

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()

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

产物浓度分析

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

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

  • 氢原子浓度按上一节中最后 指数衰减的结果来考虑,不计入中间物与反应物之间的相互影响(用相同的衰变规律表示值这一影响足矣);
  • 氧气浓度按同样的时间常数衰减,设初值为 $10$,末值为 $5$;
  • 常数 设为 $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()

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

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