202 veces más disipación y solo 13% menos temperatura — el límite de extinción de la flamelet

La extinción no es un enfriamiento gradual. La rama ardiente desaparece en el punto donde la pérdida difusiva supera la generación de calor.

La difusión aumentó, ¿por qué se apagó la llama?#

Construir una biblioteca de flamelets obliga a fijar a mano una malla: hasta dónde estirar la disipación escalar $\chi_{st}$ antes de cortar la tabla. Cortar demasiado pronto hace que el solver turbulento consulte fuera de la tabla. Cortar demasiado tarde deja a Newton girando en busca de soluciones que no existen. Para ver esa frontera se resolvieron directamente las ecuaciones de flamelet en una dimensión.

Primero el resultado. Llevar $\chi_{st}$ de 0.05 a 10.12 1/s, un factor de 202, bajó la temperatura pico de 2258.8 K a 1960.9 K: apenas 13.2%. Un paso más y la solución desapareció. La temperatura no tiende lentamente a cero; la rama que sostenía la llama se esfuma entera.

Suena al revés. La difusión es el único medio que lleva los reactivos a la superficie de llama. En una llama no premezclada el combustible y el oxidante llegan por lados opuestos y no tienen otra forma de encontrarse. Y sin embargo, al aumentar esa difusión, el fuego se apaga.

Al subir el control strain rate a la capa de mezcla se adelgaza y $\chi$ crece en el mismo $Z$ .

strain rate a = 60 1/s Z_st = 0.055 

Layer thickness sqrt(2D/a) = 0.82 mm, chi_st = 1.48 1/s against the quench limit 10.12 1/s. Raising a squeezes the layer and the sheet dims; moving Z_st slides the sheet off the chi peak, so the same strain delivers far less dissipation to the reaction zone.

Conviene mirar cómo se apaga la banda de llama. Al apartar Z_st del pico, la deformación se queda y baja solo la disipación que llega a la zona de reacción.

¿Qué hace falta para escribir una llama con un solo escalar?#

Supóngase que todas las especies difunden con el mismo coeficiente. Bajo esa hipótesis el estado de mezcla de una partícula fluida queda fijado por un único escalar conservado: la fracción de mezcla $Z$ . Aquí $Z=0$ es oxidante puro y $Z=1$ combustible puro.

Peters propuso en 1984 una transformación de coordenadas que toma $Z$ como variable independiente. Escritas sobre $Z$ en lugar de tres coordenadas espaciales, solo sobrevive la dirección normal a la llama.

\rho \frac{\partial Y_k}{\partial t} = \rho \frac{\chi}{2} \frac{\partial^2 Y_k}{\partial Z^2} + \dot{\omega}_k

$Y_k$ es la fracción másica de la especie $k$ , $\dot{\omega}_k$ su tasa de producción y $\chi$ la tasa de disipación escalar. La ecuación de energía adopta la misma forma.

\rho \frac{\partial T}{\partial t} = \rho \frac{\chi}{2} \frac{\partial^2 T}{\partial Z^2} - \frac{1}{c_p}\sum_k h_k \dot{\omega}_k

$T$ es la temperatura, $c_p$ el calor específico a presión constante y $h_k$ las entalpías de especie. Ferreira relajó después la hipótesis de $c_p$ constante y propuso una forma más general, que vuelve a esta cuando $c_p$ no varía.

Lo que desapareció es el flujo. No hay campo de velocidad, ni presión, ni malla tridimensional. Todo lo que la turbulencia le hace a la llama quedó comprimido en $\chi$ . El sitio donde un mismo $\chi$ carga dos temperaturas es donde la tabla de flamelets pliega 1759 K junto a 818 K.

¿De qué es exactamente el inverso la disipación escalar?#

La definición es simple.

\chi = 2D\,\left|\nabla Z\right|^2

$D$ es la difusividad. La dimensión es 1/s: se lee como el inverso del tiempo característico de difusión.

Cuando la turbulencia estira una flamelet, el gradiente de $Z$ se empina. Como $\chi$ va con el cuadrado de ese gradiente, crece aún más rápido. Llega más difusión molecular a la superficie de reacción y, al mismo tiempo, el calor que sale de esa superficie crece en la misma proporción. El tiempo de residencia disponible para reaccionar se acorta.

Esa es la respuesta a la paradoja. La difusión trae los reactivos, pero también se lleva el calor. La química depende exponencialmente de la temperatura, así que el lado que pierde calor gana siempre primero.

Para una llama laminar en contraflujo, la dependencia de $\chi$ con $Z$ tiene forma cerrada, la que Pitsch verificó contra DNS.

