feat: 水平卧式流化床空床流场 GPU LBM 仿真 (100 格/m)
- 3D D3Q19 LBM (numba.cuda), 240x130x265=8.27M 单元 - 恒功率风机曲线 dp*Q=P 压力进口, 工况点涌现 - 多孔布风板: 2mm 厚, 6.25% 开孔, 5 度倾斜(-x) - 4 段 DN300 侧进风 / 567->1200 渐扩 / DN800 排风 - 输出 19 张切面/矢量/压差/流量/分区图 (4x 分辨率)
This commit is contained in:
@@ -0,0 +1,4 @@
|
||||
__pycache__/
|
||||
*.pyc
|
||||
*.pyo
|
||||
.ipynb_checkpoints/
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 1.0 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 2.7 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 415 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 406 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 554 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 2.0 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 1.8 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 1.5 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 1.8 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 1.8 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 3.4 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 899 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 623 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 1.7 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 2.9 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 1.6 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 716 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 7.2 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 1.3 MiB |
Binary file not shown.
@@ -0,0 +1,400 @@
|
||||
#!/usr/bin/env python
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
水平卧式流化床 空床流场 —— 分析图 (几何/矢量/压差/流量/均匀性/衰减/倾斜/分段放大)
|
||||
读取 fbed_fields.npz。
|
||||
|
||||
python tools/fbed_analysis.py
|
||||
"""
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
import matplotlib
|
||||
matplotlib.use("Agg")
|
||||
import matplotlib.pyplot as plt
|
||||
from matplotlib.patches import Rectangle
|
||||
|
||||
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
||||
import fbed_simulation as fs
|
||||
|
||||
OUT = fs.OUT
|
||||
BG, FG, MUTED, ACCENT = fs.BG, fs.FG, fs.MUTED, fs.ACCENT
|
||||
CMAP_FLOW = fs.CMAP_FLOW
|
||||
XC, YC, ZC = fs.XC, fs.YC, fs.ZC
|
||||
X0, X1, Y0, Y1, Z0, Z1 = fs.X0, fs.X1, fs.Y0, fs.Y1, fs.Z0, fs.Z1
|
||||
DX, SCALE_U, BED_L, BED_W, N_SEC, SEC_L = fs.DX, fs.SCALE_U, fs.BED_L, fs.BED_W, fs.N_SEC, fs.SEC_L
|
||||
PPM = fs.PPM
|
||||
NZ, NY, NX = fs.NZ, fs.NY, fs.NX
|
||||
|
||||
D = np.load(os.path.join(OUT, "fbed_fields.npz"))
|
||||
UX, UY, UZ, RHO, SOLID = D["ux"], D["uy"], D["uz"], D["rho"], D["solid"]
|
||||
PLATE = D["plate"]
|
||||
FLUID = ~SOLID
|
||||
UX, UY, UZ = UX * SCALE_U, UY * SCALE_U, UZ * SCALE_U
|
||||
SPEED = np.sqrt(UX ** 2 + UY ** 2 + UZ ** 2)
|
||||
P = (RHO - 1.0) / 3.0 * (1.2 * SCALE_U ** 2)
|
||||
CELL_A = DX * DX
|
||||
IZP = fs._zi(0.0)
|
||||
|
||||
|
||||
def zi(v):
|
||||
return fs._zi(v)
|
||||
|
||||
|
||||
def yi(v):
|
||||
return fs._yi(v)
|
||||
|
||||
|
||||
def xi(v):
|
||||
return fs._xi(v)
|
||||
|
||||
|
||||
def quiver_slice(ax, which, idx, step=3, scale=None, title="", vmax=None):
|
||||
if which == "xz":
|
||||
ext = [X0, X1, Z0, Z1]; xl, yl = "x (m)", "z (m)"
|
||||
U, V, M = UX[:, idx, :], UZ[:, idx, :], SPEED[:, idx, :]
|
||||
m = SOLID[:, idx, :]
|
||||
xs, ys = XC, ZC
|
||||
elif which == "yz":
|
||||
ext = [Y0, Y1, Z0, Z1]; xl, yl = "y (m)", "z (m)"
|
||||
U, V, M = UY[:, :, idx], UZ[:, :, idx], SPEED[:, :, idx]
|
||||
m = SOLID[:, :, idx]
|
||||
xs, ys = YC, ZC
|
||||
else:
|
||||
ext = [X0, X1, Y0, Y1]; xl, yl = "x (m)", "y (m)"
|
||||
U, V, M = UX[idx, :, :], UY[idx, :, :], SPEED[idx, :, :]
|
||||
m = SOLID[idx, :, :]
|
||||
xs, ys = XC, YC
|
||||
U = np.where(m, np.nan, U); V = np.where(m, np.nan, V); M = np.where(m, np.nan, M)
|
||||
ax.imshow(np.where(m, np.nan, M), origin="lower", extent=ext, cmap=CMAP_FLOW,
|
||||
vmin=0, vmax=vmax or np.nanpercentile(M, 99), interpolation="bilinear",
|
||||
aspect="equal", alpha=0.85)
|
||||
Xs, Ys = np.meshgrid(xs[::step], ys[::step])
|
||||
Us = U[::step, ::step]; Vs = V[::step, ::step]; Ms = M[::step, ::step]
|
||||
if scale is None:
|
||||
scale = 6.0 * max(ext[1] - ext[0], ext[3] - ext[2]) / max(np.nanmax(np.hypot(Us, Vs)), 1e-6)
|
||||
ax.quiver(Xs, Ys, Us, Vs, Ms, cmap="cool", scale=scale, width=0.0022,
|
||||
clim=(0, np.nanpercentile(M, 99)))
|
||||
ax.contour(np.linspace(ext[0], ext[1], m.shape[1]),
|
||||
np.linspace(ext[2], ext[3], m.shape[0]),
|
||||
m.astype(float), levels=[0.5], colors="#adbac7", linewidths=0.8)
|
||||
ax.set_xlim(ext[0], ext[1]); ax.set_ylim(ext[2], ext[3])
|
||||
ax.set_xlabel(xl, fontsize=9); ax.set_ylabel(yl, fontsize=9)
|
||||
ax.tick_params(labelsize=8)
|
||||
if title:
|
||||
ax.set_title(title, fontsize=12)
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------
|
||||
def fig_geometry():
|
||||
fig, axes = plt.subplots(1, 2, figsize=(17.0, 7.6))
|
||||
for ax, which, idx in [(axes[0], "xz", yi(0.0)), (axes[1], "yz", xi(1.2))]:
|
||||
if which == "xz":
|
||||
ext = [X0, X1, Z0, Z1]; sol = SOLID[:, idx, :]
|
||||
xl, yl = "x (m)", "z (m)"
|
||||
else:
|
||||
ext = [Y0, Y1, Z0, Z1]; sol = SOLID[:, :, idx]
|
||||
xl, yl = "y (m)", "z (m)"
|
||||
ax.imshow(np.where(sol, 0.55, 0.08), origin="lower", extent=ext, cmap="gray",
|
||||
vmin=0, vmax=1, interpolation="nearest", aspect="equal")
|
||||
ax.contour(np.linspace(ext[0], ext[1], sol.shape[1]),
|
||||
np.linspace(ext[2], ext[3], sol.shape[0]),
|
||||
sol.astype(float), levels=[0.5], colors="#adbac7", linewidths=1.0)
|
||||
ax.set_xlabel(xl, fontsize=10); ax.set_ylabel(yl, fontsize=10)
|
||||
# 标注 (在 x-z 图)
|
||||
ax = axes[0]
|
||||
ax.annotate("4 x DN300 侧进风\n(平行床板, 2500 m³/h/段)", xy=(0.05, -0.225), xytext=(-0.35, 0.9),
|
||||
fontsize=10, color="#7ee787",
|
||||
arrowprops=dict(arrowstyle="->", color="#7ee787", lw=1.6))
|
||||
ax.annotate("布风板 2mm孔 5.6% 5°倾斜", xy=(1.2, 0.0), xytext=(1.55, -0.35),
|
||||
fontsize=10, color="#58a6ff",
|
||||
arrowprops=dict(arrowstyle="->", color="#58a6ff", lw=1.4))
|
||||
ax.annotate("顶部仅 DN800 排风\n(无进风)", xy=(1.2, 2.1), xytext=(1.45, 1.55),
|
||||
fontsize=10.5, color="#ff7b72",
|
||||
arrowprops=dict(arrowstyle="->", color="#ff7b72", lw=1.6))
|
||||
ax.set_title("几何 / 边界示意 (x-z, y=0)", fontsize=12)
|
||||
axes[1].set_title("几何 / 边界示意 (y-z, x=1.2)", fontsize=12)
|
||||
fig.suptitle("空床几何与边界条件：底部 4 段侧进风，顶部仅排风", fontsize=15, y=0.98)
|
||||
fig.tight_layout(rect=[0, 0, 1, 0.95])
|
||||
fig.savefig(os.path.join(OUT, "fbed-geometry.png"), dpi=540)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-geometry.png")
|
||||
|
||||
|
||||
def fig_vectors():
|
||||
fig, axes = plt.subplots(2, 2, figsize=(16.0, 12.5))
|
||||
quiver_slice(axes[0, 0], "xz", yi(0.0), step=3, title="(a) 垂直切面 y=0 速度矢量")
|
||||
quiver_slice(axes[0, 1], "yz", xi(1.2), step=3, title="(b) 垂直切面 x=1.2 速度矢量")
|
||||
quiver_slice(axes[1, 0], "xy", zi(0.15), step=3, title="(c) 布风板近区 z=0.15 速度矢量")
|
||||
quiver_slice(axes[1, 1], "xy", zi(0.90), step=3, title="(d) 沸腾段 z=0.90 速度矢量")
|
||||
fig.suptitle("空床流场速度矢量图（箭头方向=流向, 颜色=速度大小 m/s）", fontsize=15, y=0.99)
|
||||
fig.tight_layout(rect=[0, 0, 1, 0.97])
|
||||
fig.savefig(os.path.join(OUT, "fbed-vectors.png"), dpi=520)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-vectors.png")
|
||||
|
||||
|
||||
def fig_tilt():
|
||||
fig = plt.figure(figsize=(16.0, 10.5))
|
||||
gs = fig.add_gridspec(2, 2)
|
||||
ax = fig.add_subplot(gs[0, 0])
|
||||
iz = zi(0.05)
|
||||
quiver_slice(ax, "xy", iz, step=4, vmax=2.5,
|
||||
title="(a) 布风板正上方 z=0.05 矢量（放大看 +x 倾斜）")
|
||||
ax.set_xlim(0.6, 1.8); ax.set_ylim(-0.30, 0.30)
|
||||
|
||||
ax = fig.add_subplot(gs[0, 1])
|
||||
zs, ang, um = [], [], []
|
||||
for z in np.arange(0.03, 0.60, 0.02):
|
||||
iz = zi(z)
|
||||
m = FLUID[iz] & (np.abs(YC)[:, None] < BED_W / 2)
|
||||
if m.sum() < 10:
|
||||
continue
|
||||
ux_m, uz_m = UX[iz][m].mean(), UZ[iz][m].mean()
|
||||
zs.append(z); ang.append(np.degrees(np.arctan2(ux_m, uz_m))); um.append(ux_m)
|
||||
ax.plot(zs, ang, "-o", color=ACCENT, ms=4, label="实际偏角")
|
||||
ax.axhline(5.0, color="#ff7b72", ls="--", lw=1.2, label="理论 +5°")
|
||||
ax.axhline(0.0, color=MUTED, lw=0.8)
|
||||
ax.set_xlabel("高度 z (m)", fontsize=10); ax.set_ylabel("流向偏角 (°)", fontsize=10)
|
||||
ax.set_title("(b) 近板流向偏角随高度衰减", fontsize=12)
|
||||
ax.legend(fontsize=9); ax.grid(alpha=0.2)
|
||||
|
||||
ax = fig.add_subplot(gs[1, 0])
|
||||
zs2, uxs, uzs = [], [], []
|
||||
for z in np.arange(0.03, 0.60, 0.02):
|
||||
iz = zi(z)
|
||||
m = FLUID[iz] & (np.abs(YC)[:, None] < BED_W / 2)
|
||||
if m.sum() < 10:
|
||||
continue
|
||||
zs2.append(z); uxs.append(UX[iz][m].mean()); uzs.append(UZ[iz][m].mean())
|
||||
ax.plot(zs2, uzs, "-o", ms=4, color="#39d0d8", label="垂直分量 <uz>")
|
||||
ax.plot(zs2, uxs, "-s", ms=4, color="#ff7b72", label="水平分量 <ux>（5° 驱动）")
|
||||
ax.set_xlabel("高度 z (m)", fontsize=10); ax.set_ylabel("面平均速度 (m/s)", fontsize=10)
|
||||
ax.set_title("(c) 近板面平均速度分量", fontsize=12)
|
||||
ax.legend(fontsize=9); ax.grid(alpha=0.2)
|
||||
|
||||
ax = fig.add_subplot(gs[1, 1])
|
||||
xline = np.linspace(0.02, BED_L - 0.02, 60)
|
||||
iz = zi(0.06)
|
||||
prof = []
|
||||
for xv in xline:
|
||||
i = xi(xv)
|
||||
m = FLUID[iz, :, i]
|
||||
prof.append(UX[iz, :, i][m].mean() if m.sum() else np.nan)
|
||||
ax.plot(xline, prof, color="#ff7b72", lw=1.6)
|
||||
ax.axhline(0, color=MUTED, lw=0.8)
|
||||
ax.set_xlabel("x (m)", fontsize=10); ax.set_ylabel("<ux> (m/s)", fontsize=10)
|
||||
ax.set_title("(d) 近板水平速度沿长度分布（整体 +x 偏置）", fontsize=12)
|
||||
ax.grid(alpha=0.2)
|
||||
fig.suptitle("布风板 5° 倾斜对近板风向的影响", fontsize=15, y=0.99)
|
||||
fig.tight_layout(rect=[0, 0, 1, 0.97])
|
||||
fig.savefig(os.path.join(OUT, "fbed-tilt.png"), dpi=520)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-tilt.png")
|
||||
|
||||
|
||||
def fig_uniformity():
|
||||
fig, axes = plt.subplots(1, 3, figsize=(17.0, 5.4))
|
||||
qs = []
|
||||
for k in range(N_SEC):
|
||||
ix0, ix1 = xi(0.6 * k), xi(0.6 * (k + 1))
|
||||
q = UZ[IZP, :, ix0:ix1][FLUID[IZP, :, ix0:ix1]].sum() * CELL_A * 3600.0
|
||||
qs.append(q)
|
||||
ax = axes[0]
|
||||
bars = ax.bar([f"第{k+1}段" for k in range(N_SEC)], qs, color=ACCENT, alpha=0.85)
|
||||
ax.axhline(2500, color="#ff7b72", ls="--", lw=1.4, label="目标 2500 m³/h")
|
||||
for b, q in zip(bars, qs):
|
||||
ax.text(b.get_x() + b.get_width() / 2, q + 40, f"{q:.0f}", ha="center", fontsize=9)
|
||||
ax.set_ylabel("布风板风量 (m³/h)", fontsize=10); ax.legend(fontsize=9)
|
||||
ax.set_title("(a) 各区段风量", fontsize=12); ax.grid(alpha=0.2, axis="y")
|
||||
|
||||
ax = axes[1]
|
||||
prof = []
|
||||
xline = np.linspace(0.02, BED_L - 0.02, 120)
|
||||
for xv in xline:
|
||||
i = xi(xv)
|
||||
m = FLUID[IZP + 3, :, i]
|
||||
prof.append(UZ[IZP + 3, :, i][m].mean() if m.sum() else np.nan)
|
||||
ax.plot(xline, prof, color="#39d0d8", lw=1.8)
|
||||
ax.axhline(2500.0 / (BED_L * BED_W) / 3600.0 * 3600.0, color=MUTED, ls=":", lw=1)
|
||||
for k in range(1, N_SEC):
|
||||
ax.axvline(0.6 * k, color="#30363d", lw=0.8)
|
||||
ax.set_xlabel("x (m)", fontsize=10); ax.set_ylabel("近板 <uz> (m/s)", fontsize=10)
|
||||
ax.set_title("(b) 近板上升速度沿长度分布（均匀性）", fontsize=12); ax.grid(alpha=0.2)
|
||||
|
||||
ax = axes[2]
|
||||
yline = np.linspace(-BED_W / 2 + 0.02, BED_W / 2 - 0.02, 60)
|
||||
profy = []
|
||||
for yv in yline:
|
||||
j = yi(yv)
|
||||
m = FLUID[IZP + 3, j, :]
|
||||
profy.append(UZ[IZP + 3, j, :][m].mean() if m.sum() else np.nan)
|
||||
ax.plot(yline, profy, color="#7ee787", lw=1.8)
|
||||
ax.set_xlabel("y (m)", fontsize=10); ax.set_ylabel("近板 <uz> (m/s)", fontsize=10)
|
||||
ax.set_title("(c) 近板上升速度沿宽度分布（均匀性）", fontsize=12); ax.grid(alpha=0.2)
|
||||
fig.suptitle("布风板流量均匀性", fontsize=15, y=1.0)
|
||||
fig.tight_layout(rect=[0, 0, 1, 0.95])
|
||||
fig.savefig(os.path.join(OUT, "fbed-uniformity.png"), dpi=520)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-uniformity.png")
|
||||
|
||||
|
||||
def fig_decay():
|
||||
fig, axes = plt.subplots(1, 2, figsize=(15.0, 7.0))
|
||||
zs, qz, uz_avg, spd_avg, area = [], [], [], [], []
|
||||
for iz in range(NZ):
|
||||
m = FLUID[iz]
|
||||
if m.sum() < 20:
|
||||
continue
|
||||
zs.append(ZC[iz])
|
||||
qz.append(UZ[iz][m].sum() * CELL_A * 3600.0)
|
||||
uz_avg.append(UZ[iz][m].mean())
|
||||
spd_avg.append(SPEED[iz][m].mean())
|
||||
area.append(m.sum() * CELL_A)
|
||||
ax = axes[0]
|
||||
ax.plot(uz_avg, zs, "-o", ms=3, color="#39d0d8", label="面平均 <uz>")
|
||||
ax.plot(spd_avg, zs, "-s", ms=3, color="#ff7b72", label="面平均 |U|")
|
||||
ax.axvline(0, color=MUTED, lw=0.8)
|
||||
ax.axhline(0.0, color="#30363d", lw=0.8)
|
||||
ax.set_xlabel("速度 (m/s)", fontsize=10); ax.set_ylabel("高度 z (m)", fontsize=10)
|
||||
ax.set_title("(a) 沿高度气速衰减", fontsize=12); ax.legend(fontsize=9); ax.grid(alpha=0.2)
|
||||
|
||||
ax = axes[1]
|
||||
ax.plot(qz, zs, "-o", ms=3, color=ACCENT, label="通过水平面的净风量")
|
||||
ax.axvline(10000, color="#ff7b72", ls="--", lw=1.4, label="总风量 10000 m³/h")
|
||||
ax.set_xlabel("风量 (m³/h)", fontsize=10); ax.set_ylabel("高度 z (m)", fontsize=10)
|
||||
ax.set_title("(b) 沿高度风量守恒（连续性）", fontsize=12)
|
||||
ax.legend(fontsize=9); ax.grid(alpha=0.2)
|
||||
fig.suptitle("沿高度气速衰减与风量守恒", fontsize=15, y=0.99)
|
||||
fig.tight_layout(rect=[0, 0, 1, 0.96])
|
||||
fig.savefig(os.path.join(OUT, "fbed-decay.png"), dpi=520)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-decay.png")
|
||||
|
||||
|
||||
def fig_zoom_sections():
|
||||
fig, axes = plt.subplots(4, 1, figsize=(11.0, 16.0))
|
||||
for k, ax in enumerate(axes):
|
||||
quiver_slice(ax, "xz", yi(0.0), step=2, title=f"第 {k+1} 段 x∈[{0.6*k:.1f}, {0.6*(k+1):.1f}] m")
|
||||
ax.set_xlim(0.6 * k - 0.02, 0.6 * (k + 1) + 0.02)
|
||||
fig.suptitle("按区段剖分的 x-z 局部放大（含风箱进风/布风板/沸腾段）", fontsize=15, y=0.995)
|
||||
fig.tight_layout(rect=[0, 0, 1, 0.975])
|
||||
fig.savefig(os.path.join(OUT, "fbed-zoom-sections.png"), dpi=500)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-zoom-sections.png")
|
||||
|
||||
|
||||
def fig_pressure_drop():
|
||||
fig = plt.figure(figsize=(16.0, 10.5))
|
||||
gs = fig.add_gridspec(2, 2)
|
||||
ax = fig.add_subplot(gs[0, 0])
|
||||
ext = [X0, X1, Z0, Z1]
|
||||
Pxz = np.where(SOLID[:, yi(0), :], np.nan, P[:, yi(0), :])
|
||||
lim = np.nanpercentile(np.abs(Pxz), 99)
|
||||
im = ax.imshow(Pxz, origin="lower", extent=ext, cmap="RdBu_r", vmin=-lim, vmax=lim,
|
||||
aspect="equal", interpolation="bilinear")
|
||||
ax.axhline(0, color="#7ee787", lw=1.0, ls="--")
|
||||
ax.set_title("(a) 静压场 x-z (y=0)，绿线=布风板", fontsize=12)
|
||||
ax.set_xlabel("x (m)"); ax.set_ylabel("z (m)")
|
||||
cb = fig.colorbar(im, ax=ax, fraction=0.036, pad=0.02); cb.set_label("Pa", fontsize=9)
|
||||
|
||||
ax = fig.add_subplot(gs[0, 1])
|
||||
ext = [X0, X1, Y0, Y1]
|
||||
Pxy = np.where(SOLID[zi(0.10)], np.nan, P[zi(0.10)])
|
||||
im = ax.imshow(Pxy, origin="lower", extent=ext, cmap="RdBu_r",
|
||||
vmin=-np.nanpercentile(np.abs(Pxy), 99), vmax=np.nanpercentile(np.abs(Pxy), 99),
|
||||
aspect="equal", interpolation="bilinear")
|
||||
ax.set_title("(b) 布风板正上方 z=0.10 静压", fontsize=12)
|
||||
ax.set_xlabel("x (m)"); ax.set_ylabel("y (m)")
|
||||
cb = fig.colorbar(im, ax=ax, fraction=0.036, pad=0.02); cb.set_label("Pa", fontsize=9)
|
||||
|
||||
ax = fig.add_subplot(gs[1, 0])
|
||||
dps = []
|
||||
for k in range(N_SEC):
|
||||
ix0, ix1 = xi(0.6 * k), xi(0.6 * (k + 1))
|
||||
mp = FLUID[zi(-0.30):zi(-0.05), :, ix0:ix1]
|
||||
pp = P[zi(-0.30):zi(-0.05), :, ix0:ix1][mp].mean()
|
||||
mf = FLUID[zi(0.05):zi(0.20), :, ix0:ix1]
|
||||
pf = P[zi(0.05):zi(0.20), :, ix0:ix1][mf].mean()
|
||||
dps.append(pp - pf)
|
||||
bars = ax.bar([f"第{k+1}段" for k in range(N_SEC)], dps, color="#ff7b72", alpha=0.85)
|
||||
ax.axhline(2900, color="#7ee787", ls="--", lw=1.4, label="设计压差 2900 Pa")
|
||||
for b, v in zip(bars, dps):
|
||||
ax.text(b.get_x() + b.get_width() / 2, v + 5, f"{v:.0f}", ha="center", fontsize=9)
|
||||
ax.set_ylabel("布风板压差 (Pa)", fontsize=10); ax.legend(fontsize=9)
|
||||
ax.set_title("(c) 各段布风板压差（模型值）", fontsize=12); ax.grid(alpha=0.2, axis="y")
|
||||
|
||||
ax = fig.add_subplot(gs[1, 1])
|
||||
ext = [X0, X1, Z0, Z1]
|
||||
Pxz2 = np.where(SOLID[:, yi(0), :], np.nan, P[:, yi(0), :])
|
||||
ax.plot(XC, np.nanmean(P[:, yi(0), :], axis=0), color="#39d0d8", lw=1.4, label="沿 x 平均 (y=0)")
|
||||
for zc, col in [(zi(-0.30), "#58a6ff"), (zi(-0.05), "#7ee787"), (zi(0.10), "#ff7b72")]:
|
||||
ax.plot(XC, np.where(SOLID[zc, yi(0), :], np.nan, P[zc, yi(0), :]), lw=1.0,
|
||||
color=col, label=f"z={ZC[zc]:.2f} m")
|
||||
ax.set_xlabel("x (m)", fontsize=10); ax.set_ylabel("静压 (Pa)", fontsize=10)
|
||||
ax.set_title("(d) 静压沿长度剖面", fontsize=12); ax.legend(fontsize=8); ax.grid(alpha=0.2)
|
||||
fig.suptitle("压差分析：布风板上下压差与静压分布", fontsize=15, y=0.99)
|
||||
fig.tight_layout(rect=[0, 0, 1, 0.97])
|
||||
fig.savefig(os.path.join(OUT, "fbed-pressure-drop.png"), dpi=520)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-pressure-drop.png")
|
||||
|
||||
|
||||
def fig_flowrate():
|
||||
fig, axes = plt.subplots(1, 3, figsize=(17.0, 5.4))
|
||||
qs = []
|
||||
for k in range(N_SEC):
|
||||
ix0, ix1 = xi(0.6 * k), xi(0.6 * (k + 1))
|
||||
q = UZ[IZP, :, ix0:ix1][FLUID[IZP, :, ix0:ix1]].sum() * CELL_A * 3600.0
|
||||
qs.append(q)
|
||||
ax = axes[0]
|
||||
ax.bar([f"第{k+1}段" for k in range(N_SEC)], qs, color="#39d0d8", alpha=0.85)
|
||||
ax.axhline(2500, color="#ff7b72", ls="--", lw=1.4, label="目标 2500")
|
||||
ax.set_ylabel("风量 (m³/h)", fontsize=10); ax.legend(fontsize=9)
|
||||
ax.set_title("(a) 各区段布风板风量", fontsize=12); ax.grid(alpha=0.2, axis="y")
|
||||
|
||||
ax = axes[1]
|
||||
cum = np.cumsum([UZ[IZP, :, i][FLUID[IZP, :, i]].sum() for i in range(NX)]) * CELL_A * 3600.0
|
||||
ax.plot(XC, cum, color=ACCENT, lw=1.8)
|
||||
ax.axhline(10000, color="#ff7b72", ls="--", lw=1.2, label="总计 10000")
|
||||
ax.set_xlabel("x (m)", fontsize=10); ax.set_ylabel("累计风量 (m³/h)", fontsize=10)
|
||||
ax.set_title("(b) 沿长度累计风量", fontsize=12); ax.legend(fontsize=9); ax.grid(alpha=0.2)
|
||||
|
||||
ax = axes[2]
|
||||
iz_in = yi(fs.PIPE_Y0)
|
||||
qin = UY[:, iz_in, :][FLUID[:, iz_in, :]].sum() * CELL_A * 3600.0
|
||||
iz_out = NZ - 2
|
||||
qout = UZ[iz_out][FLUID[iz_out]].sum() * CELL_A * 3600.0
|
||||
qplate = sum(qs)
|
||||
labels = ["进风\n(DN300×4)", "布风板\n(通过)", "排风\n(DN800)"]
|
||||
vals = [qin, qplate, qout]
|
||||
cols = ["#7ee787", "#39d0d8", "#ff7b72"]
|
||||
b = ax.bar(labels, vals, color=cols, alpha=0.85)
|
||||
for bb, v in zip(b, vals):
|
||||
ax.text(bb.get_x() + bb.get_width() / 2, v + 60, f"{v:.0f}", ha="center", fontsize=9.5)
|
||||
ax.axhline(10000, color=MUTED, ls=":", lw=1)
|
||||
ax.set_ylabel("风量 (m³/h)", fontsize=10)
|
||||
ax.set_title("(c) 进风 / 布风板 / 排风 风量平衡", fontsize=12); ax.grid(alpha=0.2, axis="y")
|
||||
fig.suptitle("流量分析", fontsize=15, y=1.0)
|
||||
fig.tight_layout(rect=[0, 0, 1, 0.95])
|
||||
fig.savefig(os.path.join(OUT, "fbed-flowrate.png"), dpi=520)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-flowrate.png")
|
||||
|
||||
|
||||
def main():
|
||||
fig_geometry()
|
||||
fig_vectors()
|
||||
fig_tilt()
|
||||
fig_uniformity()
|
||||
fig_decay()
|
||||
fig_zoom_sections()
|
||||
fig_pressure_drop()
|
||||
fig_flowrate()
|
||||
print("[analysis] done")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -0,0 +1,58 @@
|
||||
#!/usr/bin/env python
|
||||
# -*- coding: utf-8 -*-
|
||||
"""空床流场诊断: 检查顶部是否有非物理进流 / 几何开口, 并采样 5 度倾斜效果。"""
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
||||
import fbed_simulation as fs
|
||||
|
||||
D = np.load(os.path.join(fs.OUT, "fbed_fields.npz"))
|
||||
UX, UY, UZ, SOLID = D["ux"], D["uy"], D["uz"], D["solid"]
|
||||
UX, UY, UZ = UX * fs.SCALE_U, UY * fs.SCALE_U, UZ * fs.SCALE_U
|
||||
fluid = ~SOLID
|
||||
NZ, NY, NX = SOLID.shape
|
||||
|
||||
print("=== 出口层 (z 顶部) 速度 ===")
|
||||
for iz in (NZ - 1, NZ - 2):
|
||||
m = fluid[iz]
|
||||
u = UZ[iz][m]
|
||||
print(f" iz={iz} z={fs.ZC[iz]:.3f} n={m.sum()} uz[min,mean,max]="
|
||||
f"[{u.min():.2f}, {u.mean():.2f}, {u.max():.2f}] m/s "
|
||||
f"inflow(<0)占比={100*np.mean(u<0):.1f}%")
|
||||
|
||||
print("\n=== 沸腾段顶部区域 (z>1.35) 是否有非物理入流 ===")
|
||||
band = fluid & (fs.ZC[:, None, None] > 1.35)
|
||||
uzb = UZ[band]
|
||||
print(f" n={uzb.size} uz mean={uzb.mean():.2f} inflow<0 占比={100*np.mean(uzb<0):.2f}%")
|
||||
|
||||
print("\n=== 采样疑似'进风口'位置 (yz 切面左上亮斑) ===")
|
||||
for (x, y, z) in [(0.15, -0.45, 1.2), (0.15, -0.35, 1.1), (2.25, -0.45, 1.2),
|
||||
(0.15, -0.50, 0.9), (1.05, -0.45, 1.2)]:
|
||||
ix, iy, iz = fs._xi(x), fs._yi(y), fs._zi(z)
|
||||
if fluid[iz, iy, ix]:
|
||||
print(f" (x={x},y={y},z={z}): |U|={np.hypot(np.hypot(UX[iz,iy,ix],UY[iz,iy,ix]),UZ[iz,iy,ix]):.2f} "
|
||||
f"ux={UX[iz,iy,ix]:+.2f} uy={UY[iz,iy,ix]:+.2f} uz={UZ[iz,iy,ix]:+.2f}")
|
||||
else:
|
||||
print(f" (x={x},y={y},z={z}): SOLID")
|
||||
|
||||
print("\n=== 5 度倾斜效果: 布风板正上方薄层水平速度 (y 平均) ===")
|
||||
for z in (0.04, 0.06, 0.10, 0.15):
|
||||
iz = fs._zi(z)
|
||||
m = fluid[iz] & (np.abs(fs.YC)[:, None] < fs.BED_W / 2)
|
||||
ux_m = UX[iz][m].mean()
|
||||
uz_m = UZ[iz][m].mean()
|
||||
print(f" z={z:.2f} <ux>={ux_m:+.3f} <uz>={uz_m:+.3f} m/s "
|
||||
f"偏角={np.degrees(np.arctan2(ux_m, uz_m)):+.2f} deg (理论 +5 deg)")
|
||||
|
||||
print("\n=== 布风板各段流量 (按 x 分 4 段) ===")
|
||||
izp = fs._zi(0.0)
|
||||
cell_a = fs.DX ** 2
|
||||
for k in range(4):
|
||||
x0, x1 = 0.6 * k, 0.6 * (k + 1)
|
||||
ix0, ix1 = fs._xi(x0), fs._xi(x1)
|
||||
sub = UZ[izp, :, ix0:ix1] * (fluid[izp, :, ix0:ix1])
|
||||
q = sub.sum() * cell_a * 3600.0
|
||||
print(f" 第{k+1}段 x[{x0:.1f},{x1:.1f}]: Q={q:.0f} m3/h")
|
||||
@@ -0,0 +1,50 @@
|
||||
#!/usr/bin/env python
|
||||
# -*- coding: utf-8 -*-
|
||||
"""沿高度净通量核查 (体积 / 质量)。"""
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
||||
import fbed_simulation as fs
|
||||
|
||||
D = np.load(os.path.join(fs.OUT, "fbed_fields.npz"))
|
||||
UX, UY, UZ, RHO, SOLID = D["ux"], D["uy"], D["uz"], D["rho"], D["solid"]
|
||||
FLUID = ~SOLID
|
||||
UX, UY, UZ = UX * fs.SCALE_U, UY * fs.SCALE_U, UZ * fs.SCALE_U
|
||||
CA = fs.DX ** 2
|
||||
|
||||
print(" z(m) 面积(m2) Q_vol(m3/h) Q_mass(m3/h) 倒流占比")
|
||||
for z in [-0.44, -0.35, -0.28, -0.22, -0.15, -0.08, -0.02, 0.02, 0.06, 0.12,
|
||||
0.25, 0.50, 0.80, 1.00, 1.20, 1.40, 1.60, 1.80, 1.95, 2.05, 2.15]:
|
||||
iz = fs._zi(z)
|
||||
m = FLUID[iz]
|
||||
if m.sum() < 5:
|
||||
continue
|
||||
area = m.sum() * CA
|
||||
qv = UZ[iz][m].sum() * CA * 3600
|
||||
qm = (RHO[iz] * UZ[iz])[m].sum() * CA * 3600
|
||||
frac = 100 * np.mean(UZ[iz][m] < 0)
|
||||
print(f" {z:+.2f} {area:7.3f} {qv:11.0f} {qm:11.0f} {frac:6.1f}%")
|
||||
|
||||
print("\n进/出总通量:")
|
||||
iy_in = fs._yi(fs.PIPE_Y0)
|
||||
m = FLUID[:, iy_in, :]
|
||||
print(f" 进风口(4xDN300) Q_vol={UY[:, iy_in, :][m].sum()*CA*3600:.0f} "
|
||||
f"Q_mass={(RHO*UY)[:, iy_in, :][m].sum()*CA*3600:.0f} m3/h")
|
||||
for iz in (fs.NZ - 1, fs.NZ - 2):
|
||||
m = FLUID[iz]
|
||||
print(f" 排风口 z={fs.ZC[iz]:.3f} Q_vol={UZ[iz][m].sum()*CA*3600:.0f} "
|
||||
f"Q_mass={(RHO[iz]*UZ[iz])[m].sum()*CA*3600:.0f} m3/h")
|
||||
|
||||
print("\n=== 压差方向核查 (布风板上下) ===")
|
||||
P = (RHO - 1.0) / 3.0 * (1.2 * fs.SCALE_U ** 2)
|
||||
for k in range(4):
|
||||
ix0, ix1 = fs._xi(0.6 * k), fs._xi(0.6 * (k + 1))
|
||||
mpl = FLUID[fs._zi(-0.35):fs._zi(-0.05), :, ix0:ix1]
|
||||
p_pl = P[fs._zi(-0.35):fs._zi(-0.05), :, ix0:ix1][mpl].mean()
|
||||
mfb = FLUID[fs._zi(0.05):fs._zi(0.20), :, ix0:ix1]
|
||||
p_fb = P[fs._zi(0.05):fs._zi(0.20), :, ix0:ix1][mfb].mean()
|
||||
print(f" 第{k+1}段: p_风箱={p_pl:7.1f} Pa p_沸腾段={p_fb:7.1f} Pa "
|
||||
f"Δp(风箱-沸腾段)={p_pl-p_fb:+7.1f} Pa {'✓正确' if p_pl>p_fb else '✗反了'}")
|
||||
@@ -0,0 +1,312 @@
|
||||
#!/usr/bin/env python
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
水平卧式流化床 空床流场 —— 批量出图
|
||||
====================================
|
||||
读取 fbed_fields.npz, 生成大量切面 / 分量 / 压力 / 分区图。
|
||||
|
||||
python tools/fbed_figures.py
|
||||
"""
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
import matplotlib
|
||||
matplotlib.use("Agg")
|
||||
import matplotlib.pyplot as plt
|
||||
from matplotlib.colors import LinearSegmentedColormap
|
||||
|
||||
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
||||
import fbed_simulation as fs
|
||||
|
||||
OUT = fs.OUT
|
||||
BG, FG, MUTED, ACCENT = fs.BG, fs.FG, fs.MUTED, fs.ACCENT
|
||||
CMAP_FLOW = fs.CMAP_FLOW
|
||||
XC, YC, ZC = fs.XC, fs.YC, fs.ZC
|
||||
X0, X1, Y0, Y1, Z0, Z1 = fs.X0, fs.X1, fs.Y0, fs.Y1, fs.Z0, fs.Z1
|
||||
PPM, DX, SCALE_U = fs.PPM, fs.DX, fs.SCALE_U
|
||||
|
||||
D = np.load(os.path.join(OUT, "fbed_fields.npz"))
|
||||
UX, UY, UZ, RHO, SOLID = D["ux"], D["uy"], D["uz"], D["rho"], D["solid"]
|
||||
SPEED = np.sqrt(UX ** 2 + UY ** 2 + UZ ** 2) * SCALE_U
|
||||
P = (RHO - 1.0) / 3.0 * (1.2 * SCALE_U ** 2)
|
||||
P[SOLID] = np.nan
|
||||
|
||||
CMAP_P = plt.get_cmap("RdBu_r")
|
||||
|
||||
|
||||
def xi(v):
|
||||
return fs._xi(v)
|
||||
|
||||
|
||||
def yi(v):
|
||||
return fs._yi(v)
|
||||
|
||||
|
||||
def zi(v):
|
||||
return fs._zi(v)
|
||||
|
||||
|
||||
def get(which, idx):
|
||||
if which == "xz":
|
||||
return (SPEED[:, idx, :], UX[:, idx, :], UZ[:, idx, :], SOLID[:, idx, :],
|
||||
[X0, X1, Z0, Z1], "x (m)", "z (m)")
|
||||
if which == "yz":
|
||||
return (SPEED[:, :, idx], UY[:, :, idx], UZ[:, :, idx], SOLID[:, :, idx],
|
||||
[Y0, Y1, Z0, Z1], "y (m)", "z (m)")
|
||||
return (SPEED[idx, :, :], UX[idx, :, :], UY[idx, :, :], SOLID[idx, :, :],
|
||||
[X0, X1, Y0, Y1], "x (m)", "y (m)")
|
||||
|
||||
|
||||
def get_scalar(field, which, idx):
|
||||
if which == "xz":
|
||||
return field[:, idx, :], SOLID[:, idx, :], [X0, X1, Z0, Z1], "x (m)", "z (m)"
|
||||
if which == "yz":
|
||||
return field[:, :, idx], SOLID[:, :, idx], [Y0, Y1, Z0, Z1], "y (m)", "z (m)"
|
||||
return field[idx, :, :], SOLID[idx, :, :], [X0, X1, Y0, Y1], "x (m)", "y (m)"
|
||||
|
||||
|
||||
def panel(ax, which, idx, vmax, stream=True, density=1.1):
|
||||
s2, u2, v2, sol2, ext, xl, yl = get(which, idx)
|
||||
ax.imshow(np.where(sol2, np.nan, s2), origin="lower", extent=ext, cmap=CMAP_FLOW,
|
||||
vmin=0, vmax=vmax, interpolation="bilinear", aspect="equal")
|
||||
if stream:
|
||||
ny, nx = sol2.shape
|
||||
xs = np.linspace(ext[0], ext[1], nx)
|
||||
ys = np.linspace(ext[2], ext[3], ny)
|
||||
U = np.where(sol2, np.nan, u2)
|
||||
V = np.where(sol2, np.nan, v2)
|
||||
with np.errstate(invalid="ignore"):
|
||||
ax.streamplot(xs, ys, U, V, density=density,
|
||||
color=(0.85, 0.88, 0.92, 0.38), linewidth=0.45,
|
||||
arrowsize=0.5, zorder=5)
|
||||
ax.contour(np.linspace(ext[0], ext[1], sol2.shape[1]),
|
||||
np.linspace(ext[2], ext[3], sol2.shape[0]),
|
||||
sol2.astype(float), levels=[0.5], colors="#adbac7", linewidths=0.7, zorder=6)
|
||||
ax.set_xlabel(xl, fontsize=8)
|
||||
ax.set_ylabel(yl, fontsize=8)
|
||||
ax.set_xlim(ext[0], ext[1])
|
||||
ax.set_ylim(ext[2], ext[3])
|
||||
ax.tick_params(labelsize=7)
|
||||
|
||||
|
||||
def panel_scalar(ax, field, which, idx, cmap, vmin=None, vmax=None, center=False, stream=True):
|
||||
f2, sol2, ext, xl, yl = get_scalar(field, which, idx)
|
||||
f2 = np.where(sol2, np.nan, f2)
|
||||
if center:
|
||||
lim = np.nanpercentile(np.abs(f2), 99) * 1.05
|
||||
im = ax.imshow(f2, origin="lower", extent=ext, cmap=cmap, vmin=-lim, vmax=lim,
|
||||
interpolation="bilinear", aspect="equal")
|
||||
else:
|
||||
im = ax.imshow(f2, origin="lower", extent=ext, cmap=cmap, vmin=vmin, vmax=vmax,
|
||||
interpolation="bilinear", aspect="equal")
|
||||
if stream:
|
||||
s2, u2, v2, _, _, _, _ = get(which, idx)
|
||||
ny, nx = sol2.shape
|
||||
xs = np.linspace(ext[0], ext[1], nx)
|
||||
ys = np.linspace(ext[2], ext[3], ny)
|
||||
with np.errstate(invalid="ignore"):
|
||||
ax.streamplot(xs, ys, np.where(sol2, np.nan, u2), np.where(sol2, np.nan, v2),
|
||||
density=1.0, color=(0.85, 0.88, 0.92, 0.3), linewidth=0.4,
|
||||
arrowsize=0.5, zorder=5)
|
||||
ax.contour(np.linspace(ext[0], ext[1], sol2.shape[1]),
|
||||
np.linspace(ext[2], ext[3], sol2.shape[0]),
|
||||
sol2.astype(float), levels=[0.5], colors="#adbac7", linewidths=0.7, zorder=6)
|
||||
ax.set_xlabel(xl, fontsize=8)
|
||||
ax.set_ylabel(yl, fontsize=8)
|
||||
ax.set_xlim(ext[0], ext[1])
|
||||
ax.set_ylim(ext[2], ext[3])
|
||||
ax.tick_params(labelsize=7)
|
||||
return im
|
||||
|
||||
|
||||
def montage(specs, ncols, suptitle, name, vmax, cbar_label="|U| (m/s)",
|
||||
figw=17.0, panel_h=3.0, stream=True, density=1.1):
|
||||
nrows = (len(specs) + ncols - 1) // ncols
|
||||
fig, axes = plt.subplots(nrows, ncols, figsize=(figw, panel_h * nrows + 1.0),
|
||||
squeeze=False)
|
||||
for ax in axes.ravel()[len(specs):]:
|
||||
ax.axis("off")
|
||||
for ax, (spec, title) in zip(axes.ravel(), specs):
|
||||
which, idx = spec
|
||||
panel(ax, which, idx, vmax, stream=stream, density=density)
|
||||
ax.set_title(title, fontsize=9.5, pad=4)
|
||||
sm = plt.cm.ScalarMappable(cmap=CMAP_FLOW, norm=plt.Normalize(0, vmax))
|
||||
cb = fig.colorbar(sm, ax=axes.ravel().tolist(), fraction=0.012, pad=0.01)
|
||||
cb.set_label(cbar_label, fontsize=10)
|
||||
cb.outline.set_edgecolor("#30363d")
|
||||
fig.suptitle(suptitle, fontsize=15, y=0.995)
|
||||
fig.savefig(os.path.join(OUT, name), dpi=500, bbox_inches="tight")
|
||||
plt.close(fig)
|
||||
print(" ->", name)
|
||||
|
||||
|
||||
def fig_xy_plenum():
|
||||
specs = [(("xy", zi(z)), f"z = {z:+.2f} m") for z in (-0.42, -0.32, -0.22, -0.12, -0.02)]
|
||||
montage(specs, 3, "风箱层水平切面 |U|（4 段侧进风分布）", "fbed-montage-plenum.png", 6.0)
|
||||
|
||||
|
||||
def fig_xy_plate():
|
||||
specs = [(("xy", zi(z)), f"z = {z:.2f} m") for z in (0.04, 0.08, 0.13, 0.20, 0.30, 0.45)]
|
||||
montage(specs, 3, "布风板近区水平切面 |U|（射流合并过程）", "fbed-montage-nearplate.png", 3.0)
|
||||
|
||||
|
||||
def fig_xy_free():
|
||||
specs = [(("xy", zi(z)), f"z = {z:.2f} m") for z in (0.60, 0.80, 1.00, 1.20, 1.45, 1.75)]
|
||||
montage(specs, 3, "沸腾段 / 渐扩 / 排风 水平切面 |U|", "fbed-montage-freeboard.png", 4.0)
|
||||
|
||||
|
||||
def fig_xz():
|
||||
specs = [(("xz", yi(y)), f"y = {y:+.2f} m") for y in (-0.26, -0.13, 0.0, 0.13, 0.26)]
|
||||
montage(specs, 3, "垂直切面 x-z（不同宽度位置）|U|", "fbed-montage-xz.png", 5.0, figw=16.0)
|
||||
|
||||
|
||||
def fig_yz():
|
||||
specs = [(("yz", xi(x)), f"x = {x:.2f} m") for x in (0.15, 0.45, 0.75, 1.05, 1.35, 1.65, 1.95, 2.25)]
|
||||
montage(specs, 4, "垂直切面 y-z（沿长度）|U|", "fbed-montage-yz.png", 5.0, figw=18.0, panel_h=3.4)
|
||||
|
||||
|
||||
def fig_components():
|
||||
fig, axes = plt.subplots(2, 2, figsize=(15.0, 12.0))
|
||||
panel(axes[0, 0], "xz", yi(0.0), 4.0)
|
||||
axes[0, 0].set_title("(a) 速度大小 |U|", fontsize=12, color=ACCENT)
|
||||
for ax, comp, title, cmap in [
|
||||
(axes[0, 1], UX, "(b) 水平分量 ux（孔板 5° 倾斜 + 向排风汇流）", "coolwarm"),
|
||||
(axes[1, 0], UY, "(c) 横向分量 uy（侧进风 + 渐扩）", "coolwarm"),
|
||||
(axes[1, 1], UZ, "(d) 垂直分量 uz（上升主流）", "coolwarm")]:
|
||||
f2, sol2, ext, xl, yl = get_scalar(comp * SCALE_U, "xz", yi(0.0))
|
||||
lim = np.nanpercentile(np.abs(np.where(sol2, np.nan, f2)), 99) * 1.05
|
||||
im = ax.imshow(np.where(sol2, np.nan, f2), origin="lower", extent=ext, cmap=cmap,
|
||||
vmin=-lim, vmax=lim, interpolation="bilinear", aspect="equal")
|
||||
ax.contour(np.linspace(ext[0], ext[1], sol2.shape[1]),
|
||||
np.linspace(ext[2], ext[3], sol2.shape[0]),
|
||||
sol2.astype(float), levels=[0.5], colors="#adbac7", linewidths=0.7)
|
||||
ax.set_title(title, fontsize=12)
|
||||
ax.set_xlabel("x (m)", fontsize=9)
|
||||
ax.set_ylabel("z (m)", fontsize=9)
|
||||
cb = fig.colorbar(im, ax=ax, fraction=0.036, pad=0.02)
|
||||
cb.set_label("m/s", fontsize=9)
|
||||
cb.outline.set_edgecolor("#30363d")
|
||||
fig.suptitle("速度分量分解（垂直切面 y=0）", fontsize=15, y=0.99)
|
||||
fig.tight_layout(rect=[0, 0, 1, 0.97])
|
||||
fig.savefig(os.path.join(OUT, "fbed-components.png"), dpi=520)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-components.png")
|
||||
|
||||
|
||||
def fig_pressure():
|
||||
fig, axes = plt.subplots(2, 2, figsize=(15.0, 12.0))
|
||||
for ax, (which, idx, title) in zip(axes.ravel(), [
|
||||
("xz", yi(0.0), "(a) 垂直切面 y=0 —— 全系统压差"),
|
||||
("yz", xi(1.2), "(b) 垂直切面 x=1.2"),
|
||||
("xy", zi(0.15), "(c) 布风板近区 z=0.15"),
|
||||
("xy", zi(0.90), "(d) 沸腾段 z=0.90")]):
|
||||
im = panel_scalar(ax, P, which, idx, CMAP_P, center=True, stream=False)
|
||||
ax.set_title(title, fontsize=12)
|
||||
cb = fig.colorbar(im, ax=ax, fraction=0.036, pad=0.02)
|
||||
cb.set_label("静压 (Pa)", fontsize=9)
|
||||
cb.outline.set_edgecolor("#30363d")
|
||||
fig.suptitle("静压场（相对出口, 单位 Pa）", fontsize=15, y=0.99)
|
||||
fig.tight_layout(rect=[0, 0, 1, 0.97])
|
||||
fig.savefig(os.path.join(OUT, "fbed-pressure.png"), dpi=520)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-pressure.png")
|
||||
|
||||
|
||||
def fig_sections():
|
||||
fig, axes = plt.subplots(1, 4, figsize=(17.0, 6.0))
|
||||
for k, ax in enumerate(axes):
|
||||
xc = 0.6 * (k + 0.5)
|
||||
panel(ax, "yz", xi(xc), 5.0)
|
||||
ax.set_title(f"第 {k+1} 段 x={xc:.1f} m", fontsize=11)
|
||||
fig.suptitle("4 个进风区段中心垂直切面 y-z（单段风量 2500 m³/h）", fontsize=14, y=1.0)
|
||||
fig.tight_layout(rect=[0, 0, 1, 0.95])
|
||||
fig.savefig(os.path.join(OUT, "fbed-sections.png"), dpi=520)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-sections.png")
|
||||
|
||||
|
||||
def fig_hero_stream():
|
||||
fig, axes = plt.subplots(1, 2, figsize=(18.0, 8.0))
|
||||
panel(axes[0], "xz", yi(0.0), 4.0, density=1.8)
|
||||
axes[0].set_title("垂直切面 y=0 —— 风箱/布风板/沸腾段/排风", fontsize=12)
|
||||
panel(axes[1], "yz", xi(1.2), 4.0, density=1.8)
|
||||
axes[1].set_title("垂直切面 x=1.2 —— 侧向渐扩", fontsize=12)
|
||||
fig.suptitle("空床流场流线（GPU D3Q19 LBM）", fontsize=15, y=0.99)
|
||||
fig.tight_layout(rect=[0, 0, 1, 0.96])
|
||||
fig.savefig(os.path.join(OUT, "fbed-streamlines.png"), dpi=560)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-streamlines.png")
|
||||
|
||||
|
||||
def fig_3d():
|
||||
from mpl_toolkits.mplot3d import Axes3D # noqa
|
||||
fig = plt.figure(figsize=(13.0, 10.0))
|
||||
ax = fig.add_subplot(111, projection="3d")
|
||||
norm = plt.Normalize(0, 4.0)
|
||||
|
||||
def surf_y(iy):
|
||||
f = SPEED[:, iy, :]
|
||||
m = SOLID[:, iy, :]
|
||||
X, Z = np.meshgrid(XC, ZC)
|
||||
Y = np.full_like(X, YC[iy])
|
||||
fc = CMAP_FLOW(norm(np.where(m, np.nan, f)))
|
||||
fc[m, 3] = 0.0
|
||||
ax.plot_surface(X, Y, Z, facecolors=fc, shade=False, rstride=1, cstride=1,
|
||||
linewidth=0, antialiased=False)
|
||||
|
||||
def surf_x(ix):
|
||||
f = SPEED[:, :, ix]
|
||||
m = SOLID[:, :, ix]
|
||||
Y, Z = np.meshgrid(YC, ZC)
|
||||
X = np.full_like(Y, XC[ix])
|
||||
fc = CMAP_FLOW(norm(np.where(m, np.nan, f)))
|
||||
fc[m, 3] = 0.0
|
||||
ax.plot_surface(X, Y, Z, facecolors=fc, shade=False, rstride=1, cstride=1,
|
||||
linewidth=0, antialiased=False)
|
||||
|
||||
def surf_z(iz):
|
||||
f = SPEED[iz, :, :]
|
||||
m = SOLID[iz, :, :]
|
||||
X, Y = np.meshgrid(XC, YC)
|
||||
Z = np.full_like(X, ZC[iz])
|
||||
fc = CMAP_FLOW(norm(np.where(m, np.nan, f)))
|
||||
fc[m, 3] = 0.0
|
||||
ax.plot_surface(X, Y, Z, facecolors=fc, shade=False, rstride=1, cstride=1,
|
||||
linewidth=0, antialiased=False)
|
||||
|
||||
surf_y(yi(0.0))
|
||||
surf_x(xi(1.2))
|
||||
surf_z(zi(0.9))
|
||||
ax.set_xlabel("x (m)")
|
||||
ax.set_ylabel("y (m)")
|
||||
ax.set_zlabel("z (m)")
|
||||
ax.set_title("三正交切面三维剖切（|U|）", fontsize=13)
|
||||
ax.view_init(elev=22, azim=-58)
|
||||
sm = plt.cm.ScalarMappable(cmap=CMAP_FLOW, norm=norm)
|
||||
cb = fig.colorbar(sm, ax=ax, fraction=0.03, pad=0.08, shrink=0.6)
|
||||
cb.set_label("|U| (m/s)", fontsize=10)
|
||||
cb.outline.set_edgecolor("#30363d")
|
||||
fig.savefig(os.path.join(OUT, "fbed-3d-cutaway.png"), dpi=520)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-3d-cutaway.png")
|
||||
|
||||
|
||||
def main():
|
||||
print("[figs] loading", os.path.join(OUT, "fbed_fields.npz"))
|
||||
fig_xy_plenum()
|
||||
fig_xy_plate()
|
||||
fig_xy_free()
|
||||
fig_xz()
|
||||
fig_yz()
|
||||
fig_components()
|
||||
fig_pressure()
|
||||
fig_sections()
|
||||
fig_hero_stream()
|
||||
fig_3d()
|
||||
print("[figs] done")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -0,0 +1,54 @@
|
||||
#!/usr/bin/env python
|
||||
# -*- coding: utf-8 -*-
|
||||
"""探针: 检查 x=0.3 的 y-z 切面是否有异常开口 / 高速射流。"""
|
||||
import os
|
||||
import sys
|
||||
|
||||
import numpy as np
|
||||
|
||||
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
||||
import fbed_simulation as fs
|
||||
|
||||
D = np.load(os.path.join(fs.OUT, "fbed_fields.npz"))
|
||||
UX, UY, UZ, SOLID = D["ux"], D["uy"], D["uz"], D["solid"]
|
||||
UX, UY, UZ = UX * fs.SCALE_U, UY * fs.SCALE_U, UZ * fs.SCALE_U
|
||||
SP = np.sqrt(UX ** 2 + UY ** 2 + UZ ** 2)
|
||||
NZ, NY, NX = SOLID.shape
|
||||
|
||||
ix = fs._xi(0.3)
|
||||
print(f"=== x=0.3 (ix={ix}) y-z 切面: 固体/流体 ASCII 图 (z 由高到低) ===")
|
||||
print("每格: '#'=固体 '.'=流体(<0.5m/s) 'o'=0.5-1.5 'O'=1.5-3 '*'=3-6 '@'=>6 m/s")
|
||||
for iz in range(NZ - 1, -1, -1):
|
||||
z = fs.ZC[iz]
|
||||
if z < -0.5 or z > 1.7:
|
||||
continue
|
||||
row = ""
|
||||
for iy in range(NY):
|
||||
if SOLID[iz, iy, ix]:
|
||||
row += "#"
|
||||
else:
|
||||
s = SP[iz, iy, ix]
|
||||
row += "." if s < 0.5 else ("o" if s < 1.5 else ("O" if s < 3 else ("*" if s < 6 else "@")))
|
||||
print(f"z={z:+.2f} {row}")
|
||||
|
||||
print("\n=== y=-0.6 墙完整性 (沿 z 扫) ===")
|
||||
iy_wall = fs._yi(-0.6)
|
||||
runs = []
|
||||
for iz in range(NZ):
|
||||
if SOLID[iz, iy_wall, ix]:
|
||||
runs.append((fs.ZC[iz], "solid"))
|
||||
else:
|
||||
runs.append((fs.ZC[iz], "FLUID"))
|
||||
open_z = [z for z, s in runs if s == "FLUID"]
|
||||
print(f" y=-0.6 处流体 z 段: {min(open_z):.2f} .. {max(open_z):.2f}" if open_z else " 全固体")
|
||||
print(f" 流体 z 计数={len(open_z)} / {NZ}")
|
||||
|
||||
print("\n=== 顶部区域 (z>1.0) 是否存在指向域内的射流 ===")
|
||||
for (x, y, z) in [(0.3, -0.5, 1.2), (0.3, -0.45, 1.3), (0.3, -0.35, 1.1),
|
||||
(2.1, -0.5, 1.2), (0.9, -0.5, 1.2), (1.5, -0.5, 1.2)]:
|
||||
i = fs._xi(x); j = fs._yi(y); k = fs._zi(z)
|
||||
if SOLID[k, j, i]:
|
||||
print(f" ({x},{y},{z}): SOLID")
|
||||
else:
|
||||
print(f" ({x},{y},{z}): ux={UX[k,j,i]:+6.2f} uy={UY[k,j,i]:+6.2f} uz={UZ[k,j,i]:+6.2f} "
|
||||
f"|U|={SP[k,j,i]:5.2f}")
|
||||
@@ -0,0 +1,401 @@
|
||||
#!/usr/bin/env python
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
水平卧式流化床 空床流场 —— RANS + 多孔跳变布风板 (交错网格投影法 + AMG)
|
||||
======================================================================
|
||||
交错网格 (MAC): p 在单元中心, u/v/w 在面心 -> 从原理上消除棋盘压力振荡。
|
||||
|
||||
预测: u* = u + dt[ -(u.grad)u + div(nu_eff grad u) + S_porous ]
|
||||
压力: div(grad p) = (rho/dt) div(u*) (pyamg 多重网格, 固体 Neumann, 出口 Dirichlet)
|
||||
修正: u = u* - (dt/rho) grad(p)
|
||||
|
||||
布风板 = 薄多孔层 (Darcy), 阻力标定 dP=2900 Pa @ 表观气速 2.04 m/s。
|
||||
进风 = 4 x DN300 速度进口 (2500 m3/h/段); 排风 = 压力出口; 其余无滑移壁。
|
||||
|
||||
运行:
|
||||
python tools/fbed_rans.py
|
||||
"""
|
||||
import os
|
||||
import time
|
||||
|
||||
import numpy as np
|
||||
import numba
|
||||
from scipy.sparse import csr_matrix
|
||||
import pyamg
|
||||
|
||||
# --------------------------------------------------------------------------
|
||||
BED_L, BED_W = 2.4, 0.567
|
||||
N_SEC = 4
|
||||
SEC_L = BED_L / N_SEC
|
||||
PLENUM_H = 0.45
|
||||
PLATE_T = 0.05
|
||||
FB_H, FB_W = 1.0, 1.2
|
||||
UPPER_H = 0.4
|
||||
REDUCER_H = 0.5
|
||||
DUCT_H = 0.3
|
||||
EXH_R = 0.4
|
||||
PIPE_R = 0.15
|
||||
PIPE_Y0 = -0.65
|
||||
|
||||
X0, X1 = 0.0, BED_L
|
||||
Y0, Y1 = PIPE_Y0, 0.65
|
||||
Z0 = -PLENUM_H
|
||||
Z1 = FB_H + UPPER_H + REDUCER_H + DUCT_H
|
||||
|
||||
PPM = int(os.environ.get("RANS_PPM", "20"))
|
||||
NX = int(round((X1 - X0) * PPM))
|
||||
NY = int(round((Y1 - Y0) * PPM))
|
||||
NZ = int(round((Z1 - Z0) * PPM))
|
||||
DX = (X1 - X0) / NX
|
||||
DY = (Y1 - Y0) / NY
|
||||
DZ = (Z1 - Z0) / NZ
|
||||
XC = X0 + (np.arange(NX) + 0.5) * DX
|
||||
YC = Y0 + (np.arange(NY) + 0.5) * DY
|
||||
ZC = Z0 + (np.arange(NZ) + 0.5) * DZ
|
||||
|
||||
Q_SEC = 2500.0 / 3600.0
|
||||
A_PIPE = np.pi * PIPE_R ** 2
|
||||
U_IN = Q_SEC / A_PIPE
|
||||
A_BED = BED_L * BED_W
|
||||
U_SUP = Q_SEC * N_SEC / A_BED
|
||||
DP_RATED = float(os.environ.get("RANS_DP", "2900.0")) # 额定点压力 (2500 m3/h/段)
|
||||
Q_RATED = Q_SEC * N_SEC # 额定总流量 m3/s
|
||||
P_FAN = DP_RATED * Q_RATED # 恒功率 (W)
|
||||
DP_SHUT = DP_RATED * float(os.environ.get("RANS_FAN_SHUT", "1.5"))
|
||||
|
||||
|
||||
def fan_dp(q_m3s):
|
||||
"""恒功率风机曲线: dp * Q = P -> dp = P/Q (Q->0 时用关死压力截断)"""
|
||||
if q_m3s < 1e-6:
|
||||
return DP_SHUT
|
||||
return min(P_FAN / q_m3s, DP_SHUT)
|
||||
RHO = 1.2
|
||||
NU = 1.5e-5
|
||||
R_POROUS = DP_RATED / (U_SUP * PLATE_T)
|
||||
NU_T = float(os.environ.get("RANS_NUT", "0.02"))
|
||||
NU_EFF = NU + NU_T
|
||||
|
||||
MAX_STEPS = int(os.environ.get("RANS_STEPS", "4000"))
|
||||
TIME_BUDGET = float(os.environ.get("RANS_TB", "16.0"))
|
||||
CFL = float(os.environ.get("RANS_CFL", "0.3"))
|
||||
|
||||
|
||||
def build_masks():
|
||||
solid = np.ones((NZ, NY, NX), dtype=bool)
|
||||
yw = BED_W / 2.0
|
||||
|
||||
def box(xa, xb, ya, yb, za, zb, val=False):
|
||||
ix0, ix1 = int(round((xa - X0) / DX)), int(round((xb - X0) / DX))
|
||||
iy0, iy1 = int(round((ya - Y0) / DY)), int(round((yb - Y0) / DY))
|
||||
iz0, iz1 = int(round((za - Z0) / DZ)), int(round((zb - Z0) / DZ))
|
||||
solid[max(iz0, 0):min(iz1, NZ), max(iy0, 0):min(iy1, NY),
|
||||
max(ix0, 0):min(ix1, NX)] = val
|
||||
|
||||
box(0.0, BED_L, -yw, yw, -PLENUM_H, 0.0)
|
||||
for k in range(1, N_SEC):
|
||||
xk = k * SEC_L
|
||||
box(xk - 0.0125, xk + 0.0125, -yw, yw, -PLENUM_H, 0.0, True)
|
||||
box(0.0, BED_L, -yw, yw, 0.0, PLATE_T)
|
||||
for iz in range(NZ):
|
||||
z = ZC[iz]
|
||||
if z <= PLATE_T or z >= FB_H:
|
||||
continue
|
||||
w = yw + (FB_W / 2.0 - yw) * (z / FB_H)
|
||||
solid[iz] = np.where(np.abs(YC[:, None]) <= w, False, solid[iz])
|
||||
box(0.0, BED_L, -FB_W / 2.0, FB_W / 2.0, FB_H, FB_H + UPPER_H)
|
||||
z_r0, z_r1 = FB_H + UPPER_H, FB_H + UPPER_H + REDUCER_H
|
||||
for iz in range(NZ):
|
||||
z = ZC[iz]
|
||||
if z < z_r0 or z >= z_r1:
|
||||
continue
|
||||
t = (z - z_r0) / REDUCER_H
|
||||
xh = BED_L / 2.0 - (BED_L / 2.0 - 0.4) * t
|
||||
yh = FB_W / 2.0 - (FB_W / 2.0 - 0.4) * t
|
||||
solid[iz] = np.where((np.abs(XC[None, :] - BED_L / 2.0) <= xh) &
|
||||
(np.abs(YC[:, None]) <= yh), False, solid[iz])
|
||||
cyl = (XC[None, :] - BED_L / 2.0) ** 2 + (YC[:, None]) ** 2 <= EXH_R ** 2
|
||||
for iz in range(NZ):
|
||||
if ZC[iz] >= z_r1:
|
||||
solid[iz] = np.where(cyl, False, solid[iz])
|
||||
|
||||
iy_p0 = int(round((PIPE_Y0 - Y0) / DY))
|
||||
iy_p1 = int(round((-yw - Y0) / DY))
|
||||
for k in range(N_SEC):
|
||||
xk = SEC_L * (k + 0.5)
|
||||
zk = -PLENUM_H / 2.0
|
||||
for iz in range(NZ):
|
||||
if abs(ZC[iz] - zk) > PIPE_R:
|
||||
continue
|
||||
m = (XC - xk) ** 2 + (ZC[iz] - zk) ** 2 <= PIPE_R ** 2
|
||||
solid[iz, iy_p0:iy_p1, m] = False
|
||||
|
||||
plate = np.zeros_like(solid)
|
||||
iz0, iz1 = int(round((0.0 - Z0) / DZ)), int(round((PLATE_T - Z0) / DZ))
|
||||
iy0, iy1 = int(round((-yw - Y0) / DY)), int(round((yw - Y0) / DY))
|
||||
plate[iz0:iz1, iy0:iy1, :] = ~solid[iz0:iz1, iy0:iy1, :]
|
||||
inlet = np.zeros_like(solid)
|
||||
inlet[:, iy_p0:iy_p0 + 2, :] = ~solid[:, iy_p0:iy_p0 + 2, :]
|
||||
outlet = np.zeros_like(solid)
|
||||
outlet[NZ - 2:, :, :] = ~solid[NZ - 2:, :, :]
|
||||
return solid, plate, inlet, outlet
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------
|
||||
# 预测步 (交错网格)
|
||||
# --------------------------------------------------------------------------
|
||||
@numba.njit(cache=True, fastmath=True)
|
||||
def predict(u, v, w, solid, plate, inlet, outlet, nu, dt, dx, dy, dz, rporous, uin):
|
||||
nz, ny, nx = solid.shape
|
||||
# u faces: (nz, ny, nx+1), u[iz,iy,ix] at west face of cell ix
|
||||
for iz in range(nz):
|
||||
for iy in range(ny):
|
||||
for ix in range(1, nx):
|
||||
if solid[iz, iy, ix] and solid[iz, iy, ix - 1]:
|
||||
u[iz, iy, ix] = 0.0
|
||||
continue
|
||||
if solid[iz, iy, ix] or solid[iz, iy, ix - 1]:
|
||||
u[iz, iy, ix] = 0.0
|
||||
continue
|
||||
uc = u[iz, iy, ix]
|
||||
ue = u[iz, iy, ix + 1] if ix + 1 <= nx else uc
|
||||
uw = u[iz, iy, ix - 1]
|
||||
un = u[iz, iy + 1, ix] if iy + 1 < ny else uc
|
||||
us = u[iz, iy - 1, ix] if iy - 1 >= 0 else uc
|
||||
ut = u[iz + 1, iy, ix] if iz + 1 < nz else uc
|
||||
ub = u[iz - 1, iy, ix] if iz - 1 >= 0 else uc
|
||||
vn = 0.5 * (v[iz, iy + 1, ix - 1] + v[iz, iy + 1, ix]) if iy + 1 < ny else 0.0
|
||||
vs = 0.5 * (v[iz, iy, ix - 1] + v[iz, iy, ix])
|
||||
wt = 0.5 * (w[iz + 1, iy, ix - 1] + w[iz + 1, iy, ix]) if iz + 1 < nz else 0.0
|
||||
wb = 0.5 * (w[iz, iy, ix - 1] + w[iz, iy, ix])
|
||||
dux = (uc - uw) / dx if uc > 0 else (ue - uc) / dx
|
||||
duy = (un - uc) / dy if vn > 0 else (uc - us) / dy
|
||||
duz = (ut - uc) / dz if wt > 0 else (uc - ub) / dz
|
||||
conv = uc * dux + 0.5 * (vn + vs) * duy + 0.5 * (wt + wb) * duz
|
||||
lap = (ue - 2.0 * uc + uw) / (dx * dx) + (un - 2.0 * uc + us) / (dy * dy) \
|
||||
+ (ut - 2.0 * uc + ub) / (dz * dz)
|
||||
unew = uc + dt * (-conv + nu * lap)
|
||||
if plate[iz, iy, ix] or plate[iz, iy, ix - 1]:
|
||||
unew /= (1.0 + dt * rporous)
|
||||
u[iz, iy, ix] = unew
|
||||
# v faces: (nz, ny+1, nx)
|
||||
for iz in range(nz):
|
||||
for iy in range(1, ny):
|
||||
for ix in range(nx):
|
||||
if solid[iz, iy, ix] or solid[iz, iy - 1, ix]:
|
||||
v[iz, iy, ix] = 0.0
|
||||
continue
|
||||
vc = v[iz, iy, ix]
|
||||
ve = v[iz, iy, ix + 1] if ix + 1 < nx else vc
|
||||
vw = v[iz, iy, ix - 1] if ix - 1 >= 0 else vc
|
||||
vn = v[iz, iy + 1, ix] if iy + 1 <= ny else vc
|
||||
vs = v[iz, iy - 1, ix]
|
||||
vt = v[iz + 1, iy, ix] if iz + 1 < nz else vc
|
||||
vb = v[iz - 1, iy, ix] if iz - 1 >= 0 else vc
|
||||
ue = 0.5 * (u[iz, iy, ix + 1] + u[iz, iy - 1, ix + 1]) if ix + 1 <= nx else 0.0
|
||||
uw = 0.5 * (u[iz, iy, ix] + u[iz, iy - 1, ix])
|
||||
wt = 0.5 * (w[iz + 1, iy, ix] + w[iz + 1, iy - 1, ix]) if iz + 1 < nz else 0.0
|
||||
wb = 0.5 * (w[iz, iy, ix] + w[iz, iy - 1, ix])
|
||||
dvy = (vc - vs) / dy if vc > 0 else (vn - vc) / dy
|
||||
dvx = (ve - vc) / dx if ue > 0 else (vc - vw) / dx
|
||||
dvz = (vt - vc) / dz if wt > 0 else (vc - vb) / dz
|
||||
conv = 0.5 * (ue + uw) * dvx + vc * dvy + 0.5 * (wt + wb) * dvz
|
||||
lap = (ve - 2.0 * vc + vw) / (dx * dx) + (vn - 2.0 * vc + vs) / (dy * dy) \
|
||||
+ (vt - 2.0 * vc + vb) / (dz * dz)
|
||||
vnew = vc + dt * (-conv + nu * lap)
|
||||
if plate[iz, iy, ix] or plate[iz, iy - 1, ix]:
|
||||
vnew /= (1.0 + dt * rporous)
|
||||
v[iz, iy, ix] = vnew
|
||||
# w faces: (nz+1, ny, nx)
|
||||
for iz in range(1, nz):
|
||||
for iy in range(ny):
|
||||
for ix in range(nx):
|
||||
if solid[iz, iy, ix] or solid[iz - 1, iy, ix]:
|
||||
w[iz, iy, ix] = 0.0
|
||||
continue
|
||||
wc = w[iz, iy, ix]
|
||||
we = w[iz, iy, ix + 1] if ix + 1 < nx else wc
|
||||
ww = w[iz, iy, ix - 1] if ix - 1 >= 0 else wc
|
||||
wn = w[iz, iy + 1, ix] if iy + 1 < ny else wc
|
||||
ws = w[iz, iy - 1, ix] if iy - 1 >= 0 else wc
|
||||
wt = w[iz + 1, iy, ix] if iz + 1 <= nz else wc
|
||||
wb = w[iz - 1, iy, ix]
|
||||
ue = 0.5 * (u[iz, iy, ix + 1] + u[iz - 1, iy, ix + 1]) if ix + 1 <= nx else 0.0
|
||||
uw = 0.5 * (u[iz, iy, ix] + u[iz - 1, iy, ix])
|
||||
vn = 0.5 * (v[iz, iy + 1, ix] + v[iz - 1, iy + 1, ix]) if iy + 1 <= ny else 0.0
|
||||
vs = 0.5 * (v[iz, iy, ix] + v[iz - 1, iy, ix])
|
||||
dwz = (wc - wb) / dz if wc > 0 else (wt - wc) / dz
|
||||
dwx = (we - wc) / dx if ue > 0 else (wc - ww) / dx
|
||||
dwy = (wn - wc) / dy if vn > 0 else (wc - ws) / dy
|
||||
conv = 0.5 * (ue + uw) * dwx + 0.5 * (vn + vs) * dwy + wc * dwz
|
||||
lap = (we - 2.0 * wc + ww) / (dx * dx) + (wn - 2.0 * wc + ws) / (dy * dy) \
|
||||
+ (wt - 2.0 * wc + wb) / (dz * dz)
|
||||
wnew = wc + dt * (-conv + nu * lap)
|
||||
if plate[iz, iy, ix] or plate[iz - 1, iy, ix]:
|
||||
wnew /= (1.0 + dt * rporous)
|
||||
w[iz, iy, ix] = wnew
|
||||
return
|
||||
|
||||
|
||||
@numba.njit(cache=True, fastmath=True)
|
||||
def apply_bc(u, v, w, solid, inlet, outlet, uin):
|
||||
nz, ny, nx = solid.shape
|
||||
# 压力进口: 边界面速度零梯度 (让气流自由进入)
|
||||
for iz in range(nz):
|
||||
for ix in range(nx):
|
||||
if inlet[iz, 0, ix]:
|
||||
v[iz, 0, ix] = v[iz, 1, ix]
|
||||
return
|
||||
|
||||
|
||||
@numba.njit(cache=True, fastmath=True)
|
||||
def divergence(u, v, w, solid, dx, dy, dz):
|
||||
nz, ny, nx = solid.shape
|
||||
d = np.zeros((nz, ny, nx))
|
||||
for iz in range(nz):
|
||||
for iy in range(ny):
|
||||
for ix in range(nx):
|
||||
if solid[iz, iy, ix]:
|
||||
continue
|
||||
d[iz, iy, ix] = (u[iz, iy, ix + 1] - u[iz, iy, ix]) / dx + \
|
||||
(v[iz, iy + 1, ix] - v[iz, iy, ix]) / dy + \
|
||||
(w[iz + 1, iy, ix] - w[iz, iy, ix]) / dz
|
||||
return d
|
||||
|
||||
|
||||
@numba.njit(cache=True, fastmath=True)
|
||||
def correct(u, v, w, p, solid, dt, rho, dx, dy, dz):
|
||||
nz, ny, nx = solid.shape
|
||||
c = dt / rho
|
||||
for iz in range(nz):
|
||||
for iy in range(ny):
|
||||
for ix in range(1, nx):
|
||||
if solid[iz, iy, ix] or solid[iz, iy, ix - 1]:
|
||||
continue
|
||||
u[iz, iy, ix] -= c * (p[iz, iy, ix] - p[iz, iy, ix - 1]) / dx
|
||||
for iz in range(nz):
|
||||
for iy in range(1, ny):
|
||||
for ix in range(nx):
|
||||
if solid[iz, iy, ix] or solid[iz, iy - 1, ix]:
|
||||
continue
|
||||
v[iz, iy, ix] -= c * (p[iz, iy, ix] - p[iz, iy - 1, ix]) / dy
|
||||
for iz in range(1, nz):
|
||||
for iy in range(ny):
|
||||
for ix in range(nx):
|
||||
if solid[iz, iy, ix] or solid[iz - 1, iy, ix]:
|
||||
continue
|
||||
w[iz, iy, ix] -= c * (p[iz, iy, ix] - p[iz - 1, iy, ix]) / dz
|
||||
return
|
||||
|
||||
|
||||
def build_poisson(solid, outlet, inlet):
|
||||
n = NX * NY * NZ
|
||||
idx = np.arange(n).reshape(NZ, NY, NX)
|
||||
fluid = ~solid
|
||||
rows, cols, vals = [], [], []
|
||||
|
||||
def add(r, c, v):
|
||||
rows.append(np.asarray(r).ravel()); cols.append(np.asarray(c).ravel())
|
||||
vals.append(np.asarray(v).ravel())
|
||||
|
||||
diag = np.zeros((NZ, NY, NX))
|
||||
m = fluid[:, :, :-1] & fluid[:, :, 1:]
|
||||
c = 1.0 / DX ** 2
|
||||
add(idx[:, :, :-1][m], idx[:, :, 1:][m], np.full(m.sum(), -c))
|
||||
add(idx[:, :, 1:][m], idx[:, :, :-1][m], np.full(m.sum(), -c))
|
||||
diag[:, :, :-1] += m * c; diag[:, :, 1:] += m * c
|
||||
m = fluid[:, :-1, :] & fluid[:, 1:, :]
|
||||
c = 1.0 / DY ** 2
|
||||
add(idx[:, :-1, :][m], idx[:, 1:, :][m], np.full(m.sum(), -c))
|
||||
add(idx[:, 1:, :][m], idx[:, :-1, :][m], np.full(m.sum(), -c))
|
||||
diag[:, :-1, :] += m * c; diag[:, 1:, :] += m * c
|
||||
m = fluid[:-1, :, :] & fluid[1:, :, :]
|
||||
c = 1.0 / DZ ** 2
|
||||
add(idx[:-1, :, :][m], idx[1:, :, :][m], np.full(m.sum(), -c))
|
||||
add(idx[1:, :, :][m], idx[:-1, :, :][m], np.full(m.sum(), -c))
|
||||
diag[:-1, :, :] += m * c; diag[1:, :, :] += m * c
|
||||
|
||||
fl = fluid.ravel()
|
||||
df = diag.ravel().copy()
|
||||
df = np.where(fl & (df == 0.0), 1.0, df)
|
||||
df[~fl] = 1.0
|
||||
add(np.arange(n), np.arange(n), df)
|
||||
A = csr_matrix((np.concatenate(vals), (np.concatenate(rows), np.concatenate(cols))),
|
||||
shape=(n, n))
|
||||
of = outlet.ravel()
|
||||
inf = inlet.ravel()
|
||||
fix = of | inf
|
||||
if fix.any():
|
||||
oi = np.where(fix)[0]
|
||||
A = A.tolil(); A[oi, :] = 0.0; A[oi, oi] = 1.0; A = A.tocsr()
|
||||
return A, idx
|
||||
|
||||
|
||||
def main():
|
||||
t_all = time.time()
|
||||
print(f"[rans] staggered grid {NX}x{NY}x{NZ}={NX*NY*NZ} PPM={PPM}")
|
||||
solid, plate, inlet, outlet = build_masks()
|
||||
print(f"[rans] fluid={int((~solid).sum())} plate={int(plate.sum())} R_porous={R_POROUS:.0f}")
|
||||
|
||||
A, idx = build_poisson(solid, outlet, inlet)
|
||||
print(f"[rans] poisson nnz={A.nnz}, AMG ...")
|
||||
t0 = time.time()
|
||||
ml = pyamg.smoothed_aggregation_solver(A, max_coarse=300)
|
||||
print(f"[rans] AMG {time.time()-t0:.2f}s levels={len(ml.levels)}")
|
||||
|
||||
u = np.zeros((NZ, NY, NX + 1))
|
||||
v = np.zeros((NZ, NY + 1, NX))
|
||||
w = np.zeros((NZ + 1, NY, NX))
|
||||
p = np.zeros((NZ, NY, NX))
|
||||
dt = CFL * min(DX, DY, DZ) / (U_IN + 1e-9)
|
||||
print(f"[rans] 恒功率风机: P={P_FAN:.0f} W 额定({Q_RATED*3600:.0f} m3/h, {DP_RATED:.0f} Pa) "
|
||||
f"关死={DP_SHUT:.0f} Pa nu_eff={NU_EFF:.4f} dt={dt:.3e}")
|
||||
|
||||
p_fan = DP_RATED
|
||||
t0 = time.time()
|
||||
for it in range(MAX_STEPS):
|
||||
predict(u, v, w, solid, plate, inlet, outlet, NU_EFF, dt, DX, DY, DZ, R_POROUS, U_IN)
|
||||
d = divergence(u, v, w, solid, DX, DY, DZ)
|
||||
rhs = (RHO / dt) * d
|
||||
rhs[outlet] = 0.0
|
||||
rhs[inlet] = p_fan
|
||||
sol = ml.solve(rhs.ravel(), x0=p.ravel(), tol=1e-8, maxiter=200)
|
||||
p = sol.reshape(NZ, NY, NX)
|
||||
correct(u, v, w, p, solid, dt, RHO, DX, DY, DZ)
|
||||
apply_bc(u, v, w, solid, inlet, outlet, U_IN)
|
||||
# 测进口流量 -> 更新风机压力
|
||||
iz_in = int(round((-0.45 - Y0) / DY))
|
||||
m_in = ~solid[:, iz_in, :]
|
||||
q_in = abs(v[:, iz_in, :][m_in].sum()) * DX * DZ
|
||||
p_fan = fan_dp(q_in)
|
||||
if it % 100 == 0:
|
||||
uc = 0.5 * (u[:, :, :-1] + u[:, :, 1:])
|
||||
vc = 0.5 * (v[:, :-1, :] + v[:, 1:, :])
|
||||
wc = 0.5 * (w[:-1, :, :] + w[1:, :, :])
|
||||
sp = np.sqrt(uc ** 2 + vc ** 2 + wc ** 2)
|
||||
print(f" it={it:5d} t={time.time()-t0:.1f}s |U|max={sp.max():.2f} "
|
||||
f"Q_in={q_in*3600:.0f} p_fan={p_fan:.0f} p=[{p.min():.0f},{p.max():.0f}]")
|
||||
if time.time() - t0 > TIME_BUDGET:
|
||||
break
|
||||
|
||||
izp = int(round((PLATE_T - Z0) / DZ))
|
||||
q_plate = w[izp][~solid[izp]].sum() * DX * DY * 3600.0
|
||||
iz_in = int(round((-0.45 - Y0) / DY))
|
||||
m_in = ~solid[:, iz_in, :]
|
||||
q_in = v[:, iz_in, :][m_in].sum() * DX * DZ * 3600.0
|
||||
iz_out = NZ - 2
|
||||
q_out = w[iz_out][~solid[iz_out]].sum() * DX * DY * 3600.0
|
||||
a0, a1 = int(round((-0.30 - Z0) / DZ)), int(round((-0.05 - Z0) / DZ))
|
||||
p_pl = p[a0:a1][~solid[a0:a1]].mean()
|
||||
b0, b1 = int(round((0.15 - Z0) / DZ)), int(round((0.35 - Z0) / DZ))
|
||||
p_fb = p[b0:b1][~solid[b0:b1]].mean()
|
||||
print(f"[rans] steps={it+1} solve={time.time()-t0:.1f}s")
|
||||
print(f"[rans] plate Q={q_plate:.0f} inlet Q={q_in:.0f} outlet Q={q_out:.0f} m3/h")
|
||||
print(f"[rans] p_风箱={p_pl:.0f} p_沸腾段={p_fb:.0f} Δp={p_pl-p_fb:+.0f} Pa")
|
||||
out = os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))), "docs", "assets")
|
||||
np.savez_compressed(os.path.join(out, "fbed_rans.npz"), u=u, v=v, w=w, p=p,
|
||||
solid=solid, plate=plate)
|
||||
print(f"[rans] total {time.time()-t_all:.1f}s")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -0,0 +1,507 @@
|
||||
#!/usr/bin/env python
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
水平卧式流化床 —— 空床气相流场快速仿真 (3D D3Q19 LBM, GPU)
|
||||
============================================================
|
||||
几何 (m)
|
||||
--------
|
||||
布风板 : 2400 (x, 长度) x 567 (y, 宽度), 水平
|
||||
孔板 : 孔径 2 mm, 开孔率 5.6%, 孔轴向 +x 倾斜 5 度 (朝出料端)
|
||||
进风分区 : 沿长度均分 4 段 (4 x 600 mm), 每段下方独立风箱
|
||||
风箱 : 每段 600 x 567 x 450 mm, DN300 进风管从 567 侧面水平居中进入
|
||||
沸腾/扩大段 : 由 2400 x 567 在 1000 mm 高度内矩形渐扩至 2400 x 1200
|
||||
上部稳压段 : 2400 x 1200 x 400 mm
|
||||
排风 : 顶部居中 DN800 圆管
|
||||
|
||||
工况
|
||||
----
|
||||
单段风量 : 2500 m^3/h (4 段共 10000 m^3/h)
|
||||
床板压差 : 2900 Pa (布风板上下)
|
||||
表观气速 : U_sup = 10000/3600 / (2.4*0.567) = 2.04 m/s
|
||||
管口气速 : DN300 单管 9.82 m/s ; 排风 DN800 5.53 m/s
|
||||
|
||||
方法
|
||||
----
|
||||
D3Q19 LBM (局部显式 stream+collide, GPU/numba.cuda), 目标 < 20 s。
|
||||
布风板用"多孔执行面"模型: 以表观气速 (含 5 度倾斜分量) 均匀喷射, 质量守恒。
|
||||
进风管 / 排风管用速度边界 (velocity inlet / outlet)。
|
||||
2900 Pa 远超声速不可压 LBM 可表达范围, 故用"给定风量"而非"给定压差"驱动。
|
||||
|
||||
输出
|
||||
----
|
||||
docs/assets/fbed-overview.png 四联总览
|
||||
docs/assets/fbed-slice-xz.png 垂直切面 y=0 (风箱-布风板-沸腾段-排风)
|
||||
docs/assets/fbed-slice-yz.png 垂直切面 x=1.2 (侧向渐扩)
|
||||
docs/assets/fbed-slice-xy-plate.png 水平切面 z=0.15 (布风板近区射流)
|
||||
docs/assets/fbed-slice-xy-free.png 水平切面 z=0.9 (沸腾段)
|
||||
|
||||
运行
|
||||
----
|
||||
python tools/fbed_simulation.py
|
||||
"""
|
||||
import os
|
||||
import time
|
||||
|
||||
import numpy as np
|
||||
from numba import cuda
|
||||
|
||||
import matplotlib
|
||||
matplotlib.use("Agg")
|
||||
import matplotlib.pyplot as plt
|
||||
from matplotlib.colors import LinearSegmentedColormap
|
||||
|
||||
# --------------------------------------------------------------------------
|
||||
# 几何参数 (m)
|
||||
# --------------------------------------------------------------------------
|
||||
BED_L, BED_W = 2.4, 0.567
|
||||
N_SEC = 4
|
||||
SEC_L = BED_L / N_SEC
|
||||
PLENUM_H = 0.45
|
||||
PLATE_T = 0.002
|
||||
|
||||
FB_H, FB_W = 1.0, 1.2
|
||||
UPPER_H = 0.4
|
||||
REDUCER_H = 0.5
|
||||
DUCT_H = 0.3
|
||||
EXH_R = 0.4
|
||||
|
||||
PIPE_R = 0.15
|
||||
PIPE_Y0 = -0.65
|
||||
|
||||
X0, X1 = 0.0, BED_L
|
||||
Y0, Y1 = PIPE_Y0, 0.65
|
||||
Z1 = FB_H + UPPER_H + REDUCER_H + DUCT_H
|
||||
Z0 = -PLENUM_H
|
||||
PPM = int(os.environ.get("FBED_PPM", "100"))
|
||||
NX = int(round((X1 - X0) * PPM))
|
||||
NY = int(round((Y1 - Y0) * PPM))
|
||||
NZ = int(round((Z1 - Z0) * PPM))
|
||||
DX = 1.0 / PPM
|
||||
|
||||
# --------------------------------------------------------------------------
|
||||
# 工况
|
||||
# --------------------------------------------------------------------------
|
||||
Q_SEC = 2500.0 / 3600.0
|
||||
A_BED = BED_L * BED_W
|
||||
U_SUP = Q_SEC * N_SEC / A_BED
|
||||
TILT = np.deg2rad(5.0)
|
||||
TILT_DIR = -1.0 # 绿色箭头方向: 孔轴向上-左 (-x)
|
||||
UP_X = TILT_DIR * U_SUP * np.sin(TILT)
|
||||
UP_Z = U_SUP * np.cos(TILT)
|
||||
TILT_TAN = np.tan(TILT)
|
||||
A_PIPE = np.pi * PIPE_R ** 2
|
||||
U_IN = Q_SEC / A_PIPE
|
||||
A_EXH = np.pi * EXH_R ** 2
|
||||
U_OUT = Q_SEC * N_SEC / A_EXH
|
||||
|
||||
U_LB = float(os.environ.get("FBED_ULB", "0.03"))
|
||||
SCALE_U = U_IN / U_LB
|
||||
U_IN_LB = U_LB
|
||||
U_OUT_LB = U_OUT / SCALE_U
|
||||
UPX_LB = UP_X / SCALE_U
|
||||
UPZ_LB = UP_Z / SCALE_U
|
||||
NU_LB = 0.04
|
||||
TAU = 0.5 + 3.0 * NU_LB
|
||||
|
||||
DP_RATED = float(os.environ.get("FBED_DP", "2900.0"))
|
||||
Q_RATED = (2500.0 * 4) / 3600.0
|
||||
P_FAN = DP_RATED * Q_RATED
|
||||
DP_SHUT = DP_RATED * 1.5
|
||||
|
||||
|
||||
def fan_dp(q_m3s):
|
||||
"""恒功率风机曲线 dp*Q=P -> dp=P/Q (关死压力截断)"""
|
||||
if q_m3s < 1e-6:
|
||||
return DP_SHUT
|
||||
return min(P_FAN / q_m3s, DP_SHUT)
|
||||
|
||||
TIME_BUDGET = float(os.environ.get("FBED_TB", "120.0"))
|
||||
MAX_STEPS = 40000
|
||||
|
||||
EX = np.array([0, 1, -1, 0, 0, 0, 0, 1, -1, 1, -1, 1, -1, 1, -1, 0, 0, 0, 0], dtype=np.int32)
|
||||
EY = np.array([0, 0, 0, 1, -1, 0, 0, 1, -1, -1, 1, 0, 0, 0, 0, 1, -1, 1, -1], dtype=np.int32)
|
||||
EZ = np.array([0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 1, -1, -1, 1, 1, -1, -1, 1], dtype=np.int32)
|
||||
WT = np.array([1/3] + [1/18]*6 + [1/36]*12, dtype=np.float64)
|
||||
OPP = np.array([0, 2, 1, 4, 3, 6, 5, 8, 7, 10, 9, 12, 11, 14, 13, 16, 15, 18, 17], dtype=np.int32)
|
||||
|
||||
XC = X0 + (np.arange(NX) + 0.5) / PPM
|
||||
YC = Y0 + (np.arange(NY) + 0.5) / PPM
|
||||
ZC = Z0 + (np.arange(NZ) + 0.5) / PPM
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------
|
||||
# GPU 核函数
|
||||
# --------------------------------------------------------------------------
|
||||
@cuda.jit
|
||||
def collide_kernel(f, rho, ux, uy, uz, solid, plate, inlet, outlet,
|
||||
ex, ey, ez, w, tau, drho, ttilt, nz, ny, nx):
|
||||
idx = cuda.grid(1)
|
||||
if idx >= nz * ny * nx:
|
||||
return
|
||||
ix = idx % nx
|
||||
iy = (idx // nx) % ny
|
||||
iz = idx // (nx * ny)
|
||||
if solid[iz, iy, ix]:
|
||||
return
|
||||
r = 0.0
|
||||
mx = 0.0
|
||||
my = 0.0
|
||||
mz = 0.0
|
||||
for k in range(19):
|
||||
fk = f[k, iz, iy, ix]
|
||||
r += fk
|
||||
mx += ex[k] * fk
|
||||
my += ey[k] * fk
|
||||
mz += ez[k] * fk
|
||||
inv = 1.0 / r
|
||||
u = mx * inv
|
||||
v = my * inv
|
||||
ww = mz * inv
|
||||
om = 1.0 / tau
|
||||
uu = u * u + v * v + ww * ww
|
||||
for k in range(19):
|
||||
eu = ex[k] * u + ey[k] * v + ez[k] * ww
|
||||
feq = w[k] * r * (1.0 + 3.0 * eu + 4.5 * eu * eu - 1.5 * uu)
|
||||
f[k, iz, iy, ix] -= (f[k, iz, iy, ix] - feq) * om
|
||||
if plate[iz, iy, ix]:
|
||||
# 布风板 5 度倾斜: 开孔处施加水平动量 (绿箭头方向 -x)
|
||||
fx = -ttilt * r * ww
|
||||
for k in range(19):
|
||||
f[k, iz, iy, ix] += 3.0 * w[k] * ex[k] * fx
|
||||
u += fx * inv
|
||||
if inlet[iz, iy, ix]:
|
||||
# 压力进口: 固定密度(= 风机压力), 速度零梯度自由
|
||||
r = 1.0 + drho
|
||||
uu = u * u + v * v + ww * ww
|
||||
for k in range(19):
|
||||
eu = ex[k] * u + ey[k] * v + ez[k] * ww
|
||||
f[k, iz, iy, ix] = w[k] * r * (1.0 + 3.0 * eu + 4.5 * eu * eu - 1.5 * uu)
|
||||
elif outlet[iz, iy, ix]:
|
||||
r = 1.0
|
||||
uu = u * u + v * v + ww * ww
|
||||
for k in range(19):
|
||||
eu = ex[k] * u + ey[k] * v + ez[k] * ww
|
||||
f[k, iz, iy, ix] = w[k] * r * (1.0 + 3.0 * eu + 4.5 * eu * eu - 1.5 * uu)
|
||||
rho[iz, iy, ix] = r
|
||||
ux[iz, iy, ix] = u
|
||||
uy[iz, iy, ix] = v
|
||||
uz[iz, iy, ix] = ww
|
||||
|
||||
|
||||
@cuda.jit
|
||||
def stream_kernel(f, ftmp, solid, ex, ey, ez, opp, nz, ny, nx):
|
||||
idx = cuda.grid(1)
|
||||
if idx >= nz * ny * nx:
|
||||
return
|
||||
ix = idx % nx
|
||||
iy = (idx // nx) % ny
|
||||
iz = idx // (nx * ny)
|
||||
if solid[iz, iy, ix]:
|
||||
return
|
||||
for k in range(19):
|
||||
sz = iz - ez[k]
|
||||
sy = iy - ey[k]
|
||||
sx = ix - ex[k]
|
||||
if 0 <= sz < nz and 0 <= sy < ny and 0 <= sx < nx and not solid[sz, sy, sx]:
|
||||
ftmp[k, iz, iy, ix] = f[k, sz, sy, sx]
|
||||
else:
|
||||
ftmp[k, iz, iy, ix] = f[opp[k], iz, iy, ix]
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------
|
||||
# 几何
|
||||
# --------------------------------------------------------------------------
|
||||
def _xi(v):
|
||||
return int(round((v - X0) * PPM))
|
||||
|
||||
|
||||
def _yi(v):
|
||||
return int(round((v - Y0) * PPM))
|
||||
|
||||
|
||||
def _zi(v):
|
||||
return int(round((v - Z0) * PPM))
|
||||
|
||||
|
||||
def _clip(i, n):
|
||||
return max(0, min(i, n))
|
||||
|
||||
|
||||
def carve_box(solid, xa, xb, ya, yb, za, zb):
|
||||
ix0, ix1 = _clip(_xi(xa), NX), _clip(_xi(xb), NX)
|
||||
iy0, iy1 = _clip(_yi(ya), NY), _clip(_yi(yb), NY)
|
||||
iz0, iz1 = _clip(_zi(za), NZ), _clip(_zi(zb), NZ)
|
||||
solid[iz0:iz1, iy0:iy1, ix0:ix1] = False
|
||||
|
||||
|
||||
def build_geometry():
|
||||
solid = np.ones((NZ, NY, NX), dtype=np.bool_)
|
||||
yw = BED_W / 2.0
|
||||
|
||||
carve_box(solid, 0.0, BED_L, -yw, yw, -PLENUM_H, 0.0)
|
||||
for k in range(1, N_SEC):
|
||||
xk = k * SEC_L
|
||||
carve_box(solid, xk - 0.0125, xk + 0.0125, -yw, yw, -PLENUM_H, 0.0)
|
||||
for k in range(1, N_SEC):
|
||||
xk = k * SEC_L
|
||||
solid[_zi(-PLENUM_H):_zi(0.0), _yi(-yw):_yi(yw), _xi(xk - 0.0125):_xi(xk + 0.0125)] = True
|
||||
|
||||
carve_box(solid, 0.0, BED_L, -yw, yw, 0.0, PLATE_T)
|
||||
|
||||
for iz in range(_zi(PLATE_T), _zi(FB_H)):
|
||||
z = ZC[iz]
|
||||
w = yw + (FB_W / 2.0 - yw) * (z / FB_H)
|
||||
solid[iz, :, :] = np.where(
|
||||
((YC[:, None] >= -w) & (YC[:, None] <= w) &
|
||||
(XC[None, :] >= 0.0) & (XC[None, :] <= BED_L)),
|
||||
False, solid[iz, :, :])
|
||||
carve_box(solid, 0.0, BED_L, -FB_W / 2.0, FB_W / 2.0, FB_H, FB_H + UPPER_H)
|
||||
|
||||
z_r0 = FB_H + UPPER_H
|
||||
z_r1 = z_r0 + REDUCER_H
|
||||
for iz in range(_zi(z_r0), _zi(z_r1)):
|
||||
t = (ZC[iz] - z_r0) / REDUCER_H
|
||||
xh = BED_L / 2.0 - (BED_L / 2.0 - 0.4) * t
|
||||
yh = FB_W / 2.0 - (FB_W / 2.0 - 0.4) * t
|
||||
solid[iz] = np.where(
|
||||
(np.abs(XC[None, :] - BED_L / 2.0) <= xh) & (np.abs(YC[:, None]) <= yh),
|
||||
False, solid[iz])
|
||||
|
||||
cyl = (XC[None, :] - BED_L / 2.0) ** 2 + (YC[:, None] - 0.0) ** 2 <= EXH_R ** 2
|
||||
iz_e0, iz_e1 = _zi(z_r1), NZ
|
||||
solid[iz_e0:iz_e1] = np.where(cyl[None, :, :], False, solid[iz_e0:iz_e1])
|
||||
|
||||
iy_p0, iy_p1 = _yi(PIPE_Y0), _yi(-yw)
|
||||
for k in range(N_SEC):
|
||||
xk = SEC_L * (k + 0.5)
|
||||
zk = -PLENUM_H / 2.0
|
||||
m = ((XC[None, :] - xk) ** 2 + (ZC[:, None] - zk) ** 2) <= PIPE_R ** 2
|
||||
m3 = np.repeat(m[:, None, :], iy_p1 - iy_p0, axis=1)
|
||||
region = solid[:, iy_p0:iy_p1, :]
|
||||
region[m3] = False
|
||||
solid[:, iy_p0:iy_p1, :] = region
|
||||
|
||||
iz_plate = _zi(0.0)
|
||||
# 多孔布风板: 固体 + 规则开孔 (孔 hs x hs 格, 间距 ss 格)
|
||||
izp0, izp1 = _zi(0.0), max(_zi(PLATE_T), _zi(0.0) + 1)
|
||||
iy0p, iy1p = _yi(-yw), _yi(yw)
|
||||
solid[izp0:izp1, iy0p:iy1p, :] = True
|
||||
hs = int(os.environ.get("FBED_HOLE", "1"))
|
||||
ss = int(os.environ.get("FBED_SPACE", "4"))
|
||||
for iz in range(izp0, izp1):
|
||||
for iy in range(iy0p, iy1p):
|
||||
for ix in range(NX):
|
||||
if ((iy - iy0p) % ss < hs) and (ix % ss < hs):
|
||||
solid[iz, iy, ix] = False
|
||||
plate = np.zeros_like(solid)
|
||||
plate[izp0:izp1, iy0p:iy1p, :] = ~solid[izp0:izp1, iy0p:iy1p, :]
|
||||
|
||||
inlet = np.zeros_like(solid)
|
||||
inlet[:, iy_p0:iy_p0 + 3, :] = ~solid[:, iy_p0:iy_p0 + 3, :]
|
||||
|
||||
outlet = np.zeros_like(solid)
|
||||
outlet[NZ - 2:, :, :] = ~solid[NZ - 2:, :, :]
|
||||
return solid, plate, inlet, outlet
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------
|
||||
# 求解
|
||||
# --------------------------------------------------------------------------
|
||||
def run_solver(solid, plate, inlet, outlet):
|
||||
f = np.zeros((19, NZ, NY, NX), dtype=np.float32)
|
||||
for k in range(19):
|
||||
f[k] = WT[k]
|
||||
ftmp = np.zeros_like(f)
|
||||
rho = np.ones((NZ, NY, NX), dtype=np.float32)
|
||||
ux = np.zeros_like(rho)
|
||||
uy = np.zeros_like(rho)
|
||||
uz = np.zeros_like(rho)
|
||||
|
||||
d = {name: cuda.to_device(a) for name, a in [
|
||||
("f", f), ("ftmp", ftmp), ("rho", rho), ("ux", ux), ("uy", uy), ("uz", uz),
|
||||
("solid", solid), ("plate", plate), ("inlet", inlet), ("outlet", outlet),
|
||||
("ex", EX), ("ey", EY), ("ez", EZ), ("w", WT), ("opp", OPP)]}
|
||||
|
||||
n = NX * NY * NZ
|
||||
tpb = 128
|
||||
bpg = (n + tpb - 1) // tpb
|
||||
|
||||
drho_fan = 3.0 * DP_RATED / (1.2 * SCALE_U ** 2)
|
||||
cuda.synchronize()
|
||||
t0 = time.time()
|
||||
step = 0
|
||||
for step in range(MAX_STEPS):
|
||||
collide_kernel[bpg, tpb](d["f"], d["rho"], d["ux"], d["uy"], d["uz"],
|
||||
d["solid"], d["plate"], d["inlet"], d["outlet"],
|
||||
d["ex"], d["ey"], d["ez"], d["w"], TAU,
|
||||
drho_fan, 0.0, NZ, NY, NX)
|
||||
stream_kernel[bpg, tpb](d["f"], d["ftmp"], d["solid"],
|
||||
d["ex"], d["ey"], d["ez"], d["opp"], NZ, NY, NX)
|
||||
d["f"], d["ftmp"] = d["ftmp"], d["f"]
|
||||
if step % 100 == 99:
|
||||
cuda.synchronize()
|
||||
vy = d["uy"].copy_to_host()
|
||||
iz_in = _yi(-0.45)
|
||||
q_in = float(vy[:, iz_in, :][~solid[:, iz_in, :]].sum()) * SCALE_U * DX * DX
|
||||
dp = fan_dp(abs(q_in))
|
||||
drho_fan = 3.0 * dp / (1.2 * SCALE_U ** 2)
|
||||
if step % 2000 == 1999:
|
||||
print(f" step {step+1:6d} Q_in={abs(q_in)*3600:.0f} m3/h "
|
||||
f"dp_fan={dp:.0f} Pa t={time.time()-t0:.1f}s")
|
||||
if time.time() - t0 > TIME_BUDGET:
|
||||
break
|
||||
cuda.synchronize()
|
||||
solve_dt = time.time() - t0
|
||||
steps = step + 1
|
||||
return (d["ux"].copy_to_host(), d["uy"].copy_to_host(), d["uz"].copy_to_host(),
|
||||
d["rho"].copy_to_host(), steps, solve_dt)
|
||||
|
||||
|
||||
def metrics(ux, uy, uz, rho, solid, plate, inlet, outlet, steps, solve_dt):
|
||||
fluid = ~solid
|
||||
speed = np.sqrt(ux ** 2 + uy ** 2 + uz ** 2) * SCALE_U
|
||||
umax = float(speed[fluid].max())
|
||||
iz_max, iy_max, ix_max = np.unravel_index(np.argmax(np.where(fluid, speed, 0)), speed.shape)
|
||||
|
||||
cell_a = DX * DX
|
||||
iz_plate = _zi(PLATE_T)
|
||||
fl_p = ~solid[iz_plate]
|
||||
q_plate = float(uz[iz_plate][fl_p].sum()) * SCALE_U * cell_a * 3600.0
|
||||
iy_in = _yi(-0.45)
|
||||
q_in = float(uy[:, iy_in, :][~solid[:, iy_in, :]].sum()) * SCALE_U * cell_a * 3600.0
|
||||
q_out = float(uz[NZ - 2][outlet[NZ - 2]].sum()) * SCALE_U * cell_a * 3600.0
|
||||
|
||||
p = (rho - 1.0) / 3.0 * (1.2 * SCALE_U ** 2)
|
||||
m_pl = ~solid[_zi(-0.30):_zi(-0.05)]
|
||||
p_pl = p[_zi(-0.30):_zi(-0.05)][m_pl].mean()
|
||||
m_fb = ~solid[_zi(PLATE_T + 0.05):_zi(PLATE_T + 0.25)]
|
||||
p_fb = p[_zi(PLATE_T + 0.05):_zi(PLATE_T + 0.25)][m_fb].mean()
|
||||
r_fluid = rho[~solid]
|
||||
print(f"[fbed] grid {NX}x{NY}x{NZ}={NX*NY*NZ}, steps={steps}, solve={solve_dt:.2f} s")
|
||||
print(f"[fbed] |U|max={umax:.2f} m/s @ (x={XC[ix_max]:.2f}, y={YC[iy_max]:.2f}, z={ZC[iz_max]:.2f})")
|
||||
print(f"[fbed] plate Q={q_plate:.0f} inlet Q={q_in:.0f} outlet Q={q_out:.0f} m3/h (target 10000)")
|
||||
print(f"[fbed] p_风箱={p_pl:.0f} Pa p_沸腾段={p_fb:.0f} Pa Δp(风箱-沸腾段)={p_pl-p_fb:+.0f} Pa")
|
||||
print(f"[fbed] rho=[{r_fluid.min():.4f},{r_fluid.max():.4f}] p99={np.percentile(r_fluid,99):.4f}")
|
||||
return speed, p
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------
|
||||
# 绘图
|
||||
# --------------------------------------------------------------------------
|
||||
BG = "#0d1117"
|
||||
FG = "#e6edf3"
|
||||
MUTED = "#8b949e"
|
||||
ACCENT = "#58a6ff"
|
||||
CMAP_FLOW = LinearSegmentedColormap.from_list(
|
||||
"flow", ["#0a1626", "#123c63", "#1f6feb", "#39d0d8", "#7ee787", "#f0f6fc"])
|
||||
plt.rcParams.update({
|
||||
"font.sans-serif": ["Microsoft YaHei", "SimHei", "DejaVu Sans"],
|
||||
"axes.unicode_minus": False,
|
||||
"figure.facecolor": BG, "savefig.facecolor": BG, "axes.facecolor": BG,
|
||||
"text.color": FG, "axes.labelcolor": FG, "xtick.color": MUTED,
|
||||
"ytick.color": MUTED, "axes.edgecolor": "#30363d",
|
||||
})
|
||||
|
||||
OUT = os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))),
|
||||
"docs", "assets")
|
||||
|
||||
|
||||
def plane(speed, ux, uy, uz, solid, which, idx):
|
||||
if which == "xz":
|
||||
return (speed[:, idx, :], ux[:, idx, :], uz[:, idx, :], solid[:, idx, :],
|
||||
[X0, X1, Z0, Z1], "x (m)", "z (m)")
|
||||
if which == "yz":
|
||||
return (speed[:, :, idx], uy[:, :, idx], uz[:, :, idx], solid[:, :, idx],
|
||||
[Y0, Y1, Z0, Z1], "y (m)", "z (m)")
|
||||
return (speed[idx, :, :], ux[idx, :, :], uy[idx, :, :], solid[idx, :, :],
|
||||
[X0, X1, Y0, Y1], "x (m)", "y (m)")
|
||||
|
||||
|
||||
def draw_plane(ax, data, vmax, density=1.1):
|
||||
s2, u2, v2, sol2, ext, xl, yl = data
|
||||
masked = np.where(sol2, np.nan, s2)
|
||||
im = ax.imshow(masked, origin="lower", extent=ext, cmap=CMAP_FLOW,
|
||||
vmin=0, vmax=vmax, interpolation="bilinear", aspect="equal")
|
||||
ny, nx = sol2.shape
|
||||
xs = np.linspace(ext[0], ext[1], nx)
|
||||
ys = np.linspace(ext[2], ext[3], ny)
|
||||
U = np.where(sol2, np.nan, u2)
|
||||
V = np.where(sol2, np.nan, v2)
|
||||
with np.errstate(invalid="ignore"):
|
||||
ax.streamplot(xs, ys, U, V, density=density, color=(0.85, 0.88, 0.92, 0.4),
|
||||
linewidth=0.5, arrowsize=0.6, zorder=5)
|
||||
ax.contour(np.linspace(ext[0], ext[1], sol2.shape[1]),
|
||||
np.linspace(ext[2], ext[3], sol2.shape[0]),
|
||||
sol2.astype(float), levels=[0.5], colors="#adbac7", linewidths=0.8, zorder=6)
|
||||
ax.set_xlabel(xl, fontsize=10)
|
||||
ax.set_ylabel(yl, fontsize=10)
|
||||
ax.set_xlim(ext[0], ext[1])
|
||||
ax.set_ylim(ext[2], ext[3])
|
||||
return im
|
||||
|
||||
|
||||
def save_slice(speed, ux, uy, uz, solid, which, idx, title, name, vmax):
|
||||
data = plane(speed, ux, uy, uz, solid, which, idx)
|
||||
fig, ax = plt.subplots(figsize=(10.5, 7.6))
|
||||
im = draw_plane(ax, data, vmax)
|
||||
cb = fig.colorbar(im, ax=ax, fraction=0.036, pad=0.02)
|
||||
cb.set_label("速度大小 |U| (m/s)", fontsize=10)
|
||||
cb.outline.set_edgecolor("#30363d")
|
||||
ax.set_title(title, fontsize=13, pad=10)
|
||||
fig.tight_layout()
|
||||
fig.savefig(os.path.join(OUT, name), dpi=540)
|
||||
plt.close(fig)
|
||||
print(" ->", name)
|
||||
|
||||
|
||||
def save_overview(speed, ux, uy, uz, solid, steps, solve_dt, umax):
|
||||
fig, axes = plt.subplots(2, 2, figsize=(16.0, 12.0))
|
||||
specs = [
|
||||
("xz", _yi(0.0), "(a) 垂直切面 y=0 —— 风箱 / 布风板 / 沸腾段 / 排风", 4.0),
|
||||
("yz", _xi(BED_L / 2.0), "(b) 垂直切面 x=1.2 —— 侧向渐扩（567→1200 mm）", 4.0),
|
||||
("xy", _zi(0.15), "(c) 水平切面 z=0.15 —— 布风板近区射流", 3.0),
|
||||
("xy", _zi(0.90), "(d) 水平切面 z=0.90 —— 沸腾段", 3.0),
|
||||
]
|
||||
for ax, (which, idx, title, vmax) in zip(axes.ravel(), specs):
|
||||
data = plane(speed, ux, uy, uz, solid, which, idx)
|
||||
im = draw_plane(ax, data, vmax)
|
||||
cb = fig.colorbar(im, ax=ax, fraction=0.036, pad=0.02)
|
||||
cb.set_label("|U| (m/s)", fontsize=9.5)
|
||||
cb.outline.set_edgecolor("#30363d")
|
||||
cb.ax.tick_params(labelsize=8.5)
|
||||
ax.set_title(title, fontsize=12, pad=8)
|
||||
fig.suptitle("水平卧式流化床 空床流场 —— 4 段侧进风 / 5.6% 开孔倾斜布风板 / 顶部 DN800 排风",
|
||||
fontsize=16, y=0.985)
|
||||
fig.text(0.5, 0.006,
|
||||
f"求解 {solve_dt:.1f} s · {steps} 步 · |U|max={umax:.1f} m/s · "
|
||||
f"总风量 10000 m³/h · D3Q19 LBM (GPU)",
|
||||
ha="center", fontsize=10.5, color=MUTED)
|
||||
fig.tight_layout(rect=[0, 0.02, 1, 0.97])
|
||||
fig.savefig(os.path.join(OUT, "fbed-overview.png"), dpi=520)
|
||||
plt.close(fig)
|
||||
print(" -> fbed-overview.png")
|
||||
|
||||
|
||||
def main():
|
||||
os.makedirs(OUT, exist_ok=True)
|
||||
print(f"[fbed] U_sup={U_SUP:.3f} m/s U_in(DN300)={U_IN:.2f} m/s "
|
||||
f"U_out(DN800)={U_OUT:.2f} m/s tilt=5 deg")
|
||||
print(f"[fbed] 恒功率风机: P={P_FAN:.0f} W 额定({Q_RATED*3600:.0f} m3/h, {DP_RATED:.0f} Pa) "
|
||||
f"关死={DP_SHUT:.0f} Pa")
|
||||
solid, plate, inlet, outlet = build_geometry()
|
||||
ux, uy, uz, rho, steps, solve_dt = run_solver(solid, plate, inlet, outlet)
|
||||
speed, p = metrics(ux, uy, uz, rho, solid, plate, inlet, outlet, steps, solve_dt)
|
||||
|
||||
np.savez_compressed(os.path.join(OUT, "fbed_fields.npz"),
|
||||
ux=ux, uy=uy, uz=uz, rho=rho, solid=solid,
|
||||
plate=plate, inlet=inlet, outlet=outlet,
|
||||
steps=steps, solve_dt=solve_dt)
|
||||
print(" -> fbed_fields.npz")
|
||||
|
||||
t0 = time.time()
|
||||
save_overview(speed, ux, uy, uz, solid, steps, solve_dt, float(speed[~solid].max()))
|
||||
print(f"[fbed] 出图 {time.time()-t0:.2f} s")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
Reference in New Issue
Block a user