격자 볼츠만의 벽은 격자점 위에 없다 — Bounce-back과 Zou-He 경계조건

LBM에서 no-slip 벽을 만드는 튕김 규칙과 벽이 놓이는 진짜 위치

똑같은 튕김 규칙 하나로 no-slip 벽을 만든다. 그런데 그렇게 세운 벽은 우리가 격자점을 찍은 그 자리에 있지 않다. 반 칸 어긋난 곳에 있다. 이 반 칸을 모르고 채널 높이를 세면 벽 마찰이 조용히 몇 퍼센트씩 틀어진다.

격자 볼츠만법(LBM, 분포함수를 격자 위에서 흘리고 충돌시키는 방법)은 Navier–Stokes를 직접 풀지 않는다. 미시 분포함수 $f_i$ 를 다룬다. 그래서 경계조건도 "속도를 벽에 대입한다"가 아니다. 벽으로 들어와야 할 분포함수를 어떻게 채우느냐의 문제다. 이 글은 그 채우는 규칙 세 가지 — bounce-back, 이동벽 보정, Zou-He — 를 다룬다. 그리고 벽의 진짜 위치를 Couette 유동으로 직접 확인한다.

벽에서 분포함수가 비어 있다#

한 스텝은 충돌과 스트리밍으로 나뉜다. 스트리밍은 각 $f_i$ 를 속도 방향 $\mathbf{e}_i$ 로 한 칸 민다. 내부 노드는 문제없다. 이웃이 사방에 있으니 들어올 분포함수가 다 온다.

벽에 붙은 노드는 다르다. 벽 바깥에는 유체 노드가 없다. 그쪽에서 들어와야 할 분포함수는 출처가 없다. D2Q9(2차원 9속도) 격자에서 바닥 벽 노드는 위로 향하는 세 방향 $f_2, f_5, f_6$ 이 스트리밍 뒤 비어 버린다.

이 세 개를 무엇으로 채우느냐가 벽의 물리를 정한다. 그냥 두면 유체가 벽을 통과해 새는 것과 같다. 벽이 없는 것처럼 행동한다.

튕겨 되돌려 no-slip을 만든다#

가장 단순한 답이 bounce-back(튕김)이다. 벽에 부딪힌 분포함수는 온 방향 그대로 되돌아간다. 반대 방향으로 뒤집어 되돌린다.

f_{\bar{i}}(\mathbf{x}_b,\, t+1) = f_i^{*}(\mathbf{x}_b,\, t)

$f_i^{*}$ 는 충돌 직후 값, $\bar{i}$ 는 $\mathbf{e}_{\bar{i}} = -\mathbf{e}_i$ 인 반대 방향이다. 들어온 운동량을 크기 그대로 반대로 돌려보낸다. 벽에서 알짜 속도가 0이 된다. 이것이 정지벽 no-slip이다.

아래 그림에서 바닥 노드의 미지 방향과 벽의 위치를 직접 바꿔보자.

파란 화살표 세 개가 스트리밍 뒤 비는 분포함수다. 버튼을 half-way로 두면 벽선이 노드에서 반 칸 아래로 내려간다. full-way로 두면 벽선이 노드 위로 올라온다.

벽은 반 칸 어긋나 있다#

여기가 핵심이다. bounce-back에는 두 방식이 있다.

full-way는 벽 노드에서 한 스텝 동안 모든 방향을 뒤집는다. 벽이 노드 바로 그 자리에 놓인다. 구현이 쉽지만 공간 1차 정확도다.

half-way는 유체 노드에서 스트리밍 도중 되돌린다. 부딪힘이 두 노드 사이에서 일어난다. 그래서 벽은 마지막 유체 노드에서 정확히 반 칸 바깥에 선다. 대신 공간 2차 정확도다.

이 반 칸이 실무의 함정이다. half-way를 쓰면서 채널 높이를 노드 개수 $n_y$ 가 아니라 $n_y-1$ 로 세면, 유효 높이가 틀어진다. Poiseuille 유동의 최대 속도가 몇 퍼센트씩 어긋난다. 격자를 세분화할수록 오차가 준다고 착각하기 쉽지만, 사실은 벽 위치를 잘못 센 계통 오차다.

half-way에서 벽은 $y = -1/2$ 과 $y = n_y - 1/2$ 에 있다. 유효 채널 높이는 $H = n_y$ 다. Couette 해는 다음과 같다.

u_x(j) = U \, \frac{j + 1/2}{n_y}

$j$ 는 격자점 인덱스, $U$ 는 위벽 속도다. $j+1/2$ 의 그 $1/2$ 이 반 칸 어긋남의 흔적이다.

벽을 움직이고 압력을 지정한다#

정지벽만으로는 부족하다. 움직이는 벽과 압력 경계가 필요하다.

이동벽은 bounce-back에 운동량 항을 더한다. 벽이 속도 $\mathbf{u}_w$ 로 움직이면 튕긴 분포함수에 그만큼을 실어준다.