\chi(Z) = \chi_{st}\,\frac{\exp\!\left(-2\left[\mathrm{erfc}^{-1}(2Z)\right]^2\right)}{\exp\!\left(-2\left[\mathrm{erfc}^{-1}(2Z_{st})\right]^2\right)}

$\chi_{st}$ es el valor en la fracción de mezcla estequiométrica $Z_{st}$ : el único número que indexa la biblioteca. Y $Z_{st}$ sale solo de la composición de las corrientes.

$Z_{st} = \frac{1}{1 + s\,Y_{F,1}/Y_{O,0}}$

$s$ es la masa de oxígeno necesaria para quemar 1 kg de combustible, $Y_{F,1}$ la fracción másica en la corriente de combustible y $Y_{O,0}$ la de oxígeno en la corriente oxidante.

¿Cómo se fija el límite de extinción en código?#

Con una reacción irreversible de un paso $F + s\,O \to (1+s)P$ la incógnita se reduce a una sola función. Como las difusividades son iguales, $Y_F + T/q$ es lineal en $Z$ .

Y_F(Z) + \frac{T(Z)}{q} = \frac{T_0}{q} + Z\left(Y_{F,1} + \frac{T_1 - T_0}{q}\right)

$q$ es el calor liberado por unidad de masa de combustible escrito en unidades de temperatura. $Y_O$ queda atado a $T$ de la misma manera. Queda entonces un solo campo por resolver, $T(Z)$ , y una sola matriz tridiagonal.

La estrategia es continuación. Se arranca con $\chi_{st}$ muy pequeño desde la solución de Burke–Schumann. Se multiplica $\chi_{st}$ por 1.35 en cada paso, se corre Newton y se pasa la solución anterior como estimación inicial. En cuanto la solución ardiente deja de existir, Newton o no converge o cae a temperatura ambiente. Cuarenta bisecciones geométricas entre esos dos valores fijan el límite.

import math
 
N = 201
H = 1.0 / (N - 1)
ZG = [i * H for i in range(N)]
T_OX, T_FU = 300.0, 300.0      # temperatura de entrada en ambos lados [K]
S_RATIO = 4.0                  # kg de oxigeno por kg de combustible
Y_OX = 0.232                   # fraccion masica de oxigeno en el aire
Q_HEAT = 36000.0               # calor de reaccion escrito como temperatura [K]
T_ACT = 24000.0                # temperatura de activacion [K]
A_PRE = 2.0e10                 # factor preexponencial [1/s]
 
def inv_erfc(y):
    lo, hi = -6.0, 6.0
    for _ in range(80):
        mid = 0.5 * (lo + hi)
        if math.erfc(mid) > y: lo = mid
        else: hi = mid
    return 0.5 * (lo + hi)
 
def chi_profile(chi_st, z_st):
    # forma de contraflujo de Pitsch. Se anula en los extremos, asi que se le pone un piso
    f = lambda z: math.exp(-2.0 * inv_erfc(2.0 * z) ** 2)
    fst = f(z_st)
    out = []
    for z in ZG:
        v = f(min(max(z, 1e-9), 1 - 1e-9)) / fst
        out.append(chi_st * max(v, 1e-4))
    return out
 
def couple_species(z, temp, y_fu):
    # funciones de acoplamiento: Y_F y Y_O se recuperan solo con T y Z
    yf = (T_OX / Q_HEAT + z * (y_fu + T_FU / Q_HEAT - T_OX / Q_HEAT)) - temp / Q_HEAT
    yo = (Y_OX + S_RATIO * T_OX / Q_HEAT
          + z * (S_RATIO * T_FU / Q_HEAT - Y_OX - S_RATIO * T_OX / Q_HEAT)) - S_RATIO * temp / Q_HEAT
    return max(yf, 0.0), max(yo, 0.0)
 
def tridiag_sweep(a, b, c, d):
    n = len(b)
    cp = [0.0] * n; dp = [0.0] * n
    cp[0] = c[0] / b[0]; dp[0] = d[0] / b[0]
    for i in range(1, n):
        m = b[i] - a[i] * cp[i - 1]
        cp[i] = c[i] / m
        dp[i] = (d[i] - a[i] * dp[i - 1]) / m
    x = [0.0] * n
    x[-1] = dp[-1]
    for i in range(n - 2, -1, -1):
        x[i] = dp[i] - cp[i] * x[i + 1]
    return x
 
