같은 χ=50에 1759 K와 818 K가 함께 걸렸다 — 화염편 라이브러리가 접히는 자리

χ는 S자 곡선에서 접히므로 테이블 축이 될 수 없다. 축을 진행변수 C로 바꾸면 같은 곡선이 단값 함수가 된다.

라이브러리를 꽂았더니 같은 셀에서 온도가 두 갈래로 갈렸다#

대향류 확산화염을 소산율만 바꿔가며 쭉 계산해 라이브러리를 만들었다. 축은 혼합분율 $Z$ 와 스칼라 소산율 $\chi$ 두 개였다. LES 쪽에서 $(Z, \chi)$ 를 키로 조회하도록 붙였다. 소염 근처 셀에서 온도가 1759 K와 818 K를 번갈아 돌려줬다. 보간 코드를 한참 뒤졌지만 버그가 아니었다.

표가 그 좌표에서 접혀 있었다. 접힘이 생기는 지점, $\chi$ 를 축으로 못 쓰는 이유, 진행변수(progress variable, 반응 진행도를 하나의 수로 압축한 양) $C$ 가 그걸 펴는 방식. 숫자는 61점 화염편 경계값 문제에서 나온다. SLF(Steady Laminar Flamelet)와 FPV(Flamelet/Progress Variable)의 차이는 물리 모델이 아니다. 같은 곡선을 어느 좌표로 읽느냐다.

axis를 chi (SLF)에서 C (FPV)로 넘기면, 한 수직선이 집어 오는 온도가 세 점에서 한 점으로 줄어든다.

table axis:  chi_st = 50.0 1/s 3 entries at one key: 300 K, 815 K, 1748 K

왼쪽은 $\chi$ – $T_{st}$ 평면의 S자 곡선이고 오른쪽은 조회가 돌려준 $T(Z)$ 다. chi에서 슬라이더를 끌면 수직선이 곡선을 두 번, 세 번 자른다. C에서는 점 하나가 접힘을 지나도 끊기지 않는다.

화염편 방정식은 χ를 어디에 들고 있나#

정상 화염편 방정식은 혼합분율 공간의 1차원 반응-확산 식이다.

\frac{\chi}{2}\frac{d^2 T}{dZ^2} + \frac{\dot\omega_T}{\rho\,c_p} = 0, \qquad T(0) = T(1) = T_u

$Z$ 는 혼합분율(연료 쪽 1, 산화제 쪽 0), $\chi = 2D|\nabla Z|^2$ 는 스칼라 소산율(단위 1/s)이다. 유동 정보는 전부 $\chi$ 하나로 압축되어 들어온다. $\chi$ 가 크면 혼합이 빨라 반응이 따라오지 못한다. 소산율을 202배 올려 소염까지 밀어본 글에서 이 항 자체를 따로 다뤘다.

단일단계 비가역 반응 $F + s\,O \to (1+s)P$ 에 Lewis 수 1을 쓰면 미지수가 $T$ 하나로 줄어든다.

Y_F = Z - \frac{T - T_u}{q}, \qquad Y_O = (1 - Z) - s\,\frac{T - T_u}{q}, \qquad \dot\omega = A\,Y_F Y_O \exp\!\left(-\frac{T_a}{T}\right)

$q = Q/c_p$ 는 연료 단위질량이 내는 온도상승이고 $s$ 는 화학양론 질량비다. $s = 4$ 면 $Z_{st} = 0.2$ , 단열 화염온도는 2100 K가 된다.

같은 χ에 해가 몇 개까지 있을 수 있나#

셋이다. 동결 해, 불안정한 중간 가지, 연소 해가 같은 $\chi$ 를 공유한다. $\chi = 50\ \text{1/s}$ 에서 $T_{st}$ 는 300.0 K, 817.8 K, 1759.3 K로 갈린다. $\chi$ 를 키로 쓰는 테이블은 이 셋 중 무엇을 돌려줄지 스스로 정할 수 없다.

접힘점은 $\chi_{ext} = 130.53\ \text{1/s}$ 에 있다. 이 값을 넘기면 연소 해와 중간 가지가 함께 사라진다. $\chi = 135$ 에서는 동결 해 하나만 남는다. 이것이 소염이다. 곡선이 그 지점에서 되돌아 나오기 때문에 $\chi$ 는 곡선을 한 번만 긁는 좌표가 아니다.

