CFDManager/docs/theory/_legacy/kbc_lbm.ipynb

Метод KBC (энтропийный Lattice Boltzmann): математика, физика и рабочий разбор

Karlin–Bösch–Chikatamarla (KBC) — энтропийная multi-relaxation-time (MRT) формулировка решёточного метода Больцмана (Lattice Boltzmann Method, LBM), которая обеспечивает безусловную нелинейную устойчивость при больших числах Рейнольдса без явной модели подсеточной турбулентности (она работает как implicit LES).

Этот ноутбук — подготовка к реализации: мы разбираем метод «на бумаге», выводим ключевые формулы, численно их самопроверяем и строим рабочий 2D-солвер (решётка D2Q9) с валидацией. В конце — перенос на целевую решётку D3Q27 и сравнение с подходом «полиномиальное равновесие + TRT + α-лимитер + Smagorinsky», который применялся ранее (это и есть «что менять для перехода к настоящему KBC»).

О источниках. Формулы по возможности взяты из оригинальных статей (Karlin, Bösch, Chikatamarla, Phys. Rev. E 90, 031302(R), 2014; Bösch, Chikatamarla, Karlin, Phys. Rev. E 92, 043309, 2015; Guo et al., Phys. Rev. E 65, 046308, 2002). Где доступ был только к вторичным источникам, формула помечена ⚠ сверить с оригиналом, а в конце собран список таких мест. Главный страж корректности здесь — численная самопроверка каждой формулы.

Сводная таблица обозначений

Везде используются решёточные единицы (lattice units): шаг сетки \Delta x = 1 (lu — lattice unit) и шаг по времени \Delta t = 1 (ts — time step). Переход к СИ описан в разделе «Единицы и обезразмеривание».

Символ	Название	Физический смысл	Единицы (lu/ts)	Единицы (СИ)	Диапазон
f_i	популяция	плотность частиц со скоростью \mathbf c_i	f_i безразм. (доля \rho)	кг·с³/м⁶	f_i>0
\mathbf c_i	дискретная скорость	направление переноса i	lu/ts	м/с	компоненты \in\{-1,0,1\}
w_i	вес квадратуры	вес направления, \sum_i w_i=1	безразм.	безразм.	\{4/9,1/9,1/36\} (D2Q9)
c_s	скорость звука решётки	c_s^2=1/3 для D2Q9/D3Q27	lu/ts	м/с	c_s=1/\sqrt3
\rho	плотность	\rho=\sum_i f_i	безразм. (опорная 1)	кг/м³	\approx 1
\mathbf u	скорость	\rho\mathbf u=\sum_i \mathbf c_i f_i	lu/ts	м/с	\lvert\mathbf u\rvert\ll c_s
\mathbf j	импульс	\mathbf j=\rho\mathbf u	—	—	—
\tau	время релаксации	задаёт вязкость	ts	с	\tau>1/2
\nu	кинематич. вязкость	\nu=c_s^2(\tau-\tfrac12)	lu²/ts	м²/с	>0
\beta	параметр сдвиговой релаксации	\beta=1/(2\tau)	безразм.	—	0<\beta<1
\gamma	энтропийный стабилизатор	релаксация высших мод	безразм.	—	\approx 2
f_i^{\mathrm{eq}}	равновесие	максимум энтропии при связях \rho,\mathbf u	как f_i	—	>0
k_i,s_i,h_i	части разложения	кинематич./сдвиг/высшие моды	как f_i	—	—
\Delta s_i,\Delta h_i	неравновесные части	s_i-s_i^{eq}, h_i-h_i^{eq}	как f_i	—	малы
H	H-функция (энтропия)	H=\sum_i f_i\ln(f_i/w_i)	безразм.	—	убывает
$\langle X	Y\rangle$	энтропийное скал. произв.	\sum_i X_iY_i/f_i^{eq}	—	—
\mathrm{Ma}	число Маха	\mathrm{Ma}=\lvert\mathbf u\rvert/c_s	безразм.	—	\lesssim 0.3
\mathrm{Re}	число Рейнольдса	\mathrm{Re}=UL/\nu	безразм.	—	задача
S_{\alpha\beta}	тензор скоростей деформации	\tfrac12(\partial_\alpha u_\beta+\partial_\beta u_\alpha)	1/ts	1/с	—
\nu_t	вихревая вязкость (LES)	\nu_t=(C_s\Delta)^2\lvert S\rvert	lu²/ts	м²/с	\ge 0
C_s	константа Смагоринского	—	безразм.	—	$0.1$–0.17
q	доля стенки	положение стенки между узлами	безразм.	—	[0,1]

0. Стиль графиков и общие настройки

Подключаем библиотеки и задаём единый «симпатичный» стиль matplotlib. Код написан так, чтобы файл работал и как ноутбук (графики встраиваются), и как обычный скрипт python kbc_lbm.py (графики сохраняются в figures/ для валидации физики).

In [1]:

import numpy as np
import matplotlib

# Определяем, выполняемся ли мы внутри Jupyter. В скрипте используем неинтерактивный
# бэкенд Agg и сохраняем PNG; в ноутбуке полагаемся на inline-вывод nbconvert.
try:
    get_ipython()  # есть только в IPython/Jupyter
    IN_NOTEBOOK = True
except NameError:
    IN_NOTEBOOK = False
if not IN_NOTEBOOK:
    matplotlib.use("Agg")
    # Консоль Windows может быть в cp1251 — выводим в UTF-8, чтобы кириллица и
    # греческие символы в print() не роняли скрипт (в ноутбуке это не нужно).
    import sys
    try:
        sys.stdout.reconfigure(encoding="utf-8")
    except Exception:
        pass

import os
import matplotlib.pyplot as plt

FIGDIR = os.path.join(os.path.dirname(os.path.abspath(__file__)) if "__file__" in globals() else ".", "figures")
os.makedirs(FIGDIR, exist_ok=True)

plt.rcParams.update({
    "figure.dpi": 110,
    "savefig.dpi": 120,
    "font.size": 11,
    "axes.titlesize": 13,
    "axes.titleweight": "bold",
    "axes.grid": True,
    "grid.alpha": 0.25,
    "axes.axisbelow": True,
    "figure.facecolor": "white",
    "axes.facecolor": "#fbfbfd",
})

def _finish(fig, name):
    """Сохранить фигуру в figures/ (для скрипта) и показать (в ноутбуке)."""
    fig.savefig(os.path.join(FIGDIR, name), bbox_inches="tight")
    if IN_NOTEBOOK:
        plt.show()           # гарантированный inline-вывод при nbconvert
    else:
        plt.close(fig)

print("matplotlib", matplotlib.__version__, "| numpy", np.__version__, "| notebook:", IN_NOTEBOOK)

matplotlib 3.10.8 | numpy 2.4.3 | notebook: True

1. Основы LBM и решётка D2Q9

Решёточный метод Больцмана отслеживает популяции f_i(\mathbf x, t) — плотность фиктивных частиц, движущихся в узле \mathbf x со скоростью \mathbf c_i. Динамика — чередование переноса (streaming) и столкновения (collision):

f_i(\mathbf x + \mathbf c_i\,\Delta t,\; t+\Delta t) = f_i(\mathbf x,t) + \Omega_i,

где \Omega_i — оператор столкновения.

Решётка D2Q9 (2D, 9 скоростей): покоящаяся частица, 4 «осевые» и 4 «диагональные». Скорость звука c_s^2 = 1/3. Веса w_i выбраны так, чтобы квадратура точно интегрировала моменты максвеллиана до нужного порядка (изотропия 4-го порядка): w_0=4/9, осевые 1/9, диагональные 1/36, \sum_i w_i = 1.

Макроскопические моменты:

\rho=\sum_i f_i,\qquad \rho\mathbf u=\sum_i \mathbf c_i f_i.

In [2]:

# --- Решётка D2Q9 -----------------------------------------------------------
# Порядок направлений: 0 — покой; 1..4 — оси (+x,+y,-x,-y); 5..8 — диагонали.
C = np.array([[0, 0], [1, 0], [0, 1], [-1, 0], [0, -1],
              [1, 1], [-1, 1], [-1, -1], [1, -1]], dtype=float)   # c_i = (c_ix, c_iy)
W = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])   # веса w_i
CS2 = 1.0 / 3.0                                                    # c_s^2
Q = 9
CX, CY = C[:, 0], C[:, 1]
OPP = np.array([0, 3, 4, 1, 2, 7, 8, 5, 6])                       # индекс противоположного c_i

# Самопроверка решёточных тождеств (нужны для корректной гидродинамики):
assert np.isclose(W.sum(), 1.0)                                   # сумма весов = 1
assert np.allclose((W[:, None] * C).sum(0), 0.0)                  # sum_i w_i c_i = 0
assert np.allclose((W[:, None, None] * C[:, :, None] * C[:, None, :]).sum(0),
                   CS2 * np.eye(2))                               # sum_i w_i c_ia c_ib = c_s^2 δ_ab
assert np.all(C[OPP] == -C)                                       # c_{opp(i)} = -c_i
print("Решётка D2Q9: тождества весов и изотропии 2-го порядка — OK")

Решётка D2Q9: тождества весов и изотропии 2-го порядка — OK

In [3]:

def make_grid(field_shape):
    """Удобный конструктор массива популяций формы (9, Ny, Nx)."""
    return np.zeros((Q,) + field_shape)


def macros(f):
    """Плотность rho (Ny,Nx) и скорость u (2,Ny,Nx) из популяций f (9,Ny,Nx)."""
    rho = f.sum(axis=0)
    jx = (CX[:, None, None] * f).sum(axis=0)
    jy = (CY[:, None, None] * f).sum(axis=0)
    u = np.stack([jx / rho, jy / rho])
    return rho, u

Рисунок: стенсиль D2Q9

Стрелки — направления \mathbf c_i, цвет/толщина — вес w_i (центральный узел имеет наибольший вес 4/9, диагонали — наименьший 1/36).

In [4]:

fig, ax = plt.subplots(figsize=(5.2, 5.2))
tier_color = {4/9: "#d1495b", 1/9: "#2e7d32", 1/36: "#1f6feb"}
for i in range(Q):
    col = tier_color[W[i]]
    if i == 0:
        ax.plot(0, 0, "o", ms=18, color=col, zorder=3)
        ax.annotate(f"$w_0={W[i]:.3f}$", (0, 0), textcoords="offset points",
                    xytext=(8, 8), fontsize=9)
    else:
        ax.annotate("", xy=(CX[i], CY[i]), xytext=(0, 0),
                    arrowprops=dict(arrowstyle="-|>", lw=1.5 + 40 * W[i], color=col))
        ax.annotate(f"$c_{{{i}}}$", (CX[i], CY[i]), textcoords="offset points",
                    xytext=(6 * np.sign(CX[i] + 0.1), 6 * np.sign(CY[i] + 0.1)), fontsize=9)
ax.set_xlim(-1.6, 1.6); ax.set_ylim(-1.6, 1.6); ax.set_aspect("equal")
ax.set_title("Решётка D2Q9: дискретные скорости $\\mathbf{c}_i$ и веса $w_i$")
ax.set_xlabel("$c_x$ [lu/ts]"); ax.set_ylabel("$c_y$ [lu/ts]")
from matplotlib.lines import Line2D
ax.legend(handles=[Line2D([0], [0], color=c, lw=3, label=f"$w={k:.4f}$")
                   for k, c in tier_color.items()], loc="upper right", fontsize=8)
_finish(fig, "01_stencil_d2q9.png")

2. Равновесие: полиномиальное и энтропийное

Полиномиальное равновесие (разложение максвеллиана по Эрмиту до O(u^2)), стандарт для BGK: $$ f_i^{\mathrm{eq,poly}} = w_i,\rho\left[1 + \frac{\mathbf c_i\cdot\mathbf u}{c_s^2}

• \frac{(\mathbf c_i\cdot\mathbf u)^2}{2c_s^4} - \frac{\mathbf u^2}{2c_s^2}\right]. $$ Коэффициенты при c_s^2=1/3: 1/c_s^2=3, 1/(2c_s^4)=4.5, 1/(2c_s^2)=1.5.

Энтропийное равновесие (Ansumali–Karlin): точный максимум дискретной H-функции при связях \sum f^{eq}=\rho, \sum \mathbf c f^{eq}=\rho\mathbf u. Для D2Q9/D3Q27 имеет замкнутую факторизованную форму ⚠ (сверить с оригиналом): $$ f_i^{\mathrm{eq}} = \rho,w_i\prod_{\alpha\in{x,y}} \underbrace{\left(2-\sqrt{1+3u_\alpha^2}\right)}{\text{префактор}} \left(\frac{2u\alpha+\sqrt{1+3u_\alpha^2}}{1-u_\alpha}\right)^{c_{i\alpha}}. $$ Здесь 3 = 1/c_s^2. Энтропийная форма точно воспроизводит \rho и \rho\mathbf u (проверим ниже) и, в отличие от полинома, гарантирует положительность и H-теорему.

In [5]:

def feq_poly(rho, u):
    """Полиномиальное равновесие O(u^2). rho:(Ny,Nx), u:(2,Ny,Nx) -> (9,Ny,Nx)."""
    cu = CX[:, None, None] * u[0] + CY[:, None, None] * u[1]   # c_i · u
    usq = u[0] ** 2 + u[1] ** 2
    return W[:, None, None] * rho * (1 + cu / CS2 + cu ** 2 / (2 * CS2 ** 2) - usq / (2 * CS2))


def feq_entropic(rho, u):
    """Энтропийное равновесие (product-form). Клип скорости для численной безопасности."""
    ux = np.clip(u[0], -0.95, 0.95)
    uy = np.clip(u[1], -0.95, 0.95)
    sx = np.sqrt(1 + 3 * ux ** 2)
    sy = np.sqrt(1 + 3 * uy ** 2)
    prefx, prefy = 2 - sx, 2 - sy                     # префактор по каждому измерению
    Bx = (2 * ux + sx) / (1 - ux)                     # основание степени по x
    By = (2 * uy + sy) / (1 - uy)
    base = rho * prefx * prefy                         # (Ny,Nx)
    Bxp = Bx[None] ** CX[:, None, None]                # (9,Ny,Nx): B_x^{c_ix}
    Byp = By[None] ** CY[:, None, None]                # (9,Ny,Nx): B_y^{c_iy}
    return W[:, None, None] * base[None] * Bxp * Byp

Самопроверка равновесий

Проверяем, что энтропийное равновесие точно даёт нужные нулевой и первый моменты (\rho, \rho\mathbf u), а второй момент при малом \mathrm{Ma} близок к \rho(c_s^2\delta_{\alpha\beta}+u_\alpha u_\beta) (тензор давления Навье–Стокса).

In [6]:

rng = np.random.default_rng(0)
rho_t = np.ones((1, 1)) * (1.0 + 0.1 * rng.standard_normal((1, 1)))
u_t = np.stack([0.06 * rng.standard_normal((1, 1)), 0.06 * rng.standard_normal((1, 1))])
fe = feq_entropic(rho_t, u_t)
rho_m, u_m = macros(fe)
Pxx = (CX[:, None, None] ** 2 * fe).sum(0)
err_rho = abs(rho_m - rho_t).max()
err_u = abs(u_m - u_t).max()
err_P = abs(Pxx - rho_t * (CS2 + u_t[0] ** 2)).max()
print(f"Энтропийное равновесие: |Σf-ρ|={err_rho:.2e}, |Σcf-ρu|={err_u:.2e} (точно),")
print(f"                       |Σc_x²f - ρ(c_s²+u_x²)|={err_P:.2e} (мало при малом Ma)")
assert err_rho < 1e-12 and err_u < 1e-12
print("Моменты равновесия — OK")

Энтропийное равновесие: |Σf-ρ|=2.22e-16, |Σcf-ρu|=5.03e-17 (точно),
                       |Σc_x²f - ρ(c_s²+u_x²)|=3.00e-09 (мало при малом Ma)
Моменты равновесия — OK

3. Разложение популяций f = k + s + h через моментные проекторы

Сердце KBC — разбиение каждого 9-вектора популяций на три части по моментным подпространствам:

• k_i — кинематическая часть, несёт только сохраняющиеся моменты (\rho, j_x, j_y);
• s_i — сдвиговая часть, несёт девиаторный тензор напряжений (N=\Pi_{xx}-\Pi_{yy} и \Pi_{xy});
• h_i — высшие («ghost») моды: след 2-го момента (объёмная вязкость), 3-й и 4-й моменты.

Строим базис из 9 независимых мономов-моментов \phi^{(a)}_i и матрицу M_{ai}=\phi^{(a)}(\mathbf c_i), тогда m = M f — вектор моментов. Проектор на группу моментов G: \;P_G = M^{-1} D_G M, где D_G — диагональная маска (единицы на моментах из G). По построению P_k+P_s+P_h = I, поэтому k+s+h=f точно. Назначение моментов в s vs h задаёт вариант KBC (ниже — KBC-N1: в s только девиаторный стресс).

In [7]:

def moment_matrix():
    """Строки — мономы-моменты, столбцы — направления. m = M @ f."""
    M = np.zeros((Q, Q))
    M[0] = 1.0                       # ρ
    M[1] = CX                        # j_x
    M[2] = CY                        # j_y
    M[3] = 3 * (CX ** 2 + CY ** 2) - 2   # след 2-го момента (энергия), центрированный
    M[4] = CX ** 2 - CY ** 2         # N — разность нормальных напряжений
    M[5] = CX * CY                   # Π_xy — сдвиговое напряжение
    M[6] = CX ** 2 * CY              # 3-й момент
    M[7] = CX * CY ** 2              # 3-й момент
    M[8] = CX ** 2 * CY ** 2         # 4-й момент
    return M

M = moment_matrix()
Minv = np.linalg.inv(M)              # существует: 9 мономов независимы на D2Q9

K_MOMENTS = [0, 1, 2]               # сохраняющиеся: ρ, j_x, j_y      -> k
S_MOMENTS = [4, 5]                  # девиаторный стресс: N, Π_xy      -> s  (вариант KBC-N1)
H_MOMENTS = [3, 6, 7, 8]           # след + 3-й + 4-й моменты          -> h

def projector(moment_idx):
    D = np.zeros((Q, Q))
    D[moment_idx, moment_idx] = 1.0
    return Minv @ D @ M

P_k, P_s, P_h = projector(K_MOMENTS), projector(S_MOMENTS), projector(H_MOMENTS)

# Самопроверка: полнота (P_k+P_s+P_h=I) и что k несёт ровно сохраняющиеся моменты.
assert np.allclose(P_k + P_s + P_h, np.eye(Q))
ftest = rng.random((Q, 3, 4))
k = np.einsum("ij,jab->iab", P_k, ftest)
s = np.einsum("ij,jab->iab", P_s, ftest)
h = np.einsum("ij,jab->iab", P_h, ftest)
assert np.allclose(k + s + h, ftest)                       # разложение полное
# у части k моменты N, Π и высшие — нулевые (несёт только ρ, j_x, j_y):
mk = np.einsum("ai,ijl->ajl", M, k)                        # моменты части k
assert np.allclose(mk[S_MOMENTS], 0.0) and np.allclose(mk[H_MOMENTS], 0.0)
print("Разложение k+s+h=f: макс. невязка =", abs(k + s + h - ftest).max())
print("Проекторы P_k+P_s+P_h=I и состав части k — OK")

Разложение k+s+h=f: макс. невязка = 1.1102230246251565e-15
Проекторы P_k+P_s+P_h=I и состав части k — OK

4. Оператор столкновений KBC и стабилизатор \gamma

Поскольку f и f^{eq} имеют одинаковые сохраняющиеся моменты, неравновесие \,f-f^{eq}\, не содержит кинематической части: \Delta k = 0, значит f-f^{eq}=\Delta s+\Delta h, где \Delta s = P_s(f-f^{eq}), \Delta h = P_h(f-f^{eq}).

Столкновение KBC:

f_i^{\text{после}} = f_i - \beta\,\bigl(2\,\Delta s_i + \gamma\,\Delta h_i\bigr).

Множитель 2 перед \Delta s фиксирует сдвиговую вязкость через \beta=1/(2\tau); высшие моды релаксируются отдельным \gamma.

Стабилизатор $\gamma$ выбирается локально из условия неувеличения дискретной энтропии (\Delta H = 0 во 2-м порядке). Результат — замкнутая формула ⚠ (сверить с оригиналом): $$ \gamma = \frac{1}{\beta} - \Bigl(2-\frac{1}{\beta}\Bigr) \frac{\langle \Delta s,|,\Delta h\rangle}{\langle \Delta h,|,\Delta h\rangle}, \qquad \langle X|Y\rangle = \sum_i \frac{X_i Y_i}{f_i^{eq}}. $$ При \gamma=2: 2\Delta s+2\Delta h = 2(f-f^{eq}), и KBC сводится к BGK с \omega=2\beta: f^{\text{после}}=f-2\beta(f-f^{eq}). Адаптивный \gamma добавляет/убавляет диссипацию высших мод там, где это нужно для устойчивости.

In [8]:

def collide_kbc(f, feq, beta, gamma_fixed=None):
    """KBC-столкновение. Возвращает (f_после, поле gamma)."""
    dfneq = f - feq
    ds = np.einsum("ij,jab->iab", P_s, dfneq)        # Δs
    dh = dfneq - ds                                  # Δh = (f-feq) - Δs, т.к. Δk=0
    if gamma_fixed is not None:
        gamma = np.full(f.shape[1:], float(gamma_fixed))
    else:
        inv = 1.0 / feq
        num = (ds * dh * inv).sum(0)                 # ⟨Δs|Δh⟩
        den = (dh * dh * inv).sum(0)                 # ⟨Δh|Δh⟩
        with np.errstate(divide="ignore", invalid="ignore"):
            gamma = np.where(den > 1e-14,
                             1.0 / beta - (2.0 - 1.0 / beta) * num / den,
                             2.0)                     # у равновесия (Δh→0) берём предел BGK
    fpost = f - beta * (2.0 * ds + gamma[None] * dh)
    return fpost, gamma


def collide_bgk(f, feq, beta):
    """Классический BGK: f - ω(f-feq), ω=2β."""
    return f - 2.0 * beta * (f - feq)

Самопроверка: предел \gamma=2 совпадает с BGK

При \gamma=2 оператор KBC обязан в точности равняться f-2\beta(f-f^{eq}).

In [9]:

rho0 = np.ones((6, 7))
u0 = np.stack([0.05 * rng.standard_normal((6, 7)), 0.05 * rng.standard_normal((6, 7))])
feq0 = feq_entropic(rho0, u0)
f0 = feq0 + 0.01 * rng.standard_normal((Q, 6, 7))      # слегка неравновесное состояние
f0 *= rho0 / f0.sum(0)                                  # сохранить ρ
beta0 = 0.7
fk2, _ = collide_kbc(f0, feq_entropic(*macros(f0)), beta0, gamma_fixed=2.0)
fbgk = collide_bgk(f0, feq_entropic(*macros(f0)), beta0)
diff = abs(fk2 - fbgk).max()
# diff — машинный ноль (масштаб |f|·eps ≈ 1e-18); важно лишь, что он ≪ порога 1e-12.
print(f"max|KBC(γ=2) - BGK| = {diff:.1e}  →  машинный ноль (порог assert 1e-12), операторы совпадают")
assert diff < 1e-12
print("Предел γ=2 → BGK — OK")

max|KBC(γ=2) - BGK| = 6.9e-18  →  машинный ноль (порог assert 1e-12), операторы совпадают
Предел γ=2 → BGK — OK

4b. Популярно: что делает каждая формула (математика + физика)

Здесь каждая формула KBC объясняется «на пальцах»: что это математически, какая за этим физика потока, и как это считается на каждом шаге.

---

Популяции f_i(\mathbf x,t) и шаг «перенос + столкновение». Математика: в каждом узле сетки храним 9 чисел — «сколько газа летит» по каждому из 9 направлений \mathbf c_i. Физика: как 9 встречных потоков на перекрёстке. Как считается: за шаг потоки сначала разъезжаются в соседние узлы (перенос, np.roll), затем в каждом узле перемешиваются (столкновение). Перенос точный и дешёвый; вся «физика жидкости» сидит в столкновении.

---

Моменты \rho=\sum_i f_i, \rho\mathbf u=\sum_i \mathbf c_i f_i. Математика: сумма всех потоков = плотность; взвешенная по направлениям = импульс. Физика: «сколько вещества» и «куда оно в среднем движется». Как считается: два суммирования по 9 числам — мгновенно. Отсюда привычные \rho и скорость \mathbf u.

---

Равновесие f_i^{\mathrm{eq}}(\rho,\mathbf u). Математика: распределение по направлениям, максимизирующее энтропию при заданных \rho,\mathbf u. Физика: «к чему стремится газ, если оставить в покое» — локальный максвелловский «штиль». Как считается: по замкнутой формуле из \rho,\mathbf u. Столкновение всегда тянет f к f^{eq} — это и есть вязкое стремление к равновесию.

---

H-функция H=\sum_i f_i\ln(f_i/w_i). Математика: дискретный аналог энтропии. Физика: мера беспорядка/необратимости; по второму началу не должна самопроизвольно убывать. Как считается: одно суммирование. В KBC H — физический компас устойчивости: из условия «энтропия не падает» выводится стабилизатор \gamma.

---

Разложение f = k + s + h. Математика: проекция 9-вектора на три группы моментов. Физика: три «слоя» движения:

• k — перенос массы и импульса (поток как целое);
• s — вязкий сдвиг (трение слоёв друг о друга — задаёт вязкость);
• h — мелкая «рябь» высших мод (мелкомасштабные флуктуации, где копится неустойчивость на грубой сетке).

Как считается: умножение на заранее посчитанные проекторы P_s,P_h. Идея KBC: релаксировать эти слои с разной скоростью.