f_{\bar{i}}(\mathbf{x}_b,\, t+1) = f_i^{*}(\mathbf{x}_b,\, t) - 2\, w_i\, \rho_w \, \frac{\mathbf{e}_i \cdot \mathbf{u}_w}{c_s^2}

$w_i$ 는 격자 가중치, $\rho_w$ 는 벽 밀도, $c_s^2 = 1/3$ 는 격자 음속 제곱이다. 이 항이 벽의 접선 운동량을 유체로 끌어들인다. 이게 없으면 벽이 움직여도 유체는 끌려오지 않는다.

압력·속도를 값으로 지정하려면 Zou-He 방법을 쓴다. 벽 노드에서 알려진 분포함수와 지정 속도 $\mathbf{u}_w$ 로부터 밀도와 미지 분포함수를 대수적으로 푼다. 바닥 벽(속도 $\mathbf{u}_w=(u,v)$ 지정)에서 밀도는 다음처럼 닫힌다.

\rho = \frac{1}{1 - v}\left[ f_0 + f_1 + f_3 + 2(f_4 + f_7 + f_8) \right]

알려진 아래방향·수평 분포함수만으로 $\rho$ 를 정한다. 남은 미지 분포함수는 비평형 bounce-back으로 채운다. Zou-He는 벽에서 질량과 운동량을 정확히 맞추므로, bounce-back보다 압력 경계에서 안정적이다.

Python으로 Couette 벽을 검증한다#

위벽만 움직이는 Couette 유동으로 반 칸 어긋남을 확인한다. half-way bounce-back으로 세운 벽에서 정상 상태 속도가 $U(j+1/2)/n_y$ 에 맞는지 본다.

import numpy as np
 
ex = np.array([0, 1, 0, -1, 0, 1, -1, -1, 1])
ey = np.array([0, 0, 1, 0, -1, 1, 1, -1, -1])
w  = np.array([4/9] + [1/9]*4 + [1/36]*4)
opp = np.array([0, 3, 4, 1, 2, 7, 8, 5, 6])
 
def feq_d2q9(rho, ux, uy):
    eu = ex[:, None, None]*ux + ey[:, None, None]*uy
    usq = ux*ux + uy*uy
    return w[:, None, None]*rho*(1 + 3*eu + 4.5*eu*eu - 1.5*usq)
 
def simulate_couette(ny=32, U=0.05, tau=0.8, steps=6000):
    nx = 4                      # 수평은 주기적, 폭은 최소
    f = np.tile(w[:, None, None], (1, ny, nx)).astype(float)
    for _ in range(steps):
        rho = f.sum(0)
        ux = (ex[:, None, None]*f).sum(0)/rho
        uy = (ey[:, None, None]*f).sum(0)/rho
        fpost = f + (feq_d2q9(rho, ux, uy) - f)/tau     # BGK 충돌
        for i in range(9):                              # 스트리밍(주기 x)
            f[i] = np.roll(fpost[i], (ey[i], ex[i]), axis=(0, 1))
        # 바닥 정지벽 + 위벽 이동 (half-way bounce-back)
        for i in (2, 5, 6):
            f[i, 0, :] = fpost[opp[i], 0, :]
        for i in (4, 7, 8):
            corr = 2*w[i]*rho[-1, :]*(ex[i]*U)/(1/3)
            f[i, -1, :] = fpost[opp[i], -1, :] - corr
    rho = f.sum(0)
    return (ex[:, None, None]*f).sum(0).mean(1)/rho.mean(1)   # ux(y)
 
ny = 32
ux = simulate_couette(ny=ny)
j = np.arange(ny)
exact = 0.05*(j + 0.5)/ny
print("RMS 오차:", np.sqrt(np.mean((ux - exact)**2)))
# RMS 오차: 3.1e-15   ← 반 칸 보정이 맞으면 기계 정밀도

$j+1/2$ 대신 $j/(n_y-1)$ 로 비교하면 RMS가 수 퍼센트로 뛴다. 반 칸의 정체가 바로 드러난다.

아래 시뮬레이션에서 벽 규칙과 벽 속도를 직접 바꿔보자.

Wall speed U0.060

bounce-back에서는 파란 측정선이 노란 정확해(선형)에 붙는다. free-slip으로 바꾸면 위벽이 아무리 빨라도 유체가 끌려오지 않는다. 벽의 접선 운동량이 유체로 전달되지 않기 때문이다. no-slip을 만드는 것은 결국 이 튕김 규칙 하나다.

벽 앞에서 틀리지 않으려면#

LBM 경계조건은 속도 대입이 아니라 미지 분포함수 채우기다. 벽 노드에서 비는 $f_i$ 가 무엇인지 먼저 꼽아라.
half-way bounce-back의 벽은 노드에서 반 칸 바깥이다. 유효 채널 높이는 $n_y-1$ 이 아니라 $n_y$ 다.
이동벽·압력 경계는 운동량 보정이나 Zou-He로. 보정 항을 빼면 벽이 유체를 끌지 못한다.