μ와 k를 같은 식에서 뽑으니 Pr이 1.98배 틀렸다 — 벽온도 123.5 K를 가른 계수 9/4

μ와 k를 같은 유도로 뽑으면 Pr은 무조건 γ가 된다#

공기를 10기압으로 압축해도 열전도도는 거의 변하지 않는다. 밀도는 10배가 되었는데 열은 같은 속도로 흐른다. 분자운동론은 이 반직관을 두 줄로 설명한다. 그런데 같은 두 줄로 점성계수까지 뽑으면 프란틀 수(Prandtl number, 점성 확산과 열 확산의 비)가 $\gamma$에 못 박힌다. 공기라면 1.400이다. 실제 측정값은 0.707이다. 1.98배 틀린 숫자를 CFD 코드에 넣으면 무엇이 어긋나는지, 그리고 1913년에 Eucken이 채워 넣은 계수 $9/4$가 정확히 무엇을 고치는지가 숫자로 갈린다. 수밀도 슬라이더 number density를 1에서 8로 올리면 crossings와 λ가 반대로 움직이고, k̂는 거의 안 움직인다.

두 막대가 반대로 가는 동안 k̂는 제자리다. n = 1에서만 $\lambda$가 벽 간격보다 커지고 k̂가 내려앉는다.

충돌면을 건너는 분자가 옮기는 두 가지#

가상의 평면 하나를 유동 안에 세운다. 이 면을 단위시간·단위면적당 건너는 분자 수는

$n$은 수밀도, $\bar{c}$는 평균 분자속력이다. 평면을 아래로 건너는 분자는 평균자유경로 (mean free path, 충돌 사이에 분자가 날아가는 평균 거리) $\lambda$만큼 위에서 마지막으로 충돌했다. 그래서 그 위치의 값을 들고 내려온다. 위로 건너는 분자는 $\lambda$만큼 아래의 값을 들고 올라간다. 운동량을 들고 오면 점성이 되고, 운동에너지를 들고 오면 열전도가 된다.

$\rho = nm$, $c_v$는 단위질량당 정적비열이다. 두 식은 앞부분이 완전히 같다. 뒤에 붙은 $c_v$ 하나만 다르다.

밀도를 8배 올려도 열전도도가 그대로인 이유#

$\lambda$는 밀도에 반비례한다. 충돌 단면적 $\sigma$에 대해

그래서 $\rho \bar{c} \lambda$ 안에서 $n$과 $1/n$이 서로를 지운다.

밀도가 사라졌다. 분자 수가 늘면 면을 건너는 분자는 많아지지만, 각 분자가 더 짧은 거리에서 출발하므로 들고 오는 온도 차가 작아진다. k̂는 노란 평면을 실제로 건넌 분자가 옮긴 온도를 센 값이다. n = 1은 연속체 가정이 깨지는 자리다.

세 경로를 한 표에 올린다#

이제 두 식을 나눈다. 프란틀 수의 정의는 $\Pr = \mu c_p / k$다.

$\rho$, $\bar{c}$, $\lambda$가 모두 약분된다. 분자 지름도, 온도도, 압력도 답에 남지 않는다. 기체의 비열비만 남는다. 이건 유도의 사고가 아니라 유도의 구조다. 운동량과 에너지를 같은 운반자에 같은 방식으로 실었기 때문이다. 점성이 분자 한 줄로 적히는 자리는 점성의 분자적 기원과 같다. 어느 한 단계의 계수를 틀린 것이 아니라, 두 개의 서로 다른 물리를 한 개의 스케일로 묶은 것이다.

회전 자유도는 운동량을 나르지 않는다#

$\Pr = \gamma$가 틀린 이유는 $c_v$를 한 덩어리로 쓴 데 있다. 이원자 분자의 $c_v$는 병진 3자유도와 회전 2자유도로 나뉜다. 운동량 전달은 병진 운동만 관여한다. 회전하는 아령이 평면을 건너도, 그 회전 에너지는 전단응력에 기여하지 않는다.

