동적 계수의 49.9%가 음수로 나왔다 — Germano 항등식과 평균 연산의 자리

동적 모델의 평균 연산은 후처리 옵션이 아니라 모델을 성립시키는 구성 요소다.

전이 영역에서 터졌고, 로그에는 음수가 찍혀 있었다#

동적 Smagorinsky 모델을 처음 켠 날 계산은 200스텝에서 발산했다. 로그를 열어보니 모델 계수가 음수였다. 코드를 의심했지만 코드는 맞았다. 음수는 버그가 아니라 모델이 데이터에서 읽어낸 값이었다. 이 글은 그 음수가 어디서 오는지, 그리고 왜 평균 연산 없이는 동적 모델이 성립하지 않는지를 24³ 난류장에서 직접 재서 보여준다.

결론부터 말하면 국소 계수의 49.9%가 음수였고, 그중 27.1%의 격자점에서 총 점성이 음수가 됐다. 같은 데이터에 평면 평균을 한 번 씌우자 그 비율은 0.0%로 떨어졌다.

상수 하나를 손으로 정하는 값#

고전 Smagorinsky 모델은 부분격자(SGS, subgrid-scale) 응력을 와점성으로 닫는다.

\nu_T = (C_s \Delta)^2 |\bar{S}|, \qquad |\bar{S}| = \sqrt{2 \bar{S}_{ij} \bar{S}_{ij}}

$\Delta$ 는 격자 필터 폭, $\bar{S}_{ij}$ 는 여과된 변형률 텐서, $C_s$ 는 상수다.

문제는 $C_s$ 가 상수가 아니라는 데 있다. 등방성 난류 감쇠에서는 0.17 근처, 채널 유동에서는 0.1 근처가 맞는다. 층류 구간에서도 $|\bar{S}| \neq 0$ 이면 와점성이 0이 되지 않는다. 벽에 붙으면 $\nu_T$ 가 $y^3$ 로 줄어야 하는데 모델은 그 극한을 모른다. PMBFS2 매뉴얼도 이 자리를 Van Driest 감쇠 함수로 덧댄다. 감쇠 함수는 벽까지의 거리를 알아야 쓸 수 있고, 복잡한 형상에서 그 거리는 잘 정의되지 않는다.

같은 응력을 두 층위에서 재면 차이가 남는다#

Germano가 1991년에 제안한 탈출구는 상수를 밖에서 주지 말고 해석된 스케일에서 읽어내자는 것이다. 격자 필터 위에 폭이 더 넓은 테스트 필터 $\widehat{(\cdot)}$ 를 하나 더 얹는다. 두 층위의 응력을 각각 쓰면

\tau_{ij} = \overline{u_i u_j} - \bar{u}_i \bar{u}_j, \qquad T_{ij} = \widehat{\overline{u_i u_j}} - \hat{\bar{u}}_i \hat{\bar{u}}_j

이고, 둘의 차이는 모르는 항이 모두 상쇄되어 계산 가능한 양만 남는다.

L_{ij} = T_{ij} - \widehat{\tau_{ij}} = \widehat{\bar{u}_i \bar{u}_j} - \hat{\bar{u}}_i \hat{\bar{u}}_j

우변에는 $\bar{u}$ 만 들어 있다. 즉 $L_{ij}$ 는 모델 없이 직접 잴 수 있다. 이것이 Germano 항등식이고, 동적 모델이 서 있는 유일한 발판이다.

아래 시뮬레이션에서 두 필터의 폭을 직접 벌려보자.

filter ratio α = Δ̂/Δ : 1.94highest wavenumber : 14

mean L = 0.0000  |  max |L| = 0.0000  |  energy kept by grid filter = 0.0 %

파란 곡선이 격자 필터, 주황 곡선이 테스트 필터다. 아래쪽 초록 띠가 $L(x)$ 다. 필터 비 $\alpha$ 를 1.2에서 4로 밀면 두 곡선의 간격이 벌어지고 $L$ 의 진폭이 함께 커진다. $\alpha \to 1$ 에서 $L$ 이 0으로 무너지는 것이 핵심이다. 잴 신호가 사라지면 계수도 정해지지 않는다.

