Skip to content
cfd-lab:~/ja/posts/2026-09-15-dynamic-smago…● online
$_
Younglin's CFD Blog
毎日の記録 · 計算流体力学
NOTE #156DAY TUE 유체역학DATE 2026.09.15READ 6 min read#Dynamic-Smagorinsky#Germano-Identity#Turbulence#Fluid-Dynamics#Numerical-Analysis

動的係数の49.9%が負になった — Germano恒等式と平均操作の位置

動的モデルの平均操作は後処理のオプションではなく、モデルを成立させる構成要素です。

遷移領域で発散し、ログには負の数が残っていた

動的 Smagorinsky モデルを初めて有効にした日、計算は 200 ステップで発散しました。ログを開くとモデル係数が負になっていました。コードを疑いましたが、コードは正しかったのです。負の値はバグではなく、モデルがデータから読み取った値でした。この記事では、その負の値がどこから来るのか、そしてなぜ平均操作なしに動的モデルが成立しないのかを、24³ の乱流場で直接測って示します。

先に結論を述べます。局所係数の 49.9% が負であり、そのうち 27.1% の格子点で総粘性が負になりました。同じデータに平面平均を一度かけると、その割合は 0.0% まで下がります。

定数を手で決めることの代償

古典的な Smagorinsky モデルは、サブグリッドスケール（SGS）応力を渦粘性で閉じます。

νT=(CsΔ)2∣Sˉ∣,∣Sˉ∣=2SˉijSˉij\nu_T = (C_s \Delta)^2 |\bar{S}|, \qquad |\bar{S}| = \sqrt{2 \bar{S}_{ij} \bar{S}_{ij}}

Δ\Delta は格子フィルタ幅、Sˉij\bar{S}_{ij} はフィルタされたひずみ速度テンソル、CsC_s が定数です。

問題は CsC_s が定数ではないことです。等方性乱流の減衰では 0.17 付近、チャネル流では 0.1 付近が適切とされます。層流域でも ∣Sˉ∣≠0|\bar{S}| \neq 0 なら渦粘性はゼロになりません。壁に近づくと νT\nu_T は y3y^3 で減衰するべきですが、モデルはその極限を知りません。PMBFS2 のマニュアルも同じ箇所を Van Driest 減衰関数で補っています。減衰関数は壁までの距離を必要とし、複雑形状ではその距離がうまく定義できません。

同じ応力を二つの階層で測ると差が残る

Germano が 1991 年に示した逃げ道は、定数を外から与えるのをやめ、解像されたスケールから読み取ることでした。格子フィルタの上に、より幅の広いテストフィルタ (⋅)^\widehat{(\cdot)} を重ねます。二つの階層の応力を書くと

τij=uiuj‾−uˉiuˉj,Tij=uiuj‾^−uˉ^iuˉ^j\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

となり、両者の差では未知項がすべて相殺して、計算可能な量だけが残ります。

Lij=Tij−τij^=uˉiuˉj^−uˉ^iuˉ^jL_{ij} = T_{ij} - \widehat{\tau_{ij}} = \widehat{\bar{u}_i \bar{u}_j} - \hat{\bar{u}}_i \hat{\bar{u}}_j

右辺には uˉ\bar{u} しか含まれません。つまり LijL_{ij} はモデルなしで直接測れます。これが Germano 恒等式であり、動的モデルが立つ唯一の足場です。

下のシミュレーションで、二つのフィルタ幅を実際に広げてみてください。

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

青い曲線が格子フィルタ、オレンジがテストフィルタ、下の緑の帯が L(x)L(x) です。フィルタ比 α\alpha を 1.2 から 4 へ動かすと、二つの曲線の間隔が開き、LL の振幅も一緒に大きくなります。注目すべきは α→1\alpha \to 1 で LL がゼロへ崩れる点です。測るべき信号が消えれば、係数も定まりません。

5 本の式をスカラー 1 個に縮める#

LijL_{ij} は偏差成分だけ数えても式が 5 本あるのに、未知数は C=Cs2C = C_s^2 の 1 個です。両階層に同じ CC を使う（スケール不変性の仮定）と整理できて

Lijd=CMij,Mij=2[Δ2∣Sˉ∣Sˉij^−Δ^2∣Sˉ^∣Sˉ^ij]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]

となります。上付きの dd は偏差部を表します。Lilly は 1992 年にこの優決定系を最小二乗で閉じました。残差 ∥Lijd−CMij∥2\|L^d_{ij} - C M_{ij}\|^2 を CC で微分してゼロと置くと

C=⟨LijdMij⟩⟨MklMkl⟩C = \frac{\langle L_{ij}^{d} M_{ij} \rangle}{\langle M_{kl} M_{kl} \rangle}

