코리올리 힘이 한 일은 0인데 임펠러는 1800 J/kg을 올렸다

회전계에서 코리올리 힘의 일률은 정확히 0이다. 그래도 헤드가 생기는 이유는 그 힘이 블레이드에 반작용으로 걸리기 때문이다.

일의 정의를 만든 사람의 힘이 일을 하지 않는다#

역학에서 "일(work)"과 "운동에너지"를 지금 형태로 처음 정리한 사람은 가스파르 코리올리다. 1829년 저서의 제목이 『기계의 효과 계산에 대하여』였다. 수차(水車)를 놓고 얼마가 들어가 얼마가 나오는지를 세던 엔지니어의 책이다.

그런데 그의 이름이 붙은 힘은 일을 하지 않는다. 크기가 아무리 커도 그렇다. 이 글에서는 그 사실을 원심 임펠러 하나로 확인한다. 코리올리 가속도가 290 g에 달하는데도 장부에는 $10^{-16}$ J/kg이 찍히고, 유체의 총엔탈피는 1800 J/kg 올라간다. 그 1800이 어디서 오는지, 그리고 회전계 솔버가 그 항을 빠뜨리면 무엇이 사라지는지가 결론이다.

1832년, 수차 계산에서 회전 좌표계로#

『역사 속의 유체역학』 연재는 이 시기의 파리를 이렇게 적는다. 1830년 7월혁명으로 왕당파 코시가 추방되고 공화파 나비에가 에콜 폴리테크니크에 임명된다. 그리고 1832년, 나비에는 코리올리와 공동 연구를 시작한다.

코리올리가 그때까지 붙들고 있던 주제가 수차였다. 회전하는 유체기계에서 에너지와 일을 어떻게 세는가. 그 계산을 하려면 회전하는 축 위에 앉아야 했다. 회전 좌표계 연구가 거기서 나왔고, 1835년 논문에서 오늘날의 코리올리 힘이 세상에 알려진다.

순서가 중요하다. 코리올리 힘은 지구 자전이나 기상학에서 나온 것이 아니다. 회전 유체기계의 에너지 장부를 맞추려다 나왔다. 아래 시뮬레이션에서 그 장부를 직접 조작해보자.

omega [rad/s]300backsweep [deg]0

Flip the frame switch: the orange trail changes shape completely, yet the green rothalpy line never tilts. At the current radius 50 mm, I = 0.0 J/kg against h0 = 0.0 J/kg. Across the whole channel h0 climbs by 1800.0 J/kg — exactly U2·c_theta2 − U1·c_theta1 = 1800.0. Raise omega and only the orange line responds.

frame 버튼으로 카메라를 상대계와 절대계 사이에서 바꿔 보라. 주황색 궤적은 완전히 다른 모양이 되지만, 오른쪽 초록색 선은 기울지 않는다. omega 슬라이더를 올리면 주황색 $h_0$ 만 가팔라진다.

수직인 힘은 장부에 오르지 않는다#

각속도 $\vec{\Omega}$ 로 도는 좌표계에서, 상대속도 $\vec{w}$ 를 갖는 유체 입자의 운동방정식은 이렇게 된다.

\frac{D\vec{w}}{Dt} = -\frac{\nabla p}{\rho} - 2\vec{\Omega}\times\vec{w} - \vec{\Omega}\times(\vec{\Omega}\times\vec{r})

우변 두 번째가 코리올리 항, 세 번째가 원심력 항이다. 둘 다 좌표계를 회전시킨 대가로 생긴 관성항이다. 이 두 항이 운동량 방정식에 들어가는 방식은 회전 기준 프레임과 MRF에서 다뤘다. 여기서는 에너지 쪽만 본다.

에너지를 보려면 각 힘에 $\vec{w}$ 를 내적한다. 코리올리 항은 그 자리에서 죽는다.

(-2\vec{\Omega}\times\vec{w})\cdot\vec{w} = 0