5개 식을 스칼라 하나로 줄이기#

$L_{ij}$ 는 편차 성분만 세어도 식이 5개인데 미지수는 $C = C_s^2$ 하나다. 두 층위에 같은 $C$ 를 쓰고(스케일 불변 가정) 정리하면

L_{ij}^{d} = C M_{ij}, \qquad M_{ij} = 2\left[\Delta^2 \widehat{|\bar{S}| \bar{S}_{ij}} - \hat{\Delta}^2 |\hat{\bar{S}}| \hat{\bar{S}}_{ij}\right]

가 된다. 윗첨자 $d$ 는 편차부(deviatoric part)를 뜻한다. Lilly는 1992년에 이 과결정계를 최소제곱으로 닫았다. 잔차 $\|L^d_{ij} - C M_{ij}\|^2$ 를 $C$ 에 대해 미분해 0으로 두면

C = \frac{\langle L_{ij}^{d} M_{ij} \rangle}{\langle M_{kl} M_{kl} \rangle}

이 나온다. $\langle \cdot \rangle$ 이 오늘의 주인공이다. 이 괄호를 어디에 씌우느냐가 모델의 성격을 바꾼다. 괄호를 빼고 격자점마다 따로 나누면 분모가 0에 가까워지는 자리가 생긴다.

Python으로 24³ 난류를 만들고 계수를 쟀다#

난수 위상만 가진 필드에는 캐스케이드가 없다. 그래서 먼저 Navier–Stokes를 40스텝 적분해 상관관계를 만든 뒤, 그 결과에 a priori 테스트를 돌렸다. 외부 라이브러리 없이 리스트만 쓴다.

import math, random
 
N, NU = 24, 0.02
NP, H = N * N * N, 2.0 * math.pi / N
GRID_W, TEST_W = 3, 5                       # 박스 필터 폭, 셀 단위
DELTA = GRID_W * H                          # 격자 필터 폭
DELTA_T = math.sqrt((GRID_W * H) ** 2 + (TEST_W * H) ** 2)   # 합성된 테스트 층위
 
_perm = {}
def shift_perm(axis, off):
    """axis 방향으로 off 칸 이동하는 주기 인덱스 맵"""
    if (axis, off) not in _perm:
        p = [0] * NP
        for i in range(N):
            for j in range(N):
                for k in range(N):
                    a, b, c = i, j, k
                    if axis == 0: a = (i + off) % N
                    elif axis == 1: b = (j + off) % N
                    else: c = (k + off) % N
                    p[(i * N + j) * N + k] = (a * N + b) * N + c
        _perm[(axis, off)] = p
    return _perm[(axis, off)]
 
def box_filter(f, w):
    r, out = w // 2, f
    for axis in (0, 1, 2):
        acc = [0.0] * NP
        for off in range(-r, r + 1):
            acc = [a + out[q] for a, q in zip(acc, shift_perm(axis, off))]
        out = [v / w for v in acc]
    return out
 
def ddx(f, axis):
    inv = 1.0 / (2.0 * H)
    return [(f[a] - f[b]) * inv for a, b in zip(shift_perm(axis, 1), shift_perm(axis, -1))]
 
def lap(f):
    out = [-6.0 * v for v in f]
    for axis in (0, 1, 2):
        for off in (1, -1):
            out = [o + f[q] for o, q in zip(out, shift_perm(axis, off))]
    return [v / (H * H) for v in out]
 
def divergence(u):
    d = ddx(u[0], 0)
    d = [a + b for a, b in zip(d, ddx(u[1], 1))]
    return [a + b for a, b in zip(d, ddx(u[2], 2))]
 