Eucken은 두 성분을 분리해 각각에 다른 확산 효율을 주었다. 병진 에너지는 속도와 상관되어 있어 빠른 분자가 더 많이 나르므로 $5/2$배의 가중을 받는다. 내부 에너지는 속도와 무관하므로 가중이 1이다. 결과는 한 줄이다.

$\gamma \to 5/3$인 단원자 기체에서 이 값은 $1 + 9/4 \cdot 2/3 = 2.5$, 즉 정확히 $5/2$가 된다. 회전이 없으니 병진 가중만 남는다. 양변을 모으면

기체 여섯 개에 세 경로를 모두 걸어봤다#

세 경로를 나란히 계산하고, 틀린 Pr이 실제 계산에서 무엇을 바꾸는지까지 이어 붙였다. 전단 가열이 있는 1차원 Couette 층이다. 한쪽 벽은 $U$로 움직이고, 반대쪽 벽은 단열이다.

$u$가 선형이면 소스항이 상수가 되고, 단열 벽의 온도 상승은 $\Pr$에만 비례한다.

import math
 
GASES = [
    ("He",  1.667, 0.68),
    ("Ar",  1.667, 0.67),
    ("N2",  1.404, 0.72),
    ("Air", 1.400, 0.707),
    ("CO2", 1.289, 0.77),
    ("H2O", 1.329, 0.98),
]
 
def pr_mean_free_path(gamma):
    """mu = (1/3) rho c_bar lam,  k = (1/3) rho c_bar lam c_v  ->  Pr = c_p/c_v."""
    return gamma
 
def pr_eucken(gamma):
    """k/(mu c_v) = 1 + (9/4)(gamma-1)  ->  Pr = 4 gamma / (9 gamma - 5)."""
    return 4.0 * gamma / (9.0 * gamma - 5.0)
 
print("gas   gamma   Pr_mfp  Pr_euck  Pr_meas   mfp/meas  euck/meas")
for name, g, meas in GASES:
    a, b = pr_mean_free_path(g), pr_eucken(g)
    print("%-5s %5.3f  %6.3f  %7.3f  %7.3f   %8.2f  %9.2f" % (name, g, a, b, meas, a / meas, b / meas))
 
def couette_wall_rise(Pr, U, cp, N=201):
    """d2T/dy2 = -Pr U^2/cp on y in [0,1]; adiabatic at y=0, T=0 at y=1. Thomas sweep."""
    h = 1.0 / (N - 1)
    src = Pr * U * U / cp
    lo = [0.0] * N; di = [0.0] * N; up = [0.0] * N; rhs = [0.0] * N
    di[0], up[0], rhs[0] = -1.0, 1.0, 0.0
    for i in range(1, N - 1):
        lo[i], di[i], up[i] = 1.0, -2.0, 1.0
        rhs[i] = -src * h * h
    lo[N - 1], di[N - 1], rhs[N - 1] = 0.0, 1.0, 0.0
    for i in range(1, N):
        w = lo[i] / di[i - 1]
        di[i] -= w * up[i - 1]
        rhs[i] -= w * rhs[i - 1]
    T = [0.0] * N
    T[N - 1] = rhs[N - 1] / di[N - 1]
    for i in range(N - 2, -1, -1):
        T[i] = (rhs[i] - up[i] * T[i + 1]) / di[i]
    return T
 
U, cp = 600.0, 1005.0
print("\nCouette shear layer, U=%.0f m/s, cp=%.0f J/kg/K, adiabatic wall at y=0" % (U, cp))
print("route            Pr    dT_wall[K]   Pr*U^2/(2cp)[K]")
rise = {}
for label, Pr in (("mean-free-path", pr_mean_free_path(1.4)),
                  ("Eucken", pr_eucken(1.4)),
                  ("measured", 0.707)):
    T = couette_wall_rise(Pr, U, cp)
    rise[label] = T[0]
    print("%-15s %5.3f  %10.1f   %14.1f" % (label, Pr, T[0], Pr * U * U / (2 * cp)))
 
