Una pequeña fluctuación de la fracción de mezcla baja 964 K la temperatura media — β-PDF presunta y el infinito en Z=0

¿Por qué falla leer la temperatura de llama en la fracción de mezcla media?#

Considérese una llama turbulenta de chorro de metano y aire. En un punto, la fracción de mezcla promediada en el tiempo es justo la estequiométrica, 0.055. ¿Cuál es la temperatura media allí? La temperatura adiabática de llama estequiométrica es 2230 K, así que la respuesta tentadora es 2230 K.

La respuesta es 1266 K, con fluctuaciones moderadas. Este artículo muestra adónde van esos 964 K. Se integra a mano la β-PDF presunta (suponer que la densidad de probabilidad tiene forma de distribución beta). Al final se sigue por qué una malla uniforme de 100 celdas se equivoca en 115 K.

El punto de partida son unas diapositivas de un seminario de combustión del verano de 2018 sobre el uso de FlameMaster. Entre la instalación, ScanMan y los comandos de FlameMaster aparecen algunas diapositivas de teoría: PDF y CDF, la función Beta, fracción de mezcla y tasa de disipación, la curva S. Solo quedan los títulos. Aquí se completa el hueco de "PDF y función Beta".

Z se encarga de la mezcla, T(Z) de la química#

La fracción de mezcla $Z$ (fracción de masa proveniente de la corriente de combustible) vale 0 en el aire y 1 en el combustible. Con química infinitamente rápida, la temperatura depende solo de $Z$. Es la solución de Burke–Schumann: un tejado con el pico en la estequiométrica $Z_{st}$ y rectas que bajan a ambos lados.

En un punto dentro de la turbulencia, $Z$ fluctúa en el tiempo. Un paquete de aire la lleva hacia 0. Un paquete de combustible la dispara hacia arriba. La temperatura media es el promedio de $T(Z)$ sobre todas esas fluctuaciones.

$\tilde T = \int_0^1 T(Z)\,P(Z)\,dZ \neq T(\tilde Z)$

$P(Z)$ es la densidad de probabilidad de $Z$ en ese punto y $\tilde Z$ es la fracción de mezcla media.

Si $T(Z)$ fuera una recta, ambos valores coincidirían. Pero el tejado se dobla en el pico. Con la media justo sobre el quiebre, cualquier fluctuación, hacia cualquier lado, baja la temperatura. Conviene manipularlo directamente en la simulación de abajo.

Cada partícula es un $Z$ muestreado de la β-PDF. La altura a la que cae sobre el tejado es su temperatura. El punto rosa es la media acumulada y pronto se pega a la línea verde $\tilde T$. Con $s=0.1$, un tercio de las partículas se amontona en $Z<0.01$. Al mover $\bar Z$ a 0.2, la pendiente derecha del tejado es suave y la brecha se reduce mucho.

β-PDF: dos números fijan la forma#

Un solver LES o RANS entrega solo dos cantidades: $\tilde Z$ y la varianza $\widetilde{Z''^2}$. Nadie resuelve la PDF completa. Por eso la forma se fija de antemano. Para una variable confinada entre 0 y 1, la distribución beta es la elección natural.

$P(Z) = \frac{Z^{a-1}(1-Z)^{b-1}}{B(a,b)}, \qquad a = \tilde Z\,\gamma,\quad b = (1-\tilde Z)\,\gamma,\quad \gamma = \frac{1}{s} - 1$

$B(a,b)$ es la función beta (constante de normalización). $s = \widetilde{Z''^2}/[\tilde Z(1-\tilde Z)]$ es la intensidad de segregación: 0 es mezcla perfecta y 1 es segregación total en 0 y 1 puros.

Cerca de la estequiométrica, $\tilde Z$ es pequeña. Así, $a$ cae por debajo de 1 enseguida. Con $\tilde Z = 0.055$, $a<1$ para todo $s > 0.052$. Entonces $Z^{a-1}$ diverge cuando $Z\to0$. La integral sigue siendo finita. Pero el valor de la PDF diverge en el extremo. Esa propiedad es la que arruina después la integral en malla uniforme.

Una recta quebrada se integra exactamente con la función beta incompleta#

Si $T(Z)$ es lineal a trozos, no hace falta cuadratura numérica. Cada tramo necesita solo dos integrales.

$\int_0^{x} P\,dZ = I_x(a,b), \qquad \int_0^{x} Z\,P\,dZ = \tilde Z\, I_x(a+1,b)$

$I_x(a,b)$ es la función beta incompleta regularizada (la CDF de la distribución beta). La segunda identidad vale porque $Z\cdot Z^{a-1} = Z^{a}$ y $B(a+1,b)/B(a,b) = a/(a+b) = \tilde Z$.