def synth_field(nmodes, kmax, seed):
    """발산이 0인 난수 Fourier 장, E(k) ~ k^(-5/3)"""
    random.seed(seed)
    u, xs = [[0.0] * NP for _ in range(3)], [i * H for i in range(N)]
    for _ in range(nmodes):
        while True:
            kv = [random.randint(-kmax, kmax) for _ in range(3)]
            km = math.sqrt(kv[0]**2 + kv[1]**2 + kv[2]**2)
            if 1.0 <= km <= kmax: break
        amp = km ** (-5.0 / 6.0)
        while True:
            r = [random.gauss(0, 1) for _ in range(3)]
            e = [r[1]*kv[2]-r[2]*kv[1], r[2]*kv[0]-r[0]*kv[2], r[0]*kv[1]-r[1]*kv[0]]
            en = math.sqrt(e[0]**2 + e[1]**2 + e[2]**2)
            if en > 1e-9: break
        e, ph = [c / en for c in e], random.uniform(0, 2 * math.pi)
        ax = [kv[0]*x for x in xs]; by = [kv[1]*x for x in xs]; cz = [kv[2]*x for x in xs]
        for i in range(N):
            for j in range(N):
                base, o = ax[i] + by[j] + ph, (i*N+j)*N
                for k in range(N):
                    c = math.cos(base + cz[k])
                    u[0][o+k] += amp*e[0]*c; u[1][o+k] += amp*e[1]*c; u[2][o+k] += amp*e[2]*c
    rms = math.sqrt(sum(v*v for comp in u for v in comp) / NP)
    return [[v / rms for v in comp] for comp in u]
 
def project(u, phi, sweeps):
    """발산 제거. 라플라시안은 2h 스텐실의 div(grad)와 일치시킨다"""
    rhs, hh = divergence(u), (2.0 * H) ** 2
    for _ in range(sweeps):
        acc = [0.0] * NP
        for axis in (0, 1, 2):
            for off in (2, -2):
                acc = [a + phi[q] for a, q in zip(acc, shift_perm(axis, off))]
        phi = [(a - hh * r) / 6.0 for a, r in zip(acc, rhs)]
    for d in range(3):
        u[d] = [v - s for v, s in zip(u[d], ddx(phi, d))]
    return u, phi
 
def rhs_ns(u):
    out = []
    for d in range(3):
        adv = [0.0] * NP
        for ax in range(3):
            g = ddx(u[d], ax)
            adv = [a + v * gg for a, v, gg in zip(adv, u[ax], g)]
        out.append([-a + NU * l for a, l in zip(adv, lap(u[d]))])
    return out
 
def advance(u, dt, nsteps, sweeps):
    """RK2 + 압력 투영: 난수 위상이 실제 캐스케이드로 바뀐다"""
    phi = [0.0] * NP
    for _ in range(nsteps):
        k1 = rhs_ns(u)
        mid = [[v + 0.5*dt*r for v, r in zip(u[d], k1[d])] for d in range(3)]
        k2 = rhs_ns(mid)
        u, phi = project([[v + dt*r for v, r in zip(u[d], k2[d])] for d in range(3)], phi, sweeps)
    return u
 
def strain_tensor(u):
    g = [[ddx(u[d], ax) for ax in range(3)] for d in range(3)]
    S = [[None]*3 for _ in range(3)]
    for a in range(3):
        for b in range(a, 3):
            S[a][b] = [0.5*(x+y) for x, y in zip(g[a][b], g[b][a])]
            S[b][a] = S[a][b]
    mag = [0.0]*NP
    for a in range(3):
        for b in range(3):
            mag = [m + 2.0*s*s for m, s in zip(mag, S[a][b])]
    return S, [math.sqrt(m) for m in mag]
 
def leonard_stress(ub, w):
    """L_ij = test(u_i u_j) - test(u_i) test(u_j), 편차부만 남긴다"""
    ut = [box_filter(c, w) for c in ub]
    L = [[None]*3 for _ in range(3)]
    for a in range(3):
        for b in range(a, 3):
            prod = box_filter([x*y for x, y in zip(ub[a], ub[b])], w)
            L[a][b] = [p - x*y for p, x, y in zip(prod, ut[a], ut[b])]
            L[b][a] = L[a][b]
    tr = [0.0]*NP
    for a in range(3):
        tr = [t + v for t, v in zip(tr, L[a][a])]
    for a in range(3):
        L[a][a] = [v - t/3.0 for v, t in zip(L[a][a], tr)]
    return L, ut
 
