湍流加速火焰的效果在哪里折弯 — Damköhler 褶皱假说与 bending

湍流火焰速度先线性上升、随后折弯、最终熄火这三个区段背后的物理

山火在风越大时蔓延越快。风速 5 m/s 时,火线跑得比人快。可是继续把风加大,过了某个风速,火反而会把自己吹灭。消防现场把这叫作 blow-off。在实验室燃烧器上测量预混湍流火焰,得到的曲线形状一模一样。

把湍流强度 $u'$ 调高,湍流火焰速度 $S_T$ 起初几乎线性上升。然后折弯。这个变平的过程称为 bending 效应。再往上推就熄火。本文用三套不同的物理来解释这条曲线的三个区段。线性段用 Damköhler 的褶皱假说,折弯与熄火用 Karlovitz 数,而火焰反过来作用于湍流的那部分用 baroclinic 力矩。

湍流为什么让火焰变快

甲烷-空气预混火焰的层流火焰速度(火焰面在静止混合气中自行传播的速度) $S_L$ 在常温常压下约为 0.40 m/s。火焰厚度 $\delta_L$ 大约 0.5 mm。靠这个数字发动机转不起来。要在 2 ms 内横穿 40 mm 缸径,需要 20 m/s,是 $S_L$ 的 50 倍。

填补这个差距的正是湍流。因此需要定义湍流火焰速度 $S_T$ 。它的定义来自反应率的积分。把平均燃料质量分数输运方程沿火焰刷法线方向积分,消去扩散项,剩下的就是这一式。

\rho_u S_T \left(Y_F^u - Y_F^b\right) = -\int_{-\infty}^{+\infty} \overline{\dot{\omega}_F}\ \mathrm{d}n

$\rho_u$ 是未燃气密度, $Y_F^u$ 与 $Y_F^b$ 分别是未燃侧与已燃侧的燃料质量分数, $\overline{\dot{\omega}_F}$ 是平均燃料反应率, $n$ 是平均火焰刷的法向坐标。

值得记住的是,这正是层流火焰所用关系式的直接推广——层流版本里把 $\overline{\dot{\omega}_F}$ 换成瞬时反应率即可。也就是说 $S_T$ 并不是按速度定义的量,而是对单位面积燃烧能力做积分得到的量。这一点解释了后面的全部内容。

错误的直觉 — 湍流把燃料混合得更好?

最常见的说法是这样的:湍流把燃料和空气混合得更好,所以反应更快。对预混火焰来说这是错的。

预混二字已经说明了一切:燃料和空气在点火之前就已经混匀。火焰面前方的气体本来就是均匀混合气,湍流没有什么可再混的了。混合决定速度的是扩散火焰(燃料与氧化剂分别供入的火焰),不是预混火焰。

还有第二个版本:湍流提高了有效扩散系数,所以火焰变快。层流火焰速度按 $S_L \sim \sqrt{D/\tau_c}$ 标度。若把 $D$ 换成湍流扩散系数 $u' l_t$ ,就得到 $S_T/S_L \sim \sqrt{u' l_t / D}$ ,是一个平方根。可实验给出的初始斜率是直线,指数对不上。

两种说法还有一处更深的偏差:它们都假定火焰不动湍流。真实情况恰恰相反。点火之后温度从 300 K 跃升到 2000 K 以上。运动粘度 $\nu$ 大致按 $T^{1.7}$ 增长,因此涨了 20 倍以上。局部 Reynolds 数 $u' l_t / \nu$ 就同比例地下降。这意味着一股充分发展的湍流有可能在点火后重新层流化。湍流在摇晃火焰的同时,火焰也在扼杀湍流。

Damköhler 的回答:关键是面积,不是速度#

Damköhler 给出了另一幅图景。湍流并不加快局部燃烧,它把火焰面折叠起来。

假设只有一条:褶皱火焰面上的每一点在局部仍以 $S_L$ 传播。化学不动,只改几何。于是质量守恒直接给出答案。穿过褶皱后真实面积 $A_T$ 的未燃气,必须等于穿过平均火焰刷投影面积 $\bar{A}$ 的未燃气。

\rho_u S_L A_T = \rho_u S_T \bar{A} \quad\Longrightarrow\quad \Sigma \equiv \frac{A_T}{\bar{A}} = \frac{S_T}{S_L}

