admin/Physics-Ai-Fluent-Analysis
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
main
Physics-Ai-Fluent-Analysis/tools/ffu_simulation.py
T

649 lines
26 KiB
Python
Raw Permalink Normal View History

初始化项目：规划文档 + FFU 洁净室流场实时仿真 Demo
2026-09-18 00:05:18 +08:00
#!/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()
Reference in New Issue Copy Permalink