# --- # jupyter: # jupytext: # text_representation: # extension: .py # format_name: percent # format_version: '1.3' # kernelspec: # display_name: Python 3 # language: python # name: python3 # --- # %% [markdown] # # Метод 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). Где доступ был только к вторичным # > источникам, формула помечена **⚠ сверить с оригиналом**, а в конце собран # > список таких мест. Главный страж корректности здесь — **численная # > самопроверка** каждой формулы. # %% [markdown] # ## Сводная таблица обозначений # # Везде используются **решёточные единицы** (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]$ | # %% [markdown] # ## 0. Стиль графиков и общие настройки # # Подключаем библиотеки и задаём единый «симпатичный» стиль matplotlib. Код # написан так, чтобы файл работал **и как ноутбук** (графики встраиваются), **и # как обычный скрипт** `python kbc_lbm.py` (графики сохраняются в `figures/` для # валидации физики). # %% 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) # %% [markdown] # ## 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. $$ # %% # --- Решётка 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") # %% 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 # %% [markdown] # ### Рисунок: стенсиль D2Q9 # # Стрелки — направления $\mathbf c_i$, цвет/толщина — вес $w_i$ (центральный узел # имеет наибольший вес $4/9$, диагонали — наименьший $1/36$). # %% 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") # %% [markdown] # ## 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-теорему. # %% 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 # %% [markdown] # ### Самопроверка равновесий # # Проверяем, что **энтропийное** равновесие точно даёт нужные нулевой и первый # моменты ($\rho$, $\rho\mathbf u$), а второй момент при малом $\mathrm{Ma}$ близок # к $\rho(c_s^2\delta_{\alpha\beta}+u_\alpha u_\beta)$ (тензор давления Навье–Стокса). # %% 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") # %% [markdown] # ## 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$ только девиаторный стресс). # %% 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") # %% [markdown] # ## 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$ # добавляет/убавляет диссипацию высших мод там, где это нужно для устойчивости. # %% 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) # %% [markdown] # ### Самопроверка: предел $\gamma=2$ совпадает с BGK # При $\gamma=2$ оператор KBC обязан в точности равняться $f-2\beta(f-f^{eq})$. # %% 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") # %% [markdown] # ## 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$). # %% [markdown] # ## 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$. # %% 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)) # %% [markdown] # ## 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) / # 5) столкновение: f <- f - beta (2 Ds + gamma Dh) # 6) перенос: f_i(x + c_i) <- f_i(x) # 7) граничные условия (bounce-back / Zou-He / sponge) # ``` # %% 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") # %% [markdown] # ## 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$. # %% 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") # %% 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") # %% [markdown] # ### Распад вихря: 10 срезов через равные интервалы времени # # Поле завихренности $\omega=\partial_x u_y-\partial_y u_x$ в 10 моментов времени, # равномерно распределённых от $t=0$ до конца прогона. **Шкала цвета единая для всех # панелей** (фиксирована по начальной завихренности $\pm|\omega|_{\max}(0)$) — поэтому # наглядно видно именно *затухание*: пространственная структура (4 противоположно # вращающихся вихря) сохраняется, а амплитуда экспоненциально гаснет. В заголовке # каждой панели — шаг $t$ и доля кинетической энергии $E/E_0$. # %% snaps = tgv["snaps"] om0 = np.abs(vorticity(snaps[0][1])).max() # масштаб цвета по t=0 KE0 = 0.5 * np.mean(snaps[0][1][0] ** 2 + snaps[0][1][1] ** 2) # начальная энергия fig, axs = plt.subplots(2, 5, figsize=(15, 6.3), constrained_layout=True) im = None for ax, (tt, uu) in zip(axs.ravel(), snaps): im = ax.imshow(vorticity(uu), origin="lower", cmap="RdBu_r", vmin=-om0, vmax=om0) KEf = 0.5 * np.mean(uu[0] ** 2 + uu[1] ** 2) / KE0 ax.set_title(f"t={tt} · $E/E_0$={KEf:.2f}", fontsize=10) ax.set_xticks([]); ax.set_yticks([]) fig.colorbar(im, ax=axs, fraction=0.025, pad=0.01, label="завихренность $\\omega$ [1/ts]") fig.suptitle("Распад вихря Тейлора–Грина: завихренность на равных интервалах времени " "(единая шкала цвета)", fontweight="bold") _finish(fig, "07b_tgv_decay_slices.png") # %% [markdown] # ## 8. Демо-2: устойчивость — двойной сдвиговый слой (BGK против KBC) # # Классический тест Minion–Brown: тонкие сдвиговые слои на **недоразрешённой** # сетке. Стандартный BGK при малой вязкости разрушается (отрицательные # популяции, NaN), а KBC за счёт адаптивного $\gamma$ остаётся устойчивым и даёт # гладкие свёрнутые вихри — это и есть «implicit LES». # %% 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") # %% 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") # %% [markdown] # ## 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 на крышке). # %% 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)") # %% [markdown] # ## 9b. Анимации: зацикленное воспроизведение тестов (видео) # # Для каждого теста собираем **зацикленный GIF** (бесконечный `loop`) из кадров, # записанных во время прогона: вихрь Тейлора–Грина (затухание), сдвиговый слой # (BGK против KBC бок о бок — видно момент разрушения BGK) и становление течения в # каверне. Файлы сохраняются в `anim/` (самостоятельные «видео»), а ниже встроены # в ноутбук и проигрываются циклически. Скорость воспроизведения — `fps` кадров/с # (физическое время между кадрами постоянно). # %% 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'')) # %% [markdown] # ### Видео 1 — затухание вихря Тейлора–Грина # %% _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)) # %% [markdown] # ### Видео 2 — сдвиговый слой: BGK (слева) против KBC (справа) # BGK замирает на последнем устойчивом кадре с пометкой «разрушился», KBC # продолжает гладко сворачивать вихри. # %% _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) # %% [markdown] # ### Видео 3 — становление течения в каверне (Re=1000) # %% _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)) # %% [markdown] # ## Демо-4: обтекание твёрдых тел разной формы (2D) # # Канал с втоком слева (Дирихле + плавный разгон), нуль-градиентным оттоком справа и # стенками сверху/снизу. Тело задаётся маской и отражает популяции (**bounce-back**). # За плохо обтекаемыми телами при Re≈150 формируется **вихревая дорожка Кармана** — # поочерёдный срыв вихрей; KBC держит течение устойчивым на грубой сетке. Тела: цилиндр, # квадрат, треугольник, профиль (под углом атаки). Ниже — анимация всех четырёх (≈10 c), # статический срез завихренности и колебания скорости в следе (мерило частоты срыва). # %% 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} (колебания → срыв вихрей)") # %% # Анимация-монтаж всех четырёх тел (завихренность, тело — тёмное). _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) # %% # Статический срез завихренности в конце прогона + колебания в следе. 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") # %% [markdown] # ## 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|$ из расчёта цилиндра — # видно, ГДЕ нужна мелкая сетка (резкие градиенты у тела и в следе). # %% 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") # %% [markdown] # ## Реализация 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×** вокруг # цилиндра (видно более тонкое разрешение у тела). # %% 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: крупное поле + наложенный мелкий блок (видно тонкое разрешение у тела). _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) # %% [markdown] # ### AMR на цилиндре на БОЛЬШОЙ области — обе реализации (как на схеме исследования) # # Эти две анимации прогнаны отдельным скриптом **`_amr_anims.py`** (большая область # 150×70, длинный прогон ~7200 шагов → 240 кадров @24fps) и встроены готовыми GIF — # чтобы тяжёлый счёт не раздувал рендер ноутбука. На большой области видна и **грубая # сетка** (уровень 0, Δx), и **учащённые зоны**: зелёный блок **2× (тело+след)** и # оранжевый **4× (у тела)**. `interpolation="nearest"` — поэтому размеры ячеек различимы. # %% from IPython.display import HTML as _HTMLbig, display as _dispbig if IN_NOTEBOOK: _dispbig(_HTMLbig( '1) Одиночный 2× AMR на цилиндре (большая область):
' '' '
2) Полный вложенный 2×+4× AMR на цилиндре (большая область):
' '')) # %% [markdown] # ## Полный 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 в миллионах. # %% 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(работа) — во сколько раз меньше работы, чем у эталона.") # %% # Графики сравнения: ошибка по методам и фронт «ошибка 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") # %% [markdown] # **Численные результаты (типичные значения; точные печатаются ячейкой выше):** # # | метод | ε (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 по времени сохраняется. # %% [markdown] # ### Анимация полного вложенного AMR (видны три разрешения: коарс / 2× / 4×) # # Тот же вложенный солвер записывает кадры. На анимации `interpolation="nearest"` — # поэтому **видны размеры ячеек**: крупные снаружи, вдвое мельче в зелёном блоке (2×), # ещё вдвое мельче в оранжевом (4×). Так наглядно работает разделение «2× для # окружения / 4× у ядра». # %% 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) # %% [markdown] # ## 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), D3Q19 (3D, экономная), **D3Q27** (3D, полная изотропия # 4-го порядка и галилеева инвариантность — предпочтительна для KBC/турбулентности). # - **Инициализация**: задать $f=f^{eq}(\rho,\mathbf u)$; для согласованности с # градиентами — добавить неравновесную часть из $\Pi^{neq}$ (схема Mei). # - **Измельчение сетки**: многосеточные блоки с пространственно-временной # интерполяцией; KBC устойчив к смене разрешения (нет явных подсеточных констант). # %% [markdown] # ## 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)$ идентичны. # %% # 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") # %% [markdown] # ### Чем «настоящий 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. # %% [markdown] # ## 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) — эталон каверны (валидация). # %% [markdown] # ## 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 для кривых стенок). # %% print("\n=== Все самопроверки и валидации пройдены ===")