def m_tensor(S, mag, ut, w, dg, dt_):
    """M_ij = 2[ D^2 test(|S|S_ij) - Dhat^2 |S_test| S_test_ij ]"""
    St, magt = strain_tensor(ut)
    M = [[None]*3 for _ in range(3)]
    for a in range(3):
        for b in range(a, 3):
            t1 = box_filter([m*s for m, s in zip(mag, S[a][b])], w)
            M[a][b] = [2.0*(dg*dg*x - dt_*dt_*mt*st) for x, mt, st in zip(t1, magt, St[a][b])]
            M[b][a] = M[a][b]
    return M
 
def contract(A, B):
    out = [0.0]*NP
    for a in range(3):
        for b in range(3):
            out = [o + x*y for o, x, y in zip(out, A[a][b], B[a][b])]
    return out
 
def plane_average(v):
    acc = [0.0]*N
    for i in range(N):
        for j in range(N):
            o = (i*N+j)*N
            for k in range(N):
                acc[k] += v[o+k]
    return [a/(N*N) for a in acc]
 
def pct(v, q):
    s = sorted(v)
    return s[min(len(s)-1, int(q*len(s)))]
 
u = advance(synth_field(40, 8, 20260915), 0.05, 40, 40)
urms = math.sqrt(sum(v*v for c in u for v in c) / NP)
dv = divergence(u)
sg = [0.0]*NP
for d in range(3):
    for a in range(3):
        sg = [x + y*y for x, y in zip(sg, ddx(u[d], a))]
print("grid %d^3  nu %.3f  t_end %.2f  u_rms %.4f" % (N, NU, 0.05*40, urms))
print("rms|div u| / rms|grad u|      : %.3f"
      % (math.sqrt(sum(v*v for v in dv)/NP) / math.sqrt(sum(sg)/NP)))
 
ub = [box_filter(c, GRID_W) for c in u]
S, mag = strain_tensor(ub)
L, ut = leonard_stress(ub, TEST_W)
 
print("--- Germano-Lilly coefficient  C = Cs^2 ---")
ref = None
for tag, dt_ in (("composed  a=%.3f" % (DELTA_T/DELTA), DELTA_T),
                 ("textbook  a=2.000", 2.0*DELTA),
                 ("test only a=%.3f" % (TEST_W/GRID_W), TEST_W*H)):
    M = m_tensor(S, mag, ut, TEST_W, DELTA, dt_)
    c = sum(contract(L, M)) / sum(contract(M, M))
    if ref is None:
        ref, Mref = c, M
        print("  %s : C = %.6f   Cs = %.4f" % (tag, c, math.sqrt(c)))
    else:
        print("  %s : C = %.6f   Cs = %.4f   (%+.1f%%)" % (tag, c, math.sqrt(c), 100*(c/ref-1)))
 
LM, MM = contract(L, Mref), contract(Mref, Mref)
Cloc = [a/b for a, b in zip(LM, MM)]
nuT = [c*DELTA*DELTA*m for c, m in zip(Cloc, mag)]
mm_mean = sum(MM) / NP
print("--- pointwise C (no averaging) ---")
print("  C < 0 fraction              : %.1f %%" % (100.0*sum(1 for c in Cloc if c < 0)/NP))
print("  C  p01 / p50 / p99          : %+.4f / %+.4f / %+.4f" % (pct(Cloc,0.01), pct(Cloc,0.5), pct(Cloc,0.99)))
print("  M:M < 1e-3 * <M:M>          : %.2f %%" % (100.0*sum(1 for m in MM if m < 1e-3*mm_mean)/NP))
print("  nu_T(global C) / nu         : %.2f" % (ref*DELTA*DELTA*(sum(mag)/NP)/NU))
print("  nu + nu_T < 0               : %.1f %%" % (100.0*sum(1 for v in nuT if NU+v < 0)/NP))
print("  worst nu_T / nu             : %.1f" % (min(nuT)/NU))
 