Al introducir las rectas pobre y rica a cada lado de $Z_{st}$, $\tilde T$ sale en forma cerrada. El infinito del extremo se trata analíticamente dentro de $I_x$. $I_x$ se evalúa con una fracción continua (método de Lentz), el mismo algoritmo que betai de Numerical Recipes.

Lo mismo ocurre con una biblioteca de flamelets. Las curvas $T(Z;\chi)$ de la tasa de disipación y el límite de extinción reemplazan al tejado. Si la malla en $Z$ de la tabla se interpola linealmente, basta aplicar la fórmula tramo a tramo.

Código: varianza de 0.01 a 0.8 en Z̄ estequiométrica#

Se fija $\tilde Z = Z_{st} = 0.055$ y solo se varía $s$. Junto a la integral exacta, se corre una integral en malla uniforme que evalúa la PDF en los centros de celda y renormaliza la suma a 1, con $N=100$ y $1000$. Las dos últimas columnas son diferencias respecto al valor exacto.

import math
 
# Temperatura de llama metano-aire de Burke-Schumann T(Z)
T_U, T_AD, Z_ST = 300.0, 2230.0, 0.055
 
def t_burke(z):
    if z <= Z_ST:
        return T_U + (T_AD - T_U) * z / Z_ST
    return T_U + (T_AD - T_U) * (1.0 - z) / (1.0 - Z_ST)
 
def betacf(a, b, x):
    # fracción continua para la beta incompleta (Lentz)
    qab, qap, qam = a + b, a + 1.0, a - 1.0
    c, d = 1.0, 1.0 - qab * x / qap
    d = 1.0 / (d if abs(d) > 1e-300 else 1e-300)
    h = d
    for m in range(1, 300):
        m2 = 2 * m
        aa = m * (b - m) * x / ((qam + m2) * (a + m2))
        d = 1.0 + aa * d; c = 1.0 + aa / c
        d = 1.0 / (d if abs(d) > 1e-300 else 1e-300); h *= d * c
        aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))
        d = 1.0 + aa * d; c = 1.0 + aa / c
        d = 1.0 / (d if abs(d) > 1e-300 else 1e-300)
        dl = d * c; h *= dl
        if abs(dl - 1.0) < 1e-14:
            break
    return h
 
def inc_beta(a, b, x):
    # beta incompleta regularizada I_x(a, b)
    if x <= 0.0: return 0.0
    if x >= 1.0: return 1.0
    lnf = (math.lgamma(a + b) - math.lgamma(a) - math.lgamma(b)
           + a * math.log(x) + b * math.log(1.0 - x))
    if x < (a + 1.0) / (a + b + 2.0):
        return math.exp(lnf) * betacf(a, b, x) / a
    return 1.0 - math.exp(lnf) * betacf(b, a, 1.0 - x) / b
 
def beta_ab(zm, s):
    # s = Z''^2 / (Zm (1 - Zm)), intensidad de segregación
    g = 1.0 / s - 1.0
    return zm * g, (1.0 - zm) * g
 
def mean_t_exact(zm, s):
    # T es lineal a trozos: cada tramo se integra con I_x
    a, b = beta_ab(zm, s)
    p_lo = inc_beta(a, b, Z_ST)                 # P(Z < Z_st)
    z_lo = zm * inc_beta(a + 1.0, b, Z_ST)      # int_0^Zst Z P dZ
    k1 = (T_AD - T_U) / Z_ST
    k2 = (T_AD - T_U) / (1.0 - Z_ST)
    lean = T_U * p_lo + k1 * z_lo
    rich = (T_U + k2) * (1.0 - p_lo) - k2 * (zm - z_lo)
    return lean + rich
 
def mean_t_grid(zm, s, n):
    # ingenuo: pdf en n centros de celda, renormalizada
    a, b = beta_ab(zm, s)
    num = den = 0.0
    for i in range(n):
        z = (i + 0.5) / n
        p = z ** (a - 1.0) * (1.0 - z) ** (b - 1.0)
        num += p * t_burke(z); den += p
    return num / den
 
zm = Z_ST
print(f"Z_mean = {zm}, T(Z_mean) = {t_burke(zm):.1f} K")
print(" s      a       b      T_exact   gap    grid100  grid1000")
for s in [0.01, 0.03, 0.05, 0.1, 0.2, 0.4, 0.8]:
    a, b = beta_ab(zm, s)
    te = mean_t_exact(zm, s)
    g1, g2 = mean_t_grid(zm, s, 100), mean_t_grid(zm, s, 1000)
    print(f"{s:4.2f} {a:7.3f} {b:8.2f} {te:9.1f} {te - t_burke(zm):7.1f}"
          f" {g1 - te:+8.1f} {g2 - te:+8.1f}")