print("\nwall-temperature gap = %.1f K   ratio = %.2f"
      % (rise["mean-free-path"] - rise["measured"], rise["mean-free-path"] / rise["measured"]))
 
print("\nk ~ sqrt(T) extrapolated from 300 K vs measured air conductivity")
k300 = 0.0263
for T_, k_meas in ((400, 0.0338), (600, 0.0469), (800, 0.0573), (1000, 0.0672)):
    k_law = k300 * math.sqrt(T_ / 300.0)
    print("T=%5d K   sqrtT=%.4f   measured=%.4f   error=%+6.1f%%"
          % (T_, k_law, k_meas, 100.0 * (k_law - k_meas) / k_meas))

gas   gamma   Pr_mfp  Pr_euck  Pr_meas   mfp/meas  euck/meas
He    1.667   1.667    0.667    0.680       2.45       0.98
Ar    1.667   1.667    0.667    0.670       2.49       0.99
N2    1.404   1.404    0.735    0.720       1.95       1.02
Air   1.400   1.400    0.737    0.707       1.98       1.04
CO2   1.289   1.289    0.781    0.770       1.67       1.01
H2O   1.329   1.329    0.764    0.980       1.36       0.78
 
Couette shear layer, U=600 m/s, cp=1005 J/kg/K, adiabatic wall at y=0
route            Pr    dT_wall[K]   Pr*U^2/(2cp)[K]
mean-free-path  1.400       249.5            250.7
Eucken          0.737       131.3            132.0
measured        0.707       126.0            126.6
 
wall-temperature gap = 123.5 K   ratio = 1.98
 
k ~ sqrt(T) extrapolated from 300 K vs measured air conductivity
T=  400 K   sqrtT=0.0304   measured=0.0338   error= -10.2%
T=  600 K   sqrtT=0.0372   measured=0.0469   error= -20.7%
T=  800 K   sqrtT=0.0429   measured=0.0573   error= -25.0%
T= 1000 K   sqrtT=0.0480   measured=0.0672   error= -28.5%

기체 선택과 Eucken 토글이 표의 세 열을 벽 온도 곡선으로 그린다. 토글을 끄면 Pr = $\gamma$ 곡선만 측정 곡선 옆에 남는다.

He를 고르면 두 포물선의 벽 끝 간격이 가장 크게 벌어진다. $\gamma$가 클수록 회전 자유도를 빠뜨린 대가가 커진다.

벽온도 123.5 K를 가른 것은 계수 9/4였다#

수치해와 해석해가 1 K 안쪽으로 맞는다. 격자 201점에서 $O(h^2)$ 오차다. 세 번째 열이 갈라 놓은 온도다. $U = 600$ m/s에서 단열 벽 온도는 Pr을 $\gamma$로 쓰면 249.5 K, 측정값으로 쓰면 126.0 K다. 123.5 K 차이다. 초음속 벽 열유속을 예측하는 계산에서 이 정도면 재료 선정이 달라진다.

Eucken 열은 모든 이원자·단원자 기체에서 측정값의 2~4% 안에 들어온다. 아령 하나를 병진과 회전으로 나눈 것만으로 2.45배 오차가 2%로 줄었다.

공기의 $\gamma=1.400$을 $4\gamma/(9\gamma-5)$에 넣으면 분자는 $4\times1.400=5.600$, 분모는 $9\times1.400-5=7.600$이다. 나누면 $0.7368$이라 표의 Eucken 값 $0.737$과 소수점 셋째 자리에서 맞다. 같은 식을 $1+(9/4)(\gamma-1)$로 풀면 $9/4\times0.400=0.900$이고, $1.400/1.900=0.737$이다. He는 $\gamma=1.667$이라 평균자유경로 경로는 $1.667$, Eucken은 $0.667$, 측정은 $0.680$이다. $1.667/0.680=2.45$라서 표의 mfp/meas와 같다. Eucken을 측정으로 나누면 $0.667/0.680=0.98$이다.

