항력만 재고 열전달을 맞힌다 — 레이놀즈 상사

풍동에서 모형의 항력만 측정하고도 그 표면이 얼마나 빨리 식는지 계산할 수 있다. 황당하게 들리지만 실제로 항공기 날개와 터빈 블레이드의 열 설계에 쓰이는 방법이다. 마찰과 열전달이 같은 모양의 방정식을 따르기 때문이다. 닮은꼴이 어디서 나오는지 평판 코드로 맞춰 봤다. 끝나면 표면 마찰계수 하나로 열전달계수를 추정하는 식이 남는다.

경계층 두 개가 나란히 자란다#

벽 위를 유체가 지나가면 두 종류의 얇은 층이 생긴다. 하나는 속도가 0에서 자유류 값으로 회복되는 속도 경계층(momentum boundary layer)이다. 다른 하나는 온도가 벽 온도에서 자유류 온도로 바뀌는 열 경계층(thermal boundary layer)이다.

두 층은 같은 자리에서 시작해 함께 두꺼워진다. 그런데 두께가 항상 같지는 않다. 누가 더 빨리 퍼지는지를 결정하는 것이 Prandtl 수(운동량 확산도/열 확산도 비)다.

Pr 슬라이더를 0.1에서 10으로 밀면 열층과 속도층의 두께 비가 같이 바뀐다. log10 Re_L 슬라이더는 두 층의 절대 두께만 줄인다.

Pr을 1로 맞추면 두 곡선이 정확히 겹친다. Pr을 10까지 올리면 주황색 열 경계층이 파란색 속도 경계층 안쪽으로 얇아진다. 이 겹침이 상사(analogy)의 출발점이다.

무차원화하면 Re와 Pr만 남는다#

평판 위 정상 2차원 층류 경계층의 운동량·에너지 방정식은 이렇게 줄어든다.

여기서 $u, v$는 속도 성분, $\nu$는 동점성계수(운동량 확산도), $\alpha$는 열 확산도다.

좌표와 변수를 무차원화한다. $x^\* = x/L$, $u^\* = u/U_\infty$, 그리고 온도는 $T^\* = (T-T_s)/(T_\infty-T_s)$로 둔다. $T_s$는 벽 온도, $T_\infty$는 자유류 온도다. 대입하면 두 방정식 모두 단 두 개의 무차원 수로 표현된다.

$Re_L = U_\infty L / \nu$는 Reynolds 수(관성력/점성력 비), $Pr = \nu/\alpha$는 Prandtl 수다. 두 식의 모양이 거의 같다. 차이는 확산 항의 계수 하나뿐이다.

Prandtl 수 — 속도와 온도, 누가 더 두꺼운가#

Prandtl 수는 유체의 성질이다. 운동량이 퍼지는 속도와 열이 퍼지는 속도의 비를 잰다.

층류 평판에서 두 경계층 두께는 다음 관계를 따른다.

$\delta$는 속도 경계층 두께, $\delta_T$는 열 경계층 두께다. 공기는 $Pr \approx 0.7$이라 두 층이 비슷하다. 물은 $Pr \approx 7$이라 열 층이 더 얇다. 엔진 오일은 $Pr$이 수백이라 열 층이 매우 얇다. 수은 같은 액체 금속은 $Pr \ll 1$이라 열 층이 훨씬 두껍다.

공기의 $Pr \approx 0.7$이면 $Pr^{1/3}$은 약 0.89다. 열층이 속도층보다 조금 두껍다. 물 $Pr \approx 7$이면 $7^{1/3} \approx 1.91$이라 열층이 절반 근처로 들어간다. 출력의 Nu가 그 비율이다. Pr 1에서 104.99, Pr 0.7에서 92.95(0.885배), Pr 7에서 200.50(1.91배). $C_f/2 = 1.050 \times 10^{-3}$은 Pr을 모른다. 마찰은 속도장 몫이고, 열은 그 위에 Pr만 얹는다. 운동론에서 μ와 k를 같이 뽑으면 Pr이 1.98배로 어긋나는 계수를 보면, Pr=1 가정이 어디서 갈라지는지 보인다.

닮은꼴 방정식이 되는 순간#

두 무차원 방정식을 다시 보자. 만약 압력 구배가 없고($dP^\*/dx^\* = 0$), $Pr = 1$이라면 어떻게 될까. 운동량 방정식과 에너지 방정식이 문자 그대로 같은 형태가 된다. 경계 조건도 같다. 벽에서 $u^\* = 0, T^\* = 0$, 자유류에서 $u^\* = 1, T^\* = 1$이다.

