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

直线段只是面积比。弯折和熄火是拟合，皱褶上限取 6，熄火 Karlovitz 数取 25。

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

把湍流强度 $u'$ 调高,湍流火焰速度 $S_T$ 起初几乎线性上升。然后折弯。这个变平的过程称为 bending 效应。再往上推就熄火。直线段是面积。折弯是褶皱的上限,熄火是最小涡的时间。

把滑块 $u'/S_L$ 往上推,褶皱变深,测得的 $\Sigma$ 会离开旁边那条直线预测。

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$ 推到 6,测得的 $\Sigma$ 停在 4 附近,直线则到 7。把 $l_t/\delta_L$ 降到 2 再把 $u'/S_L$ 推到 12,标识切到 quenching, $\Sigma$ 从 3.5 掉到 2.0。

湍流为什么让火焰变快

甲烷-空气预混火焰的层流火焰速度(火焰面在静止混合气中自行传播的速度) $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 的线性律只在其中有效。

实验显示的折弯 — 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$ 不是物理常数而是调参值,燃料、压力、燃烧器构型一变它们也跟着变。原始文献说的是同一件事:预测曲线的折弯很困难,预测熄火极限"几乎不可能"。真正从物理里长出来的只有线性段,余下两段至今仍是关联式。

拿表里的一行用手走一遍。大尺度 $l_t/\delta_L = 10$ 、 $u'/S_L = 6$ 时,Damköhler 数是湍流时间除以化学时间,所以 $10/6 \approx 1.67$ 。Karlovitz 数把同样两个比用另一套指数捆在一起。 $6^{1.5}$ 除以 $\sqrt{10}$ ,就是表上的 4.65。再除以熄火 Karlovitz 数 25,平方之后放进指数, $q \approx 0.966$ 。接近 1。这一格里,熄火因子几乎还没把曲线扳弯。扳弯的是分母。 $u'/S_L$ 除以上限 6 得到 1,分母变成 2,分子 6 缩成 3。 $1 + 0.966 \times 3 = 3.90$ 。线性律写的是 7.00。差出来的 3.1 不是化学变慢,是褶皱面积撞上了上限。

乘上层流速度 0.40 m/s,湍流火焰速度就是输出里的 1.559 m/s。用这个速度穿过 40 mm 缸径,0.040 除以 1.559,大约 0.026 秒,26 ms。前面说发动机要的 2 ms,还差十倍以上。把比值抬到 12,大尺度褶皱因子也只有 4.03,速度停在 $0.40 \times 4.03 \approx 1.61$ m/s。线性式承诺的 13 倍、5.2 m/s,在 bending 那一列里一次都没出现。

只把尺度从 10 降到 2,同样的比值 6 会让 Karlovitz 数从 4.65 跳到 10.39。 $10/2$ 开平方大约是 2.24,两个 $Ka$ 的比正好对上。小涡的时间变短了。能量沿着尺度往下传的级串里,这个时间会先掉到化学时间以下。比值 12、尺度 2 时 $Ka$ 是 29.39。除以 25 已经大于 1, $q$ 掉到大约 0.25,褶皱因子变成 2.00。不是面积又折了一层,是衰减把已经折出来的面积擦掉了。

火焰也在改变湍流 — 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" 的来源。

密度比大约 7,不只是速度的事。未燃侧 $S_L = 0.40$ m/s 时,离开火焰面的速度大约 2.8 m/s。温度从 300 K 跃过 2000 K,运动粘度按 $T^{1.7}$ 涨 20 倍以上,同样的 $u' l_t$ 下 Reynolds 数就按这个倍数下降。II 项在有膨胀的地方抽走涡量,IV 项按 $\nabla\rho$ 和 $\nabla p$ 错开的程度再补回去。

这种反馈带来一个麻烦的副作用: $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_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 项关掉。

对比两列就能看到分岔点。到 $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。湍流强度相同,尺度小的先熄。

丢掉 Ka 只留 u',熄火那一格就没了#

$u'/S_L = 12$ 这一格里, $\Sigma = 4.03$ 和 $2.00$ 同时在。分开它们的是 $l_t/\delta_L$ , $Ka$ 分别是 13.15 和 29.39。线性律两边都写成 13.00。燃烧器上 $S_T$ 上不去时先去拧反应常数,已经封顶的面积会被读成化学太慢。

RANS 在火焰附近显得怪,先看 $u'$ 在哪里取样,再动模型常数。未燃、已燃、火焰刷平均是三个数。动态系数掉成负数的 Germano 平均也是同一类平均问题。 $\Sigma_\infty = 6$ 、 $Ka_q = 25$ 是对这个燃烧器调出来的,不是带到下一种燃料上的常数。