Cpl = [a/b for a, b in zip(plane_average(LM), plane_average(MM))]
nuTp = [Cpl[n % N]*DELTA*DELTA*mag[n] for n in range(NP)]
print("--- C averaged over i-j planes ---")
print("  C range over %d planes       : %+.5f .. %+.5f" % (N, min(Cpl), max(Cpl)))
print("  negative planes             : %d / %d" % (sum(1 for c in Cpl if c < 0), N))
print("  nu + nu_T < 0               : %.1f %%" % (100.0*sum(1 for v in nuTp if NU+v < 0)/NP))

출력은 이렇게 나온다.

grid 24^3  nu 0.020  t_end 2.00  u_rms 0.5140
rms|div u| / rms|grad u|      : 0.008
--- Germano-Lilly coefficient  C = Cs^2 ---
  composed  a=1.944 : C = 0.008296   Cs = 0.0911
  textbook  a=2.000 : C = 0.008035   Cs = 0.0896   (-3.1%)
  test only a=1.667 : C = 0.001029   Cs = 0.0321   (-87.6%)
--- pointwise C (no averaging) ---
  C < 0 fraction              : 49.9 %
  C  p01 / p50 / p99          : -0.5499 / +0.0001 / +0.3246
  M:M < 1e-3 * <M:M>          : 0.07 %
  nu_T(global C) / nu         : 0.25
  nu + nu_T < 0               : 27.1 %
  worst nu_T / nu             : -78.6
--- C averaged over i-j planes ---
  C range over 24 planes       : -0.00021 .. +0.01692
  negative planes             : 1 / 24
  nu + nu_T < 0               : 0.0 %

전역 최소제곱으로 얻은 $C_s$ 는 0.0911이다. 박스 필터를 쓴 a priori 테스트에서 흔히 보고되는 0.09~0.12 구간에 들어온다.

절반이 음수인 것은 오류가 아니다#

국소 계수의 49.9%가 음수다. 중앙값은 $+0.0001$ 로 사실상 0인데, 1퍼센타일은 $-0.55$ , 99퍼센타일은 $+0.32$ 다. 평균값 0.0083의 40~66배가 양쪽으로 벌어져 있다.

이 음수에는 물리적 의미가 있다. 에너지는 큰 스케일에서 작은 스케일로만 흐르지 않는다. 국소적으로는 반대로 흐르며, 이를 역산란(backscatter)이라 부른다. 실제 난류에서 역산란은 전체 격자점의 30~50%에서 관측된다. 동적 모델은 이 방향까지 읽어낸다는 점에서 정직하다.

문제는 정직함의 대가다. $C < 0$ 이면 $\nu_T < 0$ 이고, 확산 항의 부호가 뒤집힌다. 위 계산에서 격자점의 27.1%가 $\nu + \nu_T < 0$ 상태였다. 최악의 지점은 $\nu_T / \nu = -78.6$ 이다. 분자 점성의 78배를 거꾸로 밀어넣는 셈이다. 발산은 정해진 결과다.

아래 지도에서 평균 창 크기를 직접 밀어보자.

averaging window : 1×1number of modes : 18

C < 0 : 0.0 %  |  ν + νT < 0 : 0.0 %  |  C range 0.0000 … 0.0000  |  global C = 0.00000

파란색이 음수 계수, 붉은색이 양수다. 창을 1×1에서 키우면 파란 점들이 먼저 사라지고, ν + νT < 0 비율이 0으로 내려간다. clipping 버튼은 평균 대신 아래에서 잘라내는 대안이다. 두 방식이 지도에 남기는 흔적이 어떻게 다른지 보는 것이 관찰 포인트다.

괄호를 어디에 씌우는가#

Germano와 Lilly가 남긴 $\langle \cdot \rangle$ 는 수식을 예쁘게 만드는 장치가 아니다. 위 출력의 마지막 블록이 그 증거다. 분자와 분모를 각각 $i$ – $j$ 평면에서 평균한 뒤 나누면, 24개 평면 중 음수는 1개뿐이고 그 값도 $-0.0002$ 다. $\nu + \nu_T < 0$ 인 격자점은 0.0%다. 같은 데이터, 같은 식인데 불안정 격자점이 27.1%에서 0%로 사라진다.