---

Столкновение f \leftarrow f - \beta\,(2\,\Delta s + \gamma\,\Delta h). Математика: \Delta s,\Delta h — отклонения сдвига и ряби от равновесия. Физика: вязкий сдвиг гасим фиксированным темпом 2\beta (он задаёт вязкость и Re); рябь — отдельным темпом \gamma\beta, ровно чтобы не нарушить энтропию. Как считается: две линейные комбинации 9-векторов на узел. Множитель 2 перед \Delta s воспроизводит правильную вязкость; при \gamma=2 всё сворачивается в обычный BGK.

---

Стабилизатор $\gamma = \dfrac{1}{\beta} - \Bigl(2-\dfrac{1}{\beta}\Bigr) \dfrac{\langle \Delta s|\Delta h\rangle}{\langle \Delta h|\Delta h\rangle}$. Математика: единственное \gamma, при котором изменение энтропии за шаг равно нулю (во 2-м порядке). Физика: локальный термостат потока в реальном времени — в каждой точке на каждом шаге подкручивает, насколько гасить «рябь», чтобы энтропия не росла. В спокойных зонах \gamma\approx 2 (как BGK); у резких градиентов (кромки тел, вихри) \gamma отклоняется и добавляет ровно столько диссипации, сколько нужно, не давая симуляции взорваться. Это «implicit LES» — встроенная, без подгоночных констант, модель мелкомасштабной турбулентности. Как считается: два скалярных произведения на узел → одно число \gamma.

---

Энтропийное скалярное произведение \langle X|Y\rangle=\sum_i X_iY_i/f_i^{eq}. Математика: взвешенное (на 1/f^{eq}) скалярное произведение — кривизна энтропии у равновесия. Физика: измеряет, насколько «рябь» \Delta h выровнена со сдвигом \Delta s; по этому решается, сколько диссипации добавить. Как считается: поэлементное произведение и сумма по 9 направлениям.

---

Вязкость \nu = c_s^2\bigl(\tfrac{1}{2\beta}-\tfrac12\bigr), \beta=1/(2\tau). Математика: связь темпа релаксации сдвига с вязкостью. Физика: чем медленнее гасится сдвиг (меньше \beta), тем «гуще» жидкость; через \nu задаётся \mathrm{Re}=UL/\nu. Как считается: один раз при настройке (выбрали \nu под нужный Re → получили \beta).

5. Перенос (streaming) и связь вязкости с \beta

Перенос — сдвиг каждой популяции на её вектор \mathbf c_i (здесь — периодически через np.roll). Вязкость: \nu = c_s^2\bigl(\tfrac{1}{2\beta}-\tfrac12\bigr), то есть \tau=\nu/c_s^2+\tfrac12, \beta=1/(2\tau). Число Рейнольдса \mathrm{Re}=UL/\nu; число Маха \mathrm{Ma}=U/c_s держим \lesssim 0.1.

In [10]:

def stream(f):
    """Периодический перенос популяций на c_i. f:(9,Ny,Nx)."""
    fs = np.empty_like(f)
    for i in range(Q):
        fs[i] = np.roll(f[i], shift=(int(CY[i]), int(CX[i])), axis=(0, 1))
    return fs


def nu_to_beta(nu):
    tau = nu / CS2 + 0.5
    return 1.0 / (2.0 * tau)


def H_function(f):
    """Дискретная H-функция (энтропия) H = Σ f_i ln(f_i/w_i)."""
    fp = np.maximum(f, 1e-300)
    return (fp * np.log(fp / W[:, None, None])).sum(0)


print("Пример: ν=0.01 → β =", round(nu_to_beta(0.01), 4), ", τ =", round(0.01 / CS2 + 0.5, 4))

Пример: ν=0.01 → β = 0.9434 , τ = 0.53

6. Алгоритм одного шага KBC

для каждого узла:
  1) макро-моменты:   rho = sum f_i,   rho*u = sum c_i f_i
  2) равновесие:      f_i^eq  (энтропийное product-form)
  3) неравновесие:    Ds = P_s (f - f^eq),   Dh = (f - f^eq) - Ds
  4) стабилизатор:    gamma = 1/beta - (2 - 1/beta) <Ds|Dh> / <Dh|Dh>
  5) столкновение:    f <- f - beta (2 Ds + gamma Dh)
  6) перенос:         f_i(x + c_i) <- f_i(x)
  7) граничные условия (bounce-back / Zou-He / sponge)

In [11]:

fig, ax = plt.subplots(figsize=(4.3, 6.2))
ax.set_axis_off()
steps_txt = [
    ("Макро-моменты\n$\\rho=\\sum f_i,\\ \\rho\\mathbf{u}=\\sum \\mathbf{c}_i f_i$", "#e8f0fe"),
    ("Равновесие $f_i^{eq}$\n(энтропийное)", "#e6f4ea"),
    ("Неравновесие\n$\\Delta s=P_s(f-f^{eq}),\\ \\Delta h=(f-f^{eq})-\\Delta s$", "#fff4e5"),
    ("Стабилизатор $\\gamma$\n$=1/\\beta-(2-1/\\beta)\\frac{\\langle\\Delta s|\\Delta h\\rangle}{\\langle\\Delta h|\\Delta h\\rangle}$", "#fde7ef"),
    ("Столкновение\n$f\\leftarrow f-\\beta(2\\Delta s+\\gamma\\Delta h)$", "#f3e8fd"),
    ("Перенос\n$f_i(\\mathbf{x}+\\mathbf{c}_i)\\leftarrow f_i(\\mathbf{x})$", "#e8f0fe"),
    ("Граничные условия", "#eceff1"),
]
y = 0.95
for txt, col in steps_txt:
    ax.add_patch(plt.Rectangle((0.05, y - 0.1), 0.9, 0.085, transform=ax.transAxes,
                               facecolor=col, edgecolor="#555", lw=1.2))
    ax.text(0.5, y - 0.057, txt, transform=ax.transAxes, ha="center", va="center", fontsize=8.5)
    if y > 0.2:
        ax.annotate("", xy=(0.5, y - 0.118), xytext=(0.5, y - 0.1), xycoords="axes fraction",
                    arrowprops=dict(arrowstyle="-|>", color="#555"))
    y -= 0.13
ax.set_title("Шаг KBC")
_finish(fig, "06_algorithm_flow.png")

7. Демо-1: вихрь Тейлора–Грина (валидация вязкости и H-теоремы)

Точное решение несжимаемого Навье–Стокса в периодической области:

u_x=-U_0\cos(k x)\sin(k y),\quad u_y=U_0\sin(k x)\cos(k y),\quad k=2\pi/N.

Кинетическая энергия затухает как E(t)=E_0\,e^{-4\nu k^2 t}. Прогоняем KBC и измеряем \nu по наклону \ln E(t) — он должен совпасть с заданным \nu.

In [12]:

def run_tgv(N=96, U0=0.04, nu=0.0125, steps=10000, rec=50, n_snap=10, frame_every=42):
    k = 2 * np.pi / N
    xs = np.arange(N)
    X, Yc = np.meshgrid(xs, xs)             # X — индекс x (ось 1), Yc — индекс y (ось 0)
    ux = -U0 * np.cos(k * X) * np.sin(k * Yc)
    uy = U0 * np.sin(k * X) * np.cos(k * Yc)
    rho = np.ones((N, N))
    f = feq_entropic(rho, np.stack([ux, uy]))
    beta = nu_to_beta(nu)
    # моменты времени для n_snap снимков поля скорости через равные интервалы
    snap_at = sorted(set(int(round(v)) for v in np.linspace(0, steps, n_snap)))
    snaps, frames, ts, KE, Hs = [], [], [], [], []
    for t in range(steps + 1):
        rho, u = macros(f)
        f, _ = collide_kbc(f, feq_entropic(rho, u), beta)
        f = stream(f)
        need = (t in snap_at) or (t % rec == 0) or (t % frame_every == 0)
        if need:
            _, uu = macros(f)
            if t in snap_at:
                snaps.append((t, uu.copy()))     # снимки для статической панели
            if t % frame_every == 0:
                frames.append((t, uu.copy()))    # кадры для анимации
            if t % rec == 0:
                ts.append(t)
                KE.append(0.5 * np.mean(uu[0] ** 2 + uu[1] ** 2))
                Hs.append(H_function(f).mean())
    return dict(N=N, k=k, nu=nu, U0=U0, ts=np.array(ts), KE=np.array(KE),
                H=np.array(Hs), u=u, rho=rho, snaps=snaps, frames=frames)


tgv = run_tgv()
# Измеряем ν по наклону ln(KE): slope = -4 ν k^2
mask = tgv["KE"] > 0
slope = np.polyfit(tgv["ts"][mask], np.log(tgv["KE"][mask]), 1)[0]
nu_meas = -slope / (4 * tgv["k"] ** 2)
rel_err = abs(nu_meas - tgv["nu"]) / tgv["nu"]
print(f"TGV: nu_zadano={tgv['nu']:.5f}, nu_izmereno={nu_meas:.5f}, otn.oshibka={rel_err*100:.2f}%")
assert rel_err < 0.06, "TGV: измеренная вязкость отклонилась слишком сильно"
print("Валидация TGV (вязкость по затуханию энергии) — OK")

TGV: nu_zadano=0.01250, nu_izmereno=0.01249, otn.oshibka=0.05%
Валидация TGV (вязкость по затуханию энергии) — OK

In [13]:

def vorticity(u):
    duy_dx = (np.roll(u[1], -1, 1) - np.roll(u[1], 1, 1)) * 0.5
    dux_dy = (np.roll(u[0], -1, 0) - np.roll(u[0], 1, 0)) * 0.5
    return duy_dx - dux_dy


fig, axs = plt.subplots(1, 3, figsize=(13.5, 4.0))
om = vorticity(tgv["u"])
im = axs[0].imshow(om, origin="lower", cmap="RdBu_r")
axs[0].set_title("Завихренность $\\omega$ (TGV, срез)")
axs[0].set_xlabel("$x$ [lu]"); axs[0].set_ylabel("$y$ [lu]")
plt.colorbar(im, ax=axs[0], fraction=0.046)

axs[1].semilogy(tgv["ts"], tgv["KE"], "o", ms=3, label="KBC (измерено)")
axs[1].semilogy(tgv["ts"], tgv["KE"][0] * np.exp(-4 * tgv["nu"] * tgv["k"] ** 2 * tgv["ts"]),
                "-", color="#d1495b", label="$E_0 e^{-4\\nu k^2 t}$ (аналитика)")
axs[1].set_title("Затухание кин. энергии"); axs[1].set_xlabel("$t$ [ts]")
axs[1].set_ylabel("$E(t)$"); axs[1].legend()

axs[2].plot(tgv["ts"], tgv["H"], color="#2e7d32")
axs[2].set_title("H-функция (энтропия) во времени")
axs[2].set_xlabel("$t$ [ts]"); axs[2].set_ylabel("$\\langle H\\rangle$")
fig.suptitle("Демо-1: вихрь Тейлора–Грина (KBC, D2Q9)", y=1.04, fontweight="bold")
_finish(fig, "07_taylor_green.png")

8. Демо-2: устойчивость — двойной сдвиговый слой (BGK против KBC)

Классический тест Minion–Brown: тонкие сдвиговые слои на недоразрешённой сетке. Стандартный BGK при малой вязкости разрушается (отрицательные популяции, NaN), а KBC за счёт адаптивного \gamma остаётся устойчивым и даёт гладкие свёрнутые вихри — это и есть «implicit LES».

In [15]:

def init_shear(N, U0=0.04, delta=0.05, kappa=80.0):
    xs = np.arange(N) / N
    X, Yc = np.meshgrid(xs, xs)
    ux = np.where(Yc <= 0.5, U0 * np.tanh(kappa * (Yc - 0.25)),
                  U0 * np.tanh(kappa * (0.75 - Yc)))
    uy = delta * U0 * np.sin(2 * np.pi * (X + 0.25))
    return np.stack([ux, uy])


def run_shear(scheme, N=128, U0=0.04, nu=1.7e-4, steps=8000, rec=100, frame_every=33):
    u0 = init_shear(N, U0=U0)
    rho = np.ones((N, N))
    beta = nu_to_beta(nu)
    feq0 = feq_entropic if scheme == "kbc" else feq_poly
    f = feq0(rho, u0)
    ts, umax, frames = [], [], []
    blew_up_at = None
    last_u, last_g = u0, np.full((N, N), 2.0)
    for t in range(steps + 1):
        rho, u = macros(f)
        if not np.all(np.isfinite(u)) or np.abs(u).max() > 0.5:
            blew_up_at = t
            break
        if t % frame_every == 0:
            frames.append((t, u.copy()))         # кадры (только до разрушения)
        if scheme == "kbc":
            f, g = collide_kbc(f, feq_entropic(rho, u), beta)
            last_g = g
        else:
            f = collide_bgk(f, feq_poly(rho, u), beta)
        f = stream(f)
        last_u = u
        if t % rec == 0:
            ts.append(t); umax.append(np.abs(u).max())
    return dict(scheme=scheme, ts=np.array(ts), umax=np.array(umax), frames=frames,
                u=last_u, gamma=last_g, blew_up_at=blew_up_at, N=N, U0=U0, nu=nu)