$\Sigma$ 是褶皱因子(wrinkling factor),即单位投影面积上的火焰面积。左边是面积比,右边是速度比,两者是同一个数。

褶皱的尺度由湍流决定。尺寸为 $l_t$ 的涡以速度 $u'$ 推动火焰面,而火焰面自身以 $S_L$ 往回推。这两个位移之比决定褶皱的斜率,斜率决定面积增量。Damköhler 的线性律由此而来。

\frac{S_T}{S_L} = 1 + \frac{u'}{S_L}

$u'$ 是未燃气的均方根速度脉动, $S_L$ 是层流火焰速度。 $u' = 0$ 时回到 $S_T = S_L$ 。

这幅图景成立的前提是火焰面比涡更薄。两个无量纲数用来衡量这个条件。

Da = \frac{\tau_t}{\tau_c} = \frac{l_t/u'}{\delta_L/S_L}, \qquad Ka = \frac{\tau_c}{\tau_\eta} = \left(\frac{u'}{S_L}\right)^{3/2}\left(\frac{l_t}{\delta_L}\right)^{-1/2}

$Da$ 是 Damköhler 数(湍流时间 / 化学时间), $Ka$ 是 Karlovitz 数(化学时间 / Kolmogorov 时间)。 $\tau_t$ 是涡的翻转时间, $\tau_c$ 是化学反应时间, $\tau_\eta$ 是最小涡的时间尺度。 $Da \gg 1$ 时火焰反应快于涡的翻转。 $Ka \lt 1$ 时连最小的涡也挤不进火焰厚度。两者同时成立的区域就是 flamelet 区,Damköhler 的线性律只在其中有效。

亲手做出褶皱

在下面的模拟中亲手把褶皱做出来。

u'/S_L3.0l_t/δ_L20

What to watch: the drawn front's measured arc length is forced to equal S_T/S_L, so a faster flame is a literally longer line. Raise u'/S_L and the amber linear-law ghost runs away from the real front — that gap is the bending effect. Lower l_t/δ_L to drive Ka past 25 and the sheet thickens, desaturates, then breaks into holes.

把 $u'/S_L$ 从 0 调到 6,火焰面折得越来越密,但测得的 $\Sigma$ 在 4 附近就不再上升。旁边那条直线型 Damköhler 预测已经到了 7,两条线明显分开。接着把 $l_t/\delta_L$ 降到 2,再把 $u'/S_L$ 推到 12。 $Ka$ 越过 29,标识切换到 quenching, $\Sigma$ 从 3.5 塌到 2.0。

实验显示的折弯 — bending 与熄火#

实验数据画出三个区段:先是近乎直线,然后变平,最后熄火。

变平的原因有三重叠加。其一,褶皱不只是被生成,也在被消灭。火焰面沿自身法向传播,因此不断削平凸起的尖点(cusp)。当生成与消灭平衡时, $\Sigma$ 就饱和了。其二, $Ka$ 超过 1 以后,小涡钻进预热区把火焰加厚,局部 $S_L$ 本身下降。其三,强应变率会局部熄灭火焰小片,在火焰片上打出孔洞。

本文把这三点压缩成一个式子。

\Sigma = \frac{S_T}{S_L} = 1 + q\ \frac{u'/S_L}{1 + (u'/S_L)/\Sigma_\infty}, \qquad q = \exp\left[-\left(\frac{Ka}{Ka_q}\right)^{2}\right]

$\Sigma_\infty$ 是褶皱饱和上限(此处取 6.0), $Ka_q$ 是熄火 Karlovitz 数(此处取 25.0), $q$ 是介于 0 与 1 之间的熄火衰减因子。分母中的饱和项造出 bending, $q$ 造出熄火。

这里必须坦白:这个式子不是推导出来的定律。它是为复现观测到的曲线形状而拟合的唯象形式。 $\Sigma_\infty = 6.0$ 与 $Ka_q = 25.0$ 不是物理常数而是调参值,燃料、压力、燃烧器构型一变它们也跟着变。原始文献说的是同一件事:预测曲线的折弯很困难,预测熄火极限"几乎不可能"。真正从物理里长出来的只有线性段,余下两段至今仍是关联式。

