El mismo χ=50 devolvió 1759 K y 818 K — donde se pliega una librería de flamelets

Se conectó la librería y la misma celda informó dos temperaturas#

La librería se armó barriendo una llama de difusión en contraflujo sobre la tasa de disipación escalar. Los ejes eran la fracción de mezcla $Z$ y la tasa de disipación escalar $\chi$. Del lado del LES se consultaba la tabla con $(Z, \chi)$ como clave. Las celdas cercanas a la extinción alternaban entre 1759 K y 818 K. Hubo una larga búsqueda en el código de interpolación, pero no había ningún error ahí.

La tabla estaba plegada en esa coordenada. Dónde cae el pliegue, por qué $\chi$ no sirve de eje, y cómo una variable de progreso (un único número que resume cuánto avanzó la reacción) $C$ lo despliega. Los números salen de un problema de contorno de flamelet con 61 puntos. SLF (Steady Laminar Flamelet) y FPV (Flamelet/Progress Variable) no son dos modelos físicos. Son dos maneras de leer la misma curva.

Al pasar axis de chi (SLF) a C (FPV), las temperaturas que recoge una línea vertical bajan de tres puntos a uno.

A la izquierda está la curva en S en el plano $\chi$–$T_{st}$. A la derecha, el $T(Z)$ que devolvió la consulta. En chi, arrastrar el deslizador corta la curva dos veces y luego tres. En C, un solo punto cruza el pliegue sin cortarse.

Dónde guarda χ la ecuación de flamelet#

La ecuación de flamelet estacionaria es un balance de reacción y difusión unidimensional en el espacio de fracción de mezcla.

$Z$ es la fracción de mezcla (1 del lado del combustible, 0 del lado del oxidante) y $\chi = 2D|\nabla Z|^2$ es la tasa de disipación escalar en 1/s. Toda la información del flujo entra comprimida en ese único $\chi$. Un $\chi$ grande significa que la mezcla adelanta a la química. Hay una entrada anterior que subió χ por un factor de 202 hasta la extinción y trató ese término por separado.

Con una reacción irreversible de un paso $F + s\,O \to (1+s)P$ y número de Lewis unitario, las incógnitas se reducen a $T$ sola.

Aquí $q = Q/c_p$ es el aumento de temperatura por unidad de masa de combustible quemado y $s$ es la relación de masa estequiométrica. Con $s = 4$ resulta $Z_{st} = 0.2$ y una temperatura de llama adiabática de 2100 K.

Cuántas soluciones pueden compartir un mismo χ#

Tres. Una solución congelada, una rama intermedia inestable y una solución en combustión conviven en el mismo $\chi$. En $\chi = 50\ \text{1/s}$, $T_{st}$ se parte en 300.0 K, 817.8 K y 1759.3 K. Una tabla indexada por $\chi$ no tiene forma de decidir cuál de las tres debe devolver.

El pliegue está en $\chi_{ext} = 130.53\ \text{1/s}$. Pasado ese valor, la solución en combustión y la rama intermedia desaparecen juntas. En $\chi = 135$ solo queda la congelada, y eso es la extinción. Como la curva se devuelve sobre sí misma allí, $\chi$ no es una coordenada que recorra la curva una sola vez.

Leída en horizontal, la rama en combustión da 1976.0 K en $\chi=5$, 1759.3 K en $\chi=50$, 1376.3 K en $\chi=120$ y 1245.7 K en $\chi=130$. La rama intermedia es 637.9 K, 817.8 K, 1066.9 K y 1177.2 K. En $\chi=50$ la separación entre 1759.3 y 817.8 es 941.5 K. En $\chi=120$ la separación entre 1376.3 y 1066.9 baja a 309.4 K, y en $\chi=130$ el par 1245.7 y 1177.2 queda a 68.5 K. El pliegue informado, 1212.0 K, cae dentro de esos 68.5 K. En $\chi=135$ solo queda 300.0 K. La tabla devuelve tres filas y luego se corta a una en $\chi_{ext}=130.53$ 1/s.

