#!/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()