admin
5b16b9e790
- 3D D3Q19 LBM (numba.cuda), 240x130x265=8.27M 单元 - 恒功率风机曲线 dp*Q=P 压力进口, 工况点涌现 - 多孔布风板: 2mm 厚, 6.25% 开孔, 5 度倾斜(-x) - 4 段 DN300 侧进风 / 567->1200 渐扩 / DN800 排风 - 输出 19 张切面/矢量/压差/流量/分区图 (4x 分辨率)
508 lines
19 KiB
Python
508 lines
19 KiB
Python
#!/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()
|