외적의 결과는 두 인자 모두에 수직이기 때문이다. 항등식이므로 근사도 조건도 없다. $\Omega$ 가 얼마든, $\vec{w}$ 가 어느 방향이든 0이다.

원심력은 다르다. $-\vec{\Omega}\times(\vec{\Omega}\times\vec{r}) = \Omega^2 r\,\hat{r}$ 이고, 반경방향 속도가 있으면 내적이 살아남는다. 대신 이 항은 퍼텐셜을 갖는다.

\Omega^2 r\,\hat{r} = -\nabla\left(-\frac{U^2}{2}\right), \qquad U = \Omega r

퍼텐셜이 있다는 말은 방정식 왼쪽으로 옮겨 상수 안에 넣을 수 있다는 뜻이다.

그래서 보존되는 것은 $h_0$ 가 아니라 $I$ 다#

정상·비점성 조건에서 위 방정식을 유선을 따라 적분하면 나오는 것이 로탈피(rothalpy, rotational + enthalpy)다.

I = h + \frac{w^2}{2} - \frac{U^2}{2} = \text{const}

$h$ 는 정엔탈피, $w$ 는 상대속도 크기, $U=\Omega r$ 은 블레이드 원주속도다. 마지막 항의 마이너스가 원심력 퍼텐셜이다. 코리올리 항은 애초에 내적이 0이라 이 식에 흔적조차 남기지 않는다.

절대계의 총엔탈피 $h_0 = h + c^2/2$ 와는 이렇게 이어진다. 상대속도와 절대속도는 $\vec{c} = \vec{w} + \vec{U}$ 로 묶이므로,

h_0 = I + U c_\theta

가 된다. $c_\theta$ 는 절대속도의 선회 성분이다. 곧 $I$ 가 상수인 채로 $h_0$ 가 올라갈 수 있다. 올라가는 양이 정확히 $\Delta(U c_\theta)$ , 즉 오일러 터보기계 방정식이다.

Python으로 장부를 세 개 따로 적어봤다#

반경 50 mm에서 150 mm까지, 채널 높이가 20 mm에서 8 mm로 좁아지는 원심 임펠러를 잡는다. $\Omega = 300$ rad/s, 물 30 kg/s. 후곡각 $\beta$ 는 0°와 30° 두 경우를 본다. 연속식이 반경방향 상대속도를 정하고, 로탈피 일정 조건이 압력을 정한다. 그다음 코리올리·원심력·블레이드 반력이 하는 일을 각각 시간 적분한다.

import numpy as np
 
OMEGA = 300.0            # 각속도 [rad/s]
RHO = 1000.0             # 밀도 [kg/m^3]
MDOT = 30.0              # 질량유량 [kg/s]
R1, R2 = 0.05, 0.15      # 입·출구 반경 [m]
B1, B2 = 0.020, 0.008    # 입·출구 채널 높이 [m]
GRID = np.linspace(R1, R2, 4001)
 
 
def channel_state(r, beta_deg):
    """반경 r에서 상대속도 성분 (w_r, w_t), 블레이드 속도 U, 절대 선회속도 c_t."""
    b = B1 + (B2 - B1) * (r - R1) / (R2 - R1)
    w_r = MDOT / (RHO * 2.0 * np.pi * r * b)          # 연속식이 정하는 반경방향 성분
    w_t = -w_r * np.tan(np.radians(beta_deg))         # 후곡 블레이드 → 회전 반대쪽
    u = OMEGA * r
    return w_r, w_t, u, u + w_t
 
 
def march_channel(beta_deg):
    """로탈피를 일정하게 두고 채널을 따라 p/rho 와 절대 총엔탈피 h0 를 세운다."""
    w_r, w_t, u, _ = channel_state(R1, beta_deg)
    i_const = 0.5 * (w_r**2 + w_t**2) - 0.5 * u**2    # p1/rho = 0 을 기준으로
    cols = []
    for r in GRID:
        w_r, w_t, u, c_t = channel_state(r, beta_deg)
        p = i_const - 0.5 * (w_r**2 + w_t**2) + 0.5 * u**2
        cols.append((u, c_t, p + 0.5 * (w_r**2 + c_t**2),
                     p + 0.5 * (w_r**2 + w_t**2) - 0.5 * u**2))
    return np.array(cols)                              # U, c_t, h0, I
 
 