순서가 중요하다. $C$ 를 먼저 구해 평균하면 안 된다. 분모가 0에 가까운 점의 $C$ 가 발산하기 때문이다. 위 계산에서 $M_{ij}M_{ij}$ 가 평균의 1/1000 미만인 점이 0.07% 있었다. 반드시 $\langle L{:}M \rangle$ 과 $\langle M{:}M \rangle$ 을 따로 평균한 뒤 나눈다.

평균 방향은 유동이 정한다. 채널이면 벽에 평행한 평면, 원관이면 원주 방향과 축 방향이다. 동차 방향이 없으면 유적선을 따라 평균하는 Lagrangian 동적 모델을 쓴다. 상수가 격자 폭에 어떻게 얽히는지는 LBM 격자 세밀화의 비평형 재스케일에서 다룬 적이 있다.

$\alpha$ 를 잘못 잡으면 계수의 88%가 사라진다#

$M_{ij}$ 에는 $\hat{\Delta}^2$ 이 들어간다. 이 값을 얼마로 놓느냐가 $C$ 를 통째로 움직인다. 위 코드는 폭 3셀 박스 필터 위에 폭 5셀 박스 필터를 얹었다. 테스트 층위의 실효 폭은 5셀이 아니다. 박스 필터를 두 번 걸면 2차 모멘트가 더해지므로

\hat{\Delta} = \sqrt{\Delta^2 + \Delta_{\text{test}}^2} = \sqrt{(3h)^2 + (5h)^2} = 5.83h

이고 $\alpha = \hat{\Delta}/\Delta = 1.944$ 다. 교과서 관행대로 $\alpha = 2$ 를 쓰면 $C$ 는 3.1%만 움직인다. 그러나 테스트 필터의 폭 $5h$ 를 그대로 $\hat{\Delta}$ 에 넣으면 $\alpha = 1.667$ 이 되고 $C$ 는 87.6% 줄어든다. $C_s$ 로 보면 0.0911이 0.0321로 떨어진다.

이유는 $M_{ij}$ 의 구조에 있다. 두 항의 차이로 만들어진 양이라 $\alpha^2 - 1$ 에 비례해 자란다. $\alpha$ 가 1에 가까워지면 분자와 분모가 함께 0으로 가고, 그 비율이 급격히 일그러진다. 첫 번째 시뮬레이션에서 $\alpha$ 를 1.2 쪽으로 밀 때 $L$ 이 무너지던 장면이 같은 현상이다.

PMBFS2가 동적 모델을 켜지 않은 이유#

이 글의 소재가 된 매뉴얼은 동적 모델을 구현해 두고도 실제 계산에는 대수적 Smagorinsky를 썼다고 적는다. 이유로 두 가지를 든다. 실제 유체에 최적인 SGS 모델이 아직 정해지지 않았다는 것, 그리고 이중 필터를 쓰는 동적 모델은 격자 요구가 더 엄격하다는 것이다.

두 번째가 실질적이다. 동적 절차가 성립하려면 $\Delta$ 와 $\hat{\Delta}$ 가 둘 다 관성 부영역 안에 있어야 한다. $\hat{\Delta}$ 는 $\Delta$ 의 2배다. 즉 격자를 2배 더 촘촘히 깔아야 같은 가정이 유지된다. 3차원에서 이 요구는 셀 수로 8배다. 모델 상수를 자동으로 얻는 대가를 격자로 치르는 구조다.

그래서 실전 판단은 이렇게 갈린다. 동차 방향이 하나라도 있고 전이나 이완이 중요한 문제라면 동적 모델이 값을 한다. 형상이 복잡하고 격자 예산이 빠듯하면 대수 모델에 감쇠 함수를 붙이는 쪽이 현실적이다. 물성으로 넣지 않은 상수가 계산 안에서 조용히 정해지는 사례는 열 LBM에 잠겨 있는 두 상수에서도 본 적이 있다.

로그에 음수가 찍히면 이제 코드를 먼저 의심하지 않는다. 괄호를 어디에 씌웠는지부터 본다.