출력의 연소 가지를 가로로 읽으면 $\chi=5$ 에서 1976.0 K, $\chi=50$ 에서 1759.3 K, $\chi=120$ 에서 1376.3 K, $\chi=130$ 에서 1245.7 K다. 중간 가지는 637.9 K, 817.8 K, 1066.9 K, 1177.2 K다. $\chi=50$ 에서 1759.3과 817.8의 간격은 941.5 K다. $\chi=120$ 에서는 1376.3과 1066.9의 간격이 309.4 K로 줄고, $\chi=130$ 에서는 1245.7과 1177.2가 68.5 K까지 붙는다. 접힘 보고값 1212.0 K는 그 68.5 K 안에 있다. $\chi=135$ 에서는 300.0 K 하나만 남는다. 세 줄을 돌려주다가 한 줄로 끊기는 자리가 $\chi_{ext}=130.53$ 1/s다.

기울기 줄은 접힘까지 남은 $\chi$ 로 읽힌다. 113.936에서는 16.593 떨어져 5.2 K/(1/s)다. 128.879에서는 1.650이 남고 17.9다. 130.441에서는 0.088밖에 안 남고 87.7이다. $C=909.0$ K 줄은 $\chi=130.528$ , $T_{st}=1209.0$ K로 이 접힘 바로 위에 앉는다. $C=600.0$ K는 $\chi=78.522$ , $T_{st}=900.0$ K이고 $C=1200.0$ K는 $\chi=102.079$ , $T_{st}=1500.0$ K다. $\chi$ 만 보면 78과 102가 다른 가지인데, $C$ 는 600과 1200으로 갈라 놓는다. $C=1500.0$ K는 $\chi=40.825$ , $T_{st}=1800.0$ K로 단열 2100 K 쪽에 가깝다.

중간 가지는 왜 Newton으로 잡히지 않나#

$\chi$ 를 고정하고 $T$ 만 Newton으로 풀면 야코비안 $J = \partial R / \partial T$ 를 쓴다. 접힘점에서 이 $J$ 가 특이행렬이 된다. 접힘점은 정의상 $d\chi/dC = 0$ 인 자리이고, 그 조건이 $J$ 의 영공간과 같은 것이다. 중간 가지는 이 특이점 반대편에 있다.

실제로 $\chi = 50$ 에 초기추정을 두 가지로 넣어봤다. 동결 추정은 300.000 K로, 완전연소 추정은 1759.282 K로 수렴했다. 817.8 K는 어느 쪽에서도 나오지 않았다. 중간 가지는 불안정해서 Newton의 수렴 반경 안에 들어오지 않는다. 라이브러리에 그 상태를 넣고 싶어도 넣을 방법이 없다.

χ를 축에서 내리면 무엇을 올려야 하나#

$\chi$ 를 미지수로 승격하고, 대신 진행변수 $C$ 를 고정한다. 여기서는 $C = T(Z_{st}) - T_u$ 로 잡았다. 미지수가 $N+1$ 개, 식이 $N+1$ 개가 된다.

\begin{aligned} R_i(T, \chi) &= 0, \qquad i = 1, \dots, N \\ T_{i_{st}} - T_u - C &= 0 \end{aligned}

Newton 단계는 테두리가 하나 붙은 선형계가 된다. $b_i = \partial R_i / \partial \chi$ , $e$ 는 $i_{st}$ 성분만 1인 단위벡터다.

\begin{bmatrix} J & b \\ e^{\mathsf T} & 0 \end{bmatrix} \begin{bmatrix} \delta T \\ \delta \chi \end{bmatrix} = - \begin{bmatrix} R \\ g \end{bmatrix}

