임펠러를 멈춰 세우고 푸는 법 — 회전 기준 프레임과 MRF

끌림 속도 1.2를 빼지 않으면 절대 궤적의 x가 1에서 벗어난다. 계면에서 같은 항을 빼먹으면 질량이 샌다.

회전하는 임펠러를 정지시킨 채로 해석한다. 그런데 답은 맞다. 펌프나 팬의 정상 상태 성능을 구할 때 가장 많이 쓰는 트릭이다. 회전 기준 프레임에서 운동량 방정식이 어떻게 바뀌는지, MRF가 회전을 어떻게 얼려 정상 해석으로 만드는지를 수식과 한 입자로 따라간다.

프레임 토글을 절대에서 상대(MRF)로 넘기면, 회전 영역 안쪽 화살표가 $\boldsymbol{\Omega}\times\mathbf{r}$ 강체 회전에서 0 근처로 가라앉는다.

omega1.00

Switch to the relative frame: the inner cyan arrows collapse toward zero because the rotating zone becomes steady, while the pink outer field is unchanged. The dashed circle is where the two velocity descriptions must be transformed into each other.

바깥 분홍 장은 그대로다. 인터페이스로 갈수록 절대 프레임의 화살표만 커진다. 상대 프레임은 그 끌림을 방정식 속으로 집어넣은 상태다.

절대 좌표계를 버리면 풍차가 멈춘다#

땅에 박힌 카메라로 보면 임펠러 날개는 빠르게 돈다. 시간에 따라 형상이 계속 변한다. 그래서 비정상(unsteady) 해석이 필요하다.

카메라를 날개 위에 올리면 날개는 정지해 있다. 주변 유체만 상대적으로 흐른다. 형상이 안 변하니 정상(steady) 해석이 가능해진다. 회전체 문제를 정지 문제로 바꾸는 이 발상이 회전 기준 프레임의 전부다.

대가는 공짜가 아니다. 좌표계 자체가 가속하기 때문에, 운동량 방정식에 원래 없던 항이 두 개 추가된다.

운동량 방정식에 끼어든 두 가짜 힘#

회전 각속도를 $\boldsymbol{\Omega}$ , 상대 속도를 $\mathbf{u}_r$ , 위치를 $\mathbf{r}$ 라 하자. 절대 속도와의 관계는 다음과 같다.

\mathbf{u}_a = \mathbf{u}_r + \boldsymbol{\Omega}\times\mathbf{r}

여기서 $\mathbf{u}_a$ 는 절대(관성) 좌표계 속도, $\boldsymbol{\Omega}\times\mathbf{r}$ 는 회전에 의한 끌림 속도다.

이 관계를 Navier–Stokes에 대입하면 상대 속도에 대한 정상 운동량 방정식이 나온다.

(\mathbf{u}_r\cdot\nabla)\mathbf{u}_r + 2\,\boldsymbol{\Omega}\times\mathbf{u}_r + \boldsymbol{\Omega}\times(\boldsymbol{\Omega}\times\mathbf{r}) = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u}_r

왼쪽에서 두 번째 항 $2\,\boldsymbol{\Omega}\times\mathbf{u}_r$ 가 Coriolis 가속도(운동 방향을 꺾는 항), 그 다음 항 $\boldsymbol{\Omega}\times(\boldsymbol{\Omega}\times\mathbf{r})$ 가 원심 가속도(축에서 바깥으로 미는 항)다. 두 항 모두 실제 힘이 아니다. 좌표계가 가속하기 때문에 나타나는 겉보기 항이다.

Coriolis 항은 운동량 방정식의 소스로 들어간다. 그래서 OpenFOAM의 MRFSource는 이 둘을 셀마다 더해주는 일을 한다. 코리올리가 한 일과 임펠러가 올린 에너지는 같은 항이 아니다. 그 갈림은 회전계 에너지와 로탈피에 따로 적혀 있다.

MRF — 회전을 정지로 얼리는 frozen rotor#

MRF는 영역을 둘로 나눈다. 회전체를 감싸는 회전 영역(rotating zone)과 그 바깥의 정지 영역(stationary zone).

회전 영역: 위의 상대 속도 방정식을 푼다. Coriolis·원심 소스가 켜진다.
정지 영역: 보통의 절대 좌표계 방정식을 푼다.

날개는 실제로 움직이지 않는다. 메쉬도 고정이다. 단지 회전 영역의 방정식에 회전 효과를 소스로 주입할 뿐이다. 그래서 이 방식을 frozen rotor(얼어붙은 로터)라 부른다. 날개와 주변 유체의 상대 위치가 한 순간에 동결된 채 정상 해를 구한다.

앞에서 토글을 넘긴 화면이 그 동결이다. 절대 프레임에서는 안쪽 유체가 $\boldsymbol{\Omega}\times\mathbf{r}$ 만큼 강체 회전한다. 화살표가 인터페이스로 갈수록 커진다. 상대(MRF) 프레임으로 바꾸면 안쪽 화살표가 0 근처로 주저앉는다. 회전 영역의 흐름이 정지한 것처럼 보이기 때문이다. 바깥 분홍 장은 그대로다.

