노드 하나 때문에 오차가 0.21%에서 76%로 뛰었다 — 멱법칙 유체 LBM의 국소 완화시간

멱법칙 LBM에서 τ 상한은 안정성 장치가 아니라 정상상태에 도달하기 위한 조건이다.

수렴하지 않는 로그에 찍힌 inf 한 줄#

혈액을 멱법칙(power-law) 유체로 두고 격자 볼츠만(LBM) 채널 유동을 돌렸다. 속도장은 그럴듯했다. 잔차만 내려가지 않았다. 2만 스텝을 돌려도 $10^{-3}$ 근처에서 멈췄다.

완화시간 $\tau$ 를 노드별로 찍어보니 채널 중심선 한 노드에 inf가 있었다. 오타도 아니고 0으로 나눈 것도 아니었다. 그 노드는 정말로 무한대의 점성을 요구하고 있었다.

발산은 버그가 아니다. 멱법칙 자체의 성질이다. 멱지수 슬라이더를 1.0에서 0.4로 내리면 중심 노드의 $\tau$ 가 빨간 띠로 밀린다.

power-law index n = 0.60 tau cap = 5.00 nodes  cap tau

Uncapped, the centre node reads tau = infinity — the shear rate vanishes there, so a shear-thinning fluid asks for infinite viscosity. Green band 0.55 < tau < 2 is the usable range; red on the left is the stability edge, red on the right is where bounce-back walls start leaking.

왼쪽 속도분포는 중심에서 평평해진다. cap tau를 끄면 그 노드는 축 밖으로 사라진다. 0.21%와 76%를 가른 자리가 거기 있다.

LBM이 변형률을 공짜로 주는 자리#

유한체적법에서 비뉴턴 점성을 쓰려면 먼저 속도 기울기가 필요하다. 이웃 셀을 읽고, 최소자승이든 Green–Gauss든 재구성을 돌리고, 그다음에야 $\dot\gamma$ 가 나온다. 병렬 계산에서는 여기에 통신이 한 번 더 붙는다.

LBM은 그 단계를 건너뛴다. Chapman–Enskog 전개(Knudsen 수의 멱급수로 분포함수를 쪼개는 다중스케일 기법)를 한 차수만 따라가면, 비평형 분포의 2차 모멘트가 변형률 텐서와 직접 묶인다.

S_{\alpha\beta} = -\frac{1}{2\rho c_s^2 \tau \Delta t} \sum_i c_{i\alpha} c_{i\beta}\, f_i^{neq}

$f_i^{neq} = f_i - f_i^{eq}$ 는 그 노드가 가진 평형으로부터의 이탈, $c_s^2 = 1/3$ 은 격자 음속의 제곱, $\tau$ 는 완화시간이다. 오른쪽은 전부 한 노드 안의 값이다. 이웃을 읽지 않는다.

같은 2차 모멘트는 다른 문제에서는 결손으로 얼굴을 내밀기도 한다. 여기서는 반대로 선물이다. 점성 응력이 이미 분포함수 안에 들어 있으니, 꺼내 쓰기만 하면 된다.

멱법칙은 이 $\dot\gamma$ 하나만 있으면 닫힌다.

\dot\gamma = \sqrt{2 S_{\alpha\beta} S_{\alpha\beta}}, \qquad \nu = K \dot\gamma^{\,n-1}

$K$ 는 점조도(consistency index), $n$ 은 멱지수다. $n<1$ 이면 전단박화(shear-thinning) — 빠르게 비빌수록 묽어진다. 혈액은 $n \approx 0.6$ 부근이다. $n=1$ 이면 뉴턴 유체로 돌아간다.

점성이 τ 안에 있고 τ가 점성 안에 있다#

문제는 식 두 개를 나란히 놓으면 보인다. BGK 충돌에서 점성은 완화시간이 정한다.

\nu = c_s^2 \left( \tau - \frac{1}{2} \right) \Delta t