火焰也在改变湍流 — baroclinic 力矩#

气体穿过火焰面时会膨胀。质量流量必须守恒,于是

\rho_u S_L = \rho_b u_b \quad\Longrightarrow\quad \frac{u_b}{S_L} = \frac{\rho_u}{\rho_b} \approx 7

$\rho_b$ 是已燃气密度, $u_b$ 是离开火焰面的气体速度。火焰面是一台把流动加速七倍的装置。这个加速会直接作用到湍流上。

涡量输运方程的四项展示了作用路径。

\frac{\mathrm{D}\boldsymbol{\omega}}{\mathrm{D}t} = \underbrace{(\boldsymbol{\omega}\cdot\nabla)\mathbf{u}}_{\text{I}} - \underbrace{\boldsymbol{\omega}(\nabla\cdot\mathbf{u})}_{\text{II}} + \underbrace{\nu\nabla^{2}\boldsymbol{\omega}}_{\text{III}} + \underbrace{\frac{1}{\rho^{2}}\nabla\rho\times\nabla p}_{\text{IV}}

$\boldsymbol{\omega}$ 是涡量矢量, $\mathbf{u}$ 是速度, $\rho$ 是密度, $p$ 是压力。I 是涡拉伸,II 是膨胀(dilatation),III 是粘性耗散,IV 是 baroclinic 力矩(密度梯度与压力梯度不共线时产生的旋转力矩)。

在冷态湍流中 II 与 IV 都是零。有火焰以后两者都活了。在火焰面上 $\nabla \cdot \mathbf{u} \gt 0$ ,所以 II 会削减涡量。而 IV 则制造新的涡量。横跨火焰面的 $\nabla \rho$ 沿法向而且很大,加速产生的 $\nabla p$ 方向不同,两者的叉积不为零。这就是所谓 "flame-generated turbulence" 的来源。

在下面的模拟中把一个涡扔向火焰面。

u_θ/S_L2.0r/δ_L10ρ_u/ρ_b6.0

What to watch: at the defaults (u_θ/S_L = 2, r/δ_L = 10) the vortex dies on the front and a violet counter-rotating pair is all that is left. Push u_θ/S_L past ~9 and it punches through, dragging a cusp that stays. The arrows get ρ_u/ρ_b times longer on the burnt side — that jump is the dilatation term; set ρ_u/ρ_b = 1 and the baroclinic pair disappears entirely.

把 $u_\theta/S_L$ 放在 2 附近再按 Release,涡会在火焰面上被撕碎,并在原处生出反向旋转的涡量。把同一个值提到 10 附近,涡几乎原封不动地穿过去,只留下一道很深的尖点。把 $\rho_u/\rho_b$ 降到 1, $\nabla\rho$ 消失,baroclinic 生成整体被关掉,这一点同样能看到。

