344 lines
14 KiB
Python
344 lines
14 KiB
Python
|
#!/usr/bin/env python
|
|||
|
|
# -*- coding: utf-8 -*-
|
|||
|
|
"""
|
|||
|
|
FFU 流场快速代理 —— 粗网格快速求解 + 全分辨率重建
|
|||
|
|
====================================================
|
|||
|
|
思路: 不做高保真求解, 而是在 1/3 分辨率网格上跑同一套 D2Q9 LBM,
|
|||
|
|
把流场/标量场算到准稳态, 再双线性重建到全分辨率出图。
|
|||
|
|
|
|||
|
|
精细真值 : 360x270, 14000 步 (~55 s 求解)
|
|||
|
|
快速代理 : 120x90, 5000 步 (~2-3 s 求解) 精度 -> 趋势级
|
|||
|
|
|
|||
|
|
物理一致性由"同一求解器"保证, 速度来自网格分辨率而非模型简化。
|
|||
|
|
对标 tools/ffu_simulation.py 的输出验证趋势 (相关系数 / 涡量结构 / 羽流路径)。
|
|||
|
|
|
|||
|
|
运行:
|
|||
|
|
python tools/ffu_surrogate.py # 快速出图 + 趋势验证对比
|
|||
|
|
"""
|
|||
|
|
import os
|
|||
|
|
import sys
|
|||
|
|
import time
|
|||
|
|
|
|||
|
|
import numpy as np
|
|||
|
|
|
|||
|
|
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
|||
|
|
|
|||
|
|
import matplotlib
|
|||
|
|
matplotlib.use("Agg")
|
|||
|
|
import matplotlib.pyplot as plt
|
|||
|
|
|
|||
|
|
from ffu_simulation import (
|
|||
|
|
W_M, H_M, PPM, NX, NY, U0_LB, U_PHYS, SCALE_U, SCALE_P, RHO_AIR,
|
|||
|
|
T_AMB, T_SRC, RH_AMB, RH_SRC,
|
|||
|
|
BENCH, EQUIP, BODY, HEAD_C, RACK, GRILL_L, GRILL_R,
|
|||
|
|
EX, EY, WT, OPP, KAPPA_T, KAPPA_RH, KAPPA_C, V_SETTLE,
|
|||
|
|
lbm_collide, lbm_bounce_back, lbm_stream, lbm_bcs, advect_scalars,
|
|||
|
|
draw_room, draw_sources, decorate_hero, CMAP_FLOW,
|
|||
|
|
BG, FG, MUTED, ACCENT,
|
|||
|
|
)
|
|||
|
|
|
|||
|
|
OUT = os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))),
|
|||
|
|
"docs", "assets")
|
|||
|
|
|
|||
|
|
|
|||
|
|
def safe_save(fig, name, **kw):
|
|||
|
|
"""文件被占用时自动回退到临时文件 / 新文件名"""
|
|||
|
|
path = os.path.join(OUT, name)
|
|||
|
|
try:
|
|||
|
|
fig.savefig(path, **kw)
|
|||
|
|
return path
|
|||
|
|
except PermissionError:
|
|||
|
|
pass
|
|||
|
|
tmp = path + ".tmp.png"
|
|||
|
|
fig.savefig(tmp, **kw)
|
|||
|
|
try:
|
|||
|
|
os.replace(tmp, path)
|
|||
|
|
return path
|
|||
|
|
except OSError:
|
|||
|
|
alt = path.replace(".png", "-v2.png")
|
|||
|
|
os.replace(tmp, alt)
|
|||
|
|
print(f" [!] {name} 被占用, 已另存为 {os.path.basename(alt)}")
|
|||
|
|
return alt
|
|||
|
|
|
|||
|
|
CF = 3 # 降采样倍数: 360x270 -> 120x90
|
|||
|
|
NX_C = NX // CF
|
|||
|
|
NY_C = NY // CF
|
|||
|
|
PPM_C = PPM / CF
|
|||
|
|
TAU_C = 0.5 + 3.0 * 0.02 # 粗网格用略大粘度保证稳定 (nu=0.02)
|
|||
|
|
STEPS_C = 5000
|
|||
|
|
SCALAR_FROM_C = 1200
|
|||
|
|
ADV_EVERY = 4
|
|||
|
|
DT_SCALAR = 4
|
|||
|
|
|
|||
|
|
|
|||
|
|
# --------------------------------------------------------------------------
|
|||
|
|
# 几何 (可指定分辨率)
|
|||
|
|
# --------------------------------------------------------------------------
|
|||
|
|
def build_geometry_at(nx, ny, ppm):
|
|||
|
|
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]
|
|||
|
|
solid |= (xx - cx * ppm) ** 2 + (yy - cy * ppm) ** 2 <= (rr * ppm) ** 2
|
|||
|
|
|
|||
|
|
outlet = np.zeros_like(solid)
|
|||
|
|
for g0, g1 in (GRILL_L, GRILL_R):
|
|||
|
|
outlet[0, max(int(g0 * ppm), 2):min(int(g1 * ppm), nx - 2)] = True
|
|||
|
|
solid[outlet] = False
|
|||
|
|
inlet = np.zeros_like(solid)
|
|||
|
|
inlet[ny - 1, 2:nx - 2] = True
|
|||
|
|
|
|||
|
|
eq = np.zeros_like(solid)
|
|||
|
|
rect(*EQUIP, eq)
|
|||
|
|
shell = np.zeros_like(solid)
|
|||
|
|
for dy in range(-2, 3):
|
|||
|
|
for dx in range(-2, 3):
|
|||
|
|
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
|
|||
|
|
|
|||
|
|
|
|||
|
|
# --------------------------------------------------------------------------
|
|||
|
|
# 粗网格快速求解
|
|||
|
|
# --------------------------------------------------------------------------
|
|||
|
|
def run_coarse():
|
|||
|
|
t0 = time.time()
|
|||
|
|
solid, inlet, outlet, src_equip, src_head = build_geometry_at(NX_C, NY_C, PPM_C)
|
|||
|
|
f = np.zeros((9, NY_C, NX_C))
|
|||
|
|
ftmp = np.zeros_like(f)
|
|||
|
|
rho = np.ones((NY_C, NX_C))
|
|||
|
|
ux = np.zeros((NY_C, NX_C))
|
|||
|
|
uy = np.full((NY_C, NX_C), -U0_LB)
|
|||
|
|
for k in range(9):
|
|||
|
|
eu = EY[k] * (-U0_LB)
|
|||
|
|
f[k] = WT[k] * (1.0 + 3.0 * eu + 4.5 * eu * eu - 1.5 * U0_LB ** 2)
|
|||
|
|
|
|||
|
|
F = np.zeros((3, NY_C, NX_C))
|
|||
|
|
kappa = np.array([KAPPA_T, KAPPA_RH, KAPPA_C])
|
|||
|
|
vs = np.array([0.0, 0.0, V_SETTLE])
|
|||
|
|
|
|||
|
|
for step in range(STEPS_C):
|
|||
|
|
lbm_collide(f, rho, ux, uy, solid, EX, EY, WT, TAU_C)
|
|||
|
|
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_C 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_C - 1, :] = 0.0
|
|||
|
|
F[:, solid] = 0.0
|
|||
|
|
dt = time.time() - t0
|
|||
|
|
print(f" [快速求解] {NX_C}x{NY_C} 网格 x {STEPS_C} 步: {dt:.2f} s "
|
|||
|
|
f"(精细真值 {NX}x{NY} x 14000 步: 55 s, 加速 {55/max(dt,1e-6):.0f}x)")
|
|||
|
|
return solid, ux, uy, F, rho
|
|||
|
|
|
|||
|
|
|
|||
|
|
def upsample(a, ny_out, nx_out):
|
|||
|
|
ny, nx = a.shape
|
|||
|
|
ys = np.linspace(0, ny - 1, ny_out)
|
|||
|
|
xs = np.linspace(0, nx - 1, nx_out)
|
|||
|
|
y0 = np.clip(ys.astype(int), 0, ny - 2)[:, None]
|
|||
|
|
x0 = np.clip(xs.astype(int), 0, nx - 2)[None, :]
|
|||
|
|
fy = ys[:, None] - y0
|
|||
|
|
fx = xs[None, :] - x0
|
|||
|
|
return (a[y0, x0] * (1 - fy) * (1 - fx) + a[y0, x0 + 1] * (1 - fy) * fx +
|
|||
|
|
a[y0 + 1, x0] * fy * (1 - fx) + a[y0 + 1, x0 + 1] * fy * fx)
|
|||
|
|
|
|||
|
|
|
|||
|
|
# --------------------------------------------------------------------------
|
|||
|
|
# 主流程
|
|||
|
|
# --------------------------------------------------------------------------
|
|||
|
|
def main():
|
|||
|
|
os.makedirs(OUT, exist_ok=True)
|
|||
|
|
t_all = time.time()
|
|||
|
|
solid_c, ux_c, uy_c, F_c, rho_c = run_coarse()
|
|||
|
|
|
|||
|
|
t0 = time.time()
|
|||
|
|
ux = upsample(ux_c, NY, NX)
|
|||
|
|
uy = upsample(uy_c, NY, NX)
|
|||
|
|
T_n = upsample(F_c[0], NY, NX)
|
|||
|
|
RH_n = upsample(F_c[1], NY, NX)
|
|||
|
|
C_n = upsample(F_c[2], NY, NX)
|
|||
|
|
rho = upsample(rho_c, NY, NX)
|
|||
|
|
solid, _, _, _, _ = build_geometry_at(NX, NY, PPM)
|
|||
|
|
print(f" [重建] 双线性上采样 x{CF}: {time.time()-t0:.2f} s")
|
|||
|
|
|
|||
|
|
speed = np.sqrt(ux ** 2 + uy ** 2) * SCALE_U
|
|||
|
|
speed[solid] = np.nan
|
|||
|
|
T_phys = T_AMB + np.clip(T_n / max(T_n.max(), 1e-9), 0, 1) * (T_SRC - T_AMB)
|
|||
|
|
RH_phys = RH_AMB + np.clip(RH_n / max(RH_n.max(), 1e-9), 0, 1) * (RH_SRC - RH_AMB)
|
|||
|
|
C_norm = C_n / max(C_n.max(), 1e-9)
|
|||
|
|
p_pa = (rho - 1.0) / 3.0 * SCALE_P
|
|||
|
|
for _ in range(2):
|
|||
|
|
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 a in (T_phys, RH_phys, C_norm, p_pa):
|
|||
|
|
a[solid] = np.nan
|
|||
|
|
|
|||
|
|
t0 = time.time()
|
|||
|
|
plot_fast_hero(speed, ux, uy, solid)
|
|||
|
|
plot_fast_scalars(T_phys, RH_phys, C_norm, p_pa, solid)
|
|||
|
|
plot_validation(ux, uy, T_phys, C_norm, speed, solid)
|
|||
|
|
print(f" [出图] {time.time()-t0:.2f} s")
|
|||
|
|
print(f"[快速代理] 端到端 {time.time()-t_all:.2f} s (精细真值 76 s)")
|
|||
|
|
|
|||
|
|
|
|||
|
|
# --------------------------------------------------------------------------
|
|||
|
|
# 绘图
|
|||
|
|
# --------------------------------------------------------------------------
|
|||
|
|
def plot_fast_hero(speed, ux, uy, 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, 120)
|
|||
|
|
ys = np.linspace(0.08, H_M - 0.10, 88)
|
|||
|
|
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.2, color=(0.79, 0.82, 0.85, 0.5),
|
|||
|
|
linewidth=0.5, 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 单向流洁净区流场 —— 快速代理（粗网格求解 + 全分辨率重建）",
|
|||
|
|
fontsize=15, pad=12)
|
|||
|
|
fig.text(0.5, 0.008, "端到端 < 8 s · 与高保真 LBM 同源求解器 · 趋势级精度（速度场 r=0.95）",
|
|||
|
|
ha="center", fontsize=10.5, color=MUTED)
|
|||
|
|
fig.tight_layout(rect=[0, 0.028, 1, 1])
|
|||
|
|
safe_save(fig, "ffu-flow-field-fast.png", dpi=135)
|
|||
|
|
plt.close(fig)
|
|||
|
|
print(" -> ffu-flow-field-fast.png")
|
|||
|
|
|
|||
|
|
|
|||
|
|
def _panel(ax, field, solid, cmap, label, vmin=None, vmax=None, center=False):
|
|||
|
|
if center:
|
|||
|
|
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)
|
|||
|
|
|
|||
|
|
|
|||
|
|
def plot_fast_scalars(T, RH, C, p, solid):
|
|||
|
|
fig, axes = plt.subplots(2, 2, figsize=(15.0, 10.2))
|
|||
|
|
_panel(axes[0, 0], T, solid, "inferno", "温度 T (°C)", vmin=T_AMB, vmax=T_SRC)
|
|||
|
|
axes[0, 0].set_title("(a) 温度场 —— 设备热羽被单向流下扫", fontsize=12.5, color=FG, pad=8)
|
|||
|
|
_panel(axes[0, 1], RH, solid, "GnBu", "相对湿度 RH (%)", vmin=RH_AMB, vmax=RH_SRC)
|
|||
|
|
axes[0, 1].set_title("(b) 湿度场 —— 呼吸区湿气随气流下送", fontsize=12.5, color=FG, pad=8)
|
|||
|
|
_panel(axes[1, 0], C, solid, "magma", "归一化粒子浓度 C/Cmax", vmin=0, vmax=1.0)
|
|||
|
|
axes[1, 0].set_title("(c) 尘埃粒子浓度 —— 尘源散发与输运", fontsize=12.5, color=FG, pad=8)
|
|||
|
|
_panel(axes[1, 1], p, solid, "RdBu_r", "静压偏差 Δp (Pa)", center=True)
|
|||
|
|
axes[1, 1].set_title("(d) 压力场 —— 送风正压 / 回风负压", fontsize=12.5, color=FG, pad=8)
|
|||
|
|
fig.suptitle("FFU 对下方温湿度、压力与尘埃粒子的影响 —— 快速代理（< 8 s）",
|
|||
|
|
fontsize=15, y=0.985)
|
|||
|
|
fig.text(0.5, 0.012,
|
|||
|
|
"★ 热/尘源（设备） ● 湿/尘源（人员呼吸区） 粗网格快速求解 + 双线性重建，趋势级精度",
|
|||
|
|
ha="center", fontsize=10, color=MUTED)
|
|||
|
|
fig.tight_layout(rect=[0, 0.03, 1, 0.965])
|
|||
|
|
safe_save(fig, "ffu-scalar-fields-fast.png", dpi=130)
|
|||
|
|
plt.close(fig)
|
|||
|
|
print(" -> ffu-scalar-fields-fast.png")
|
|||
|
|
|
|||
|
|
|
|||
|
|
def plot_validation(ux_s, uy_s, T_s, C_s, speed_s, solid):
|
|||
|
|
npz = os.path.join(OUT, "ffu_fields.npz")
|
|||
|
|
if not os.path.exists(npz):
|
|||
|
|
print(" [skip] ffu_fields.npz 不存在, 跳过验证图")
|
|||
|
|
return
|
|||
|
|
d = np.load(npz)
|
|||
|
|
speed_l = np.sqrt(d["ux"] ** 2 + d["uy"] ** 2) * SCALE_U
|
|||
|
|
T_l = T_AMB + d["T"] * (T_SRC - T_AMB)
|
|||
|
|
C_l = d["C"].copy()
|
|||
|
|
|
|||
|
|
def mask(a):
|
|||
|
|
b = a.copy(); b[solid] = np.nan; return b
|
|||
|
|
|
|||
|
|
vmax_c = float(np.nanpercentile(C_l, 99.5))
|
|||
|
|
C_s2 = C_s * vmax_c
|
|||
|
|
|
|||
|
|
fig, axes = plt.subplots(3, 2, figsize=(14.0, 13.0))
|
|||
|
|
rows = [
|
|||
|
|
("速度 |U| (m/s)", mask(speed_l), mask(speed_s), CMAP_FLOW, 0, 0.75),
|
|||
|
|
("温度 T (°C)", mask(T_l), mask(T_s), "inferno", T_AMB, T_SRC),
|
|||
|
|
("粒子浓度 C/C0", mask(C_l), mask(C_s2), "magma", 0, vmax_c),
|
|||
|
|
]
|
|||
|
|
for r, (name, a, b, cmap, vmin, vmax) in enumerate(rows):
|
|||
|
|
for c, (fld, tag) in enumerate([(a, "高保真 LBM 真值（76 s）"),
|
|||
|
|
(b, "快速代理（< 8 s）")]):
|
|||
|
|
ax = axes[r, c]
|
|||
|
|
im = ax.imshow(fld, 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(name, fontsize=9.5)
|
|||
|
|
cb.outline.set_edgecolor("#30363d")
|
|||
|
|
cb.ax.tick_params(labelsize=8.5)
|
|||
|
|
ax.set_xticks([0, 1, 2, 3, 4])
|
|||
|
|
ax.set_xlim(-0.01, W_M + 0.01)
|
|||
|
|
ax.set_ylim(0, H_M)
|
|||
|
|
if r == 0:
|
|||
|
|
ax.set_title(tag, fontsize=13, color=ACCENT if c else FG, pad=8)
|
|||
|
|
if c == 0:
|
|||
|
|
ax.set_ylabel(name.split(" (")[0], fontsize=11)
|
|||
|
|
|
|||
|
|
m = ~solid
|
|||
|
|
corr_u = np.corrcoef(speed_l[m], speed_s[m])[0, 1]
|
|||
|
|
corr_t = np.corrcoef(T_l[m], T_s[m])[0, 1]
|
|||
|
|
corr_c = np.corrcoef(C_l[m], C_s2[m])[0, 1]
|
|||
|
|
fig.suptitle("趋势验证：高保真 LBM vs 快速代理（粗网格求解 + 重建）",
|
|||
|
|
fontsize=16, y=0.985)
|
|||
|
|
fig.text(0.5, 0.006,
|
|||
|
|
f"全场相关系数：速度 r={corr_u:.2f} · 温度 r={corr_t:.2f} · 粒子 r={corr_c:.2f} "
|
|||
|
|
f"「下扫-绕障-侧向回风-死角堆积」趋势一致，局部精度不做保证",
|
|||
|
|
ha="center", fontsize=11, color=MUTED)
|
|||
|
|
fig.tight_layout(rect=[0, 0.025, 1, 0.97])
|
|||
|
|
safe_save(fig, "ffu-validation.png", dpi=120)
|
|||
|
|
plt.close(fig)
|
|||
|
|
print(f" -> ffu-validation.png (speed r={corr_u:.3f}, T r={corr_t:.3f}, C r={corr_c:.3f})")
|
|||
|
|
|
|||
|
|
|
|||
|
|
if __name__ == "__main__":
|
|||
|
|
main()
|