Las filas de pendiente se leen como la distancia que queda hasta el pliegue. En 113.936 el pliegue está a 16.593 y la pendiente es 5.2 K/(1/s). En 128.879 quedan 1.650 y la pendiente es 17.9. En 130.441 el resto es 0.088 y la pendiente es 87.7. La línea $C=909.0$ K se sienta sobre ese pliegue: $\chi=130.528$, $T_{st}=1209.0$ K. $C=600.0$ K da $\chi=78.522$ y $T_{st}=900.0$ K, mientras $C=1200.0$ K da $\chi=102.079$ y $T_{st}=1500.0$ K. Con solo $\chi$, 78 y 102 son ramas distintas. $C$ las separa en 600 y 1200. $C=1500.0$ K es $\chi=40.825$, $T_{st}=1800.0$ K, cerca de los 2100 K adiabáticos.

Por qué Newton nunca encuentra la rama intermedia#

Fijar $\chi$ y resolver solo para $T$ con Newton usa el jacobiano $J = \partial R / \partial T$. En el pliegue ese $J$ es singular. El pliegue es por definición el punto donde $d\chi/dC = 0$, y esa condición coincide con el núcleo de $J$. La rama intermedia vive al otro lado de esa singularidad.

Se probó $\chi = 50$ desde dos estimaciones iniciales. La congelada convergió a 300.000 K y la de combustión completa a 1759.282 K. Ninguna de las dos produjo 817.8 K. La rama intermedia es inestable, así que nunca cae dentro del radio de convergencia de Newton. Aunque se quiera meter ese estado en la librería, no hay manera de hacerlo.

Si χ sale del eje, qué se pone en su lugar#

Se promueve $\chi$ a incógnita y se fija en su lugar una variable de progreso $C$. Aquí se tomó $C = T(Z_{st}) - T_u$. Eso deja $N+1$ incógnitas y $N+1$ ecuaciones.

El paso de Newton se vuelve un sistema lineal con un borde añadido. Aquí $b_i = \partial R_i / \partial \chi$ y $e$ es el vector unitario cuya única componente no nula está en $i_{st}$.

La matriz con borde sigue siendo no singular donde $J$ no lo es, de modo que atraviesa el pliegue sin detenerse. Por eliminación por bloques se resuelve con dos barridos tridiagonales. Se calculan $u = J^{-1}(-R)$ y $v = J^{-1}(-b)$, y luego se lee $\delta\chi = (-g - u_{i_{st}}) / v_{i_{st}}$.

import math
 
T0, Q, S, TA, A_PRE, N = 300.0, 9000.0, 4.0, 16000.0, 4.0e8, 61
ZST = 1.0 / (1.0 + S)
H = 1.0 / (N + 1)
ZS = [(i + 1) * H for i in range(N)]
IST = min(range(N), key=lambda i: abs(ZS[i] - ZST))
 
def rate(z, t):
    r = (t - T0) / Q
    yf, yo = z - r, (1.0 - z) - S * r
    if yf <= 0.0 or yo <= 0.0:
        return 0.0, 0.0
    e = math.exp(-TA / t)
    return (A_PRE * yf * yo * e,
            A_PRE * e * (-(yo + S * yf) / Q + yf * yo * TA / (t * t)))
 
def flamelet_residual(T, chi):
    R, diag, b = [0.0]*N, [0.0]*N, [0.0]*N
    for i in range(N):
        tl = T0 if i == 0 else T[i-1]
        tr = T0 if i == N-1 else T[i+1]
        lap = (tl - 2*T[i] + tr) / (H*H)
        w, dw = rate(ZS[i], T[i])
        R[i], diag[i], b[i] = 0.5*chi*lap + Q*w, -chi/(H*H) + Q*dw, 0.5*lap
    return R, diag, chi/(2*H*H), b
 
def tridiag_solve(diag, off, rhs):
    n = len(rhs); c = [0.0]*n; d = [0.0]*n
    c[0], d[0] = off/diag[0], rhs[0]/diag[0]
    for i in range(1, n):
        m = diag[i] - off*c[i-1]
        c[i], d[i] = off/m, (rhs[i] - off*d[i-1])/m
    x = [0.0]*n; x[-1] = d[-1]
    for i in range(n-2, -1, -1):
        x[i] = d[i] - c[i]*x[i+1]
    return x
 