Re_shear = 0.04 * 128 / 1.7e-4
print(f"Sdvigovyj sloj: Re ~ {Re_shear:.0f}, beta ~ {nu_to_beta(1.7e-4):.4f}")
res_bgk = run_shear("bgk")
res_kbc = run_shear("kbc")
print("BGK:", ("raz=" + str(res_bgk["blew_up_at"])) if res_bgk["blew_up_at"] else "ustojchiv")
print("KBC:", ("raz=" + str(res_kbc["blew_up_at"])) if res_kbc["blew_up_at"] else "ustojchiv")
assert res_kbc["blew_up_at"] is None, "KBC должен оставаться устойчивым"
print("Демонстрация устойчивости BGK vs KBC — OK")

Sdvigovyj sloj: Re ~ 30118, beta ~ 0.9990

BGK: raz=2011
KBC: ustojchiv
Демонстрация устойчивости BGK vs KBC — OK

In [16]:

fig, axs = plt.subplots(1, 3, figsize=(13.5, 4.0))
om_k = vorticity(res_kbc["u"])
im0 = axs[0].imshow(om_k, origin="lower", cmap="RdBu_r")
axs[0].set_title("KBC: завихренность (устойчиво)")
axs[0].set_xlabel("$x$ [lu]"); axs[0].set_ylabel("$y$ [lu]")
plt.colorbar(im0, ax=axs[0], fraction=0.046)

im1 = axs[1].imshow(res_kbc["gamma"] - 2.0, origin="lower", cmap="magma")
axs[1].set_title("Поле стабилизатора $\\gamma-2$")
axs[1].set_xlabel("$x$ [lu]"); axs[1].set_ylabel("$y$ [lu]")
plt.colorbar(im1, ax=axs[1], fraction=0.046)

axs[2].plot(res_kbc["ts"], res_kbc["umax"], color="#2e7d32", label="KBC")
axs[2].plot(res_bgk["ts"], res_bgk["umax"], color="#d1495b", label="BGK")
if res_bgk["blew_up_at"]:
    axs[2].axvline(res_bgk["blew_up_at"], color="#d1495b", ls="--", alpha=0.6)
axs[2].set_title("$\\max|\\mathbf{u}|$ во времени")
axs[2].set_xlabel("$t$ [ts]"); axs[2].set_ylabel("$\\max|\\mathbf{u}|$ [lu/ts]"); axs[2].legend()
fig.suptitle(f"Демо-2: сдвиговый слой, Re~{Re_shear:.0f} (BGK разрушается, KBC держит)",
             y=1.04, fontweight="bold")
_finish(fig, "08_shear_stability.png")

9. Демо-3: каверна с движущейся крышкой против эталона Ghia (Re=1000)

Квадратная каверна, верхняя стенка движется со скоростью U. На неподвижных стенках — bounce-back (отскок), на крышке — bounce-back с поправкой импульса подвижной стенки -2 w_i \rho\,(\mathbf c_i\cdot\mathbf u_w)/c_s^2. Сравниваем профили скорости по центральным линиям с данными Ghia, Ghia & Shin (1982). Здесь сетка намеренно грубая (для быстрого прогона), а bounce-back — первого порядка, поэтому совпадение качественное; оно улучшается измельчением сетки, числом шагов и переходом к ГУ второго порядка (Zou–He / Grad на крышке).

In [17]:

UP = [2, 5, 6]      # c_iy = +1  (вверх)
DOWN = [4, 7, 8]    # c_iy = -1  (вниз)
RIGHT = [1, 5, 8]   # c_ix = +1
LEFT = [3, 6, 7]    # c_ix = -1


def run_cavity(N=110, Re=1000.0, U=0.05, steps=14000, tol=1e-7, check=500, frame_every=58):
    nu = U * N / Re
    beta = nu_to_beta(nu)
    rho = np.ones((N, N))
    u = np.zeros((2, N, N))
    f = feq_entropic(rho, u)
    uw = np.array([U, 0.0])
    prev = None
    frames = []
    t = 0
    for t in range(steps):
        rho, u = macros(f)
        if t % frame_every == 0:
            frames.append((t, u.copy()))         # кадры становления течения
        fcol, _ = collide_kbc(f, feq_entropic(rho, u), beta)
        f = stream(fcol)
        # bounce-back: на каждой стенке неизвестны направления, входящие в область.
        for i in UP:        # низ (y=0)
            f[i, 0, :] = fcol[OPP[i], 0, :]
        for i in RIGHT:     # лево (x=0)
            f[i, :, 0] = fcol[OPP[i], :, 0]
        for i in LEFT:      # право (x=N-1)
            f[i, :, -1] = fcol[OPP[i], :, -1]
        for i in DOWN:      # верх (крышка) + поправка подвижной стенки
            io = OPP[i]     # io ударяет в крышку (идёт вверх)
            corr = 2.0 * W[io] * (CX[io] * uw[0] + CY[io] * uw[1]) / CS2
            f[i, -1, :] = fcol[io, -1, :] - corr
        if t % check == 0 and t > 0:
            cur = u.copy()
            if prev is not None and np.abs(cur - prev).max() < tol:
                break
            prev = cur
    rho, u = macros(f)
    return dict(N=N, Re=Re, U=U, u=u, steps=t, frames=frames)


cav = run_cavity()
print(f"Kaverna: Re={cav['Re']:.0f}, N={cav['N']}, ostanovleno na shage {cav['steps']}")

# Эталон Ghia, Ghia & Shin (1982), Re=1000.
ghia_y = np.array([0.0000, 0.0547, 0.0625, 0.0703, 0.1016, 0.1719, 0.2813, 0.4531,
                   0.5000, 0.6172, 0.7344, 0.8516, 0.9531, 0.9609, 0.9688, 0.9766, 1.0000])
ghia_u = np.array([0.0000, -0.18109, -0.20196, -0.22220, -0.29730, -0.38289, -0.27805,
                   -0.10648, -0.06080, 0.05702, 0.18719, 0.33304, 0.46604, 0.51117,
                   0.57492, 0.65928, 1.0000])
ghia_x = np.array([0.0000, 0.0625, 0.0703, 0.0781, 0.0938, 0.1563, 0.2266, 0.2344,
                   0.5000, 0.8047, 0.8594, 0.9063, 0.9453, 0.9531, 0.9609, 0.9688, 1.0000])
ghia_v = np.array([0.0000, 0.27485, 0.29012, 0.30353, 0.32627, 0.37095, 0.33075, 0.32235,
                   0.02526, -0.31966, -0.42665, -0.51550, -0.39188, -0.33714, -0.27669,
                   -0.21388, 0.0000])