Z_mean = 0.055, T(Z_mean) = 2230.0 K
 s      a       b      T_exact   gap    grid100  grid1000
0.01   5.445    93.55    1895.7  -334.3     +5.4     -0.0
0.03   1.778    30.55    1663.1  -566.9     -5.9     -0.2
0.05   1.045    17.95    1513.3  -716.7     +3.0     -0.2
0.10   0.495     8.50    1266.2  -963.8   +114.6    +32.3
0.20   0.220     3.78     980.9 -1249.1   +335.4   +167.4
0.40   0.083     1.42     682.6 -1547.4   +585.9   +382.0
0.80   0.014     0.24     391.7 -1838.3   +664.2   +506.9

El código de arriba puede ejecutarse directamente en el navegador.

Al pulsar Play, $s$ recorre de 0.01 a 0.8 en escala logarítmica y se trazan tres curvas. Con los valores por defecto, la tabla de la derecha reproduce la salida anterior. La banda roja es la región $a<1$, y la curva azul de error se dispara justo al cruzarla. Al subir $\bar Z$ a 0.2, la banda roja se corre a la derecha y el error de $N_1=100$ también disminuye.

Lo que dicen 964 K y +115 K#

El primer número es física. Con $s=0.1$, la fluctuación rms es 0.072. Es mayor que la media 0.055 y es habitual en el borde de una capa de mezcla. Allí $\tilde T$ vale 1266 K, 964 K por debajo de $T(\tilde Z)$. Incluso una fluctuación mínima, $s=0.01$, quita 334 K. Una media sobre el pico es el punto más sensible.

Si $s$ sube a 0.8, $\tilde T$ cae a 392 K. En el límite $s\to1$, la PDF se divide en dos deltas en $Z=0$ y $Z=1$. Ambos extremos están a 300 K. Que la temperatura media vuelva a la ambiente es justo lo correcto.

El segundo número es numérico. Para $s\le0.05$, el error con 100 celdas no supera 6 K, porque $a$ está por encima de 1 o cerca. En $s=0.1$, $a=0.495$, y el error salta a +115 K. Con 1000 celdas aún quedan +32 K. En $s=0.8$, incluso 1000 celdas se pasan por +507 K. Refinar la malla 10 veces reduce el error solo 3.6 veces.

Con $a<1$, la malla uniforme pierde el aire frío#

El signo siempre es positivo, y esa es la pista. La integral en malla calienta la media. Está contando de menos la masa del lado frío.

Con $s=0.1$, la probabilidad real en la primera celda $[0, 0.01]$ es $I_{0.01}(a,b) = 0.320$. Un tercio de los paquetes de aire está en esa sola celda. Evaluar la PDF en el centro $Z=0.005$ le asigna solo 0.247. $Z^{a-1}$ se curva muy por encima de su valor central dentro de la celda. El 0.073 que falta pasa a celdas más calientes a través de la renormalización.

La regla del punto medio suele converger como $h^2$. Pero con un integrando singular como $Z^{a-1}$, el error relativo de la primera celda es fijo, sin importar la malla. Como la primera celda contiene una masa proporcional a $h^a$, el error total también baja como $h^a$. Con $a=0.495$, una malla 10 veces más fina reduce el error a $10^{-0.495}\approx 0.32$. Eso casi coincide con 115 → 32 (0.28×) de la salida.

Hay dos arreglos. Uno es asignar la masa exacta a cada celda mediante diferencias de $I_x$, como arriba. El otro es concentrar la malla geométricamente hacia $Z=0$. El primero es más barato e infalible. Además, cuando $s$ es diminuto y $\gamma$ es enorme, calcular $B(a,b)$ directamente desborda, así que se trabaja en espacio logarítmico con lgamma.

Dónde revisar la malla en Z al construir una tabla de flamelets#

La última diapositiva de FlameMaster es una sola línea: usar ListTool para volcar archivos .kg y construir una biblioteca. Cuando esa biblioteca alimenta una tabla SLF/FPV, la integración con β-PDF es el paso inmediato. Las soluciones de flamelet suelen concentrar puntos cerca de $Z_{st}$. Sin embargo, la PDF diverge en $Z=0$.

Tras construir la tabla, basta una comprobación. En $\tilde Z = Z_{st}$ y $s=0.1$, se calcula $\tilde T$ con cuadratura en malla y con $I_x$. Si difieren en decenas de kelvin, esa diferencia llega tal cual al campo de temperatura media del LES. Y en cantidades exponencialmente sensibles a la temperatura, como tasas de reacción o formación de NO, crece varias veces.