def solve_at_fixed_chi(T, chi):
    """SLF: se fija chi como coordenada y solo se resuelve T."""
    T = T[:]
    for _ in range(200):
        R, diag, off, _ = flamelet_residual(T, chi)
        dT = tridiag_solve(diag, off, [-r for r in R])
        big = max(abs(v) for v in dT)
        lam = min(1.0, 150.0/max(1e-30, big))
        T = [max(T0, T[i] + lam*dT[i]) for i in range(N)]
        if big < 1e-9:
            return T, True
    return T, False
 
def trace_scurve(c_list):
    """FPV: se fija la variable de progreso C y se resuelve (T, chi) a la vez."""
    T, chi, arc = [T0 + 1e-3]*N, 1.0e-7, []
    for C in c_list:
        for _ in range(100):
            R, diag, off, b = flamelet_residual(T, chi)
            u = tridiag_solve(diag, off, [-r for r in R])
            v = tridiag_solve(diag, off, [-x for x in b])
            dchi = (-(T[IST] - T0 - C) - u[IST]) / v[IST]
            dT = [u[i] + v[i]*dchi for i in range(N)]
            big = max(abs(x) for x in dT)
            lam = min(1.0, 100.0/max(1e-30, big),
                      0.3*abs(chi)/max(1e-30, abs(dchi)))
            T = [T[i] + lam*dT[i] for i in range(N)]
            chi += lam*dchi
            if big < 1e-9:
                break
        arc.append((C, chi, T[:]))
    return arc
 
CMAX = Q*ZST
arc = trace_scurve([CMAX*k/600.0 for k in range(1, 600)])
chis = [p[1] for p in arc]
tsts = [p[2][IST] for p in arc]
fold = max(range(len(chis)), key=lambda i: chis[i])
print("Z_st = %.3f   T_ad(Z_st) = %.1f K   grid N = %d" % (ZST, T0+CMAX, N))
print("extinction fold: chi_ext = %.2f 1/s,  T_st = %.1f K" % (chis[fold], tsts[fold]))
 
def branches_at_chi(target):
    hits = [T0]                                    # la solucion congelada, presente para cualquier chi
    for i in range(len(arc)-1):
        a, b = chis[i], chis[i+1]
        if (a-target)*(b-target) <= 0 and a != b:
            f = (target-a)/(b-a)
            hits.append(tsts[i] + f*(tsts[i+1]-tsts[i]))
    return sorted(hits)
 
print("\n[SLF axis] how many T_st does one chi_st key return?")
for target in (5.0, 50.0, 120.0, 130.0, 135.0):
    hits = branches_at_chi(target)
    print("  chi=%6.1f -> %d entries:  %s" % (target, len(hits), ",  ".join("%.1f K" % v for v in hits)))
 
print("\n[SLF solve] chi=50 fixed, Newton from two different guesses")
for name, guess in (("frozen ", [T0]*N),
                    ("burning", [T0 + Q*min(z, (1-z)/S) for z in ZS])):
    T, ok = solve_at_fixed_chi(guess, 50.0)
    print("  init=%s -> T_st = %8.3f K  (converged=%s)" % (name, T[IST], ok))
print("  -> the middle branch at %.1f K is unreachable from any guess" % branches_at_chi(50.0)[1])
 
print("\n[FPV axis] one C gives exactly one (chi, T_st)")
for Ct in (300.0, 600.0, 909.0, 1200.0, 1500.0):
    i = min(range(len(arc)), key=lambda k: abs(arc[k][0]-Ct))
    print("  C = %7.1f K -> chi = %8.3f 1/s,  T_st = %7.1f K" % (arc[i][0], chis[i], tsts[i]))
 
print("\n[sensitivity] dT_st/dchi as the fold is approached  [K/(1/s)]")
for off in (60, 20, 5):
    i = fold - off
    print("  chi = %7.3f (fold is %6.3f away) -> %10.1f"
          % (chis[i], chis[fold]-chis[i], (tsts[i+1]-tsts[i])/(chis[i+1]-chis[i])))
 