$J$ 가 특이해도 이 테두리 행렬은 정칙이다. 그래서 접힘점을 그냥 통과한다. 블록 소거로 풀면 삼중대각 해법 두 번으로 끝난다. $u = J^{-1}(-R)$ , $v = J^{-1}(-b)$ 를 구한 뒤 $\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: chi를 좌표로 고정하고 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: 진행변수 C를 좌표로 고정하고 (T, chi)를 함께 푼다."""
    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]                                    # 어느 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)))

출력:

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

axis:  Z_ref = 0.190 cells = 24 single-valued, cell stretch 4:1

가로축은 S자 곡선을 따라간 호길이 $s$ 다. axis를 chi로 바꾸면 곡선이 접히면서 빨간 구간이 뜨고, 균일 셀 경계가 한 레벨에서 두 번 찍힌다. Z_ref 슬라이더를 끌면 셀 경계가 아래 띠에서 뭉쳤다 풀린다.

C는 어디서 읽어야 테이블이 고르게 펴지나#

$C$ 가 곡선을 따라 단조이기만 하면 단값 테이블이 된다. 그런데 단조성만으로는 부족하다. $C$ 를 균일 격자로 잘랐을 때 곡선 위의 셀 간격이 얼마나 벌어지는지가 보간 정확도를 정한다.

$C = T(Z_{ref}) - T_u$ 에서 $Z_{ref}$ 를 옮겨보면 차이가 크다. $Z_{ref} = 0.194$ (화학양론 근처)는 최대/최소 셀 간격비가 1.8 : 1로 거의 균일하다. $Z_{ref} = 0.048$ 에서는 717.5 : 1, $Z_{ref} = 0.726$ 에서는 733.9 : 1까지 벌어진다. 단조는 전부 만족한다. 그래도 화염에서 먼 지점을 쓰면 테이블 한쪽이 셀 하나로 뭉개진다.

24칸으로 자른 그 비가 1.8인 줄은 $Z_{ref}=0.194$ 뿐이다. 0.048은 717.5, 0.387은 734.2, 0.726과 0.952는 733.9다. 출력의 yes 다섯 줄이 간격까지 같다는 뜻은 아니다. 717.5와 1.8을 나란히 두면 같은 24칸인데 호길이의 대부분이 칸 한두 개에 쌓인다. 그 칸 안에서 1759.3 K와 817.8 K가 한 셀에 섞일 수 있다.

접힘 근처에서 테이블 민감도는 얼마나 커지나#

접힘점에서 $dT_{st}/d\chi$ 가 발산한다. 곡선이 $\chi_{ext}$ 근처에서 $\chi \approx \chi_{ext} - k(C - C_{ext})^2$ 꼴이므로 역함수가 제곱근 특이성을 갖는다. 숫자로 재보면 접힘까지 16.593 떨어진 지점에서 5.2 K/(1/s), 1.650 떨어진 지점에서 17.9, 0.088 떨어진 지점에서 87.7이다.

LES 셀의 $\chi$ 에는 아화격자 모델 오차가 늘 섞여 있다. 접힘 근처에서 $\chi$ 를 1% 틀리게 추정하면 온도가 수십 K 틀어진다. $C$ 축에서는 같은 구간이 그냥 평탄하다. 접힘은 물리가 아니라 좌표 선택이 만든 특이성이다.

FlameMaster로 이 표를 실제로 뽑는 순서는#

강의노트의 순서를 그대로 적으면 네 단계다. 먼저 CreateBinFile로 열역학·수송 데이터를 바이너리로 묶는다(thermo.dat, trans.dat → a.bin). 다음으로 ScanMan이 메커니즘 파일을 읽어 전처리 결과를 낸다. 그 다음 FlameMan이 입력파일의 MaxIter, TolRes, TolDy와 maxGridPoints(100~200)로 정상해 하나를 푼다.

마지막이 libgen이다. 화염편 하나하나를 모아 $(Z, \chi)$ 또는 $(Z, C)$ 로 색인한 라이브러리를 만든다. 여기서 datalist 파일에 어떤 화염편을 넣을지, config 파일에 어떤 화학종을 테이블에 남길지 정한다. 이 색인 축을 $\chi$ 로 두면 위에서 본 접힘이 그대로 파일에 들어간다. 중간 가지를 담으려면 FlameMan을 $\chi$ 고정으로 돌려서는 안 된다.

그래도 SLF로 충분한 경우는 언제인가#

$\chi$ 가 $\chi_{ext}$ 에서 멀고 국소소염이 없는 화염이면 연소 가지만 쓰면 된다. 그 영역에서는 $\chi \mapsto T_{st}$ 가 단값이고 기울기도 완만하다. 부상화염이나 점화·소염을 보지 않는 설계 계산에서는 SLF 테이블이 더 싸고 축도 직관적이다.

반대로 국소소염과 재점화가 섞이는 조건이면 중간 가지가 통계에 들어온다. 난류가 화염을 흔드는 정도와 화학시간의 비가 이 경계를 정한다. 난류 화염속도가 Damköhler 수에서 꺾이는 지점과 같은 비율이 여기서도 기준이 된다. 그 영역에서는 축을 $C$ 로 바꾸는 비용이 테이블 용량 몇 배보다 싸다.

표 하나를 두 번 접지 않으려면#

라이브러리를 만들 때 정하는 것은 물리가 아니라 색인이다. $\chi$ 는 화염편 방정식의 계수로는 자연스럽지만, 해집합의 좌표로는 접힌다. 접힘이 있는 좌표에 보간을 걸면 그 보간은 어떤 값을 돌려줘도 틀릴 수 있다.

$C$ 를 쓸 때도 공짜는 아니다. $C$ 는 LES에서 수송방정식 하나를 더 풀어야 하는 양이고, 그 방정식의 소스항은 다시 테이블에서 읽어온다. 그래서 $C$ 의 정의는 단조성과 셀 간격 두 가지를 동시에 봐야 한다. 위 실험에서 $Z_{ref}$ 하나를 옮겼을 때 간격비가 1.8에서 734로 뛴 것이 그 이야기다. 축을 고르고 나면 다음에 재볼 숫자는 정해진다. 균일 격자로 자른 셀 간격비, 그리고 접힘까지의 거리다.