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

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

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

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

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

아래 시뮬레이션에서 직접 조작해보자.

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$이라 열 층이 훨씬 두껍다.

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

두 무차원 방정식을 다시 보자. 만약 압력 구배가 없고($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$이면 원래 식으로 돌아간다. 아래 그래프에서 직접 조작해보자.

주황색 곡선(생짜 $St$)은 $Pr$에 따라 위아래로 움직인다. 그러나 청록색 곡선($St \cdot Pr^{2/3}$)은 $Pr$이 무엇이든 점선($C_f/2$) 위에 딱 붙어 있다. 보정 지수 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$다.

기억할 점#

• 속도 경계층과 열 경계층은 같은 모양의 무차원 방정식을 따른다. 차이는 확산 계수 $Re$와 $Re \cdot Pr$ 뿐이다.
• $Pr = 1$, 압력 구배 0이면 속도장과 온도장이 포개지고 $C_f/2 = St$다. 이것이 레이놀즈 상사다.
• 일반적인 유체는 $C_f/2 = St \cdot Pr^{2/3}$의 Chilton-Colburn 상사로 마찰에서 열전달을 추정한다.