인터페이스에서 속도를 다시 꿰매다#

두 영역은 서로 다른 좌표계로 푼다. 그러면 경계에서 속도가 안 맞는다. 회전 영역은 $\mathbf{u}_r$ , 정지 영역은 $\mathbf{u}_a$ 로 저장하니까.

인터페이스에서는 위의 변환식을 그대로 적용한다. 회전 영역에서 정지 영역으로 넘어갈 때 $\mathbf{u}_a = \mathbf{u}_r + \boldsymbol{\Omega}\times\mathbf{r}$ 로 더해주고, 반대 방향이면 빼준다. 압력과 난류량 같은 스칼라는 그대로 통과한다. 벡터만 끌림 속도를 보정한다.

이 변환이 빠지면 인터페이스에서 질량과 운동량이 새기 시작한다. 잔차가 떨어지지 않거나, 인터페이스를 가로지르는 비물리적 제트가 생긴다. MRF 설정에서 가장 흔한 실수가 회전 영역 경계를 날개에 너무 바짝 붙이는 것이다. 인터페이스는 유동이 거의 축대칭이 되는 자리에 둬야 변환 오차가 작다.

MRF와 sliding mesh, 무엇을 언제#

MRF는 정상 상태 근사다. 날개와 볼류트(volute) 같은 정지 구조의 상대 위치가 시간에 따라 변하는 효과(rotor–stator 상호작용)를 잡지 못한다. 한 순간의 위치로 동결하기 때문이다.

진짜 비정상 효과가 필요하면 sliding mesh(미끄럼 격자)로 간다. 메쉬가 실제로 회전하고, 매 스텝 인터페이스를 다시 보간한다. 비싸지만 그 상대 운동을 남긴다.

항목	MRF (frozen rotor)	Sliding mesh
메쉬 모션	없음 (고정)	실제 회전
시간	정상 상태	비정상
비용	낮음	높음
rotor–stator 상호작용	근사 못 함	잡아냄
용도	설계 점검, 성능 곡선	소음, 압력 맥동

실무 흐름은 보통 이렇다. MRF로 빠르게 작동점을 잡고, 그 해를 초기장으로 sliding mesh를 돌려 비정상 디테일을 본다.

Python — 회전 좌표계 속 자유 입자의 궤적#

아무 힘도 안 받는 자유 입자로 소스 항을 확인한다. 절대 좌표계에서는 직선으로 날아가야 한다. 회전 좌표계에서 Coriolis·원심 항만으로 적분한 뒤, 다시 절대 좌표계로 되돌렸을 때 직선이 복원되면 소스 항의 부호가 맞다.

import numpy as np
 
OMEGA = np.array([0.0, 0.0, 1.2])   # 회전 각속도 벡터 (rad/s), z축 둘레
 
def rotating_frame_acceleration(r, v_rel):
    """자유 입자(실제 힘 0)가 회전 좌표계에서 받는 겉보기 가속도.
    Coriolis + 원심 항만 남는다."""
    coriolis = -2.0 * np.cross(OMEGA, v_rel)              # -2 Ω × u_r
    centrifugal = -np.cross(OMEGA, np.cross(OMEGA, r))    # -Ω × (Ω × r)
    return coriolis + centrifugal
 
def integrate_rotating_particle(r0, v0, dt, steps):
    """velocity-Verlet 변형으로 상대 좌표계 궤적을 적분."""
    r = np.array(r0, float); v = np.array(v0, float)
    traj = [r.copy()]
    a = rotating_frame_acceleration(r, v)
    for _ in range(steps):
        r = r + v * dt + 0.5 * a * dt * dt
        a_half = rotating_frame_acceleration(r, v)        # 속도 의존 -> 1회 보정
        v = v + 0.5 * (a + a_half) * dt
        a = rotating_frame_acceleration(r, v)
        traj.append(r.copy())
    return np.array(traj)
 
def relative_to_absolute(traj_rel, dt):
    """상대 좌표계 궤적을 절대 좌표계로 되돌린다: 각 시각 Ωt 만큼 역회전."""
    w = OMEGA[2]
    out = []
    for k, r in enumerate(traj_rel):
        t = k * dt
        c, s = np.cos(w * t), np.sin(w * t)
        out.append([c * r[0] - s * r[1], s * r[0] + c * r[1], r[2]])
    return np.array(out)
 
# 절대 좌표계에서 직선 운동하도록 초기조건을 잡는다
V = 0.8
r0 = np.array([1.0, 0.0, 0.0])
v_abs0 = np.array([0.0, V, 0.0])                # 절대 속도: +y 방향 직선
v_rel0 = v_abs0 - np.cross(OMEGA, r0)           # u_r = u_a - Ω × r
 