그런데 위의 $S_{\alpha\beta}$ 식에도 $\tau$ 가 분모에 있다. $\tau$ 를 알아야 $\dot\gamma$ 를 재고, $\dot\gamma$ 를 알아야 $\nu$ 가 나오고, $\nu$ 가 나와야 $\tau$ 가 정해진다. 노드마다 작은 비선형 방정식이 하나씩 놓인 셈이다. $\nu = c_s^2(\tau - 1/2)\Delta t$ 에서 반 스텝이 어디서 붙는지는 이산 LBE의 2차 전개에 따로 있다. 아래에서는 그 식을 노드마다 다시 푸는 쪽만 다룬다.

많은 구현이 여기서 지난 스텝의 $\tau$ 를 그대로 쓴다. 이 지연(lagging)은 정상상태에서는 맞지만, 계면이나 박동처럼 전단이 빠르게 변하는 곳에서는 한 스텝만큼 틀린 점성으로 충돌을 돌린다.

지연하지 않고 풀 수도 있다. 열쇠는 $\tau$ 에 의존하지 않는 양을 먼저 분리하는 것이다. 비평형 2차 모멘트의 편차(deviatoric) 부분을 $\Pi'_{\alpha\beta}$ 라 하고

Q \equiv \sqrt{2\, \Pi'_{\alpha\beta} \Pi'_{\alpha\beta}}

로 두면, $Q$ 는 그 순간의 $f_i$ 만으로 계산된다. $\tau$ 가 들어갈 자리가 없다. 그리고 $\dot\gamma = Q / (2\rho c_s^2 \tau)$ 다.

고정점 하나로 닫히는 국소 방정식#

$m \equiv c_s^2 \tau$ 로 변수를 바꾸면 세 식이 한 줄로 합쳐진다.

m = \frac{c_s^2}{2} + B\, m^{1-n}, \qquad B \equiv K \left( \frac{Q}{2\rho} \right)^{n-1}

미지수는 $m$ 하나, 계수 $B$ 는 상수다. 이 형태가 좋은 이유는 따로 있다. 우변을 $g(m)$ 이라 하면

g'(m) = B(1-n)\, m^{-n}

이고, 해에서는 $m = c_s^2/2 + B m^{1-n} > B m^{1-n}$ 이므로 $B m^{-n} < 1$ 이다. 따라서 $|g'| < 1-n < 1$ 이 $n<1$ 에서 항상 성립한다. 전단박화 유체라면 단순 대입 반복이 무조건 수축한다. Newton도, 안전장치도, 초기추정도 필요 없다.

D2Q9에서 $\Pi_{xy}$ 에 기여하는 것은 대각 링크 네 개뿐이다. $c_{ix}c_{iy}$ 가 0인 축 방향 링크는 이 모멘트에 안 들어온다.

같은 채널에 상한 다섯 개를 걸어봤다#

물체력으로 구동하는 완전발달 채널이다. 벽은 halfway bounce-back, 폭은 31노드. 멱법칙 채널 유동에는 해석해가 있으니 오차를 그대로 잴 수 있다.

u(y) = \frac{n}{n+1} \left( \frac{g}{K} \right)^{1/n} \left( h^{\frac{n+1}{n}} - |y|^{\frac{n+1}{n}} \right)

$g$ 는 단위질량당 물체력, $h$ 는 반폭이다. $n=1$ 을 넣으면 익숙한 포물선으로 돌아간다.

import numpy as np
 
CX = np.array([0, 1, 0, -1, 0, 1, -1, -1, 1])
CY = np.array([0, 0, 1, 0, -1, 1, 1, -1, -1])
W  = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])
OPP = np.array([0, 3, 4, 1, 2, 7, 8, 5, 6])
CS2, G, NY = 1/3, 2e-6, 31
 
def eq(rho, ux):
    cu = CX[:, None] * ux
    return W[:, None] * rho * (1 + cu/CS2 + cu*cu/(2*CS2**2) - ux*ux/(2*CS2))
 
def relax_time(Q, rho, n, K, sweeps=40):
    """m = cs2*tau. 이웃 노드를 한 번도 읽지 않는 고정점 반복."""
    with np.errstate(divide='ignore', over='ignore', invalid='ignore'):
        B = K * (Q / (2*rho))**(n - 1)
        m = np.ones_like(Q)
        for _ in range(sweeps):
            m = CS2/2 + B * m**(1 - n)
    return m / CS2
 