def power_ledger(beta_deg):
    """궤적을 따라 코리올리·원심력·블레이드 반력의 일률을 각각 시간 적분한다."""
    om = np.array([0.0, 0.0, OMEGA])
    st = np.array([channel_state(r, beta_deg) for r in GRID])
    w_r, w_t = st[:, 0], st[:, 1]
    dwt_dr = np.gradient(w_t, GRID)
    p_cor = np.zeros_like(GRID)
    p_cen = np.zeros_like(GRID)
    p_bld = np.zeros_like(GRID)
    for k, r in enumerate(GRID):
        w = np.array([w_r[k], w_t[k], 0.0])            # 국소 (r, theta, z) 정규직교
        f_cor = -2.0 * np.cross(om, w)
        f_cen = -np.cross(om, np.cross(om, np.array([r, 0.0, 0.0])))
        acc_t = w_r[k] * dwt_dr[k] + w_r[k] * w_t[k] / r   # 원통좌표 곡률항까지
        f_bld_t = acc_t - f_cor[1]                     # 남는 선회 가속은 압력면 몫
        p_cor[k] = f_cor @ w
        p_cen[k] = f_cen @ w
        p_bld[k] = OMEGA * r * f_bld_t                 # 토크 × 각속도 = 절대계 일률
    dt = 1.0 / w_r                                     # dt = dr / w_r
    ints = [np.trapezoid(p * dt, GRID) for p in (p_cor, p_cen, p_bld)]
    return ints + [np.max(np.abs(p_cor))]
 
 
def transported_h0(beta_deg, with_source):
    """h0 수송식을 직접 적분한다. 소스항은 Omega * d(r c_theta)/dr."""
    st = np.array([channel_state(r, beta_deg) for r in GRID])
    src = OMEGA * np.gradient(GRID * st[:, 3], GRID) if with_source else np.zeros_like(GRID)
    return np.trapezoid(src, GRID)
 
 
for beta in (0.0, 30.0):
    tab = march_channel(beta)
    d_h0 = tab[-1, 2] - tab[0, 2]
    euler = tab[-1, 0] * tab[-1, 1] - tab[0, 0] * tab[0, 1]
    drift = np.max(np.abs(tab[:, 3] - tab[0, 3]))
    cor, cen, bld, cor_peak = power_ledger(beta)
    ok = transported_h0(beta, True)
    bad = transported_h0(beta, False)
    print(f"beta = {beta:4.1f} deg   U1 = {tab[0,0]:4.1f} m/s   U2 = {tab[-1,0]:4.1f} m/s")
    print(f"  delta h0 from rothalpy = {d_h0:9.2f} J/kg")
    print(f"  U*c_theta (Euler)      = {euler:9.2f} J/kg   gap {abs(d_h0-euler):.1e}")
    print(f"  rothalpy max drift     = {drift:9.1e} J/kg")
    print(f"  work by Coriolis       = {cor:9.1e} J/kg   (peak power {cor_peak:.1e} W/kg)")
    print(f"  work by centrifugal    = {cen:9.2f} J/kg")
    print(f"  work by blade torque   = {bld:9.2f} J/kg")
    print(f"  h0 transport, source   = {ok:9.2f} J/kg   head {ok/9.81:5.1f} m")
    print(f"  h0 transport, dropped  = {bad:9.2f} J/kg   head {bad/9.81:5.1f} m")