print("\n[progress variable] choosing Z_ref in C = T(Z_ref) - T_u")
print("  Z_ref    monotone?   max/min cell width (24 cells)")
for zr in (0.048, 0.194, 0.387, 0.726, 0.952):
    j = min(range(N), key=lambda k: abs(ZS[k]-zr))
    vals = [p[2][j] - T0 for p in arc]
    mono = all(vals[i+1] >= vals[i] - 1e-9 for i in range(len(vals)-1))
    lo, hi = min(vals), max(vals)
    pos = []
    for k in range(25):
        lv = lo + (hi-lo)*k/24.0
        for i in range(len(vals)-1):
            if (vals[i]-lv)*(vals[i+1]-lv) <= 0 and vals[i] != vals[i+1]:
                f = (lv-vals[i])/(vals[i+1]-vals[i])
                pos.append((i+f)/(len(vals)-1))
    pos.sort()
    g = [pos[i+1]-pos[i] for i in range(len(pos)-1) if pos[i+1]-pos[i] > 1e-9]
    print("  %.3f    %-8s  %8.1f : 1" % (ZS[j], "yes" if mono else "FOLDS", max(g)/min(g)))

Salida:

Z_st = 0.200   T_ad(Z_st) = 2100.0 K   grid N = 61
extinction fold: chi_ext = 130.53 1/s,  T_st = 1212.0 K
 
[SLF axis] how many T_st does one chi_st key return?
  chi=   5.0 -> 3 entries:  300.0 K,  637.9 K,  1976.0 K
  chi=  50.0 -> 3 entries:  300.0 K,  817.8 K,  1759.3 K
  chi= 120.0 -> 3 entries:  300.0 K,  1066.9 K,  1376.3 K
  chi= 130.0 -> 3 entries:  300.0 K,  1177.2 K,  1245.7 K
  chi= 135.0 -> 1 entries:  300.0 K
 
[SLF solve] chi=50 fixed, Newton from two different guesses
  init=frozen  -> T_st =  300.000 K  (converged=True)
  init=burning -> T_st = 1759.282 K  (converged=True)
  -> the middle branch at 817.8 K is unreachable from any guess
 
[FPV axis] one C gives exactly one (chi, T_st)
  C =   300.0 K -> chi =    2.012 1/s,  T_st =   600.0 K
  C =   600.0 K -> chi =   78.522 1/s,  T_st =   900.0 K
  C =   909.0 K -> chi =  130.528 1/s,  T_st =  1209.0 K
  C =  1200.0 K -> chi =  102.079 1/s,  T_st =  1500.0 K
  C =  1500.0 K -> chi =   40.825 1/s,  T_st =  1800.0 K
 
[sensitivity] dT_st/dchi as the fold is approached  [K/(1/s)]
  chi = 113.936 (fold is 16.593 away) ->        5.2
  chi = 128.879 (fold is  1.650 away) ->       17.9
  chi = 130.441 (fold is  0.088 away) ->       87.7
 
[progress variable] choosing Z_ref in C = T(Z_ref) - T_u
  Z_ref    monotone?   max/min cell width (24 cells)
  0.048    yes          717.5 : 1
  0.194    yes            1.8 : 1
  0.387    yes          734.2 : 1
  0.726    yes          733.9 : 1
  0.952    yes          733.9 : 1

El eje horizontal es la longitud de arco $s$ recorrida sobre la curva en S. Al pasar axis a chi, la curva se pliega, aparece un tramo rojo y los bordes de celda uniformes caen dos veces en un mismo nivel. Al arrastrar $Z_{ref}$, esos bordes se juntan y se separan en la franja de abajo.

Dónde conviene leer C para que la tabla se estire de forma pareja#

Mientras $C$ sea monótona a lo largo de la curva, la tabla queda univaluada. Pero la monotonía sola no alcanza. Lo que fija la precisión de la interpolación es cuán desigual queda el espaciado de celdas sobre la curva una vez que $C$ se corta en una malla uniforme.

Mover $Z_{ref}$ en $C = T(Z_{ref}) - T_u$ cambia mucho ese reparto. En $Z_{ref} = 0.194$, cerca del estequiométrico, la relación entre el ancho máximo y el mínimo de celda es 1.8 : 1, casi uniforme. En $Z_{ref} = 0.048$ llega a 717.5 : 1 y en $Z_{ref} = 0.726$ a 733.9 : 1. Todas cumplen la monotonía. Aun así, al elegir un punto lejano a la llama, un extremo de la tabla se colapsa en una sola celda.