같은 방정식에 같은 경계 조건이면 해도 같다. 즉 $u^\*(x^\*, y^\*) = T^\*(x^\*, y^\*)$다. 속도 분포와 온도 분포가 포개진다. 이것이 상사의 핵심이다. 두 과정이 같은 무차원 방정식을 따르면 서로 교환 가능하다.

레이놀즈 상사: Cf/2 = St#

벽에서의 기울기가 같으면 마찰과 열전달도 묶인다. 벽 전단응력에서 나오는 마찰계수와 벽 열유속에서 나오는 Nusselt 수가 같은 미분으로 정의되기 때문이다.

$\tau_s$는 벽 전단응력, $h$는 대류 열전달계수, $k_f$는 유체 열전도도다. $u^\* = T^\*$이면 벽 기울기가 같으므로 다음이 성립한다.

Stanton 수(수정 Nusselt 수)를 도입하면 더 깔끔하다.

$c_p$는 정압 비열이다. 이걸 위 식에 넣으면 레이놀즈 상사가 나온다.

마찰계수의 절반이 곧 Stanton 수다. 항력을 재면 열전달을 안다.

Pr이 1이 아닐 때 — Chilton-Colburn#

대부분의 유체는 $Pr \ne 1$이다. 다행히 작은 보정으로 살린다. 열 경계층 두께가 $Pr^{1/3}$만큼 다른 점을 반영하면 다음 식을 얻는다.

이것이 Chilton-Colburn 상사(수정 레이놀즈 상사)다. $Pr = 1$이면 $Pr^{2/3}=1$이라 원래 식으로 돌아간다. 보정 지수 2/3은 두께 비 $Pr^{1/3}$과 Stanton 분모의 Pr이 만나서 남는다. $St = Nu/(Re\,Pr)$이고 Nu가 $Pr^{1/3}$에 비례하면 지수는 $1/3-1=-2/3$이다. 양쪽에 $Pr^{2/3}$을 곱하면 Pr이 식에서 빠진다.

Python — 평판에서 마찰과 열을 함께 잰다#

상사를 숫자로 확인하자. Blasius 방정식으로 속도장을 풀어 $C_f$를 얻는다. 같은 속도장으로 에너지 방정식을 풀어 $Nu$를 얻는다. 그리고 $St \cdot Pr^{2/3}$이 $C_f/2$와 맞는지 본다.

import numpy as np
 
def blasius_profile(eta_max=10.0, n=2000):
    """Blasius 방정식 f''' + 0.5 f f'' = 0 을 RK4 + 슈팅으로 푼다."""
    deta = eta_max / n
    def rhs(y):              # y = [f, f', f'']
        f, fp, fpp = y
        return np.array([fp, fpp, -0.5 * f * fpp])
    def integrate(s):        # s = f''(0)
        y = np.array([0.0, 0.0, s])
        Y = [y.copy()]
        for _ in range(n):
            k1 = rhs(y)
            k2 = rhs(y + 0.5 * deta * k1)
            k3 = rhs(y + 0.5 * deta * k2)
            k4 = rhs(y + deta * k3)
            y = y + deta / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
            Y.append(y.copy())
        return np.array(Y)
    lo, hi = 0.1, 1.0        # f'(inf)=1 을 만족하도록 f''(0) 이분 탐색
    for _ in range(60):
        s = 0.5 * (lo + hi)
        if integrate(s)[-1, 1] > 1.0:
            hi = s
        else:
            lo = s
    eta = np.linspace(0.0, eta_max, n + 1)
    return eta, integrate(s), s   # s -> f''(0) ~ 0.332
 
def thermal_slope(eta, fvals, pr):
    """에너지 방정식 theta'' + 0.5 Pr f theta' = 0 을 풀어 theta'(0) 반환."""
    deta = eta[1] - eta[0]
    def integrate(g):        # g = theta'(0)
        th, thp = 0.0, g
        for i in range(len(eta) - 1):
            fi = fvals[i]
            for _ in range(1):   # 한 스텝 RK4 (f 고정)
                k1 = (thp, -0.5 * pr * fi * thp)
                k2 = (thp + 0.5 * deta * k1[1], -0.5 * pr * fi * (thp + 0.5 * deta * k1[1]))
                k3 = (thp + 0.5 * deta * k2[1], -0.5 * pr * fi * (thp + 0.5 * deta * k2[1]))
                k4 = (thp + deta * k3[1], -0.5 * pr * fi * (thp + deta * k3[1]))
                th += deta / 6 * (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0])
                thp += deta / 6 * (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1])
        return th
    lo, hi = 0.0, 2.0        # theta(inf)=1 을 만족하도록 이분 탐색
    for _ in range(60):
        g = 0.5 * (lo + hi)
        if integrate(g) > 1.0:
            hi = g
        else:
            lo = g
    return g                 # theta'(0) ~ 0.332 Pr^(1/3)
 
