혼합분율이 조금만 흔들려도 평균 온도가 964 K 내려간다 — presumed β-PDF와 Z=0의 무한대

평균 혼합분율에서 화염 온도를 읽으면 왜 틀릴까#

메탄-공기 난류 제트 화염을 하나 떠올려 보자. 어느 한 지점에서 시간 평균한 혼합분율이 정확히 화학양론값 0.055다. 그 지점의 평균 온도는 몇 K일까. 화학양론에서 단열 화염 온도가 2230 K이니 2230 K라고 답하고 싶어진다.

답은 1266 K다. 혼합분율의 흔들림이 평균의 0.07배 수준일 때 그렇다. 이 글은 그 964 K가 어디서 오는지 보여준다. 그리고 그 평균을 계산하는 presumed β-PDF(확률밀도함수의 모양을 베타 분포로 가정하는 방법)를 직접 적분해 본다. 마지막으로 균일 격자 100칸으로 적분하면 왜 115 K가 틀리는지까지 따라간다.

출발점은 2018년 여름 연소 세미나에서 FlameMaster 사용법을 정리한 슬라이드다. 설치·ScanMan·FlameMaster 실행법 사이에 이론 슬라이드 몇 장이 끼어 있다. PDF와 CDF, Beta 함수, 혼합분율과 소산율, S-curve다. 슬라이드는 제목만 남아 있다. 여기서는 그 빈칸 중 "PDF와 Beta 함수"를 채운다.

Z는 섞임을, T(Z)는 화학을 맡는다#

혼합분율 $Z$(연료 흐름에서 온 질량의 비율)는 공기에서 0, 연료에서 1이다. 반응 속도가 무한히 빠르다고 보면 온도는 $Z$만의 함수가 된다. 이것이 Burke–Schumann 해다. 화학양론 $Z_{st}$에서 꼭짓점을 찍고 양쪽으로 직선을 내린 지붕 모양이다.

난류 속의 한 점에서 $Z$는 시간에 따라 흔들린다. 공기 덩어리가 지나가면 0 근처로, 연료 덩어리가 지나가면 큰 값으로 튄다. 평균 온도는 이 흔들림 전체에 대한 $T(Z)$의 평균이다.

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

$P(Z)$는 그 점에서 $Z$의 확률밀도, $\tilde Z$는 평균 혼합분율이다.

$T(Z)$가 직선이면 두 값이 같다. 지붕 모양은 꼭짓점에서 꺾인다. 꺾인 점 바로 위에 평균을 두면 흔들림은 어느 쪽으로든 온도를 내린다. 아래 시뮬레이션에서 직접 조작해보자.

입자 하나하나가 β-PDF에서 뽑은 $Z$다. 지붕 곡선에 떨어진 높이가 그 입자의 온도다. 분홍 점이 입자들의 누적 평균이고, 곧 초록 선 $\tilde T$에 붙는다. $s$를 0.1로 두면 입자의 3분의 1이 $Z<0.01$에 쏟아진다. $\bar Z$를 0.2로 옮겨 보자. 지붕의 오른쪽 경사가 완만해서 간격이 훨씬 줄어든다.

β-PDF: 평균과 분산 두 개로 모양을 정한다#

LES·RANS 솔버가 넘겨주는 것은 $\tilde Z$와 분산 $\widetilde{Z''^2}$ 두 개뿐이다. PDF 전체를 풀지 않는다. 그래서 모양을 미리 정해 둔다. 0과 1 사이에 갇힌 변수에는 베타 분포가 자연스럽다.

$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)$는 베타 함수(정규화 상수)다. $s = \widetilde{Z''^2}/[\tilde Z(1-\tilde Z)]$는 분리 강도(segregation intensity)로, 0이면 완전 혼합, 1이면 0과 1만 있는 완전 분리다.

화학양론 근처에서는 $\tilde Z$가 작다. 그래서 $a$가 금방 1 밑으로 내려간다. $\tilde Z = 0.055$이면 $s > 0.052$에서 $a<1$이다. 이때 $Z^{a-1}$은 $Z\to0$에서 무한대로 간다. 적분은 유한하다. 그러나 PDF 값 자체는 끝점에서 발산한다. 이 성질이 뒤에서 균일 격자 적분을 망친다.

꺾인 직선은 불완전 베타 함수로 정확히 적분된다#

$T(Z)$가 구간별 직선이면 수치 적분이 필요 없다. 각 구간에서 필요한 것은 두 적분뿐이다.

$\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)$는 정규화된 불완전 베타 함수(베타 분포의 누적분포함수)다. 두 번째 식은 $Z\cdot Z^{a-1} = Z^{a}$이고 $B(a+1,b)/B(a,b) = a/(a+b) = \tilde Z$이기 때문에 성립한다.

$Z_{st}$를 경계로 희박 쪽과 과농 쪽 직선을 각각 넣으면 $\tilde T$가 닫힌 꼴로 나온다. 끝점의 무한대는 $I_x$ 안에서 해석적으로 처리된다. $I_x$는 연분수(Lentz 방법)로 계산한다. Numerical Recipes의 betai와 같은 알고리즘이다.

화염편 라이브러리에서도 사정은 같다. 소산율과 소염 한계에서 본 $T(Z;\chi)$ 곡선이 지붕 대신 들어간다. 테이블의 $Z$ 격자가 구간별 직선 보간이라면 위 공식을 구간마다 쓰면 된다.