这种反馈带来一个麻烦的副作用: $S_T(u')$ 曲线里说的" $u'$ "到底指哪里的 $u'$ 变得含糊。在未燃气中测、在已燃气中测、在整个火焰刷上测,是三个不同的数。不做密度加权的火焰刷平均,会把未燃与已燃两态之间来回切换的间歇性也算成脉动,数值被抬高。模型很少区分这一点,这正是 RANS 代码在火焰附近出现非物理行为的原因之一。

用 Python 画出 bending 曲线#

把上面的模型原样搬进代码,在两种湍流尺度下并排打印线性预测与 bending 预测。

import math
 
SIGMA_INF = 6.0     # 褶皱饱和上限
KA_QUENCH = 25.0    # 熄火 Karlovitz 数
 
 
def damkohler_number(u_ratio, l_ratio):
    """Da = 湍流时间 / 化学时间。"""
    return l_ratio / u_ratio
 
 
def karlovitz_number(u_ratio, l_ratio):
    """Ka = 化学时间 / Kolmogorov 时间。"""
    return u_ratio**1.5 * l_ratio**-0.5
 
 
def quench_factor(ka):
    """Ka 增大时趋于 0 的熄火衰减因子。"""
    return math.exp(-(ka / KA_QUENCH) ** 2)
 
 
def wrinkling_factor(u_ratio, l_ratio):
    """Sigma = S_T / S_L(含饱和与熄火)。"""
    ka = karlovitz_number(u_ratio, l_ratio)
    return 1.0 + quench_factor(ka) * u_ratio / (1.0 + u_ratio / SIGMA_INF)
 
 
def turbulent_flame_speed(s_laminar, u_ratio, l_ratio):
    """湍流火焰速度 [m/s]。"""
    return s_laminar * wrinkling_factor(u_ratio, l_ratio)
 
 
U_RATIOS = [0.5, 1, 2, 3, 4, 6, 8, 10, 12]
 
for l_ratio in (10.0, 2.0):
    print(f"l_t/delta_L = {l_ratio:.0f}")
    print(f"{'u/SL':>6}{'Da':>9}{'Ka':>9}{'linear':>9}{'bending':>9}")
    for u_ratio in U_RATIOS:
        da = damkohler_number(u_ratio, l_ratio)
        ka = karlovitz_number(u_ratio, l_ratio)
        lin = 1.0 + u_ratio                      # Damkohler 线性律
        ben = wrinkling_factor(u_ratio, l_ratio)  # 饱和 + 熄火
        print(f"{u_ratio:6.1f}{da:9.2f}{ka:9.2f}{lin:9.2f}{ben:9.2f}")
    print()
 
print("S_T at S_L = 0.40 m/s, u'/S_L = 6, l_t/delta_L = 10 :",
      round(turbulent_flame_speed(0.40, 6.0, 10.0), 3), "m/s")

输出如下。

l_t/delta_L = 10
  u/SL       Da       Ka   linear  bending
   0.5    20.00     0.11     1.50     1.46
   1.0    10.00     0.32     2.00     1.86
   2.0     5.00     0.89     3.00     2.50
   3.0     3.33     1.64     4.00     2.99
   4.0     2.50     2.53     5.00     3.38
   6.0     1.67     4.65     7.00     3.90
   8.0     1.25     7.16     9.00     4.16
  10.0     1.00    10.00    11.00     4.20
  12.0     0.83    13.15    13.00     4.03
 
l_t/delta_L = 2
  u/SL       Da       Ka   linear  bending
   0.5     4.00     0.25     1.50     1.46
   1.0     2.00     0.71     2.00     1.86
   2.0     1.00     2.00     3.00     2.49
   3.0     0.67     3.67     4.00     2.96
   4.0     0.50     5.66     5.00     3.28
   6.0     0.33    10.39     7.00     3.52
   8.0     0.25    16.00     9.00     3.28
  10.0     0.20    22.36    11.00     2.68
  12.0     0.17    29.39    13.00     2.00
 
S_T at S_L = 0.40 m/s, u'/S_L = 6, l_t/delta_L = 10 : 1.559 m/s

对比两列就能看到分岔点。到 $u'/S_L = 1$ 为止,线性与 bending 在 15% 以内重合。到了 $u'/S_L = 6$ 已经是 7.00 对 3.90。而在大尺度( $l_t/\delta_L = 10$ )下,即使推到 $u'/S_L = 12$ , $\Sigma$ 仍守在 4 附近。在小尺度( $l_t/\delta_L = 2$ )下,同样的 $u'$ 让 $Ka$ 越过 29, $\Sigma$ 掉到 2.0。湍流强度相同,尺度小的先熄。

下次遇到湍流火焰时

把 $S_T$ 读成面积而不是速度,问题就解决了一半。湍流并没有让化学变快,它只是把以同样 $S_L$ 燃烧的面铺得更宽。实验里 $S_T$ 上不去时,先看火焰面积为什么饱和,再去怀疑反应速率。

不用单一的 $u'$ 下判断,是另一半。 $u'/S_L$ 相同而 $l_t/\delta_L$ 不同, $Ka$ 会差十倍。上面的表里 $u'/S_L = 12$ 时 $\Sigma$ 有时是 4.03,有时是 2.00。区分它们的是湍流尺度而非湍流强度。不给出 $Ka$ 的 $S_T$ 关联式只说了一半。

还有,不要把计算器吐出的 bending 曲线当成定律。线性段来自几何,折弯与熄火是拟合数据得到的。RANS 结果在火焰附近显得奇怪时,先确认 $u'$ 是在哪里测的,比动湍流模型常数要快。