#!/usr/bin/env python # -*- coding: utf-8 -*- """ FFU 单向流洁净区流场仿真 Demo ================================ 用 D2Q9 格子玻尔兹曼方法 (LBM) 模拟 FFU (Fan Filter Unit, 风机过滤单元) 下方洁净区的流场, 并用半拉格朗日输运方程求解: - 温度场 T (设备发热源) - 相对湿度 RH (人员呼吸/发湿源) - 尘埃粒子浓度 C (设备 + 人员散发, 含重力沉降) - 静压场 p (LBM 直接输出, 转换为 Pa) 物理设置 -------- 计算域 : 4.0 m (宽) x 3.0 m (高), 90 px/m FFU 送风 : 顶棚满布, 面风速 0.45 m/s (ISO 14644 单向流推荐值) 回风口 : 底部两侧 (0~0.5 m, 3.5~4.0 m), 高架地板/下侧回风 障碍物 : 工作台 + 设备(热/尘源) + 操作人员(湿/尘源) + 物料架 雷诺数 : Re = U*H/nu ~ 540 (层流单向流) 网格 : 360 x 270 = 97200 格点, D2Q9 BGK, tau = 0.62 输出 ---- docs/assets/ffu-flow-field.png 速度场 + 流线 (主图) docs/assets/ffu-scalar-fields.png 温度 / 湿度 / 粒子浓度 / 压力 四联图 docs/assets/ffu-particle-tracks.png 粒子拉格朗日轨迹 docs/assets/ffu-flow-field.svg 主图矢量版 运行 ---- python tools/ffu_simulation.py """ import os import numpy as np from numba import njit import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from matplotlib.patches import Rectangle, Circle, FancyArrow from matplotlib.colors import LinearSegmentedColormap from matplotlib.collections import LineCollection # -------------------------------------------------------------------------- # 全局样式 # -------------------------------------------------------------------------- BG = "#0d1117" PANEL = "#0d1117" FG = "#e6edf3" MUTED = "#8b949e" ACCENT = "#58a6ff" WARM = "#ff7b72" COOL = "#39d0d8" GREEN = "#7ee787" plt.rcParams.update({ "font.sans-serif": ["Microsoft YaHei", "SimHei", "DejaVu Sans"], "axes.unicode_minus": False, "svg.fonttype": "path", "figure.facecolor": BG, "savefig.facecolor": BG, "axes.facecolor": BG, "text.color": FG, "axes.labelcolor": FG, "xtick.color": MUTED, "ytick.color": MUTED, "axes.edgecolor": "#30363d", }) # -------------------------------------------------------------------------- # 物理 / 数值参数 # -------------------------------------------------------------------------- W_M, H_M = 4.0, 3.0 # 计算域尺寸 [m] PPM = 90 # 像素/米 NX, NY = int(W_M * PPM), int(H_M * PPM) # 360 x 270 U0_LB = 0.08 # FFU 送风速度 [lattice] U_PHYS = 0.45 # FFU 送风速度 [m/s] SCALE_U = U_PHYS / U0_LB # lattice -> m/s NU_LB = 0.04 # 运动粘度 [lattice] TAU = 0.5 + 3.0 * NU_LB RE = U0_LB * NY / NU_LB T_AMB, T_SRC = 22.0, 35.0 # 环境温度 / 设备表面温度 [degC] RH_AMB, RH_SRC = 45.0, 70.0 # 送风湿度 / 呼吸区湿度 [%] RHO_AIR = 1.2 # 空气密度 [kg/m3] SCALE_P = RHO_AIR * SCALE_U ** 2 # lattice -> Pa DT_SCALAR = 4 # 标量输运时间步 (每 4 个 LBM 步输运一次) ADV_EVERY = 4 KAPPA_T = 0.0020 # 温度扩散系数 [lattice] KAPPA_RH = 0.0020 KAPPA_C = 0.0015 V_SETTLE = 0.004 # 粒子沉降速度 [lattice] (~2 cm/s) N_STEPS = 14000 # LBM 总步数 SCALAR_FROM = 6000 # 从该步开始输运标量 TRACK_FROM = 9000 # 从该步开始记录粒子轨迹 # D2Q9 EX = np.array([0, 1, 0, -1, 0, 1, -1, -1, 1], dtype=np.int32) EY = np.array([0, 0, 1, 0, -1, 1, 1, -1, -1], dtype=np.int32) WT = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]) OPP = np.array([0, 3, 4, 1, 2, 7, 8, 5, 6], dtype=np.int32) # -------------------------------------------------------------------------- # 几何: 单位 m # -------------------------------------------------------------------------- BENCH = (0.60, 0.00, 1.80, 0.90) # 工作台 (x0,y0,x1,y1) EQUIP = (0.85, 0.90, 1.50, 1.25) # 设备 (热源 + 尘源) BODY = (2.35, 0.00, 2.75, 1.45) # 人员躯干 HEAD_C = (2.55, 1.62, 0.15) # 人员头部 (圆心, 半径) RACK = (3.00, 0.00, 3.50, 1.15) # 物料架 GRILL_L = (0.00, 0.50) # 左侧回风口 x 范围 GRILL_R = (3.50, 4.00) # 右侧回风口 x 范围 def build_geometry(): """构建固体掩膜 + 源项掩膜""" solid = np.zeros((NY, NX), dtype=np.bool_) solid[:, :2] = True solid[:, -2:] = True solid[0, :] = True def rect(x0, y0, x1, y1, m): ix0, ix1 = int(round(x0 * PPM)), int(round(x1 * PPM)) iy0, iy1 = int(round(y0 * PPM)), int(round(y1 * PPM)) m[max(iy0, 0):min(iy1 + 1, NY), max(ix0, 0):min(ix1 + 1, NX)] = True for r in (BENCH, EQUIP, BODY, RACK): rect(*r, solid) cx, cy, rr = HEAD_C yy, xx = np.mgrid[0:NY, 0:NX] head = (xx - cx * PPM) ** 2 + (yy - cy * PPM) ** 2 <= (rr * PPM) ** 2 solid |= head # 回风口 (底部, 非固体 => 压力出口) outlet = np.zeros_like(solid) iy = 0 for g0, g1 in (GRILL_L, GRILL_R): ix0, ix1 = int(g0 * PPM), int(g1 * PPM) outlet[iy, max(ix0, 2):min(ix1, NX - 2)] = True solid[outlet] = False inlet = np.zeros_like(solid) inlet[NY - 1, 2:NX - 2] = True # 设备热/尘源: 设备外扩 3 格的流体单元 eq = np.zeros_like(solid) rect(*EQUIP, eq) shell = np.zeros_like(solid) for dy in range(-3, 4): for dx in range(-3, 4): shell |= np.roll(np.roll(eq, dy, 0), dx, 1) src_equip = shell & ~solid # 人员呼吸/发湿区: 头部前方圆形区域 br = (xx - (cx + 0.20) * PPM) ** 2 + (yy - cy * PPM) ** 2 <= (0.22 * PPM) ** 2 src_head = br & ~solid return solid, inlet, outlet, src_equip, src_head # -------------------------------------------------------------------------- # LBM 求解器 (numba) # -------------------------------------------------------------------------- @njit(cache=True, fastmath=True) def lbm_collide(f, rho, ux, uy, solid, ex, ey, w, tau): ny, nx = rho.shape for iy in range(ny): for ix in range(nx): if solid[iy, ix]: continue r = 0.0; mx = 0.0; my = 0.0 for k in range(9): fk = f[k, iy, ix] r += fk mx += ex[k] * fk my += ey[k] * fk u = mx / r v = my / r rho[iy, ix] = r ux[iy, ix] = u uy[iy, ix] = v uu = u * u + v * v for k in range(9): eu = ex[k] * u + ey[k] * v feq = w[k] * r * (1.0 + 3.0 * eu + 4.5 * eu * eu - 1.5 * uu) f[k, iy, ix] -= (f[k, iy, ix] - feq) / tau @njit(cache=True, fastmath=True) def lbm_bounce_back(f, solid, opp): ny, nx = solid.shape for iy in range(ny): for ix in range(nx): if solid[iy, ix]: for k in range(9): ko = opp[k] if k < ko: tmp = f[k, iy, ix] f[k, iy, ix] = f[ko, iy, ix] f[ko, iy, ix] = tmp @njit(cache=True, fastmath=True) def lbm_stream(f, ftmp, ex, ey): ny, nx = f.shape[1], f.shape[2] for k in range(9): dx = ex[k]; dy = ey[k] for iy in range(ny): sy = iy - dy for ix in range(nx): sx = ix - dx if 0 <= sy < ny and 0 <= sx < nx: ftmp[k, iy, ix] = f[k, sy, sx] else: ftmp[k, iy, ix] = f[k, iy, ix] @njit(cache=True, fastmath=True) def lbm_bcs(f, ux, uy, solid, inlet, outlet, ex, ey, w, u0): ny, nx = solid.shape for ix in range(nx): if inlet[ny - 1, ix]: u = 0.0; v = -u0; r = 1.0 uu = u * u + v * v for k in range(9): eu = ex[k] * u + ey[k] * v f[k, ny - 1, ix] = w[k] * r * (1.0 + 3.0 * eu + 4.5 * eu * eu - 1.5 * uu) if outlet[0, ix]: u = ux[1, ix]; v = uy[1, ix]; r = 1.0 uu = u * u + v * v for k in range(9): eu = ex[k] * u + ey[k] * v f[k, 0, ix] = w[k] * r * (1.0 + 3.0 * eu + 4.5 * eu * eu - 1.5 * uu) # -------------------------------------------------------------------------- # 标量输运 (半拉格朗日 + 扩散) # -------------------------------------------------------------------------- @njit(cache=True, fastmath=True) def advect_scalars(F, ux, uy, solid, kappa, dt, vs): nc, ny, nx = F.shape out = np.empty_like(F) for c in range(nc): k = kappa[c] s = vs[c] for iy in range(ny): for ix in range(nx): fc = F[c, iy, ix] if solid[iy, ix]: out[c, iy, ix] = 0.0 continue xd = ix - ux[iy, ix] * dt yd = iy - (uy[iy, ix] + s) * dt if xd < 0.0: xd = 0.0 if xd > nx - 1.001: xd = nx - 1.001 if yd < 0.0: yd = 0.0 if yd > ny - 1.001: yd = ny - 1.001 x0 = int(xd); y0 = int(yd) if x0 > nx - 2: x0 = nx - 2 if y0 > ny - 2: y0 = ny - 2 fx = xd - x0; fy = yd - y0 v00 = F[c, y0, x0]; v01 = F[c, y0, x0 + 1] v10 = F[c, y0 + 1, x0]; v11 = F[c, y0 + 1, x0 + 1] if solid[y0, x0]: v00 = fc if solid[y0, x0 + 1]: v01 = fc if solid[y0 + 1, x0]: v10 = fc if solid[y0 + 1, x0 + 1]: v11 = fc val = (v00 * (1 - fx) * (1 - fy) + v01 * fx * (1 - fy) + v10 * (1 - fx) * fy + v11 * fx * fy) l = fc; r = fc; d = fc; u = fc if ix > 0 and not solid[iy, ix - 1]: l = F[c, iy, ix - 1] if ix < nx - 1 and not solid[iy, ix + 1]: r = F[c, iy, ix + 1] if iy > 0 and not solid[iy - 1, ix]: d = F[c, iy - 1, ix] if iy < ny - 1 and not solid[iy + 1, ix]: u = F[c, iy + 1, ix] out[c, iy, ix] = val + k * dt * (l + r + d + u - 4.0 * fc) return out # -------------------------------------------------------------------------- # 粒子轨迹 (拉格朗日) # -------------------------------------------------------------------------- @njit(cache=True, fastmath=True) def sample_uv(ux, uy, x, y): ny, nx = ux.shape if x < 0: x = 0.0 if x > nx - 1.001: x = nx - 1.001 if y < 0: y = 0.0 if y > ny - 1.001: y = ny - 1.001 x0 = int(x); y0 = int(y) if x0 > nx - 2: x0 = nx - 2 if y0 > ny - 2: y0 = ny - 2 fx = x - x0; fy = y - y0 u = (ux[y0, x0] * (1 - fx) * (1 - fy) + ux[y0, x0 + 1] * fx * (1 - fy) + ux[y0 + 1, x0] * (1 - fx) * fy + ux[y0 + 1, x0 + 1] * fx * fy) v = (uy[y0, x0] * (1 - fx) * (1 - fy) + uy[y0, x0 + 1] * fx * (1 - fy) + uy[y0 + 1, x0] * (1 - fx) * fy + uy[y0 + 1, x0 + 1] * fx * fy) return u, v @njit(cache=True, fastmath=True) def step_particles(px, py, alive, age, ux, uy, solid, vs): n = px.shape[0] ny, nx = solid.shape for i in range(n): if not alive[i]: continue x, y = px[i], py[i] u1, v1 = sample_uv(ux, uy, x, y) xm, ym = x + 0.5 * u1, y + 0.5 * (v1 - vs) u2, v2 = sample_uv(ux, uy, xm, ym) xn, yn = x + u2, y + (v2 - vs) ix, iy = int(xn), int(yn) if 0 <= ix < nx and 0 <= iy < ny and not solid[iy, ix]: px[i] = xn; py[i] = yn else: ix, iy = int(xn), int(y) if 0 <= ix < nx and 0 <= iy < ny and not solid[iy, ix]: px[i] = xn else: ix, iy = int(x), int(yn) if 0 <= ix < nx and 0 <= iy < ny and not solid[iy, ix]: py[i] = yn age[i] += 1.0 if py[i] < 1.5 or px[i] < 2.0 or px[i] > nx - 3.0 or age[i] > 2500: alive[i] = False # -------------------------------------------------------------------------- # 主循环 # -------------------------------------------------------------------------- def main(): out_dir = os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))), "docs", "assets") os.makedirs(out_dir, exist_ok=True) print(f"[FFU] grid {NX}x{NY} Re={RE:.0f} tau={TAU:.3f} steps={N_STEPS}") solid, inlet, outlet, src_equip, src_head = build_geometry() f = np.zeros((9, NY, NX)) ftmp = np.zeros_like(f) rho = np.ones((NY, NX)) ux = np.zeros((NY, NX)) uy = np.zeros((NY, NX)) uy[:] = -U0_LB for k in range(9): eu = EX[k] * 0.0 + EY[k] * (-U0_LB) f[k] = WT[k] * (1.0 + 3.0 * eu + 4.5 * eu * eu - 1.5 * U0_LB ** 2) # 标量: [T, RH, C] F = np.zeros((3, NY, NX)) kappa = np.array([KAPPA_T, KAPPA_RH, KAPPA_C]) vs = np.array([0.0, 0.0, V_SETTLE]) # 粒子 N_PART = 2400 px = np.zeros(N_PART); py = np.zeros(N_PART) alive = np.zeros(N_PART, dtype=np.bool_) age = np.zeros(N_PART) gen = np.zeros(N_PART) gen_counter = 0 next_p = 0 max_pts = 400000 traj = np.zeros((max_pts, 4)) n_pts = 0 src_px = np.concatenate([np.where(src_equip)[1], np.where(src_head)[1]]) src_py = np.concatenate([np.where(src_equip)[0], np.where(src_head)[0]]) rng = np.random.default_rng(7) for step in range(N_STEPS): lbm_collide(f, rho, ux, uy, solid, EX, EY, WT, TAU) lbm_bounce_back(f, solid, OPP) lbm_stream(f, ftmp, EX, EY) f, ftmp = ftmp, f lbm_bcs(f, ux, uy, solid, inlet, outlet, EX, EY, WT, U0_LB) if step >= SCALAR_FROM and step % ADV_EVERY == 0: F = advect_scalars(F, ux, uy, solid, kappa, DT_SCALAR, vs) F[0][src_equip] = 1.0 F[1][src_head] = 1.0 F[2][src_equip] = 1.0 F[2][src_head] = 1.0 F[:, NY - 1, :] = 0.0 F[:, solid] = 0.0 if step >= TRACK_FROM: if step % 40 == 0: for _ in range(18): i = next_p % N_PART next_p += 1 gen_counter += 1 j = rng.integers(0, src_px.size) px[i] = src_px[j] + rng.uniform(-1, 1) py[i] = src_py[j] + rng.uniform(-1, 1) alive[i] = True age[i] = 0.0 gen[i] = gen_counter step_particles(px, py, alive, age, ux, uy, solid, V_SETTLE) if step % 4 == 0: idx = np.where(alive)[0] m = min(idx.size, max_pts - n_pts) if m > 0: sel = idx[:m] traj[n_pts:n_pts + m, 0] = px[sel] traj[n_pts:n_pts + m, 1] = py[sel] traj[n_pts:n_pts + m, 2] = age[sel] traj[n_pts:n_pts + m, 3] = gen[sel] n_pts += m if step % 2000 == 0: umax = np.sqrt(ux ** 2 + uy ** 2).max() * SCALE_U print(f" step {step:6d} |U|max={umax:6.3f} m/s rho=[{rho.min():.4f},{rho.max():.4f}]") print(f"[FFU] simulation done, trajectory points: {n_pts}") traj = traj[:n_pts] np.savez_compressed(os.path.join(out_dir, "ffu_fields.npz"), ux=ux, uy=uy, T=F[0], RH=F[1], C=F[2], rho=rho, traj=traj) make_figures(out_dir, ux, uy, F, rho, solid, src_equip, src_head, traj) print("[FFU] figures written to", out_dir) def make_figures(out_dir, ux, uy, F, rho, solid, src_equip, src_head, traj): speed = np.sqrt(ux ** 2 + uy ** 2) * SCALE_U speed[solid] = np.nan T_phys = T_AMB + F[0] * (T_SRC - T_AMB) RH_phys = RH_AMB + F[1] * (RH_SRC - RH_AMB) C_norm = F[2].copy() p_pa = (rho - 1.0) / 3.0 * SCALE_P for _ in range(3): p_s = np.zeros_like(p_pa) for dy in (-1, 0, 1): for dx in (-1, 0, 1): p_s += np.roll(np.roll(p_pa, dy, 0), dx, 1) p_pa = p_s / 9.0 for arr in (T_phys, RH_phys, C_norm, p_pa): arr[solid] = np.nan plot_hero(out_dir, ux, uy, speed, solid) plot_scalars(out_dir, T_phys, RH_phys, C_norm, p_pa, solid, src_equip, src_head) plot_tracks(out_dir, speed, traj, solid) def replot(): """仅重绘图形, 不重跑仿真""" out_dir = os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))), "docs", "assets") d = np.load(os.path.join(out_dir, "ffu_fields.npz")) solid, inlet, outlet, src_equip, src_head = build_geometry() F = np.stack([d["T"], d["RH"], d["C"]]) make_figures(out_dir, d["ux"], d["uy"], F, d["rho"], solid, src_equip, src_head, d["traj"]) print("[FFU] replotted from ffu_fields.npz") # -------------------------------------------------------------------------- # 绘图 # -------------------------------------------------------------------------- CMAP_FLOW = LinearSegmentedColormap.from_list( "flow", ["#0a1626", "#123c63", "#1f6feb", "#39d0d8", "#7ee787", "#f0f6fc"]) def draw_room(ax, solid): """绘制房间轮廓与障碍物""" for (x0, y0, x1, y1) in (BENCH, EQUIP, BODY, RACK): ax.add_patch(Rectangle((x0, y0), x1 - x0, y1 - y0, facecolor="#2d333b", edgecolor="#adbac7", linewidth=1.1, zorder=6)) cx, cy, rr = HEAD_C ax.add_patch(Circle((cx, cy), rr, facecolor="#2d333b", edgecolor="#adbac7", linewidth=1.1, zorder=6)) # FFU 顶棚 ax.add_patch(Rectangle((0.06, H_M - 0.06), W_M - 0.12, 0.045, facecolor="#1f6feb", edgecolor="#58a6ff", linewidth=0.8, alpha=0.85, zorder=6)) # 回风口 for g0, g1 in (GRILL_L, GRILL_R): ax.add_patch(Rectangle((g0 + 0.02, 0.0), g1 - g0 - 0.04, 0.035, facecolor="#7ee787", edgecolor="#7ee787", alpha=0.65, zorder=6)) def draw_sources(ax): ax.plot(1.18, 1.25, marker="*", ms=13, color="#ff7b72", zorder=7) ax.plot(2.75, 1.62, marker="o", ms=8, color="#39d0d8", zorder=7) ax.plot(1.18, 1.35, marker="^", ms=7, color="#f0e68c", zorder=7) def decorate_hero(ax): """FFU 送风箭头 + 回风箭头 + 标注""" for x in np.linspace(0.45, W_M - 0.45, 9): ax.annotate("", xy=(x, H_M - 0.30), xytext=(x, H_M - 0.075), arrowprops=dict(arrowstyle="-|>", color="#58a6ff", lw=1.4, alpha=0.9), zorder=7) ax.text(W_M / 2, H_M - 0.145, "FFU 送风口 0.45 m/s（满布率 100%）", ha="center", va="bottom", fontsize=11.5, color="#58a6ff", zorder=8) for g0, g1 in (GRILL_L, GRILL_R): xm = (g0 + g1) / 2 ax.annotate("", xy=(xm, -0.02), xytext=(xm, 0.16), arrowprops=dict(arrowstyle="-|>", color="#7ee787", lw=1.6), annotation_clip=False, zorder=7) ax.text(xm, -0.09, "回风口", ha="center", va="top", fontsize=10, color="#7ee787", clip_on=False, zorder=8) ax.annotate("工作台", xy=(1.20, 0.92), xytext=(1.55, 1.62), fontsize=10.5, color="#e6edf3", zorder=8, arrowprops=dict(arrowstyle="-", color="#8b949e", lw=0.8)) ax.annotate("设备（热源+尘源）", xy=(1.18, 1.25), xytext=(0.42, 2.02), fontsize=10.5, color="#ff7b72", zorder=8, arrowprops=dict(arrowstyle="-", color="#8b949e", lw=0.8)) ax.annotate("操作人员（湿源+尘源）", xy=(2.58, 1.66), xytext=(2.62, 2.30), fontsize=10.5, color="#39d0d8", zorder=8, arrowprops=dict(arrowstyle="-", color="#8b949e", lw=0.8)) ax.annotate("物料架", xy=(3.25, 1.18), xytext=(3.10, 1.72), fontsize=10.5, color="#e6edf3", zorder=8, arrowprops=dict(arrowstyle="-", color="#8b949e", lw=0.8)) def plot_hero(out_dir, ux, uy, speed, solid): fig, ax = plt.subplots(figsize=(15.5, 9.2)) im = ax.imshow(speed, origin="lower", extent=[0, W_M, 0, H_M], cmap=CMAP_FLOW, vmin=0, vmax=0.75, interpolation="bilinear") xs = np.linspace(0.10, W_M - 0.10, 150) ys = np.linspace(0.08, H_M - 0.10, 110) XS, YS = np.meshgrid(xs, ys) ix = np.clip((XS * PPM).astype(int), 0, NX - 1) iy = np.clip((YS * PPM).astype(int), 0, NY - 1) UU = np.where(solid[iy, ix], np.nan, ux[iy, ix] * SCALE_U) VV = np.where(solid[iy, ix], np.nan, uy[iy, ix] * SCALE_U) with np.errstate(invalid="ignore"): ax.streamplot(xs, ys, UU, VV, density=1.5, color=(0.79, 0.82, 0.85, 0.5), linewidth=0.55, arrowsize=0.7, zorder=5) draw_room(ax, solid) draw_sources(ax) decorate_hero(ax) cb = fig.colorbar(im, ax=ax, fraction=0.032, pad=0.015) cb.set_label("速度大小 |U| (m/s)", fontsize=11) cb.outline.set_edgecolor("#30363d") ax.set_xlim(-0.01, W_M + 0.01) ax.set_ylim(-0.16, H_M + 0.01) ax.set_xlabel("x (m)", fontsize=11) ax.set_ylabel("y (m)", fontsize=11) ax.set_title("FFU 单向流洁净区流场仿真（D2Q9 LBM，Re≈540，网格 360×270）", fontsize=15, pad=14) ax.text(0.5, 1.012, "送风 0.45 m/s · 顶棚满布 · 底部两侧回风 · 流线为气流轨迹", transform=ax.transAxes, ha="center", va="bottom", fontsize=10.5, color=MUTED) fig.tight_layout() fig.savefig(os.path.join(out_dir, "ffu-flow-field.png"), dpi=155) fig.savefig(os.path.join(out_dir, "ffu-flow-field.svg")) plt.close(fig) print(" -> ffu-flow-field.png / .svg") def _panel(ax, field, solid, cmap, label, vmin=None, vmax=None, center=None): if center is not None: lim = np.nanmax(np.abs(field)) * 1.05 im = ax.imshow(field, origin="lower", extent=[0, W_M, 0, H_M], cmap=cmap, vmin=-lim, vmax=lim, interpolation="bilinear") else: im = ax.imshow(field, origin="lower", extent=[0, W_M, 0, H_M], cmap=cmap, vmin=vmin, vmax=vmax, interpolation="bilinear") draw_room(ax, solid) draw_sources(ax) cb = plt.colorbar(im, ax=ax, fraction=0.036, pad=0.02) cb.set_label(label, fontsize=10) cb.outline.set_edgecolor("#30363d") cb.ax.tick_params(labelsize=9) ax.set_xlim(-0.01, W_M + 0.01) ax.set_ylim(0, H_M) ax.set_xticks([0, 1, 2, 3, 4]) ax.set_xlabel("x (m)", fontsize=10) return im def plot_scalars(out_dir, T, RH, C, p, solid, src_equip, src_head): fig, axes = plt.subplots(2, 2, figsize=(15.5, 10.6)) ax = axes[0, 0] _panel(ax, T, solid, "inferno", "温度 T (°C)", vmin=T_AMB, vmax=T_SRC) ax.set_title("(a) 温度场 —— 设备热羽被单向流下扫", fontsize=12.5, color=FG, pad=8) ax = axes[0, 1] _panel(ax, RH, solid, "GnBu", "相对湿度 RH (%)", vmin=RH_AMB, vmax=RH_SRC) ax.set_title("(b) 湿度场 —— 呼吸区湿气随气流下送", fontsize=12.5, color=FG, pad=8) ax = axes[1, 0] _panel(ax, C, solid, "magma", "归一化粒子浓度 C/C₀", vmin=0, vmax=0.6) ax.set_title("(c) 尘埃粒子浓度 —— 尘源散发与输运", fontsize=12.5, color=FG, pad=8) ax = axes[1, 1] _panel(ax, p, solid, "RdBu_r", "静压偏差 Δp (Pa)", center=0) ax.set_title("(d) 压力场 —— FFU 送风正压与回风负压", fontsize=12.5, color=FG, pad=8) fig.suptitle("FFU 对下方温湿度、压力与尘埃粒子浓度的影响（稳态，D2Q9 LBM + 标量输运）", fontsize=15, y=0.985) fig.text(0.5, 0.012, "★ 热/尘源（设备） ● 湿/尘源（人员呼吸区） ▲ 尘源（设备上方） " "送风 0.45 m/s，设备表面 35 °C，呼吸区 RH 70%", ha="center", fontsize=10, color=MUTED) fig.tight_layout(rect=[0, 0.03, 1, 0.965]) fig.savefig(os.path.join(out_dir, "ffu-scalar-fields.png"), dpi=150) plt.close(fig) print(" -> ffu-scalar-fields.png") def plot_tracks(out_dir, speed, traj, solid): fig, ax = plt.subplots(figsize=(15.0, 8.9)) ax.imshow(speed, origin="lower", extent=[0, W_M, 0, H_M], cmap=CMAP_FLOW, vmin=0, vmax=0.75, interpolation="bilinear", alpha=0.9) if traj.shape[0] > 1: x = traj[:, 0] / PPM y = traj[:, 1] / PPM pid = traj[:, 3] age = traj[:, 2] / 60.0 m = pid[1:] == pid[:-1] segs = np.stack([np.column_stack([x[:-1][m], y[:-1][m]]), np.column_stack([x[1:][m], y[1:][m]])], axis=1) lc = LineCollection(segs, cmap="cool", linewidths=0.7, alpha=0.6) lc.set_array(age[1:][m]) ax.add_collection(lc) cb = plt.colorbar(lc, ax=ax, fraction=0.032, pad=0.015) cb.set_label("粒子在洁净区内的滞留时间 (s)", fontsize=10.5) cb.outline.set_edgecolor("#30363d") draw_room(ax, solid) draw_sources(ax) decorate_hero(ax) ax.set_xlim(-0.01, W_M + 0.01) ax.set_ylim(-0.16, H_M + 0.01) ax.set_xlabel("x (m)", fontsize=11) ax.set_ylabel("y (m)", fontsize=11) ax.set_title("尘埃粒子拉格朗日轨迹 —— 输运路径与滞留时间", fontsize=15, pad=14) ax.text(0.5, 1.012, "粒子释放自设备与人员源 · 含重力沉降 · 颜色越暖表示滞留越久（易沉积/死角）", transform=ax.transAxes, ha="center", va="bottom", fontsize=10.5, color=MUTED) fig.tight_layout() fig.savefig(os.path.join(out_dir, "ffu-particle-tracks.png"), dpi=150) plt.close(fig) print(" -> ffu-particle-tracks.png") if __name__ == "__main__": import sys if "--replot" in sys.argv: replot() else: main()