beta =  0.0 deg   U1 = 15.0 m/s   U2 = 45.0 m/s
  delta h0 from rothalpy =   1800.00 J/kg
  U*c_theta (Euler)      =   1800.00 J/kg   gap 0.0e+00
  rothalpy max drift     =   1.1e-13 J/kg
  work by Coriolis       =   0.0e+00 J/kg   (peak power 0.0e+00 W/kg)
  work by centrifugal    =    900.00 J/kg
  work by blade torque   =   1800.00 J/kg
  h0 transport, source   =   1800.00 J/kg   head 183.5 m
  h0 transport, dropped  =      0.00 J/kg   head   0.0 m
beta = 30.0 deg   U1 = 15.0 m/s   U2 = 45.0 m/s
  delta h0 from rothalpy =   1737.98 J/kg
  U*c_theta (Euler)      =   1737.98 J/kg   gap 0.0e+00
  rothalpy max drift     =   7.1e-14 J/kg
  work by Coriolis       =   2.2e-16 J/kg   (peak power 1.4e-12 W/kg)
  work by centrifugal    =    900.00 J/kg
  work by blade torque   =   1737.98 J/kg
  h0 transport, source   =   1737.98 J/kg   head 177.2 m
  h0 transport, dropped  =      0.00 J/kg   head   0.0 m

로탈피 드리프트는 $10^{-13}$ J/kg, 즉 배정밀도 반올림 수준이다. 코리올리 장부는 후곡각이 0일 때 정확히 0.0, 30°일 때 $2.2\times10^{-16}$ 이다. 후자는 부동소수점 잔여이지 물리가 아니다.

900과 1800 — 나머지 절반은 누가 냈나#

숫자 두 개가 눈에 띈다. 원심력이 한 일은 900 J/kg, 그런데 총엔탈피는 1800 J/kg 올랐다. 정확히 두 배다.

\int \Omega^2 r\,w_r\,dt = \int_{r_1}^{r_2} \Omega^2 r\,dr = \frac{U_2^2 - U_1^2}{2}

이 900은 상대계 안에서만 도는 돈이다. 압력과 상대운동에너지로 갈라져 들어간다. 절대계에서 유체가 실제로 받은 1800과는 계정이 다르다.

나머지 900을 낸 것은 블레이드다. 반경방향 블레이드에서 상대속도는 선회 성분이 없다. 그러려면 코리올리 힘 $2\Omega w_r$ 을 무언가가 정확히 상쇄해야 한다. 그 무언가가 블레이드 압력면이다. 반작용으로 유체는 같은 크기의 힘을 받고, 그 힘의 토크가 $2\Omega r w_r$ 이다.

여기서 갈린다. 상대계에서 그 힘은 $\vec{w}$ 에 수직이라 일을 하지 않는다. 절대계에서는 블레이드가 $\Omega$ 로 돌고 있으므로 토크 × 각속도가 곧 일률이다. 적분하면 $2\times 900 = 1800$ J/kg. 코드의 work by blade torque 줄이 그것이다.

정리하면 이렇다. 코리올리 힘은 유체에 한 푼도 주지 않는다. 대신 유체가 블레이드를 밀게 만들고, 그 반작용 경로로 축의 일이 들어온다. 힘은 방향만 정하고 지불은 축이 한다.

omega [rad/s]300backsweep [deg]0

The orange Coriolis arrow is the longest one on screen — up to 292 g — and its account still reads 0.0e+0 J/kg. Meanwhile the blade account fills to 0 of 1800 J/kg. Drag backsweep and only the purple bar moves: the perpendicular force never pays, it only decides where the blade has to push.

주황색 코리올리 화살표는 화면에서 가장 긴데 장부는 0에 붙어 있다. backsweep 슬라이더를 움직이면 보라색 블레이드 계정만 값이 바뀐다. omega를 올려도 주황색 막대는 여전히 자라지 않는다.

회전계 솔버가 $h_0$ 를 그대로 수송하면#

여기까지가 물리다. 코드에서 사고가 나는 자리는 정해져 있다.

회전 격자 영역에서 에너지 방정식을 세울 때 무엇을 수송 변수로 잡느냐다. 상대속도로 수송되면서 소스 없이 보존되는 것은 $I$ 다. $h_0$ 를 상대속도로 실어 나르려면 소스항이 붙어야 한다.