La razón 1.8 en un corte de 24 celdas aparece solo en $Z_{ref}=0.194$. La fila 0.048 da 717.5, 0.387 da 734.2, y tanto 0.726 como 0.952 dan 733.9. Cinco líneas yes no dicen que el espaciado coincida. Junto a 1.8, 717.5 mete casi toda la longitud de arco en una o dos de esas 24 celdas. Dentro de esa celda, 1759.3 K y 817.8 K pueden compartir un mismo intervalo.

Cuánto crece la sensibilidad de la tabla junto al pliegue#

En el pliegue, $dT_{st}/d\chi$ diverge. Cerca de $\chi_{ext}$ la curva se comporta como $\chi \approx \chi_{ext} - k(C - C_{ext})^2$, de modo que su inversa arrastra una singularidad de raíz cuadrada. Medido directamente: 5.2 K/(1/s) a 16.593 del pliegue, 17.9 a 1.650 y 87.7 a 0.088.

El $\chi$ de una celda LES siempre carga el error del modelo de submalla. Errar $\chi$ en un uno por ciento junto al pliegue desplaza la temperatura en decenas de kelvin. Sobre el eje $C$, ese mismo tramo es simplemente suave. El pliegue es una singularidad que creó la elección de coordenada, no la física.

Cómo es en la práctica la secuencia de FlameMaster para esta tabla#

Escrita en el orden que dan los apuntes, son cuatro pasos. Primero CreateBinFile empaqueta los datos termodinámicos y de transporte en un binario (thermo.dat, trans.dat → a.bin). Después ScanMan lee el archivo de mecanismo y emite un resultado preprocesado. Luego FlameMan resuelve un flamelet estacionario con MaxIter, TolRes, TolDy y maxGridPoints (de 100 a 200) tomados del archivo de entrada.

El último paso es libgen. Reúne los flamelets individuales en una librería indexada por $(Z, \chi)$ o $(Z, C)$. El archivo datalist decide qué flamelets entran y el archivo config decide qué especies sobreviven en la tabla. Si se indexa esa librería por $\chi$, el pliegue de arriba pasa tal cual al archivo. Capturar la rama intermedia implica que FlameMan no puede correrse a $\chi$ fijo.

Cuándo sigue alcanzando con SLF#

Si $\chi$ se mantiene lejos de $\chi_{ext}$ y la llama no tiene extinción local, basta con la rama en combustión. En ese rango $\chi \mapsto T_{st}$ es univaluada y su pendiente es suave. Para llamas levantadas que no se van a resolver, o para cálculos de diseño que nunca miran ignición ni extinción, una tabla SLF sale más barata y su eje es más intuitivo.

En cambio, si se mezclan extinción local y reignición, la rama intermedia entra en la estadística. La relación entre cuánto sacude la turbulencia a la llama y el tiempo químico fija ese límite. La misma relación decide las cosas en el punto donde la velocidad de llama turbulenta se dobla con el número de Damköhler. En ese régimen, pagar por mover el eje a $C$ cuesta menos que unas pocas veces el tamaño de la tabla.

Para no plegar dos veces una misma tabla#

Lo que se fija al construir una librería es el índice, no la física. $\chi$ resulta natural como coeficiente dentro de la ecuación de flamelet, pero como coordenada sobre el conjunto de soluciones se pliega. Al interpolar sobre una coordenada plegada, esa interpolación puede estar equivocada sea cual sea el valor que devuelva.

$C$ tampoco es gratis. Es una cantidad que el LES tiene que transportar con una ecuación más, y el término fuente de esa ecuación se vuelve a leer de la tabla. Por eso la definición de $C$ debe responder a la vez por la monotonía y por el espaciado de celdas. Mover un solo $Z_{ref}$ en el experimento de arriba llevó la relación de anchos de 1.8 a 734, y ese es justamente el punto. Una vez elegido el eje, los números a volver a medir quedan fijos: la relación de anchos de celda bajo malla uniforme y la distancia que resta hasta el pliegue.