dt = 0.002
traj_rel = integrate_rotating_particle(r0, v_rel0, dt, 1500)
traj_abs = relative_to_absolute(traj_rel, dt)
 
x_dev = np.abs(traj_abs[:, 0] - 1.0).max()      # 직선이면 x ≈ 1.0 유지
print(f"절대 좌표계 x 최대 이탈: {x_dev:.4e}")
print(f"절대 좌표계 끝점: {traj_abs[-1].round(3)}")

초기 조건은 코드에 박혀 있다. $\boldsymbol{\Omega}=(0,0,1.2)$ rad/s, $\mathbf{r}_0=(1,0,0)$ . 끌림 속도 $\boldsymbol{\Omega}\times\mathbf{r}$ 의 크기는 $\Omega r=1.2$ 다. 절대 속도는 $V=0.8$ 로 $+y$ 방향이라 $\mathbf{u}_a=(0,0.8,0)$ . 상대 속도는 그 차다. $y$ 성분만 $0.8-1.2=-0.4$ . v_rel0 한 줄이 이 뺄셈이다.

$t=0$ 에서 코리올리 가속도의 크기를 식에서 읽으면 $|2\boldsymbol{\Omega}\times\mathbf{u}_r|=2\times 1.2\times 0.4=0.96$ 이다. 방향은 $-x$ . 원심 가속도의 크기는 $\Omega^2 r=1.44$ 이고, 코드의 부호는 축 바깥인 $+x$ 다. 운동량 방정식 좌변의 두 항을 우변으로 넘기면 코드의 마이너스와 맞는다.

시간 간격은 $dt=0.002$ , 걸음은 1500. 적분 시간은 3초다. 그 사이 회전판은 $\Omega t=1.2\times 3=3.6$ rad를 돈다. 절대 궤적의 $x$ 가 1 근처에 남으면 이 부호가 맞다는 뜻이다. 부호를 하나라도 뒤집으면 입자가 바깥으로 튕기며 $x$ 이탈이 커진다.

적분은 velocity-Verlet 변형이다. 위치 갱신에 $0.5\,a\,dt^2$ 가 들어가고, 가속도가 속도에 의존하므로 속도를 한 번 더 고쳐 넣는다. 되돌리기 회전각은 각 시각 $wt$ 이고, $w$ 는 $\Omega$ 의 $z$ 성분 1.2다. 한 걸음의 각은 $1.2\times 0.002=0.0024$ rad다. 1500걸음이 3.6 rad를 채운다. 절대 좌표에서 $x$ 가 1에서 떨어지는지는, 그 회전을 되돌린 뒤의 최댓값으로 본다.

Ω 슬라이더를 올리면 회색 직선은 그대로고, 시안 궤적만 나선으로 감긴다.

omega0.90

At omega = 0 both paths coincide. Increase omega and the cyan trail curls into a spiral — the Coriolis deflection that an observer riding the turntable mistakes for a force.

$\Omega=0$ 이면 둘이 포개진다. 슬라이더가 코드의 1.2 rad/s에 오면, 3초 동안 3.6 rad를 돈 관측자가 보는 휘어짐이 시안이다. 회색이 곧게 남는지가 소스 부호의 판정이다.

코드로 옮기기 전 짚을 함정#

회전축과 $\boldsymbol{\Omega}$ 방향을 도면과 한 번 더 대조한다. 부호가 반대면 추력이 정반대로 나온다. 발산하지 않고 그럴듯하게 틀린 답이 나오기 때문에 잡기 어렵다.

인터페이스는 유동이 부드러운 곳에 둔다. 날개 끝 와류 한가운데를 가르면 변환 오차가 커진다. 회전 영역은 날개보다 충분히 크게 잡는다.

MRF 해를 마지막 답으로 두지 않는다. 정상 상태 근사다. rotor–stator 맥동이 중요한 문제면 sliding mesh로 한 번 더 확인한다.

부호가 뒤집히면 직선이 먼저 죽는다#

자유 입자는 힘이 없다. 절대 좌표에서 $x=1$ 을 유지해야 한다. v_rel0에서 $\boldsymbol{\Omega}\times\mathbf{r}$ 를 빼지 않거나, Coriolis와 원심의 부호를 하나만 뒤집으면 그 직선이 먼저 무너진다. 3초, 3.6 rad를 도는 동안 $x$ 이탈이 커지면 소스 항을 의심한다. 추력 곡선은 그 다음에 틀린다.

임펠러 계산으로 옮길 때도 순서는 같다. 축 방향, 계면의 $\mathbf{u}_a=\mathbf{u}_r+\boldsymbol{\Omega}\times\mathbf{r}$ , 그 다음 성능 곡선이다. 입자에서 직선이 복원되지 않은 부호를 셀 소스에 넣으면, 잔차는 내려도 추력만 반대로 나온다. 그 잔차가 압력 격자 무늬로 남으면 SIMPLE의 압력-속도 연계와 같은 자리를 의심한다.