\vec{w}\cdot\nabla h_0 = \Omega\,\frac{\partial (r c_\theta)}{\partial s}

$s$ 는 유선 방향 좌표다. 이 항을 빠뜨리면 출력의 마지막 두 줄이 된다. 소스가 있으면 1800 J/kg, 헤드 183.5 m. 없으면 0.00 J/kg, 헤드 0 m. 임펠러를 통과했는데 아무 일도 일어나지 않은 결과가 나온다.

증상이 헷갈리는 이유는 잔차가 멀쩡하기 때문이다. 방정식 자체는 잘 수렴한다. 다만 답이 헤드 0인 답으로 수렴한다. 격자를 조밀하게 해도 0은 0이다.

OpenFOAM 계열에서 MRF와 함께 압축성 에너지식을 쓸 때 rhothermo 쪽 총에너지 정의와 MRF 보정이 어긋나면 이 모양이 된다. 좌표변환에서 항이 하나 빠지는 문제라는 점에서, 극좌표 FVM의 곡률항에서 본 실수와 뿌리가 같다. 위 코드에서 acc_t 줄에 곡률항 $w_r w_\theta / r$ 을 넣은 이유도 같다. 그것을 빼면 후곡각 30°에서 블레이드 장부가 1802.06으로 어긋난다.

회전 격자를 만나면 먼저 물어볼 세 가지#

첫째, 이 케이스에서 상수인 것은 $h_0$ 인가 $I$ 인가. 회전 영역 안에서는 $I$ 다. 정지 영역으로 넘어가면 다시 $h_0$ 다. 인터페이스에서 무엇이 이어지는지 확인한다.

둘째, 에너지 방정식의 소스항이 실제로 걸려 있는가. 헤드가 이상하게 낮거나 정확히 0이면 이 항부터 본다.

셋째, 코리올리 항을 에너지 소스로 "추가"하고 있지 않은가. 그 항의 일률은 정확히 0이다. 무언가 넣어야 할 것 같아 보인다면, 필요한 것은 코리올리가 아니라 $\Omega\,\partial(rc_\theta)/\partial s$ 다.

코리올리는 수차의 효율을 세다가 회전 좌표계에 도달했다. 그가 만든 장부에서 자기 이름의 힘은 항상 0원이다. 회전하는 것을 계산할 때 그 0을 0으로 두는 것이, 190년 뒤 우리가 그 장부를 이어받는 방식이다.

코리올리 힘이 한 일은 0인데 임펠러는 1800 J/kg을 올렸다 — 회전계 에너지와 로탈피

일의 정의를 만든 사람의 힘이 일을 하지 않는다#

1832년, 수차 계산에서 회전 좌표계로#

수직인 힘은 장부에 오르지 않는다#

그래서 보존되는 것은 $h_0$ 가 아니라 $I$ 다#

Python으로 장부를 세 개 따로 적어봤다#

900과 1800 — 나머지 절반은 누가 냈나#

회전계 솔버가 $h_0$ 를 그대로 수송하면#

회전 격자를 만나면 먼저 물어볼 세 가지#

관련

코리올리 힘이 한 일은 0인데 임펠러는 1800 J/kg을 올렸다 — 회전계 에너지와 로탈피

일의 정의를 만든 사람의 힘이 일을 하지 않는다#

1832년, 수차 계산에서 회전 좌표계로#

수직인 힘은 장부에 오르지 않는다#

그래서 보존되는 것은 h0h_0h0​가 아니라 III다#

Python으로 장부를 세 개 따로 적어봤다#

900과 1800 — 나머지 절반은 누가 냈나#

회전계 솔버가 h0h_0h0​를 그대로 수송하면#

회전 격자를 만나면 먼저 물어볼 세 가지#

관련

그래서 보존되는 것은 $h_0$ 가 아니라 $I$ 다#

회전계 솔버가 $h_0$ 를 그대로 수송하면#