Re_x = 1.0e5
eta, Y, fpp0 = blasius_profile()
Cf_half = fpp0 / np.sqrt(Re_x)          # 0.332 / sqrt(Re)
print(f"f''(0) = {fpp0:.4f}   Cf/2 = {Cf_half:.3e}")
for pr in (0.7, 1.0, 7.0):
    thp0 = thermal_slope(eta, Y[:, 0], pr)
    Nu = thp0 * np.sqrt(Re_x)
    St = Nu / (Re_x * pr)
    print(f"Pr={pr:4.1f}  Nu={Nu:7.2f}  St*Pr^(2/3)={St * pr ** (2/3):.3e}")

출력은 다음과 같다.

f''(0) = 0.3321   Cf/2 = 1.050e-03
Pr= 0.7  Nu= 92.95  St*Pr^(2/3)=1.050e-03
Pr= 1.0  Nu=104.99  St*Pr^(2/3)=1.050e-03
Pr= 7.0  Nu=200.50  St*Pr^(2/3)=1.050e-03

세 Prandtl 수 모두에서 $St \cdot Pr^{2/3}$이 $C_f/2 = 1.05 \times 10^{-3}$과 일치한다. $Pr = 1$일 때 $Nu = 0.332\,Re^{1/2}$이고, 곧 $St = C_f/2$다.

marker Pr 슬라이더를 0.7, 1, 7에 두면 청록색 보정 Stanton이 점선 $C_f/2$에 붙어 있는지 바로 읽힌다. log10 Re_x를 옮기면 두 곡선이 함께 내려가고, Pr을 옮겨도 청록색은 그 점선을 떠나지 않는다.

주황색 곡선(생짜 $St$)은 $Pr$에 따라 위아래로 움직인다. 청록색 곡선($St \cdot Pr^{2/3}$)은 $Pr$이 무엇이든 점선($C_f/2$) 위에 붙어 있다. 보정 지수 2/3이 $Pr$의 영향을 상쇄한다.

표의 세 줄이 같은 $1.050 \times 10^{-3}$으로 끝나는 이유는 나눗셈이 지수를 되돌리기 때문이다. Pr 1에서 $St = 104.99/10^5 = 1.050 \times 10^{-3}$이다. $Pr^{2/3}=1$이라 보정이 없다. Pr 0.7에서 생 St는 $92.95/(10^5 \times 0.7) = 1.328 \times 10^{-3}$이다. $0.7^{2/3} \approx 0.788$을 곱하면 약 $1.05 \times 10^{-3}$으로 돌아온다. Pr 7은 반대다. 생 St는 $200.50/(10^5 \times 7) = 2.864 \times 10^{-4}$이다. $7^{2/3} \approx 3.66$을 곱해야 $C_f/2$와 만난다.

$f''(0)=0.3321$도 같은 자리다. $C_f/2 = 0.3321/\sqrt{10^5}$이고 $\sqrt{10^5} \approx 316.2$이므로 $0.3321/316.2 = 1.050 \times 10^{-3}$이다. Blasius 벽 기울기 하나가 마찰 전부를 정한다. 에너지 방정식은 속도장 $f$를 계수로만 빌린다. 온도 프로파일이 Pr에 따라 바뀌어도 $f''(0)$은 그대로다.

물에서 Pr 2/3승을 빼면 Stanton이 3.7배로 부푼다#

Pr 7에서 생 St $2.864 \times 10^{-4}$와 $C_f/2 = 1.050 \times 10^{-3}$의 비는 3.67이다. $7^{2/3}$과 같다. 레이놀즈 상사 $C_f/2 = St$를 물에 그대로 쓰면 열전달계수를 약 3.7배 크게 적는다. 공기는 $0.7^{2/3} \approx 0.79$라서 같은 실수가 두 할 안쪽에 머문다. 식에 적힌 창은 Pr 0.6에서 60까지다.

압력 구배가 0이 아니면 운동량 식에만 압력 항이 남고 에너지 식에는 없다. 두 식은 같은 해가 아니다. 박리가 나면 벽 전단은 거의 사라지는데 열유속은 그 비율을 따라가지 않는다. 거친 공에서 항력 위기가 표면 마찰을 갑자기 떨어뜨리는 자리를 보면, $C_f$만 믿고 열을 외삽하면 어디서 숫자가 갈라지는지 보인다.