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

508 lines
19 KiB
Python
Raw Permalink Normal View History

feat: 水平卧式流化床空床流场 GPU LBM 仿真 (100 格/m)
2026-09-19 13:37:51 +08:00
#!/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()
Reference in New Issue Copy Permalink