def newton_flamelet(chi, guess, y_fu, itmax=60):
    T = list(guess)
    n = N - 2
    for _ in range(itmax):
        a = [0.0] * n; b = [0.0] * n; c = [0.0] * n; r = [0.0] * n
        for k in range(n):
            i = k + 1
            yf, yo = couple_species(ZG[i], T[i], y_fu)
            ex = math.exp(-T_ACT / T[i])
            w = A_PRE * yf * yo * ex
            dw = A_PRE * ex * (-yo / Q_HEAT - yf * S_RATIO / Q_HEAT + yf * yo * T_ACT / T[i] ** 2)
            if yf <= 0.0 or yo <= 0.0: w = 0.0; dw = 0.0
            cf = chi[i] / (2.0 * H * H)
            r[k] = -(cf * (T[i + 1] - 2 * T[i] + T[i - 1]) + Q_HEAT * w)
            a[k] = cf if k > 0 else 0.0
            b[k] = -2.0 * cf + Q_HEAT * dw
            c[k] = cf if k < n - 1 else 0.0
        d = tridiag_sweep(a, b, c, r)
        lim = max(abs(v) for v in d)
        scale = min(1.0, 200.0 / lim) if lim > 200.0 else 1.0   # limita la actualizacion de temperatura
        for k in range(n):
            T[k + 1] = max(T[k + 1] + scale * d[k], 280.0)
        if lim < 1e-6:
            return T, True
    return T, False
 
def burke_schumann(y_fu):
    # limite de velocidad infinita: el punto de partida de la continuacion
    z_st = 1.0 / (1.0 + S_RATIO * y_fu / Y_OX)
    T = []
    for z in ZG:
        if z <= z_st: T.append(T_OX + (T_FU - T_OX) * z + Q_HEAT * y_fu * z)
        else: T.append(T_OX + (T_FU - T_OX) * z + Q_HEAT * y_fu * z_st * (1 - z) / (1 - z_st))
    return z_st, T
 
def extinction_scan(y_fu, chi0=0.05, growth=1.35, steps=60):
    z_st, T = burke_schumann(y_fu)
    chi_st = chi0
    hist = []
    last_ok = None
    while steps > 0:
        chi = chi_profile(chi_st, z_st)
        Tn, ok = newton_flamelet(chi, T, y_fu)
        tmax = max(Tn)
        if not ok or tmax < 1000.0:
            break
        T = Tn; last_ok = (chi_st, tmax, list(Tn))
        hist.append((chi_st, tmax))
        chi_st *= growth; steps -= 1
    lo, hi = last_ok[0], chi_st
    base = last_ok[2]
    for _ in range(40):                      # ajusta la ultima solucion ardiente con una biseccion geometrica
        mid = math.sqrt(lo * hi)
        Tn, ok = newton_flamelet(chi_profile(mid, z_st), base, y_fu)
        if ok and max(Tn) > 1000.0: lo, base = mid, Tn
        else: hi = mid
    return z_st, hist, lo, base
 
if __name__ == "__main__":
    z_st, hist, chi_q, Tq = extinction_scan(1.0)
    print("Z_st = %.4f   T_ad(Burke-Schumann) = %.1f K" % (z_st, max(burke_schumann(1.0)[1])))
    print("  chi_st [1/s]    T_max [K]   drop [K]")
    tad = hist[0][1]
    for c, t in hist[::3]:
        print("    %9.4f    %8.1f    %7.1f" % (c, t, tad - t))
    print("    %9.4f    %8.1f    %7.1f   <- quench" % (chi_q, max(Tq), tad - max(Tq)))
    print("chi range %.4f -> %.4f = x%.1f,  T_max drop = %.1f %%"
          % (hist[0][0], chi_q, chi_q / hist[0][0], 100.0 * (tad - max(Tq)) / tad))
    print()
    print("  Y_F1     Z_st    T_ad [K]   chi_q [1/s]")
    for yf in (1.0, 0.5, 0.25, 0.125):
        zs, h, cq, _ = extinction_scan(yf)
        print("   %5.3f   %6.4f    %7.1f    %9.3f" % (yf, zs, max(burke_schumann(yf)[1]), cq))

Z_st = 0.0548   T_ad(Burke-Schumann) = 2273.2 K
  chi_st [1/s]    T_max [K]   drop [K]
       0.0500      2258.8        0.0
       0.1230      2251.1        7.7
       0.3027      2245.4       13.4
       0.7447      2233.0       25.8
       1.8322      2210.9       47.9
       4.5079      2164.0       94.8
      10.1232      1960.9      297.9   <- quench
chi range 0.0500 -> 10.1232 = x202.5,  T_max drop = 13.2 %
 
  Y_F1     Z_st    T_ad [K]   chi_q [1/s]
   1.000   0.0548     2273.2       10.123
   0.500   0.1039     2168.8        8.662
   0.250   0.1883     1991.3        3.942
   0.125   0.3169     1719.8        0.561