Ncav = cav["N"]
yc = (np.arange(Ncav) + 0.5) / Ncav
xc = (np.arange(Ncav) + 0.5) / Ncav
u_centerline = cav["u"][0, :, Ncav // 2] / cav["U"]     # u_x(y) на вертикали x=0.5
v_centerline = cav["u"][1, Ncav // 2, :] / cav["U"]     # u_y(x) на горизонтали y=0.5

fig, axs = plt.subplots(1, 3, figsize=(14, 4.4))
spd = np.sqrt(cav["u"][0] ** 2 + cav["u"][1] ** 2) / cav["U"]
axs[0].imshow(spd, origin="lower", cmap="viridis", extent=[0, 1, 0, 1])
sx = np.linspace(0, 1, Ncav)
axs[0].streamplot(sx, sx, cav["u"][0], cav["u"][1], density=1.1, color="white", linewidth=0.6)
axs[0].set_title("Линии тока и $|\\mathbf{u}|/U$"); axs[0].set_xlabel("$x/L$"); axs[0].set_ylabel("$y/L$")

axs[1].plot(u_centerline, yc, "-", color="#1f6feb", label="KBC")
axs[1].plot(ghia_u, ghia_y, "o", ms=5, mfc="none", color="#d1495b", label="Ghia 1982")
axs[1].set_title("$u_x/U$ по вертикали $x=0.5$"); axs[1].set_xlabel("$u_x/U$")
axs[1].set_ylabel("$y/L$"); axs[1].legend()

axs[2].plot(xc, v_centerline, "-", color="#1f6feb", label="KBC")
axs[2].plot(ghia_x, ghia_v, "o", ms=5, mfc="none", color="#d1495b", label="Ghia 1982")
axs[2].set_title("$u_y/U$ по горизонтали $y=0.5$"); axs[2].set_xlabel("$x/L$")
axs[2].set_ylabel("$u_y/U$"); axs[2].legend()
fig.suptitle("Демо-3: каверна с движущейся крышкой, Re=1000 (KBC vs Ghia)",
             y=1.04, fontweight="bold")
_finish(fig, "09_cavity_ghia.png")

u_interp = np.interp(ghia_y, yc, u_centerline)
rms = np.sqrt(np.mean((u_interp - ghia_u) ** 2))
print(f"Kaverna: RMS-otklonenie u_x ot Ghia ~ {rms:.3f} (gruboe BB -> dopustimo)")

Kaverna: Re=1000, N=110, ostanovleno na shage 13999

Kaverna: RMS-otklonenie u_x ot Ghia ~ 0.167 (gruboe BB -> dopustimo)

9b. Анимации: зацикленное воспроизведение тестов (видео)

Для каждого теста собираем зацикленный GIF (бесконечный loop) из кадров, записанных во время прогона: вихрь Тейлора–Грина (затухание), сдвиговый слой (BGK против KBC бок о бок — видно момент разрушения BGK) и становление течения в каверне. Файлы сохраняются в anim/ (самостоятельные «видео»), а ниже встроены в ноутбук и проигрываются циклически. Скорость воспроизведения — fps кадров/с (физическое время между кадрами постоянно).

In [18]:

import io as _io
from PIL import Image
from IPython.display import HTML, display

ANIMDIR = os.path.join(os.path.dirname(os.path.abspath(__file__)) if "__file__" in globals() else ".", "anim")
os.makedirs(ANIMDIR, exist_ok=True)


def render_gif(n_frames, draw_fn, name, figsize, fps=18, dpi=80):
    """Отрисовать n_frames кадров (draw_fn(fig, i)) и сохранить зацикленный GIF."""
    path = os.path.join(ANIMDIR, name)
    fig = plt.figure(figsize=figsize, dpi=dpi)
    pil = []
    for i in range(n_frames):
        fig.clf()
        draw_fn(fig, i)
        buf = _io.BytesIO()
        fig.savefig(buf, format="png")        # фикс. размер кадра (figsize×dpi)
        buf.seek(0)
        pil.append(Image.open(buf).convert("RGB").convert("P", palette=Image.ADAPTIVE))
    plt.close(fig)
    pil[0].save(path, save_all=True, append_images=pil[1:],
                duration=int(1000 / fps), loop=0, optimize=True)   # loop=0 → бесконечно
    print("GIF:", os.path.relpath(path), f"({n_frames} кадров)")
    return path


def show_gif(path, width=460):
    """Встроить GIF в ноутбук (ссылкой на файл — авто-цикл в браузере/Jupyter)."""
    if IN_NOTEBOOK:
        display(HTML(f'<img src="anim/{os.path.basename(path)}" width="{width}"/>'))

Видео 1 — затухание вихря Тейлора–Грина

In [19]:

_tgv_fr = tgv["frames"]
_om0 = np.abs(vorticity(_tgv_fr[0][1])).max()
_KE0 = 0.5 * np.mean(_tgv_fr[0][1][0] ** 2 + _tgv_fr[0][1][1] ** 2)

def _draw_tgv(fig, i):
    t, u = _tgv_fr[i]
    ax = fig.add_subplot(111)
    ax.imshow(vorticity(u), origin="lower", cmap="RdBu_r", vmin=-_om0, vmax=_om0)
    ke = 0.5 * np.mean(u[0] ** 2 + u[1] ** 2) / _KE0
    ax.set_title(f"Тейлор–Грин · t={t} · $E/E_0$={ke:.2f}", fontsize=10)
    ax.set_xticks([]); ax.set_yticks([])
    fig.tight_layout()

show_gif(render_gif(len(_tgv_fr), _draw_tgv, "tgv_decay.gif", (4.6, 4.4), fps=24))

GIF: anim\tgv_decay.gif (239 кадров)

Видео 2 — сдвиговый слой: BGK (слева) против KBC (справа)

BGK замирает на последнем устойчивом кадре с пометкой «разрушился», KBC продолжает гладко сворачивать вихри.

In [20]:

_kf, _bf = res_kbc["frames"], res_bgk["frames"]
_om0s = np.percentile(np.abs(vorticity(_kf[len(_kf) // 2][1])), 99.5)

def _draw_shear(fig, i):
    ax1 = fig.add_subplot(1, 2, 1)
    ax2 = fig.add_subplot(1, 2, 2)
    dead = i >= len(_bf)
    tb, ub = _bf[min(i, len(_bf) - 1)]
    ax1.imshow(vorticity(ub), origin="lower", cmap="RdBu_r", vmin=-_om0s, vmax=_om0s)
    ax1.set_title("BGK — разрушился" if dead else f"BGK · t={tb}",
                  color="red" if dead else "black", fontsize=10)
    tk, uk = _kf[i]
    ax2.imshow(vorticity(uk), origin="lower", cmap="RdBu_r", vmin=-_om0s, vmax=_om0s)
    ax2.set_title(f"KBC · t={tk}", fontsize=10)
    for a in (ax1, ax2):
        a.set_xticks([]); a.set_yticks([])
    fig.tight_layout()

show_gif(render_gif(len(_kf), _draw_shear, "shear_bgk_vs_kbc.gif", (8.4, 4.3), fps=24), width=720)

GIF: anim\shear_bgk_vs_kbc.gif (243 кадров)

Видео 3 — становление течения в каверне (Re=1000)

In [21]:

_cf = cav["frames"]
_vmax = max(np.sqrt(u[0] ** 2 + u[1] ** 2).max() for _, u in _cf) / cav["U"]

def _draw_cav(fig, i):
    t, u = _cf[i]
    ax = fig.add_subplot(111)
    sp = np.sqrt(u[0] ** 2 + u[1] ** 2) / cav["U"]
    ax.imshow(sp, origin="lower", cmap="viridis", vmin=0, vmax=_vmax, extent=[0, 1, 0, 1])
    ax.set_title(f"Каверна Re=1000 · t={t} · $|\\mathbf{{u}}|/U$", fontsize=10)
    ax.set_xlabel("$x/L$"); ax.set_ylabel("$y/L$")
    fig.tight_layout()

show_gif(render_gif(len(_cf), _draw_cav, "cavity_transient.gif", (4.9, 4.5), fps=24))

GIF: anim\cavity_transient.gif (242 кадров)

Демо-4: обтекание твёрдых тел разной формы (2D)

Канал с втоком слева (Дирихле + плавный разгон), нуль-градиентным оттоком справа и стенками сверху/снизу. Тело задаётся маской и отражает популяции (bounce-back). За плохо обтекаемыми телами при Re≈150 формируется вихревая дорожка Кармана — поочерёдный срыв вихрей; KBC держит течение устойчивым на грубой сетке. Тела: цилиндр, квадрат, треугольник, профиль (под углом атаки). Ниже — анимация всех четырёх (≈10 c), статический срез завихренности и колебания скорости в следе (мерило частоты срыва).

In [22]:

def smoothstep(x):
    x = min(max(x, 0.0), 1.0)
    return x * x * (3.0 - 2.0 * x)


def solid_mask(Nx, Ny, shape, cx, cy, D):
    Y, X = np.mgrid[0:Ny, 0:Nx]
    if shape == "circle":
        m = (X - cx) ** 2 + (Y - cy) ** 2 <= (D / 2) ** 2
    elif shape == "square":
        m = (np.abs(X - cx) <= D / 2) & (np.abs(Y - cy) <= D / 2)
    elif shape == "triangle":          # вершиной против потока
        tt = (X - (cx - D / 2)) / D
        m = (tt >= 0) & (tt <= 1) & (np.abs(Y - cy) <= tt * (D / 2))
    elif shape == "airfoil":           # наклонный эллипс (угол атаки 18°)
        aoa = np.deg2rad(18.0); a = D * 0.95; b = D * 0.26
        xr = (X - cx) * np.cos(aoa) + (Y - cy) * np.sin(aoa)
        yr = -(X - cx) * np.sin(aoa) + (Y - cy) * np.cos(aoa)
        m = (xr / a) ** 2 + (yr / b) ** 2 <= 1
    else:
        m = np.zeros((Ny, Nx), bool)
    m[0, :] = True; m[-1, :] = True    # стенки канала (верх/низ)
    return m


def run_obstacle(shape, Nx=150, Ny=72, D=14, U=0.06, Re=150.0, steps=9600,
                 ramp=1000, frame_every=40):
    cx, cy = Nx // 4, Ny // 2
    solid = solid_mask(Nx, Ny, shape, cx, cy, D)
    beta = nu_to_beta(U * D / Re)
    f = feq_entropic(np.ones((Ny, Nx)), np.zeros((2, Ny, Nx)))
    frames, probe = [], []
    px, py = cx + 3 * D, cy           # зонд в следе
    ones_col = np.ones((Ny, 1))
    for t in range(steps + 1):
        rho, u = macros(f)
        fcol, _ = collide_kbc(f, feq_entropic(rho, u), beta)
        for i in range(Q):
            fcol[i, solid] = f[OPP[i], solid]              # bounce-back (тело + стенки)
        f = stream(fcol)
        Uin = U * smoothstep(t / ramp)                      # плавный разгон входа
        uin = np.zeros((2, Ny, 1)); uin[0, :, 0] = Uin
        f[:, 1:-1, 0] = feq_entropic(ones_col, uin)[:, 1:-1, 0]   # вток (Дирихле)
        f[:, 1:-1, -1] = f[:, 1:-1, -2]                            # отток (нуль-градиент)
        if t % frame_every == 0:
            _, uu = macros(f)
            frames.append((t, uu.copy()))
            probe.append(uu[1, py, px])
    return dict(shape=shape, solid=solid, frames=frames, probe=np.array(probe),
                Nx=Nx, Ny=Ny, U=U, Re=Re, D=D, cx=cx, cy=cy)


SHAPES = [("circle", "цилиндр"), ("square", "квадрат"),
          ("triangle", "треугольник"), ("airfoil", "профиль")]
obs = {key: run_obstacle(key) for key, _ in SHAPES}
for key, name in SHAPES:
    p = obs[key]["probe"]; half = p[len(p) // 2:]
    print(f"{name:12}: max|u_y(зонд)|={np.abs(half).max():.4f} (колебания → срыв вихрей)")

цилиндр     : max|u_y(зонд)|=0.0296 (колебания → срыв вихрей)
квадрат     : max|u_y(зонд)|=0.0140 (колебания → срыв вихрей)
треугольник : max|u_y(зонд)|=0.0348 (колебания → срыв вихрей)
профиль     : max|u_y(зонд)|=0.0076 (колебания → срыв вихрей)

In [23]:

# Анимация-монтаж всех четырёх тел (завихренность, тело — тёмное).
_keys = [k for k, _ in SHAPES]
_titles = {k: n for k, n in SHAPES}
_nf = min(len(obs[k]["frames"]) for k in _keys)
_vmax = max(np.percentile(np.abs(vorticity(obs[k]["frames"][_nf // 2][1])), 99.0) for k in _keys)
_ocmap = plt.get_cmap("RdBu_r").copy(); _ocmap.set_bad("#1b1b1b")

def _draw_obst(fig, i):
    for j, k in enumerate(_keys):
        ax = fig.add_subplot(2, 2, j + 1)
        t, u = obs[k]["frames"][i]
        vort = np.where(obs[k]["solid"], np.nan, vorticity(u))
        ax.imshow(vort, origin="lower", cmap=_ocmap, vmin=-_vmax, vmax=_vmax)
        ax.set_title(f"{_titles[k]} · t={t}", fontsize=9)
        ax.set_xticks([]); ax.set_yticks([])
    fig.suptitle("Обтекание тел (KBC, Re≈150): завихренность", fontsize=11)
    fig.tight_layout(rect=(0, 0, 1, 0.96))

show_gif(render_gif(_nf, _draw_obst, "obstacle_shapes.gif", (11, 6.2), fps=24), width=780)

GIF: anim\obstacle_shapes.gif (241 кадров)

In [24]:

# Статический срез завихренности в конце прогона + колебания в следе.
fig, axs = plt.subplots(2, 2, figsize=(12, 5.6), constrained_layout=True)
for ax, k in zip(axs.ravel(), _keys):
    t, u = obs[k]["frames"][-1]
    vort = np.where(obs[k]["solid"], np.nan, vorticity(u))
    ax.imshow(vort, origin="lower", cmap=_ocmap, vmin=-_vmax, vmax=_vmax)
    ax.set_title(f"{_titles[k]} (t={t})", fontsize=10)
    ax.set_xticks([]); ax.set_yticks([])
fig.suptitle("Обтекание тел: поле завихренности в конце прогона", fontweight="bold")
_finish(fig, "10_obstacle_shapes.png")

fig, ax = plt.subplots(figsize=(9.5, 3.4))
for k in _keys:
    ax.plot(np.arange(len(obs[k]["probe"])), obs[k]["probe"], label=_titles[k])
ax.set_title("Поперечная скорость в следе (зонд за телом): колебания = срыв вихрей")
ax.set_xlabel("кадр"); ax.set_ylabel("$u_y$ в зонде [lu/ts]"); ax.legend(fontsize=8, ncol=4)
_finish(fig, "10b_obstacle_probe.png")

KBC на сетках с измельчением у границ (grid refinement)

Может ли KBC работать с неоднородным разрешением? Да — но важно, как именно:

1. Чистый LBM живёт на РАВНОМЕРНОЙ решётке. Перенос f_i(\mathbf x+\mathbf c_i) требует, чтобы сосед был ровно на расстоянии \mathbf c_i — шаг сетки постоянен. Поэтому одна «кривая» неравномерная сетка для классического LBM не подходит.
2. Неоднородное разрешение делают блочным измельчением (multi-block / AMR): область делят на вложенные РАВНОМЕРНЫЕ блоки разного шага — крупный вдали, мелкий у тел и в следе. На стыке блоков поля интерполируют и перемасштабируют: при переходе на вдвое более мелкий блок \Delta x\!\to\!\Delta x/2, \Delta t\!\to\!\Delta t/2 (акустическое масштабирование, c_s сохраняется), время релаксации меняется как \tau_f = 2\tau_c - \tfrac12 (вязкость \nu постоянна), а неравновесная часть f^{neq} масштабируется множителем отношения шагов. ⚠ детали сверить с оригиналом.
3. KBC полностью совместим с измельчением. Это локальный оператор столкновения: на каждом блоке своё \beta, \gamma — поузловой как обычно. Энтропийная устойчивость особенно помогает на стыках крупный↔мелкий блок и в пограничном слое у стенок, где градиенты резкие. KBC применяли с измельчением и на телесно-подогнанных сетках (Dorschner, Chikatamarla, Karlin, JFM 2018; grid refinement — arXiv:1608.06915).
4. Альтернатива у кривых стенок — интерполированный bounce-back (Bouzidi) на равномерной решётке (доля стенки q). Измельчение добавляют там, где надо разрешить тонкий пограничный слой/след.

Полноценный AMR-решатель — отдельная большая работа. Ниже схема многоуровневого измельчения вокруг тела, наложенная на реальное поле |\omega| из расчёта цилиндра — видно, ГДЕ нужна мелкая сетка (резкие градиенты у тела и в следе).

In [25]:

import matplotlib.patches as mpatches

rc = obs["circle"]
_t, _u = rc["frames"][-1]
vmag = np.where(rc["solid"], np.nan, np.abs(vorticity(_u)))
fig, ax = plt.subplots(figsize=(11, 5.6))
ax.imshow(vmag, origin="lower", cmap="inferno",
          vmax=np.nanpercentile(vmag, 99))
cx, cy, D = rc["cx"], rc["cy"], rc["D"]

def _grid_block(x0, y0, w, h, step, color, label):
    ax.add_patch(mpatches.Rectangle((x0, y0), w, h, fill=False, edgecolor=color, lw=2.2))
    for xx in np.arange(x0, x0 + w + 0.1, step):
        ax.plot([xx, xx], [y0, y0 + h], color=color, lw=0.4, alpha=0.45)
    for yy in np.arange(y0, y0 + h + 0.1, step):
        ax.plot([x0, x0 + w], [yy, yy], color=color, lw=0.4, alpha=0.45)
    ax.text(x0 + 1, y0 + h - 4, label, color=color, fontsize=8, weight="bold")

_grid_block(0, 0, rc["Nx"] - 1, rc["Ny"] - 1, 16, "#4fc3f7", "уровень 0: Δx (крупный)")
_grid_block(cx - 1.5 * D, cy - 1.9 * D, 6.5 * D, 3.8 * D, 8, "#aed581", "уровень 1: Δx/2 (тело+след)")
_grid_block(cx - 1.1 * D, cy - 1.3 * D, 2.6 * D, 2.6 * D, 4, "#ff8a65", "уровень 2: Δx/4 (у тела)")
ax.set_xlim(0, rc["Nx"] - 1); ax.set_ylim(0, rc["Ny"] - 1)
ax.set_title("Схема блочного измельчения вокруг тела (поверх $|\\omega|$): "
             "сетка мельче там, где градиенты резче")
ax.set_xlabel("x [lu]"); ax.set_ylabel("y [lu]")
_finish(fig, "11c_refinement_schematic.png")

Реализация AMR-солвера (2 уровня) и демонстрация на цилиндре

По исследованию выше реализуем рабочий 2-уровневый блочный refinement:

• крупная сетка покрывает весь канал (с цилиндром и стенками);
• мелкий блок 2× охватывает цилиндр (то же тело bounce-back, но вдвое мельче);
• связка: \tau_f = 2\tau_c-\tfrac12 (вязкость сохраняется); на каждый крупный шаг — 2 мелких подшага; граница блока заполняется из крупной сетки (билинейная интерполяция \rho,\mathbf u + равновесие + перемасштаб f^{neq} на \tau_f/(2\tau_c)), а внутренние крупные узлы обновляются из мелкой (обратный масштаб); рестрикция — только по жидким узлам.

Проверка связки на отдельном гладком тесте (Тейлора–Грина с внутренним мелким блоком): схема устойчива, отн. отклонение AMR-крупной сетки от равномерной крупной ≈ 0.7 %, «шов» на границе блока ≈ 1 % — мелкий блок не вносит артефактов, KBC-γ работает на обоих уровнях. ⚠ Это минимальный AMR (линейная пространственная интерполяция, без временно́й интерполяции границы между подшагами); детали продакшн-схемы сверить с оригиналом (Lagrava et al. 2012; Dorschner et al., JFM 2018).

Ниже — анимация: крупное поле завихренности с наложенным мелким блоком 2× вокруг цилиндра (видно более тонкое разрешение у тела).

In [26]:

def run_amr_cylinder(Nx=120, Ny=60, D=10, U=0.06, Re=120.0, steps=5400,
                     ramp=1000, frame_every=22):
    cx, cy = Nx // 4, Ny // 2
    nu = U * D / Re
    tau_c = nu / CS2 + 0.5; beta_c = nu_to_beta(nu)
    tau_f = 2 * tau_c - 0.5; beta_f = 1.0 / (2 * tau_f)
    Rcf = tau_f / (2 * tau_c); Rfc = 1.0 / Rcf
    # мелкий блок (в крупных клетках), не касается стенок канала
    bx0 = max(cx - 2 * D, 2); bx1 = min(cx + 3 * D, Nx - 2)
    by0 = max(cy - 2 * D, 2); by1 = min(cy + 2 * D, Ny - 2)
    Wx, Wy = bx1 - bx0, by1 - by0
    Nfx, Nfy = 2 * Wx + 1, 2 * Wy + 1
    # маски тела
    Yc_, Xc_ = np.mgrid[0:Ny, 0:Nx]
    solid_c = (Xc_ - cx) ** 2 + (Yc_ - cy) ** 2 <= (D / 2) ** 2
    solid_c[0, :] = True; solid_c[-1, :] = True
    fcx, fcy = 2 * (cx - bx0), 2 * (cy - by0)
    Yf_, Xf_ = np.mgrid[0:Nfy, 0:Nfx]
    solid_f = (Xf_ - fcx) ** 2 + (Yf_ - fcy) ** 2 <= float(D) ** 2
    # билинейная интерполяция крупная → мелкая граница (индексы фиксированы)
    FX, FY = np.meshgrid(bx0 + 0.5 * np.arange(Nfx), by0 + 0.5 * np.arange(Nfy))
    ix0 = np.floor(FX).astype(int); iy0 = np.floor(FY).astype(int)
    ix1 = np.minimum(ix0 + 1, Nx - 1); iy1 = np.minimum(iy0 + 1, Ny - 1)
    wx = FX - ix0; wy = FY - iy0

    def interp(fld):
        return (fld[..., iy0, ix0] * (1 - wx) * (1 - wy) + fld[..., iy0, ix1] * wx * (1 - wy)
                + fld[..., iy1, ix0] * (1 - wx) * wy + fld[..., iy1, ix1] * wx * wy)

    slx, sly = slice(2, 2 * Wx, 2), slice(2, 2 * Wy, 2)
    block_fluid = ~solid_c[by0 + 1:by1, bx0 + 1:bx1]
    onecol = np.ones((Ny, 1))
    fc = feq_entropic(np.ones((Ny, Nx)), np.zeros((2, Ny, Nx)))
    ff = feq_entropic(np.ones((Nfy, Nfx)), np.zeros((2, Nfy, Nfx)))
    frames = []
    for t in range(steps + 1):
        rho, u = macros(fc)
        pre = fc.copy()
        fc = collide_kbc(fc, feq_entropic(rho, u), beta_c)[0]
        for i in range(Q):
            fc[i, solid_c] = pre[OPP[i], solid_c]
        fc = stream(fc)
        uin = np.zeros((2, Ny, 1)); uin[0, :, 0] = U * smoothstep(t / ramp)
        fc[:, 1:-1, 0] = feq_entropic(onecol, uin)[:, 1:-1, 0]
        fc[:, 1:-1, -1] = fc[:, 1:-1, -2]
        # крупная → мелкая (ghost)
        rcd, ucd = macros(fc); neqc = fc - feq_entropic(rcd, ucd)
        gh = (feq_entropic(interp(rcd), np.stack([interp(ucd[0]), interp(ucd[1])]))
              + Rcf * interp(neqc))
        for _ in range(2):
            rf, uf = macros(ff); pref = ff.copy()
            ff = collide_kbc(ff, feq_entropic(rf, uf), beta_f)[0]
            for i in range(Q):
                ff[i, solid_f] = pref[OPP[i], solid_f]
            ff = stream(ff)
            ff[:, 0, :] = gh[:, 0, :]; ff[:, -1, :] = gh[:, -1, :]
            ff[:, :, 0] = gh[:, :, 0]; ff[:, :, -1] = gh[:, :, -1]
        # мелкая → крупная (рестрикция, только жидкие узлы)
        rf, uf = macros(ff); neqf = ff - feq_entropic(rf, uf)
        newv = feq_entropic(rf[sly, slx], uf[:, sly, slx]) + Rfc * neqf[:, sly, slx]
        cur = fc[:, by0 + 1:by1, bx0 + 1:bx1]
        fc[:, by0 + 1:by1, bx0 + 1:bx1] = np.where(block_fluid[None], newv, cur)
        if t % frame_every == 0:
            _, uuc = macros(fc); _, uuf = macros(ff)
            frames.append((t, uuc.copy(), uuf.copy()))
    return dict(frames=frames, solid_c=solid_c, solid_f=solid_f, Nx=Nx, Ny=Ny,
                bx0=bx0, bx1=bx1, by0=by0, by1=by1, tau_c=tau_c, tau_f=tau_f, Rcf=Rcf)


amr = run_amr_cylinder()
print(f"AMR-цилиндр: крупная {amr['Nx']}x{amr['Ny']} + мелкий блок 2x, "
      f"tau_c={amr['tau_c']:.3f} tau_f={amr['tau_f']:.3f} Rcf={amr['Rcf']:.3f}, "
      f"кадров={len(amr['frames'])} (устойчиво)")

AMR-цилиндр: крупная 120x60 + мелкий блок 2x, tau_c=0.515 tau_f=0.530 Rcf=0.515, кадров=246 (устойчиво)

In [27]:

# Анимация AMR: крупное поле + наложенный мелкий блок (видно тонкое разрешение у тела).
_af = amr["frames"]
_vmaxA = np.nanpercentile(np.abs(np.where(amr["solid_c"], np.nan,
                                          vorticity(_af[len(_af) // 2][1]))), 99.0)
_cmapA = plt.get_cmap("RdBu_r").copy(); _cmapA.set_bad("#1b1b1b")

def _draw_amr(fig, i):
    t, uc, uf = _af[i]
    ax = fig.add_subplot(111)
    # interpolation="nearest" — чтобы были ВИДНЫ ячейки: крупные снаружи, вдвое
    # мельче в блоке (иначе imshow сглаживает и разница незаметна).
    ax.imshow(np.where(amr["solid_c"], np.nan, vorticity(uc)), origin="lower",
              cmap=_cmapA, vmin=-_vmaxA, vmax=_vmaxA, extent=(0, amr["Nx"], 0, amr["Ny"]),
              interpolation="nearest")
    omf = np.where(amr["solid_f"], np.nan, vorticity(uf))[1:-1, 1:-1]   # без ghost-рамки
    ax.imshow(omf, origin="lower", cmap=_cmapA, vmin=-_vmaxA, vmax=_vmaxA,
              extent=(amr["bx0"] + 0.5, amr["bx1"] - 0.5, amr["by0"] + 0.5, amr["by1"] - 0.5),
              interpolation="nearest")
    ax.add_patch(mpatches.Rectangle((amr["bx0"], amr["by0"]),
                 amr["bx1"] - amr["bx0"], amr["by1"] - amr["by0"],
                 fill=False, edgecolor="#ffe600", lw=3.0))
    ax.text(amr["bx0"] + 1, amr["by1"] - 4, "мелкий блок 2×", color="#ffe600", fontsize=9, weight="bold")
    ax.set_title(f"AMR: крупная сетка + мелкий блок 2× (жёлтая рамка) · t={t}", fontsize=10)
    ax.set_xticks([]); ax.set_yticks([])
    fig.tight_layout()

show_gif(render_gif(len(_af), _draw_amr, "amr_cylinder.gif", (9.2, 4.7), fps=24), width=780)

GIF: anim\amr_cylinder.gif (246 кадров)

AMR на цилиндре на БОЛЬШОЙ области — обе реализации (как на схеме исследования)

Эти две анимации прогнаны отдельным скриптом _amr_anims.py (большая область 150×70, длинный прогон ~7200 шагов → 240 кадров @24fps) и встроены готовыми GIF — чтобы тяжёлый счёт не раздувал рендер ноутбука. На большой области видна и грубая сетка (уровень 0, Δx), и учащённые зоны: зелёный блок 2× (тело+след) и оранжевый 4× (у тела). interpolation="nearest" — поэтому размеры ячеек различимы.

In [28]:

from IPython.display import HTML as _HTMLbig, display as _dispbig
if IN_NOTEBOOK:
    _dispbig(_HTMLbig(
        '<b>1) Одиночный 2× AMR на цилиндре (большая область):</b><br>'
        '<img src="anim/amr2x_cyl.gif" width="900"/>'
        '<br><b>2) Полный вложенный 2×+4× AMR на цилиндре (большая область):</b><br>'
        '<img src="anim/amr_nested_cyl.gif" width="900"/>'))

1) Одиночный 2× AMR на цилиндре (большая область):

2) Полный вложенный 2×+4× AMR на цилиндре (большая область):

Полный AMR: вложенный 2×+4× с временно́й интерполяцией и численное сравнение

Расширяем минимальный AMR до полного вложенного (как в исследовании — зона 2× для окружения/хвоста и 4× для завихрений у тела):

• три уровня: коарс → 2× (большой блок) → 4× (вложенный блок у тела);
• рекурсивный sub-cycling 1:2:4 — на 1 крупный шаг 2 подшага уровня-2× и 4 подшага уровня-4× (шаг по времени делится на 2 при каждом измельчении);
• временна́я интерполяция: граница каждого блока между подшагами линейно интерполируется во времени между состояниями родителя в моменты t и t+\Delta t (минимальный AMR держал границу постоянной — это вносило ошибку);
• перемасштаб на каждой границе: \tau_{child}=2\tau_{parent}-\tfrac12, f^{neq}\!\cdot\!\tau_{child}/(2\tau_{parent}).

Методика сравнения. Контролируемый тест — локализованный резкий вихрь (его ядро крупная сетка разрешает плохо). Эталон точности — равномерная мелкая 4× сетка. Сравниваем: равномерную крупную, минимальный 2× (без врем.интерп.), одиночный 4×, вложенный 2×+4× без и с врем.интерп. Производительность меряем в cell-updates (объём работы — аппаратно-независимо); wall-time в чистом numpy искажён Python- оверхедом мелких блоков, поэтому приводим его лишь справочно. Коэффициент производительности/точности: \eta = 1/(\varepsilon\cdot W_6), где \varepsilon — отн. L2-ошибка к эталону, W_6 — cell-updates в миллионах.

In [29]:

import time as _time

_Nc = 48; _L = 48; _xc = _yc = 24.0; _sig = 1.4; _A = 0.115
_nu_c = 0.008; _tau_c = _nu_c / CS2 + 0.5; _BSTEPS = 400

def _vfield(X, Y):
    rr = (X - _xc) ** 2 + (Y - _yc) ** 2; e = np.exp(-rr / (2 * _sig ** 2))
    return np.stack([-_A * (Y - _yc) / _sig ** 2 * e, _A * (X - _xc) / _sig ** 2 * e])

def _grid_v(N, h):
    xs = np.arange(N) * h; Xg, Yg = np.meshgrid(xs, xs); return _vfield(Xg, Yg)

def _adv(f, beta):
    r, u = macros(f); return stream(collide_kbc(f, feq_entropic(r, u), beta)[0])

def _patch(pN, a, b, r):
    Wc = b - a; Nf = r * Wc + 1; pos = a + np.arange(Nf) / r; FX, FY = np.meshgrid(pos, pos)
    x0 = np.floor(FX).astype(int); y0 = np.floor(FY).astype(int)
    x1 = np.minimum(x0 + 1, pN - 1); y1 = np.minimum(y0 + 1, pN - 1)
    return dict(Nf=Nf, a=a, b=b, Wc=Wc, r=r, x0=x0, y0=y0, x1=x1, y1=y1,
                tx=FX - x0, ty=FY - y0, sl=slice(r, r * Wc, r))

def _pint(fld, P):
    return (fld[..., P["y0"], P["x0"]] * (1 - P["tx"]) * (1 - P["ty"]) + fld[..., P["y0"], P["x1"]] * P["tx"] * (1 - P["ty"])
            + fld[..., P["y1"], P["x0"]] * (1 - P["tx"]) * P["ty"] + fld[..., P["y1"], P["x1"]] * P["tx"] * P["ty"])

def _ghost(pf, P, Rcf):
    r, u = macros(pf); neq = pf - feq_entropic(r, u)
    return feq_entropic(_pint(r, P), np.stack([_pint(u[0], P), _pint(u[1], P)])) + Rcf * _pint(neq, P)

def _fill(cf, gh):
    cf[:, 0, :] = gh[:, 0, :]; cf[:, -1, :] = gh[:, -1, :]; cf[:, :, 0] = gh[:, :, 0]; cf[:, :, -1] = gh[:, :, -1]; return cf

def _restrict(cf, pf, P, Rfc):
    r, u = macros(cf); neq = cf - feq_entropic(r, u); sl = P["sl"]
    pf[:, P["a"] + 1:P["b"], P["a"] + 1:P["b"]] = feq_entropic(r[sl, sl], u[:, sl, sl]) + Rfc * neq[:, sl, sl]; return pf

def _tauc(tp, r): return r * tp - (r - 1) / 2

def _run_uniform(N, h, tau, nsteps):
    beta = 1 / (2 * tau); f = feq_entropic(np.ones((N, N)), _grid_v(N, h)); t0 = _time.perf_counter()
    for _ in range(nsteps): f = _adv(f, beta)
    return macros(f)[1], _time.perf_counter() - t0, N * N * nsteps

def _init_patch(P, a_cc, h_par):
    pos = a_cc + np.arange(P["Nf"]) * (h_par / P["r"]); X, Y = np.meshgrid(pos, pos)
    return feq_entropic(np.ones((P["Nf"], P["Nf"])), _vfield(X, Y))

def _run_single(r, a, b, temporal):
    P = _patch(_Nc, a, b, r); tau_f = _tauc(_tau_c, r)
    b0, bf = 1 / (2 * _tau_c), 1 / (2 * tau_f); Rcf = tau_f / (r * _tau_c); Rfc = 1 / Rcf
    f0 = feq_entropic(np.ones((_Nc, _Nc)), _grid_v(_Nc, 1.0)); ff = _init_patch(P, a, 1.0)
    t0 = _time.perf_counter()
    for t in range(_BSTEPS):
        f0o = f0.copy(); f0 = _adv(f0, b0)
        gho = _ghost(f0o, P, Rcf); ghn = _ghost(f0, P, Rcf)
        for s in range(r):
            ff = _adv(ff, bf); ff = _fill(ff, (1 - (s + 1) / r) * gho + ((s + 1) / r) * ghn if temporal else ghn)
        f0 = _restrict(ff, f0, P, Rfc)
    return macros(f0)[1], _time.perf_counter() - t0, _Nc * _Nc * _BSTEPS + P["Nf"] ** 2 * r * _BSTEPS

def _run_nested(temporal):
    P1 = _patch(_Nc, 8, 40, 2); P2 = _patch(P1["Nf"], 16, 48, 2)
    t1 = _tauc(_tau_c, 2); t2 = _tauc(t1, 2)
    b0, b1, b2 = 1 / (2 * _tau_c), 1 / (2 * t1), 1 / (2 * t2)
    R01 = t1 / (2 * _tau_c); Rf01 = 1 / R01; R12 = t2 / (2 * t1); Rf12 = 1 / R12
    f0 = feq_entropic(np.ones((_Nc, _Nc)), _grid_v(_Nc, 1.0))
    f1 = _init_patch(P1, 8, 1.0); f2 = _init_patch(P2, 16.0, 0.5)
    t0 = _time.perf_counter()
    for t in range(_BSTEPS):
        f0o = f0.copy(); f0 = _adv(f0, b0)
        g1o = _ghost(f0o, P1, R01); g1n = _ghost(f0, P1, R01)
        for s1 in range(2):
            f1o = f1.copy(); f1 = _adv(f1, b1)
            f1 = _fill(f1, (1 - (s1 + 1) / 2) * g1o + ((s1 + 1) / 2) * g1n if temporal else g1n)
            g2o = _ghost(f1o, P2, R12); g2n = _ghost(f1, P2, R12)
            for s2 in range(2):
                f2 = _adv(f2, b2); f2 = _fill(f2, (1 - (s2 + 1) / 2) * g2o + ((s2 + 1) / 2) * g2n if temporal else g2n)
            f1 = _restrict(f2, f1, P2, Rf12)
        f0 = _restrict(f1, f0, P1, Rf01)
    cells = _Nc * _Nc * _BSTEPS + P1["Nf"] ** 2 * 2 * _BSTEPS + P2["Nf"] ** 2 * 4 * _BSTEPS
    return macros(f0)[1], _time.perf_counter() - t0, cells

# эталон: равномерная мелкая 4× (то же физ. время → 4×шагов)
_uref, _tref, _cref = _run_uniform(_Nc * 4, 0.25, _tauc(_tau_c, 4), 4 * _BSTEPS)
_uref = _uref[:, ::4, ::4]; _rms = np.sqrt(np.mean(_uref[0] ** 2 + _uref[1] ** 2))
def _err(u): return np.sqrt(np.mean((u[0] - _uref[0]) ** 2 + (u[1] - _uref[1]) ** 2)) / _rms

_bench = []
u, t, c = _run_uniform(_Nc, 1.0, _tau_c, _BSTEPS);  _bench.append(("uniform-coarse", _err(u), t, c))
u, t, c = _run_single(2, 12, 36, False);            _bench.append(("AMR 2× (мин., без вр.инт.)", _err(u), t, c))
u, t, c = _run_single(4, 12, 36, True);             _bench.append(("AMR 4× single (вр.инт.)", _err(u), t, c))
u, t, c = _run_nested(False);                       _bench.append(("AMR 2×+4× nested (без вр.инт.)", _err(u), t, c))
u, t, c = _run_nested(True);                        _bench.append(("AMR 2×+4× nested ПОЛНЫЙ", _err(u), t, c))

print(f"Эталон uniform-fine 4× (192²): время={_tref:.1f}c, cell-updates={_cref:.2e}")
print(f"{'метод':32}{'ε,%':>8}{'cell-upd':>11}{'η=1/(εW)':>10}{'×fine(работа)':>14}")
for name, e, t, c in _bench:
    eta = 1.0 / (e * c / 1e6)
    print(f"{name:32}{e*100:>8.2f}{c:>11.2e}{eta:>10.2f}{_cref/c:>14.1f}")
print("\nη = 1/(ε·W6), W6 — cell-updates в млн. ×fine(работа) — во сколько раз меньше работы, чем у эталона.")

Эталон uniform-fine 4× (192²): время=33.0c, cell-updates=5.90e+07
метод                                ε,%   cell-upd  η=1/(εW) ×fine(работа)
uniform-coarse                      5.28   9.22e+05     20.55          64.0
AMR 2× (мин., без вр.инт.)          2.56   2.84e+06     13.74          20.8
AMR 4× single (вр.инт.)             2.23   1.60e+07      2.81           3.7
AMR 2×+4× nested (без вр.инт.)      2.26   1.11e+07      4.00           5.3
AMR 2×+4× nested ПОЛНЫЙ             2.15   1.11e+07      4.20           5.3

η = 1/(ε·W6), W6 — cell-updates в млн. ×fine(работа) — во сколько раз меньше работы, чем у эталона.

In [30]:

# Графики сравнения: ошибка по методам и фронт «ошибка vs работа».
names = [b[0] for b in _bench]; errs = np.array([b[1] * 100 for b in _bench])
cells = np.array([b[3] for b in _bench]); cols = ["#9e9e9e", "#1f6feb", "#2e7d32", "#fb8c00", "#d1495b"]
fig, axs = plt.subplots(1, 2, figsize=(13.5, 4.4))
axs[0].barh(range(len(names)), errs, color=cols)
axs[0].set_yticks(range(len(names))); axs[0].set_yticklabels(names, fontsize=8); axs[0].invert_yaxis()
axs[0].set_xlabel("отн. L2-ошибка к эталону, %"); axs[0].set_title("Точность (меньше — лучше)")
for i, e in enumerate(errs): axs[0].text(e + 0.05, i, f"{e:.2f}", va="center", fontsize=8)
axs[1].scatter(cells, errs, c=cols, s=90, zorder=3, edgecolors="k")
for n, c2, e in zip(names, cells, errs):
    axs[1].annotate(n, (c2, e), fontsize=7, xytext=(5, 4), textcoords="offset points")
axs[1].axvline(_cref, color="#d1495b", ls="--", alpha=0.6); axs[1].text(_cref, errs.max()*0.9, " работа эталона", color="#d1495b", fontsize=8)
axs[1].set_xscale("log"); axs[1].set_xlabel("cell-updates (работа, лог)"); axs[1].set_ylabel("ошибка, %")
axs[1].set_title("Фронт «точность–работа»: AMR ближе к низ-лево (точно и дёшево)")
fig.suptitle("AMR: численное сравнение точности и производительности (KBC, тест-вихрь)", fontweight="bold")
_finish(fig, "11d_amr_comparison.png")

Численные результаты (типичные значения; точные печатаются ячейкой выше):

метод	ε (L2)	cell-updates	работа vs эталон
uniform-coarse	~5.3 %	~9.2·10⁵	×64 меньше
AMR 2× (мин., без вр.инт.)	~2.6 %	~2.8·10⁶	×21 меньше
AMR 4× single (вр.инт.)	~2.2 %	~1.6·10⁷	×3.7 меньше
AMR 2×+4× nested (без вр.инт.)	~2.3 %	~1.1·10⁷	×5.3 меньше
AMR 2×+4× nested ПОЛНЫЙ	~2.15 %	~1.1·10⁷	×5.3 меньше
uniform-fine 4× (эталон)	0 (опорная)	~5.9·10⁷	—

Выводы:

• Вложенность эффективна: полный 2×+4× точнее одиночного 4× (≈2.15 % против ≈2.2 %) при в ~1.45× меньшей работе — 4× только у ядра, 2× вокруг.
• Временна́я интерполяция помогает: при той же работе ошибка падает (≈2.26 %→ ≈2.15 %, на ~5 % относительно).
• AMR vs эталон: полный AMR даёт точность ~2 % за в ~5× меньше работы, чем равномерная мелкая сетка.
• Крупная сетка дёшева, но её точность (~5 %) принципиально хуже и не улучшается без измельчения.
• ⚠ wall-time в чистом numpy искажён Python-оверхедом мелких блоков; на GPU/в оптимизированном коде время ∝ cell-updates, и преимущество AMR по времени сохраняется.

Анимация полного вложенного AMR (видны три разрешения: коарс / 2× / 4×)

Тот же вложенный солвер записывает кадры. На анимации interpolation="nearest" — поэтому видны размеры ячеек: крупные снаружи, вдвое мельче в зелёном блоке (2×), ещё вдвое мельче в оранжевом (4×). Так наглядно работает разделение «2× для окружения / 4× у ядра».

In [31]:

def _run_nested_anim(steps=750, frame_every=3):
    P1 = _patch(_Nc, 8, 40, 2); P2 = _patch(P1["Nf"], 16, 48, 2)
    t1 = _tauc(_tau_c, 2); t2 = _tauc(t1, 2)
    b0, b1, b2 = 1 / (2 * _tau_c), 1 / (2 * t1), 1 / (2 * t2)
    R01 = t1 / (2 * _tau_c); Rf01 = 1 / R01; R12 = t2 / (2 * t1); Rf12 = 1 / R12
    f0 = feq_entropic(np.ones((_Nc, _Nc)), _grid_v(_Nc, 1.0))
    f1 = _init_patch(P1, 8, 1.0); f2 = _init_patch(P2, 16.0, 0.5)
    frames = []
    for t in range(steps + 1):
        f0o = f0.copy(); f0 = _adv(f0, b0)
        g1o = _ghost(f0o, P1, R01); g1n = _ghost(f0, P1, R01)
        for s1 in range(2):
            f1o = f1.copy(); f1 = _adv(f1, b1); f1 = _fill(f1, (1 - (s1 + 1) / 2) * g1o + ((s1 + 1) / 2) * g1n)
            g2o = _ghost(f1o, P2, R12); g2n = _ghost(f1, P2, R12)
            for s2 in range(2):
                f2 = _adv(f2, b2); f2 = _fill(f2, (1 - (s2 + 1) / 2) * g2o + ((s2 + 1) / 2) * g2n)
            f1 = _restrict(f2, f1, P2, Rf12)
        f0 = _restrict(f1, f0, P1, Rf01)
        if t % frame_every == 0:
            frames.append((t, macros(f0)[1], macros(f1)[1], macros(f2)[1]))
    return frames

_naf = _run_nested_anim()
_vmaxN = np.nanpercentile(np.abs(vorticity(_naf[0][1])), 99.0)
_cN = plt.get_cmap("RdBu_r").copy(); _cN.set_bad("#1b1b1b")

def _draw_nested(fig, i):
    t, u0, u1, u2 = _naf[i]
    ax = fig.add_subplot(111)
    ax.imshow(vorticity(u0), origin="lower", cmap=_cN, vmin=-_vmaxN, vmax=_vmaxN,
              extent=(0, _Nc, 0, _Nc), interpolation="nearest")               # коарс: крупные ячейки
    ax.imshow(vorticity(u1)[1:-1, 1:-1], origin="lower", cmap=_cN, vmin=-_vmaxN, vmax=_vmaxN,
              extent=(8.5, 39.5, 8.5, 39.5), interpolation="nearest")          # 2×: вдвое мельче
    ax.imshow(vorticity(u2)[1:-1, 1:-1], origin="lower", cmap=_cN, vmin=-_vmaxN, vmax=_vmaxN,
              extent=(16.25, 31.75, 16.25, 31.75), interpolation="nearest")    # 4×: ещё вдвое мельче
    ax.add_patch(mpatches.Rectangle((8, 8), 32, 32, fill=False, edgecolor="#7CFC00", lw=2.5))
    ax.add_patch(mpatches.Rectangle((16, 16), 16, 16, fill=False, edgecolor="#ff8a00", lw=2.5))
    ax.text(8.3, 40.4, "2× блок", color="#7CFC00", fontsize=9, weight="bold")
    ax.text(16.3, 32.4, "4× блок", color="#ff8a00", fontsize=9, weight="bold")
    ax.set_title(f"Полный вложенный AMR: коарс + 2× + 4× (видны размеры ячеек) · t={t}", fontsize=10)
    ax.set_xticks([]); ax.set_yticks([]); fig.tight_layout()

show_gif(render_gif(len(_naf), _draw_nested, "amr_nested.gif", (6.2, 6.2), fps=24), width=620)

GIF: anim\amr_nested.gif (251 кадров)

10. Сопутствующие технологии для полной симуляции

KBC задаёт только оператор столкновения. Полный решатель требует ещё ряда компонентов.

10.1 Граничные условия

• Bounce-back (отскок) — стенка без проскальзывания (no-slip): входящая популяция отражается, f_{\bar i}=f_i. Первый порядок (стенка на полпути между узлами). Для движущейся стенки добавляется импульс -2 w_i\rho\,(\mathbf c_i\cdot\mathbf u_w)/c_s^2 (использован в каверне).
• Zou–He — заданные скорость/давление на входе/выходе: недостающие популяции восстанавливаются из условий на \rho,\mathbf u и «bounce-back неравновесной части». Второй порядок.
• Grad / регуляризованные — реконструкция неравновесной части через тензор напряжений \Pi^{neq}; устойчивее при больших градиентах (рекомендуются с энтропийными моделями).
• Интерполированный bounce-back (Bouzidi, Yu–Mei–Luo–Shyy) — для кривых стенок: используется доля стенки q\in[0,1] из signed-distance field (SDF), две ветви q\!\ge\!0.5 и q\!<\!0.5.
• Sponge / absorbing-слой — у границ плавно подмешивается целевое равновесие, чтобы гасить отражения волн (использовался в прежней реализации на 6 гранях).

10.2 Объёмные силы (Guo)

Схема Guo et al. (2002): к столкновению добавляется $F_i = \bigl(1-\tfrac{\beta}{...}\bigr)w_i\Bigl[\tfrac{\mathbf c_i-\mathbf u}{c_s^2} +\tfrac{(\mathbf c_i\cdot\mathbf u)}{c_s^4}\mathbf c_i\Bigr]\cdot\mathbf F$, а скорость берётся со сдвигом \mathbf u\to\mathbf u+\tfrac{\mathbf F}{2\rho}. Даёт второй порядок и точное сохранение импульса.

10.3 Турбулентность: implicit LES vs Smagorinsky

• KBC — implicit LES: адаптивный \gamma добавляет диссипацию высших мод ровно там, где растут неразрешённые масштабы — параметр-free.
• Smagorinsky (явная LES): вихревая вязкость \nu_t=(C_s\Delta)^2|S|, |S|=\sqrt{2 S_{\alpha\beta}S_{\alpha\beta}} из неравновесного тензора напряжений; эффективное \tau_{\mathrm{eff}}=3(\nu+\nu_t)+\tfrac12. $C_s!\approx!0.1$–0.17. Иногда комбинируют с KBC у стенок.

10.4 Единицы и обезразмеривание (lu/ts ↔ СИ)

Выбираем масштабы \Delta x (м/lu), \Delta t (с/ts), \rho_0 (кг/м³). Тогда $$ \nu_{\text{СИ}} = \nu_{\text{lu}},\frac{\Delta x^2}{\Delta t},\quad \mathbf u_{\text{СИ}} = \mathbf u_{\text{lu}},\frac{\Delta x}{\Delta t},\quad p_{\text{СИ}} = c_s^2\rho_{\text{lu}},\rho_0\frac{\Delta x^2}{\Delta t^2}. $$ В решёточных единицах \Delta x=\Delta t=1 (узел и шаг — за единицу), поэтому скорости и вязкости — безразмерные «решёточные». Ограничение сжимаемости: $\mathrm{Ma}=|\mathbf u|/c_s\lesssim 0.1$–0.3 (ошибка \propto\mathrm{Ma}^2); отсюда Mach-лимитер шага \Delta t \le \mathrm{Ma}\,c_s\,\Delta x/|\mathbf u|.

10.5 Решётки, инициализация, измельчение сетки

• Решётки: D2Q9 (2D), D3Q27 (3D, полная изотропия 4-го порядка и галилеева инвариантность — предпочтительна для KBC/турбулентности).
• Инициализация: задать f=f^{eq}(\rho,\mathbf u); для согласованности с градиентами — добавить неравновесную часть из \Pi^{neq} (схема Mei).
• Измельчение сетки: многосеточные блоки с пространственно-временной интерполяцией; KBC устойчив к смене разрешения (нет явных подсеточных констант).

11. Перенос на D3Q27 и связь с целевой реализацией

Вся математика KBC переносится на 3D без изменений по сути: меняются решётка (27 скоростей, веса 8/27,\,2/27,\,1/54,\,1/216), число моментов (моментный базис 27×27) и состав сдвиговой части s (девиаторный тензор напряжений 3D — 5 независимых компонент). Формула \gamma и оператор f\leftarrow f-\beta(2\Delta s+\gamma\Delta h) идентичны.

In [32]:

# D3Q27: скорости и веса.
C3 = np.array([[x, y, z] for x in (-1, 0, 1) for y in (-1, 0, 1) for z in (-1, 0, 1)], dtype=float)
nz = (C3 != 0).sum(1)
W3 = np.select([nz == 0, nz == 1, nz == 2, nz == 3], [8/27, 2/27, 1/54, 1/216])
assert np.isclose(W3.sum(), 1.0)
assert np.allclose((W3[:, None, None] * C3[:, :, None] * C3[:, None, :]).sum(0), CS2 * np.eye(3))
print("D3Q27: Σw=1 и Σw c c = c_s² I — OK; число направлений:", len(C3))

fig = plt.figure(figsize=(6.4, 5.6))
ax = fig.add_subplot(111, projection="3d")
tier3 = {8/27: ("#d1495b", "покой"), 2/27: ("#2e7d32", "оси (×6)"),
         1/54: ("#1f6feb", "рёбра (×12)"), 1/216: ("#9c27b0", "углы (×8)")}
for wv, (col, lab) in tier3.items():
    m = np.isclose(W3, wv)
    ax.scatter(C3[m, 0], C3[m, 1], C3[m, 2], s=70 + 5000 * wv, color=col, label=f"{lab}, w={wv:.4f}",
               depthshade=True, edgecolors="k", linewidths=0.4)
ax.set_xlabel("$c_x$"); ax.set_ylabel("$c_y$"); ax.set_zlabel("$c_z$")
ax.set_title("Решётка D3Q27 (цель реализации)")
ax.legend(loc="upper left", fontsize=7, bbox_to_anchor=(0.0, 1.0))
_finish(fig, "11_d3q27.png")

D3Q27: Σw=1 и Σw c c = c_s² I — OK; число направлений: 27

Чем «настоящий KBC» отличается от прежнего подхода

Прежняя реализация (GPU-решатель) фактически использовала полиномиальное равновесие + TRT (две скорости релаксации) + α-лимитер положительности (ELBM-типа, по симметричной части) + явный Smagorinsky. Это устойчивая и рабочая схема, но не KBC. Переход к настоящему KBC:

Аспект	Настоящий KBC	Прежний подход (TRT + α)
Равновесие	энтропийное product-form	полиномиальное O(u^2)
Столкновение	энтропийный MRT, k/s/h	TRT (симм./антисимм.)
Стабилизация	\gamma из энтропии (высшие моды)	α-лимитер положительности (симм. часть)
H-теорема	встроена в вывод \gamma	контролируется отдельно (диагностика)
Турбулентность	implicit LES (параметр-free)	явный Smagorinsky
Решётка	D3Q27	D3Q27 (совпадает)

Дорожная карта реализации на GPU: (1) энтропийное f^{eq} (product-form); (2) моментные проекторы P_s,P_h для D3Q27 (предвычислить матрицы); (3) per-node \gamma по формуле; (4) коллизия f-\beta(2\Delta s+\gamma\Delta h); (5) опционально — комбинировать с пристеночным Smagorinsky.

12. Формулы под сверку с оригиналом (⚠) и литература

Под сверку с первоисточниками (доступ к части PDF был ограничен — проверить коэффициенты по оригиналам):

1. Энтропийное равновесие product-form: префактор (2-\sqrt{1+3u_\alpha^2}) и основание \frac{2u_\alpha+\sqrt{1+3u_\alpha^2}}{1-u_\alpha} — сверить вид и область применимости (Ansumali–Karlin; Chikatamarla–Karlin). Здесь проверено численно: даёт \rho,\rho\mathbf u точно.
2. Формула $\gamma$ =\tfrac1\beta-(2-\tfrac1\beta)\tfrac{\langle\Delta s|\Delta h\rangle}{\langle\Delta h|\Delta h\rangle} и порядок разложения \Delta H, из которого она получена (Bösch et al. 2015). Здесь проверено: предел \gamma=2 совпадает с BGK.
3. Назначение моментов в s vs h для вариантов KBC-N1…N4 (D2Q9/D3Q27) — сверить точный состав по таблицам оригинала.
4. Множитель 2 перед \Delta s и связь \beta=1/(2\tau) — сверить нормировку.

Литература (первоисточники):

• I. V. Karlin, F. Bösch, S. S. Chikatamarla, Gibbs' principle for the lattice-kinetic theory of fluid dynamics, Phys. Rev. E 90, 031302(R) (2014).
• F. Bösch, S. S. Chikatamarla, I. V. Karlin, Entropic multirelaxation lattice Boltzmann models for turbulent flows, Phys. Rev. E 92, 043309 (2015).
• S. S. Chikatamarla, I. V. Karlin, Entropy and Galilean invariance of lattice Boltzmann theories, Phys. Rev. Lett. (2006) — энтропийное равновесие.
• B. Dorschner, S. S. Chikatamarla, I. V. Karlin и др. — применения KBC (каналы, обтекание тел), измельчение сетки, подвижные границы.
• Z. Guo, C. Zheng, B. Shi, Discrete lattice effects on the forcing term in the lattice Boltzmann method, Phys. Rev. E 65, 046308 (2002) — форсинг.
• U. Ghia, K. N. Ghia, C. T. Shin (1982) — эталон каверны (валидация).

13. Итоги

• KBC = энтропийный MRT-LBM: сохраняющиеся моды не трогаем, сдвиг релаксируем фиксированным 2\beta (задаёт вязкость), высшие моды — адаптивным \gamma из условия неувеличения энтропии. При \gamma=2 — это в точности BGK.
• Зачем: безусловная нелинейная устойчивость и параметр-free implicit LES — считает турбулентность на грубых сетках там, где BGK разрушается (продемонстрировано на сдвиговом слое).
• Проверено численно: точность равновесия и разложения (до 10^{-15}), правильная вязкость (TGV, ошибка <0.1\%), устойчивость (сдвиговый слой), разумное совпадение с эталоном Ghia (каверна).
• Дальше: перенос на D3Q27 на GPU (Vulkan), сверка коэффициентов с оригинальными статьями, при необходимости — комбинация с Smagorinsky и точными ГУ (Grad/Zou–He, Bouzidi для кривых стенок).

In [33]:

print("\n=== Все самопроверки и валидации пройдены ===")

=== Все самопроверки и валидации пройдены ===