def exact(n, K):
    h = NY/2.0
    y = np.arange(NY) - (NY - 1)/2.0
    return n/(n+1) * (G/K)**(1/n) * (h**((n+1)/n) - np.abs(y)**((n+1)/n))
 
def consistency_K(n, umax=8e-3):
    h = NY/2.0
    return G / (umax*(n+1)/(n*h**((n+1)/n)))**n
 
def channel(n, K, cap, steps=20000, tol=1e-12):
    ue = exact(n, K)
    f = eq(np.ones(NY), ue)          # 해석해에서 출발해 정상상태 보정만 이완시킨다
    prev, tau = ue.copy(), np.full(NY, 0.8)
    for it in range(steps):
        rho = f.sum(0)
        ux = (CX[:, None]*f).sum(0)/rho + 0.5*G
        fe = eq(rho, ux)
        Pxy = (CX[:, None]*CY[:, None]*(f - fe)).sum(0)
        raw = relax_time(2*np.abs(Pxy), rho, n, K)
        tau = np.clip(np.nan_to_num(raw, nan=cap, posinf=cap), 0.5005, cap)
        cu = CX[:, None]*ux
        src = W[:, None]*(1 - 1/(2*tau))*rho*G*((CX[:, None]-ux)/CS2 + cu*CX[:, None]/CS2**2)
        post = f - (f - fe)/tau + src
        for i in range(9):
            f[i] = np.roll(post[i], CY[i])
        for i in range(9):                                  # halfway bounce-back
            if CY[i] == 1:  f[i, 0]  = post[OPP[i], 0]
            if CY[i] == -1: f[i, -1] = post[OPP[i], -1]
        if it % 200 == 199:
            res = np.abs(ux - prev).max()/ux.max()
            if res < tol: break
            prev = ux.copy()
    rho = f.sum(0)
    ux = (CX[:, None]*f).sum(0)/rho + 0.5*G
    err = np.linalg.norm(ux - ue)/np.linalg.norm(ue)
    return err, res, it + 1, tau, raw
 
K1 = consistency_K(1.0)
e, r, it, tau, _ = channel(1.0, K1, 1e6)
print(f"Newtonian n=1.0  L2={e:.4%}  tau={tau.min():.4f}..{tau.max():.4f}  steps={it}")
 
n, K = 0.6, consistency_K(0.6)
_, _, _, _, raw = channel(n, K, 4.0, steps=200)
print(f"n=0.6  tau at wall {raw[0]:.4f}   tau at centre {raw[NY//2]:.4e}")
print(f"{'cap':>9} {'L2 err':>9} {'residual':>10} {'steps':>7} {'capped':>8}")
for cap in [0.8, 1.5, 4.0, 40.0, 1e6]:
    e, r, it, tau, _ = channel(n, K, cap)
    print(f"{cap:9.1f} {e:9.4%} {r:10.2e} {it:7d} {int((tau >= cap - 1e-9).sum()):5d}/{NY}")

Newtonian n=1.0  L2=0.1362%  tau=0.5901..0.5901  steps=20000
n=0.6  tau at wall 0.5690   tau at centre inf
      cap    L2 err   residual   steps   capped
      0.8   0.2112%   1.08e-06   20000     3/31
      1.5   0.2121%   1.14e-06   20000     1/31
      4.0   0.2134%   1.07e-06   20000     1/31
     40.0   0.4850%   1.02e-08   20000     1/31
1000000.0  75.9752%   7.69e-03   20000     1/31

Pi_xy를 올리면 초록 곡선이 내려가고, 고정점 계단이 대각선에 붙는 칸 수가 줄어든다.

Pi_xy = 4.00e-4 n = 0.60 K = 2.5e-3  run

The staircase is the whole solver: 0 steps of m → cs2/2 + B m^(1-n) land on tau = 0.5894 (gamma = 2.04e-3, nu = 2.98e-2). Freezing tau at the previous step instead would answer 0.6011. No neighbour node is read at any point.