Al girar el dial de $\chi_{st}$ hacia la derecha, el punto de la rama naranja se mueve con la disipación.

chi_st = 0.500 1/s fuel Y_F1 = 1.000 

Z_st = 0.055, T_ad = 2274 K, quench limit chi_q = 10.20 1/s, current T_max = 2273 K (burning). Push chi_st past the red line and the dot falls off the orange branch. Pull it back down — it stays on the floor, because the frozen branch is also a steady solution. Only re-ignite lifts it.

Al llevar el punto hacia la derecha sobre la rama naranja y cruzar la línea roja, cae al suelo. Y al tirar de $\chi_{st}$ otra vez hacia la izquierda, no vuelve a subir: la solución congelada también es una solución estacionaria. La histéresis es la esencia de esta curva.

Vale la pena mirar lo planas que son las seis primeras filas. Con $\chi_{st}$ multiplicado por 90, la temperatura cedió 94.8 K. El último paso se lleva de golpe los 203 K restantes. Esa es la rodilla de la curva en S.

Con $s=4$ y $Y_{O,0}=0.232$ dentro de $Z_{st}=1/(1+s Y_{F,1}/Y_{O,0})$ , un $Y_{F,1}=1$ deja el denominador en $1+4/0.232$ y el recíproco es el $Z_{st}=0.0548$ impreso. Con $Y_{F,1}=0.125$ la misma fórmula devuelve el 0.3169 de la segunda tabla. Burke–Schumann a 2273.2 K y la solución con $\chi_{st}=0.05$ a 2258.8 K ya distan 14.4 K. Ese hueco separa el límite de velocidad infinita de la primera solución ardiente.

De $\chi_{st}=0.0500$ a 4.5079 la temperatura solo baja de 2258.8 K a 2164.0 K, los 94.8 K ya dichos. El paso a 10.1232 se lleva el resto, hasta 1960.9 K, y la caída de 297.9 K de la tabla es la suma de ambos tramos. En 0.1230, 0.3027, 0.7447 y 1.8322 las caídas son 7.7 K, 13.4 K, 25.8 K y 47.9 K. Casi una recta hasta la rodilla. Multiplicar por 1.35 cada vez cruza ese último tramo con paso grueso. Cuarenta bisecciones geométricas son lo que fija $\chi_q=10.123$ .

¿Por qué diluir el combustible rompe el límite con mucha más fuerza?#

La segunda tabla de la salida lo muestra. Diluir la fracción másica de combustible de 1.0 a 0.125, un factor de 8, bajó la temperatura adiabática de llama de 2273.2 K a 1719.8 K, un 24%. El límite de extinción $\chi_q$ pasó de 10.123 a 0.561: una caída de 18 veces.

Una caída del 24% en temperatura compró una caída del 1800% en el límite. El término de Arrhenius lo explica. La velocidad de reacción cuelga de $\exp(-T_a/T)$ , con $T_a = 24000$ K en este caso. Mover $T$ de 2273 a 1720 mueve el exponente de $-10.6$ a $-14.0$ . Una diferencia de 3.4 en el exponente son 30 veces en el producto.

El desplazamiento de $Z_{st}$ de 0.055 a 0.317 empuja en el mismo sentido. Se ve moviendo Z_st hacia la derecha en la primera simulación. El pico de $\chi(Z)$ está cerca de $Z=0.5$ , así que un $Z_{st}$ mayor acerca la zona de reacción a ese pico. La misma deformación la golpea ahora con más disipación.

La recirculación de gases de escape (EGR) es donde esta combinación aparece en la práctica. Al diluir combustible u oxidante con un inerte, $T_{ad}$ baja con suavidad mientras el límite de extinción se desploma. Primero aparece la extinción local, luego el despegue de llama (lift-off) y al final el soplado (blow-off).

Consultar el lado derecho de la rodilla devuelve la temperatura de entrada#

Llenar la malla de $\chi_{st}$ hasta $\chi_q=10.123$ deja la tabla densa, pero el paso de 4.5079 a 10.1232 baja la temperatura de 2164.0 K a 1960.9 K, 203 K en un solo tramo. Cuántos puntos se gastan en esa rodilla es el trueque real. Una celda con $\chi_{st} > \chi_q$ no es una consulta. El modelo elige un estado congelado, una probabilidad de extinción o una flamelet no estacionaria.

$\chi_{st}$ sale del modelo de turbulencia. Un índice desviado por un factor de dos invierte el lado de la rodilla. Un cálculo que cae a la temperatura de entrada se revisa en la entrada de $\chi_{st}$ antes que en el mecanismo. Un barrido tridiagonal de esa tabla cuesta el mismo tipo de trabajo que la eliminación de Thomas para difusión implícita.