が得られます。この ⟨⋅⟩\langle \cdot \rangle が本題です。この括弧をどこにかけるかでモデルの性格が変わります。括弧を外して格子点ごとに割ると、分母がゼロに近づく場所が生まれます。

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):
    """発散ゼロのランダム 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 %

全体最小二乗で得た CsC_s は 0.0911 です。ボックスフィルタを使った a priori テストでよく報告される 0.09〜0.12 の範囲に入ります。

半分が負であることは誤りではない

局所係数の 49.9% が負です。中央値は +0.0001+0.0001 で実質ゼロなのに、1 パーセンタイルは −0.55-0.55、99 パーセンタイルは +0.32+0.32 です。平均値 0.0083 の 40〜66 倍が両側に広がっています。

この負の値には物理的な意味があります。エネルギーは大スケールから小スケールへ一方向に流れるだけではありません。局所的には逆向きにも流れ、これを逆散乱（バックスキャッター）と呼びます。実際の乱流では格子点の 30〜50% で観測されます。動的モデルはその向きまで読み取る点で正直です。

問題はその正直さの代償です。C<0C < 0 なら νT<0\nu_T < 0 となり、拡散項の符号が反転します。上の計算では格子点の 27.1% が ν+νT<0\nu + \nu_T < 0 でした。最悪の点は νT/ν=−78.6\nu_T / \nu = -78.6 です。分子粘性の 78 倍を逆向きに押し込むことになります。発散は必然の結果です。

下の地図で平均窓の大きさを動かしてみてください。

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

青が負の係数、赤が正です。窓を 1×1 から広げると青い点が先に消え、ν + νT < 0 の割合がゼロへ下がります。clipping ボタンは平均の代わりに下側で切り落とす代替案です。二つの方法が地図に残す痕跡の違いが観察ポイントです。

括弧をどこにかけるか

Germano と Lilly が残した ⟨⋅⟩\langle \cdot \rangle は、式を整えるための装飾ではありません。出力の最後のブロックがその証拠です。分子と分母をそれぞれ ii–jj 平面で平均してから割ると、24 枚の平面のうち負は 1 枚だけで、その値も −0.0002-0.0002 です。ν+νT<0\nu + \nu_T < 0 の格子点は 0.0% です。同じデータ、同じ式なのに、不安定な格子点が 27.1% からゼロへ消えます。

順序が重要です。CC を先に求めてから平均してはいけません。分母がゼロに近い点で CC が発散するからです。上の計算では MijMijM_{ij}M_{ij} が平均の 1/1000 未満の点が 0.07% ありました。必ず ⟨L:M⟩\langle L{:}M \rangle と ⟨M:M⟩\langle M{:}M \rangle を別々に平均してから割ります。

平均方向は流れが決めます。チャネルなら壁に平行な平面、円管なら周方向と軸方向です。一様方向がまったくない場合は、流跡線に沿って平均する Lagrangian 動的モデルを使います。定数が格子幅とどう絡むかはLBM 格子細分化の非平衡リスケールでも扱いました。

α\alpha を取り違えると係数の 88% が消える#

MijM_{ij} には Δ^2\hat{\Delta}^2 が入ります。この値をいくつに置くかが CC を丸ごと動かします。上のコードは幅 3 セルのボックスフィルタの上に幅 5 セルを重ねています。テスト階層の実効幅は 5 セルではありません。ボックスフィルタを二度かけると二次モーメントが加算されるので

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

となり、α=Δ^/Δ=1.944\alpha = \hat{\Delta}/\Delta = 1.944 です。教科書どおり α=2\alpha = 2 を使っても CC は 3.1% しか動きません。しかしテストフィルタの幅 5h5h をそのまま Δ^\hat{\Delta} に入れると α=1.667\alpha = 1.667 になり、CC は 87.6% 減ります。CsC_s で見れば 0.0911 が 0.0321 に落ちます。

理由は MijM_{ij} の構造にあります。二つの項の差として作られる量なので、α2−1\alpha^2 - 1 に比例して育ちます。α\alpha が 1 に近づくと分子と分母が一緒にゼロへ向かい、その比が急激に歪みます。最初のシミュレーションで α\alpha を 1.2 側へ押したときに LL が崩れた場面と同じ現象です。

PMBFS2 が動的モデルを有効にしなかった理由#

この記事の題材になったマニュアルは、動的モデルを実装しておきながら、実際の計算には代数的 Smagorinsky モデルを使ったと書いています。理由は二つ挙げられています。実在流体に最適な SGS モデルがまだ定まっていないこと、そして二重フィルタを使う動的モデルは格子要求がより厳しいことです。

実務的なのは二つ目です。動的手続きが成立するには Δ\Delta と Δ^\hat{\Delta} の両方が慣性小領域の中になければならず、Δ^\hat{\Delta} は Δ\Delta の 2 倍です。つまり同じ仮定を保つには格子を 2 倍細かく刻む必要があり、三次元ではセル数で 8 倍になります。モデル定数を自動で得る代償を格子で払う構造です。

そこで実務判断はこう分かれます。一様方向が一つでもあり、遷移や緩和が重要な問題なら動的モデルは値打ちがあります。形状が複雑で格子予算が厳しければ、代数モデルに減衰関数を付ける方が現実的です。物性として与えていない定数が計算の中で静かに決まる例は、熱 LBM に潜む二つの定数でも見ました。

ログに負の数が出たら、もうコードを最初に疑いません。括弧をどこにかけたかから見ます。

関連記事

水流が乱れる瞬間 — レイノルズ臨界数と乱流境界層の多層構造層流はどこで崩れ、乱流境界層はどう積み重なるのか
May 5, 2026·6 min read
【OpenFOAMチュートリアル】後方ステップ流の解析 — simpleFoamで再循環領域を捉えるOpenFOAMのsimpleFoamを使用して、後方ステップ（Backward Facing Step）における非圧縮性流れを段階的に解析するガイド
Mar 26, 2026·2 min read
乱流が火炎を速くする効果が折れ曲がる点 — Damköhler のしわ仮説と bending乱流火炎速度が直線的に上昇し、折れ曲がり、消炎に至る三つの領域の物理
Jul 28, 2026·8 min read

役に立ったらシェアしてください。

GitHub でソースを見る ↗誤字・改善提案 ↗