$n$ 을 1 위로 올리면 기울기가 음수가 되어 계단이 나선을 그린다. $n<2$ 면 그래도 수축한다. 로제트의 축 방향 링크가 회색인 이유는 $c_{ix}c_{iy} = 0$ 이라 $\Pi_{xy}$ 에 안 실리기 때문이다.

0.21%와 76%를 가른 것은 노드 하나였다#

뉴턴 유체 기준선이 0.1362%다. 솔버 자체는 맞게 돌아간다. 같은 격자에 $n=0.6$ 을 걸고 상한만 바꿨더니 오차가 0.2112%에서 75.9752%까지 움직였다.

capped 열을 보자. 상한이 1.5 이상일 때 걸린 노드는 31개 중 정확히 1개다. 채널 중심선 노드 하나다. 그 한 노드의 $\tau$ 를 $10^6$ 으로 두는 것과 4.0으로 두는 것 사이에 오차 75.76%p가 놓여 있다.

이유는 충돌식에 있다. $f_i \leftarrow f_i - (f_i - f_i^{eq})/\tau$ 에서 $\tau \to \infty$ 면 괄호 앞 계수가 0이 된다. 그 노드는 평형으로 되돌아가지 않는다. 초기 분포에 실려 있던 오차가 그대로 남고, 주변 노드가 스트리밍으로 가져간다. 잔차가 $7.69 \times 10^{-3}$ 에서 멈춘 것이 그 결과다.

왜 중심선에서 $\dot\gamma$ 가 정확히 0인가. 대칭 때문이다. 채널 중심에서 $\partial u / \partial y$ 가 부호를 바꾸므로 그 자리에서 0을 지난다. 이산적으로도 $\Pi_{xy}$ 가 상쇄되어 정확히 0이 된다. 그러면 $\nu = K \dot\gamma^{\,n-1}$ 에서 $n-1 = -0.4$ 이므로 $\nu \to \infty$ 다. 연속체에서도 같다. 멱법칙은 전단률 0에서 정의되지 않는다.

벽의 $\tau = 0.5690$ 을 점성으로 되돌리면 $\nu = c_s^2(\tau - 1/2) = (1/3)(0.0690) = 0.0230$ 이다. 뉴턴 줄의 $\tau = 0.5901$ 은 $(1/3)(0.0901) = 0.030033$ 이다. 같은 $G = 2\times 10^{-6}$ , 폭 31노드인데 벽 점성만 이미 다르다. 코드의 하한 클립은 0.5005다. 거기서는 $\nu = (1/3)(0.0005) = 1.67\times 10^{-4}$ 까지 내려가, 벽의 0.0230보다 두 자리 이상 작다. $n = 0.6$ 이면 $1 - n = 0.4$ 라 고정점 기울기는 $|g'| < 0.4$ 안에 있다. 수축은 된다. 다만 중심에서 $Q = 0$ 이면 지수 $n - 1 = -0.4$ 가 $B$ 를 터뜨린다.

상한 0.8과 4.0의 L2는 0.2112%와 0.2134%다. 차이는 0.0022퍼센트포인트다. 1.5는 0.2121%로, 0.8보다 0.0009퍼센트포인트 크다. 걸린 노드는 3개에서 1개로 준다. 상한을 $10^6$ 으로 열면 오차 75.9752%, 잔차 $7.69\times 10^{-3}$ 이다. 상한 4.0의 잔차 $1.07\times 10^{-6}$ 보다 세 자리 이상 크다. 2만 스텝은 양쪽 다 20000으로 같다. 상한을 닫으면 $1.07\times 10^{-6}$ 까지 내려가고, 열면 $7.69\times 10^{-3}$ 에서 멈춘다. 도입부의 $10^{-3}$ 은 그 열린 쪽이다. 걸린 노드는 여전히 31개 중 1개다. 뉴턴 기준 L2는 0.1362%이고 $\tau$ 는 0.5901로 채널 전체가 같다. $n = 0.6$ 으로 바꾸자마자 중심 한 칸만 inf가 된다.

상한값을 어디에 두든 거의 같은 답이 나오는 이유#

여기서 실무적으로 중요한 사실이 나온다. 0.8, 1.5, 4.0의 오차가 0.2112%, 0.2121%, 0.2134%다. 세 자리까지 같다.

