目录
-
氢原子浓度分析
- 解法
- 程序源码与结果
- 对规律的修正
- 产物浓度分析
这是一次燃烧学课程作业的内容。
在氢气燃烧过程中，有三种成分的浓度值得关注：
- 氧气（反应物的代表）
- 氢原子（中间物的代表，决定了反应速率）
- 水（产物）
根据对各步基元反应的分析，可以得出：反应的整体速率，由最慢的一步反应
H+O2→OH+O
所决定，因此氢原子浓度至关重要。考虑到氢原子的增殖与终止，可以列出一个一阶线性常微分方程表征氢原子浓度的变化：
dcHdτ=wH+fcH−gcH
其中 wH 表征了氢原子的初始生成速率，f 表示了增殖效应，而 g 表征了销毁效应。
氢原子浓度分析
本问题提法如下：已知 ωH=0.01，f=3。
- 分别求 g=2.9,3,3.1 的时候，cH 随时间变化的曲线。
- f=3，而 t<30 时，g=2.9；t≥30 时，g=0.1(50−cH)。求 cH 随时间变化的曲线。
此处将两部分合而为一，将总计 4 条曲线绘制在同一幅图中。
解法
课上建议采用最为简便的 Euler 法数值求解这一常微分方程问题。此处打算对这一做法进行最简便的一种优化：将 Euler 法与后退 Euler 法结合，组成一个预估-校正公式。对一个一般的一阶线性方程边值问题
{y′(x)=f(x,y)y(x0)=y0
列出如下所示的迭代格式：
{pi=yi+f(xi,yi)⋅Δxyi+1=yi+f(xi+1,pi)⋅Δx
其中 pi 相当于是用 Euler 法对 yi+1 做的事先预估，再将其代入到后退 Euler 法中以避免隐式方程的求解。
以下便遵循这个步骤来求解问题。
程序源码与结果
-
import numpy as np
-
import 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: (t, c)
-
else: (t, c, linewidth=1, linestyle='--')
-
-
-
('$c_H$ versus $t$')
-
('$t$')
-
('$c_H$')
-
(['$g=2.9$', '$g=3.0$', '$g=3.1$', '$g=0.1(50-c_H)$'])
-
()
从结果可以看出，当 g=2.9 时生成速率较高，氢原子浓度呈指数上升；当 g=3.0 时增殖与销毁持平，氢原子仅按初始生成速率线性生成；当 g=3.1 时销毁速率高于增殖速率，氢原子浓度几乎没有变化。
特别地，若服从问题 (2) 中提出的销毁因子变化规律，氢原子的浓度将在 g 的突变点处发生转折，迅速降至接近于 0 的水平。推敲此种规律的提出背景，可能来自于反应物的输入突然停止、反应温度骤降等诸多可能。作为体现总包反应速率的中间产物，氢原子浓度的骤降直接意味着反应的终结。
对规律的修正
仔细推敲氢原子从生成、增殖到全部销毁（尽管从数学上与实在上来说，其都不可能全部销毁）的过程，可以看出：在反应物没有持续供应的情况下，导致氢原子最终散尽的不是 g 值的上升，而应该是 f 值的下降——更进一步的说，是反应物的耗尽，导致 f 从 3 很快的掉到 1 左右——与消耗持平，变为不分支的链反应。
为了模拟 f 值的下降，考虑以下两种规律：
- f 值在 10 秒内从 3 线性的降低到 0；
- f 值以 τ=20e 的时间常数指数下降，趋向于 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]
-
(t, c)
-
-
# 指数下降规律
-
init()
-
tau = 20 /
-
for i in range(1000):
-
t = t + [t[i] + dt]
-
f = 3 * (- t[i] / tau)
-
p = c[i] + (w + (f-g) * c[i]) * dt
-
c = c + [c[i] + (w + (f - g) * p) * dt]
-
(t, c)
-
-
('$c_H$ versus $t$ (with $f$ decreasing)')
-
('$t$')
-
('$c_H$')
-
(['$f=0.3(10-t)$', '$f=3\exp(-t\cdot\mathrm{e}/20)$'])
-
()
从绘图结果可以看出，经过精心取定的两个规律大体上产生了相近的结果：在 t=10 以内的时间之中，氢原子浓度很快地攀升，又很快的由于反应物的耗尽而回落到接近于 0 的水平。
产物浓度分析
以上仅分析了中间物——氢原子的浓度变化规律。除此以外，产物——水的浓度变化，也是非常重要的一个指标，值得研究。
根据对各基元反应情况的分析，水的生成速率可按
dcH2Odt=2kcHcO2
计算，其中反应常数 k 在等温条件下维持不变；而 cO2 为反应物之一的消耗情况，若我们认为反应时氧气过量，且反应物均按指数规律衰减——与之前氢原子浓度分析时的最后一种情景吻合，则其将很快以指数规律衰减到一个正数（将氢气耗尽后的剩余部分）cO2,res。以下不妨对反应做这样的假设：
- 氢原子浓度按上一节中最后 f 按 τ=20/e 指数衰减的结果来考虑，不计入中间物与反应物之间的相互影响（用相同的衰变规律表示值这一影响足矣）；
- 氧气浓度按同样的时间常数衰减，设初值为 10，末值为 5；
- 常数 2k=K 设为 0.002。
接下来据以上假设进行模拟——这种模拟只能推算出产物浓度的总体变化规律，具体的数值并无参考价值。
-
# 部分参数直接用上面已有的结果
-
-
c_w = [0]
-
for i in range(1000):
-
c_o = 5 + 5 * (-t[i] / tau)
-
c_w = c_w + [c_w[i] + .002 * c[i] * c_o]
-
-
('Concentration of water over time')
-
('$t$')
-
('$c_{H_2O}$')
-
(t, c_w)
-
()
可以看到，在这种假设下，水的浓度呈现「由慢转快——稳定生成——由快转慢——停止生成」的增长速率变化规律，这与我们的预期相符合。以上也可以说明，之前对氢原子浓度变化的假设基本合理，应该与某种情形之下的反应条件相吻合——只是具体的细节需要修正与改进。