코드: 화학양론 Z̄에서 분산을 0.01부터 0.8까지#

$\tilde Z = Z_{st} = 0.055$에 고정하고 $s$만 바꾼다. 정확한 적분과 함께 "PDF를 셀 중심에서 찍고 합을 1로 다시 맞추는" 균일 격자 적분을 $N=100$, $1000$으로 돌린다. 마지막 두 열은 정확한 값과의 차이다.

import math
 
# 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):
    # 불완전 베타 함수의 연분수 (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):
    # 정규화된 불완전 베타 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)), 분리 강도
    g = 1.0 / s - 1.0
    return zm * g, (1.0 - zm) * g
 
def mean_t_exact(zm, s):
    # T가 구간별 직선이라 각 구간을 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):
    # 단순 방식: 셀 중심 n개에서 pdf를 찍고 재정규화
    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

위 코드를 브라우저에서 직접 돌려보자.

Play를 누르면 $s$를 0.01에서 0.8까지 로그 간격으로 쓸면서 세 곡선을 그린다. 기본값에서 오른쪽 표가 위 출력과 같은 값을 낸다. 빨간 영역이 $a<1$ 구간이고, 파란 오차 곡선은 그 경계를 넘자마자 치솟는다. $\bar Z$를 0.2로 올리면 빨간 영역이 오른쪽으로 밀리고 $N_1=100$의 오차도 함께 줄어든다.

964 K와 +115 K가 말하는 것#

첫째 숫자는 물리다. $s=0.1$은 rms 흔들림이 0.072다. 평균 0.055보다 크다. 이 정도면 혼합층 가장자리에서 흔한 값이다. 그 점에서 $\tilde T$는 1266 K이고 $T(\tilde Z)$보다 964 K 낮다. $s=0.01$로 아주 작게 흔들려도 334 K가 빠진다. 꼭짓점 위의 평균은 흔들림에 가장 민감하다.

$s$를 0.8까지 올리면 $\tilde T$는 392 K까지 내려간다. 극한 $s\to1$에서 PDF는 $Z=0$과 $Z=1$의 두 델타로 갈라진다. 두 끝은 모두 300 K다. 평균 온도가 상온으로 돌아가는 것이 맞다.

둘째 숫자는 수치다. $s\le0.05$에서 100칸 격자 오차는 6 K 이하다. $a$가 1을 넘거나 1 근처이기 때문이다. $s=0.1$에서 $a=0.495$가 되면 오차가 +115 K로 뛴다. 1000칸으로 늘려도 +32 K가 남는다. $s=0.8$에서는 1000칸에서도 +507 K다. 격자를 10배 촘촘히 해도 오차는 3.6분의 1밖에 줄지 않는다.

$a<1$이면 균일 격자는 차가운 공기를 놓친다#

부호가 늘 플러스라는 점이 원인을 알려준다. 격자 적분은 평균을 덥게 만든다. 차가운 쪽 질량을 덜 세고 있다는 뜻이다.

$s=0.1$에서 첫 셀 $[0, 0.01]$에 실제로 들어 있는 확률은 $I_{0.01}(a,b) = 0.320$이다. 공기 덩어리의 3분의 1이 이 한 칸에 있다. 셀 중심 $Z=0.005$에서 PDF를 찍으면 이 칸의 몫이 0.247로 잡힌다. $Z^{a-1}$이 셀 안에서 중심값보다 훨씬 크게 휘기 때문이다. 빠진 0.073은 재정규화를 거쳐 더 뜨거운 셀들로 흘러간다.

셀 중심 규칙의 오차는 원래 $h^2$로 준다. 하지만 피적분함수가 $Z^{a-1}$처럼 특이하면 첫 셀의 상대 오차가 격자와 무관하게 일정하다. 첫 셀 질량이 $h^a$에 비례하니 전체 오차도 $h^a$로 준다. $a=0.495$이면 격자를 10배 촘촘히 할 때 오차가 $10^{-0.495}\approx 0.32$배다. 출력의 115 → 32(0.28배)와 거의 맞는다.

고치는 방법은 둘이다. 하나는 위처럼 셀마다 $I_x$의 차이로 질량을 정확히 넣는 것이다. 다른 하나는 $Z=0$ 근처로 격자를 기하급수적으로 몰아주는 것이다. 첫 번째가 싸고 확실하다. 그리고 $s$가 아주 작아 $\gamma$가 커지면 $\ln B(a,b)$ 계산에서 넘침이 생기니, lgamma로 로그 공간에서 계산한다.

화염편 테이블을 만들 때 Z 격자를 다시 볼 자리#

FlameMaster 슬라이드의 마지막은 ListTool로 .kg 파일을 뽑아 라이브러리를 만든다는 한 줄이다. 그 라이브러리를 SLF·FPV 테이블로 쓸 때 β-PDF 적분이 바로 다음 단계다. 화염편 해는 보통 $Z_{st}$ 근처에 점을 몰아 둔다. 정작 PDF가 발산하는 곳은 $Z=0$이다.

테이블을 만든 뒤 한 번만 확인해 보자. $\tilde Z = Z_{st}$, $s=0.1$에서 $\tilde T$를 격자 적분과 $I_x$ 적분으로 각각 계산한다. 두 값이 수십 K 이상 다르면 그 차이가 그대로 LES의 평균 온도장에 실린다. 그리고 반응률·NO 생성처럼 온도에 지수적으로 민감한 양에서는 몇 배로 커진다.