그 노드는 전단응력을 거의 나르지 않기 때문이다. 완전발달 채널에서 운동량 평형이 정하는 전단응력은 $\sigma(y) = \rho g y$ 다. 중심에서 $y=0$ 이니 $\sigma=0$ 이다. 점성을 틀리게 넣어도 곱해질 전단률이 0이라 응력 기여가 없다.

즉 발산하는 것은 $\nu$ 지, $\nu\dot\gamma$ 가 아니다. 유동이 실제로 느끼는 것은 후자다. 그래서 상한은 정확도 손실 없이 걸 수 있다. 상한이 40으로 커지면 오차가 0.4850%로 두 배가 되는데, 이때는 $\tau$ 가 커서 그 노드의 이완이 느려진 탓이다. 정확도가 아니라 수렴 속도 문제다.

반대로 상한을 너무 낮추면 다른 대가가 있다. 0.8에서는 3개 노드가 걸린다. 중심 근처 두 노드가 실제보다 묽어지고, 그 오차는 응력에 실린다. 표에서는 아직 0.2112%로 작지만, $n$ 이 더 작아지면 걸리는 노드 수가 빠르게 는다. 첫 번째 시뮬레이션의 nodes 선택을 81로 올리고 $n$ 을 0.4로 내려보면 붉은 띠에 들어가는 노드가 눈에 띄게 늘어난다.

혈액을 풀 때 이 이야기가 달라지는 지점#

멱법칙이 전단률 0에서 발산한다는 것은 유변학에서도 알려진 결함이다. 실제 혈액은 전단률이 0으로 가도 점성이 유한한 값(영전단 점성)에 머문다. Carreau–Yasuda 모델이 그 부분을 고친다.

\nu = \nu_\infty + (\nu_0 - \nu_\infty)\left[ 1 + (\lambda \dot\gamma)^a \right]^{\frac{n-1}{a}}

$\nu_0$ 와 $\nu_\infty$ 가 양 끝을 막고, $\lambda$ 가 전단박화가 시작되는 전단률의 역수다. 이 모델을 쓰면 $\dot\gamma \to 0$ 에서 $\nu \to \nu_0$ 로 유한하다. 상한을 손으로 걸 필요가 없다. $\nu_0$ 가 바로 상한이다.

그러면 왜 멱법칙을 쓰는가. 관 중앙을 제외한 대부분의 영역에서 두 모델의 차이가 작고, 멱법칙은 해석해가 있어 검증이 되기 때문이다. 게다가 위의 고정점은 멱법칙일 때만 한 줄로 닫힌다. Carreau–Yasuda에서는 같은 반복이 되지만 수축 증명이 그렇게 깔끔하지 않다.

LBM 쪽에서 $\tau$ 가 커질 때 조심할 것이 하나 더 있다. bounce-back 벽의 유효 위치가 $\tau$ 에 따라 움직인다는 점이다. 전단박화 유체는 벽 근처에서 $\tau$ 가 작으므로 이 오차는 오히려 작아진다. 대신 $\tau$ 가 $0.5$ 에 가까워지면 안정성 여유가 줄어든다. 위 계산에서 벽 노드는 0.5690이었다. 여유는 있지만 넉넉하지는 않다. $\tau - 0.5 = 0.069$ 라 하한까지의 간격은 그 정도다.

발산하는 것과 못 쓰는 것은 다르다#

inf는 0으로 나눈 버그가 아니었다. 전단률 0에서 $n-1 = -0.4$ 이면 $\nu \to \infty$ 가 맞다. 계산이 멈춘 이유는 BGK가 $\tau^{-1}$ 을 계수로 쓰기 때문이다. 점성이 무한대라는 말과, 그 노드를 영원히 이완하지 않겠다는 말이 같은 식이다.

상한은 후자를 막는 이산화다. 표에서 1.5와 4.0의 L2는 0.2121%와 0.2134%로 거의 같다. $n=1$ 이면 같은 채널 식이 포물선으로 돌아가고, 그 기준선은 포아즈이유 해석해 쪽이다. 로그에 걸린 노드가 1을 크게 넘으면 손볼 것은 상한이 아니라 $K$ 나 격자다.