Couette에서 $U=600$, $c_p=1005$이면 $U^2/(2c_p)=360000/2010=179.1$ K다. 여기에 $1.400$을 곱하면 $250.7$ K, $0.707$을 곱하면 $126.6$ K, $0.737$을 곱하면 $132.0$ K다. 격자 201점의 수치해는 $249.5$, $126.0$, $131.3$이다. 해석해와 각각 $1.2$ K, $0.6$ K, $0.7$ K 차이라 1 K 안쪽이다. $249.5-126.0=123.5$이고 $249.5/126.0=1.98$이다. $300$ K의 $k=0.0263$을 $\sqrt{1000/300}$으로 밀면 $\sqrt{3.333}\times0.0263$이 $0.0480$ 근처다. 측정 $0.0672$와의 상대차는 $(0.0480-0.0672)/0.0672$로 약 $-28.6\%$라, 출력의 $-28.5\%$와 반올림 안에서 같다. H₂O는 $0.764/0.980=0.780$이라 표의 euck/meas = 0.78이고, 약 22% 낮다.

√T 법칙이 1000 K에서 28.5% 모자란 이유#

같은 유도는 온도 의존성도 준다. $\bar{c} \propto \sqrt{T}$이고 $\sigma$가 상수라면

마지막 출력 블록이 이 법칙을 300 K에서 외삽한 결과다. 400 K에서 10.2%, 1000 K에서 28.5% 낮게 나온다. 부호가 한쪽으로만 몰린다는 것이 힌트다. $\sigma$를 상수로 본 가정이 깨진 것이다. 빠른 분자는 서로를 더 가까이 스쳐 지나가므로 유효 단면적이 줄어든다. Chapman–Enskog 이론이 $\Omega^{(2,2)*}$ 충돌적분으로 채우는 자리이고, 실무에서는 Sutherland 식의 상수 $S$가 대신 흡수한다. 벽에서 마찰과 열을 한 비로 묶는 일은 레이놀즈 상사가 항력만으로 열전달을 맞히는 자리와 같다. 그 비에 넣는 Pr이 1.98배 틀리면 열유속이 같이 틀린다.

H₂O에서 Eucken이 22% 어긋나는 자리#

표에서 H₂O만 euck/meas = 0.78이다. 22% 낮다. 물 분자는 극성이고 쌍극자 상호작용이 충돌을 지배한다. 회전 완화가 병진보다 훨씬 느려서, 내부 에너지가 "속도와 무관하게 실려 간다"는 가정이 성립하지 않는다. 수정 Eucken 식은 내부 자유도에 $\rho D_{\text{int}}/\mu$를 따로 곱해 이 효과를 넣는다. 수증기 비율이 높은 연소 계산에서 Pr을 상수로 두면 이 차이가 벽 열유속에 그대로 실린다.

코드에 Pr = 0.71을 적을 때 실제로 서명하는 것#

Pr = 0.71 한 줄은 상수 선언이 아니다. 운동량과 에너지가 같은 분자에 실려 오는데도 서로 다른 속도로 확산한다는 주장이고, 그 비율이 회전 자유도 두 개에서 나온다는 주장이다. 그 주장이 깨지는 조건은 셋이다. 단원자(0.67로 내려간다), 극성 다원자(0.98까지 올라간다), 그리고 $\lambda$가 격자보다 커지는 희박 기체(첫 시뮬레이션의 n = 1)다. 이 셋 중 하나에 들어가는 계산이면, 상수 하나를 바꾸는 것으로 끝나지 않는다.