소산율을 202배 올려도 온도는 13%만 내려갔다 — 층류 화염편의 소염 한계
소염은 온도가 천천히 내려가다 일어나지 않는다. 확산 손실이 반응 생성을 앞지르는 한 점에서 해가 사라진다.
확산이 늘어났는데 왜 화염이 꺼지는가#
화염편 라이브러리를 만들 때 격자 하나를 손으로 정해야 한다. 스칼라 소산율 를 어디까지 늘려놓고 표를 끊을지다. 너무 짧게 끊으면 난류 계산 중에 표 밖을 조회하고, 너무 길게 잡으면 존재하지도 않는 해를 찾느라 Newton이 헛돈다. 그 경계를 눈으로 보려고 1차원 화염편 방정식을 직접 풀어봤다.
결과부터 적는다. 를 0.05에서 10.12 1/s까지, 202배 키우는 동안 최고 온도는 2258.8 K에서 1960.9 K로 13.2%만 내려갔다. 그리고 한 걸음 더 갔더니 해가 없어졌다. 온도가 서서히 0으로 가는 것이 아니라, 버티던 해가 통째로 사라진다.
역설처럼 들린다. 확산은 반응물을 화염면으로 실어 나르는 유일한 수단이다. 비예혼합 화염에서 연료와 산화제는 서로 다른 쪽에서 오고, 분자확산 말고는 만날 길이 없다. 그런데 그 확산을 키웠더니 불이 꺼졌다.
strain rate a 슬라이더를 올리면 혼합층이 얇아지고, 같은 의 가 커진다.
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.
화염면 띠가 흐려지는 쪽을 본다. Z_st를 봉우리에서 밀면 신장률은 그대로인데, 반응대가 받는 소산만 빠진다.
화염을 스칼라 하나로 적을 수 있는 조건은 무엇인가#
모든 화학종의 확산계수가 같다고 두자. 이 가정을 쓰면 유동 입자의 혼합 상태는 보존 스칼라 하나로 결정된다. 혼합분율(mixture fraction) 다. 이 순수 산화제, 이 순수 연료다.
Peters는 1984년에 를 독립변수로 삼는 좌표 변환을 제안했다. 3차원 공간 좌표 대신 위에서 방정식을 적으면, 화염면에 수직인 방향 하나만 남는다.
는 화학종 의 질량분율, 는 생성률, 는 스칼라 소산율이다. 에너지 방정식도 같은 꼴을 갖는다.
는 온도, 는 정압비열, 는 화학종 엔탈피다. Ferreira는 뒤에 가 일정하다는 조건을 풀어 더 일반적인 형태를 제시했지만, 가 상수면 위 식으로 돌아온다.
여기서 눈에 걸리는 것은 유동이 사라졌다는 점이다. 속도장도, 압력도, 3차원 격자도 없다. 난류가 화염에 하는 일이 전부 하나에 압축되어 들어갔다. 같은 에서 온도가 두 갈래로 접히는 자리는 화염편 표가 1759 K와 818 K를 같이 거는 곳이다.
스칼라 소산율은 무엇의 역수인가#
의 정의는 단순하다.
는 확산계수다. 차원을 보면 1/s다. 확산에 걸리는 특성 시간의 역수로 읽으면 된다.
난류가 화염편을 잡아 늘리면 의 구배가 커진다. 는 구배의 제곱에 붙으므로 더 빠르게 커진다. 반응면으로 들어오는 분자확산이 늘어나고, 동시에 반응면에서 빠져나가는 열전달도 같은 비율로 늘어난다. 반응물이 반응할 수 있는 체류시간은 줄어든다.
이것이 앞의 역설의 답이다. 확산은 반응물을 데려오지만 열도 가져간다. 화학반응은 온도에 지수적으로 붙어 있으므로, 열을 잃는 쪽이 항상 먼저 이긴다.
대향류 층류 화염에서 의 의존성은 닫힌 형태로 적을 수 있다. Pitsch가 DNS로 검증한 형태다.
는 화학당량 혼합분율 에서의 값이다. 라이브러리의 인덱스가 되는 단 하나의 숫자다. 그리고 는 연료와 공기의 조성만으로 정해진다.
는 연료 1 kg을 태우는 데 필요한 산소 질량, 은 연료 쪽 유입 질량분율, 은 산화제 쪽 산소 질량분율이다.
소염 한계는 코드로 어떻게 잡는가#
한 단계 비가역 반응 를 쓰면 미지 함수가 하나로 줄어든다. 확산계수가 같으므로 는 에 대해 선형이다.
는 연료 단위 질량당 발열량을 온도 단위로 적은 것이다. 도 같은 방식으로 에 묶인다. 그래서 풀 것은 하나, 3중대각 행렬 하나다.
전략은 연속화다. 아주 작은 에서 Burke–Schumann 해로 출발한다. 를 1.35배씩 올리며 Newton을 돌리고, 이전 해를 다음 초기추정으로 넘긴다. 연소 해가 사라지는 순간 Newton은 수렴하지 못하거나 상온으로 떨어진다. 그 사이를 기하평균으로 40번 이분하면 한계값이 나온다.
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 # 양쪽 유입 온도 [K]
S_RATIO = 4.0 # 연료 1 kg당 산소 [kg]
Y_OX = 0.232 # 공기 중 산소 질량분율
Q_HEAT = 36000.0 # 발열량을 온도로 환산 [K]
T_ACT = 24000.0 # 활성화 온도 [K]
A_PRE = 2.0e10 # 전지수 인자 [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):
# Pitsch의 대향류 형태. 양끝에서 0이 되므로 바닥을 깔아준다
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):
# 결합 함수: Y_F, Y_O를 T와 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 # 온도 갱신 제한
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):
# 무한 반응속도 극한 — 연속화의 출발점
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): # 마지막 연소 해를 기하 이분으로 조인다
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다이얼을 오른쪽으로 밀면 주황 가지 위의 점이 소산율과 함께 움직인다.
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.
점을 주황색 가지 위에서 오른쪽으로 밀다가 빨간 선을 넘기면 바닥으로 떨어진다. 이때 를 다시 왼쪽으로 당겨도 올라오지 않는다. 동결 해도 정상해이기 때문이다 — 히스테리시스가 이 곡선의 본질이다.
표의 앞 6줄이 평평한 것에 주목하자. 가 90배 늘어나는 동안 온도는 94.8 K 떨어졌다. 마지막 한 칸에서 나머지 203 K가 한꺼번에 빠진다. 이것이 S-곡선의 무릎이다.
에 , 를 넣는다. 이면 분모가 이고, 그 역수가 출력의 이다. 이면 같은 식이 둘째 표의 0.3169를 되돌린다. Burke–Schumann 2273.2 K와 의 2258.8 K는 이미 14.4 K 떨어져 있다. 무한 반응속도 극한과 첫 연소 해가 그 간격이다.
그 다음 가 0.0500에서 4.5079로 가는 동안 온도는 2258.8에서 2164.0으로, 94.8 K만 빠진다. 10.1232의 한 칸이 1960.9 K까지 나머지 203 K를 가져간다. 표의 drop 297.9 K와 같다. 0.1230, 0.3027, 0.7447, 1.8322를 지나는 동안 하락은 7.7 K, 13.4 K, 25.8 K, 47.9 K다. 무릎 전까지는 거의 직선이다. 1.35배씩 올리는 연속화는 마지막 칸을 성기게 지나친다. 40번의 기하 이분이 을 조이는 이유다.
연료를 희석하면 한계가 왜 더 크게 무너지는가#
출력의 두 번째 표가 그것을 보여준다. 연료 질량분율을 1.0에서 0.125로, 8배 희석했다. 단열 화염온도는 2273.2 K에서 1719.8 K로 24% 내려갔다. 그런데 소염 한계 는 10.123에서 0.561로, 18배 떨어졌다.
온도가 24% 내려갈 때 한계가 1800% 내려간 셈이다. 이유는 Arrhenius 항에 있다. 반응률은 에 붙어 있고, 여기서 K다. 가 2273에서 1720으로 가면 지수는 에서 으로 움직인다. 지수 차이 3.4는 곱으로 30배다.
가 0.055에서 0.317로 옮겨간 것도 같은 방향으로 작용한다. 앞의 첫 번째 시뮬레이션에서 Z_st를 오른쪽으로 밀어보면 보인다. 의 봉우리는 근처에 있으므로, 가 커질수록 반응대가 봉우리에 가까워진다. 같은 신장률로 더 많은 소산을 맞는다.
실무에서 이 조합이 나오는 자리가 배기가스 재순환(EGR)이다. 연료나 산화제를 불활성 기체로 희석하면 는 완만하게 내려가는데 소염 한계는 급격히 내려간다. 국소 소염이 먼저 나타나고, 그다음에 화염 부상(lift-off)이, 마지막에 날아감(blow-off)이 온다.
무릎 오른쪽을 조회하면 상온이 나온다#
격자를 까지 채우면 표는 촘촘해지지만, 4.5079에서 2164.0 K이던 온도가 10.1232에서 1960.9 K로 한 칸에 203 K가 빠진다. 그 무릎을 점 몇 개로 적을지가 흥정이다. 인 셀의 반환값은 표가 아니라 모델이 정한다. 동결 해로 보내거나, 소염 확률로 가중하거나, 비정상 화염편으로 넘긴다.
자체는 난류 모델의 출력이다. 인덱스가 2배 어긋나면 무릎의 좌우가 뒤집힌다. 상온으로 떨어진 계산은 반응 메커니즘보다 입력부터 본다. 그 표를 만드는 3중대각 한 번의 비용은 음해 확산의 Thomas 소거와 같은 자리다.
관련
도움이 됐다면 공유해주세요.