#!/usr/bin/env python # -*- coding: utf-8 -*- """ 水平卧式流化床 空床流场 —— RANS + 多孔跳变布风板 (交错网格投影法 + AMG) ====================================================================== 交错网格 (MAC): p 在单元中心, u/v/w 在面心 -> 从原理上消除棋盘压力振荡。 预测: u* = u + dt[ -(u.grad)u + div(nu_eff grad u) + S_porous ] 压力: div(grad p) = (rho/dt) div(u*) (pyamg 多重网格, 固体 Neumann, 出口 Dirichlet) 修正: u = u* - (dt/rho) grad(p) 布风板 = 薄多孔层 (Darcy), 阻力标定 dP=2900 Pa @ 表观气速 2.04 m/s。 进风 = 4 x DN300 速度进口 (2500 m3/h/段); 排风 = 压力出口; 其余无滑移壁。 运行: python tools/fbed_rans.py """ import os import time import numpy as np import numba from scipy.sparse import csr_matrix import pyamg # -------------------------------------------------------------------------- BED_L, BED_W = 2.4, 0.567 N_SEC = 4 SEC_L = BED_L / N_SEC PLENUM_H = 0.45 PLATE_T = 0.05 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 Z0 = -PLENUM_H Z1 = FB_H + UPPER_H + REDUCER_H + DUCT_H PPM = int(os.environ.get("RANS_PPM", "20")) NX = int(round((X1 - X0) * PPM)) NY = int(round((Y1 - Y0) * PPM)) NZ = int(round((Z1 - Z0) * PPM)) DX = (X1 - X0) / NX DY = (Y1 - Y0) / NY DZ = (Z1 - Z0) / NZ XC = X0 + (np.arange(NX) + 0.5) * DX YC = Y0 + (np.arange(NY) + 0.5) * DY ZC = Z0 + (np.arange(NZ) + 0.5) * DZ Q_SEC = 2500.0 / 3600.0 A_PIPE = np.pi * PIPE_R ** 2 U_IN = Q_SEC / A_PIPE A_BED = BED_L * BED_W U_SUP = Q_SEC * N_SEC / A_BED DP_RATED = float(os.environ.get("RANS_DP", "2900.0")) # 额定点压力 (2500 m3/h/段) Q_RATED = Q_SEC * N_SEC # 额定总流量 m3/s P_FAN = DP_RATED * Q_RATED # 恒功率 (W) DP_SHUT = DP_RATED * float(os.environ.get("RANS_FAN_SHUT", "1.5")) def fan_dp(q_m3s): """恒功率风机曲线: dp * Q = P -> dp = P/Q (Q->0 时用关死压力截断)""" if q_m3s < 1e-6: return DP_SHUT return min(P_FAN / q_m3s, DP_SHUT) RHO = 1.2 NU = 1.5e-5 R_POROUS = DP_RATED / (U_SUP * PLATE_T) NU_T = float(os.environ.get("RANS_NUT", "0.02")) NU_EFF = NU + NU_T MAX_STEPS = int(os.environ.get("RANS_STEPS", "4000")) TIME_BUDGET = float(os.environ.get("RANS_TB", "16.0")) CFL = float(os.environ.get("RANS_CFL", "0.3")) def build_masks(): solid = np.ones((NZ, NY, NX), dtype=bool) yw = BED_W / 2.0 def box(xa, xb, ya, yb, za, zb, val=False): ix0, ix1 = int(round((xa - X0) / DX)), int(round((xb - X0) / DX)) iy0, iy1 = int(round((ya - Y0) / DY)), int(round((yb - Y0) / DY)) iz0, iz1 = int(round((za - Z0) / DZ)), int(round((zb - Z0) / DZ)) solid[max(iz0, 0):min(iz1, NZ), max(iy0, 0):min(iy1, NY), max(ix0, 0):min(ix1, NX)] = val box(0.0, BED_L, -yw, yw, -PLENUM_H, 0.0) for k in range(1, N_SEC): xk = k * SEC_L box(xk - 0.0125, xk + 0.0125, -yw, yw, -PLENUM_H, 0.0, True) box(0.0, BED_L, -yw, yw, 0.0, PLATE_T) for iz in range(NZ): z = ZC[iz] if z <= PLATE_T or z >= FB_H: continue w = yw + (FB_W / 2.0 - yw) * (z / FB_H) solid[iz] = np.where(np.abs(YC[:, None]) <= w, False, solid[iz]) box(0.0, BED_L, -FB_W / 2.0, FB_W / 2.0, FB_H, FB_H + UPPER_H) z_r0, z_r1 = FB_H + UPPER_H, FB_H + UPPER_H + REDUCER_H for iz in range(NZ): z = ZC[iz] if z < z_r0 or z >= z_r1: continue t = (z - 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]) ** 2 <= EXH_R ** 2 for iz in range(NZ): if ZC[iz] >= z_r1: solid[iz] = np.where(cyl, False, solid[iz]) iy_p0 = int(round((PIPE_Y0 - Y0) / DY)) iy_p1 = int(round((-yw - Y0) / DY)) for k in range(N_SEC): xk = SEC_L * (k + 0.5) zk = -PLENUM_H / 2.0 for iz in range(NZ): if abs(ZC[iz] - zk) > PIPE_R: continue m = (XC - xk) ** 2 + (ZC[iz] - zk) ** 2 <= PIPE_R ** 2 solid[iz, iy_p0:iy_p1, m] = False plate = np.zeros_like(solid) iz0, iz1 = int(round((0.0 - Z0) / DZ)), int(round((PLATE_T - Z0) / DZ)) iy0, iy1 = int(round((-yw - Y0) / DY)), int(round((yw - Y0) / DY)) plate[iz0:iz1, iy0:iy1, :] = ~solid[iz0:iz1, iy0:iy1, :] inlet = np.zeros_like(solid) inlet[:, iy_p0:iy_p0 + 2, :] = ~solid[:, iy_p0:iy_p0 + 2, :] outlet = np.zeros_like(solid) outlet[NZ - 2:, :, :] = ~solid[NZ - 2:, :, :] return solid, plate, inlet, outlet # -------------------------------------------------------------------------- # 预测步 (交错网格) # -------------------------------------------------------------------------- @numba.njit(cache=True, fastmath=True) def predict(u, v, w, solid, plate, inlet, outlet, nu, dt, dx, dy, dz, rporous, uin): nz, ny, nx = solid.shape # u faces: (nz, ny, nx+1), u[iz,iy,ix] at west face of cell ix for iz in range(nz): for iy in range(ny): for ix in range(1, nx): if solid[iz, iy, ix] and solid[iz, iy, ix - 1]: u[iz, iy, ix] = 0.0 continue if solid[iz, iy, ix] or solid[iz, iy, ix - 1]: u[iz, iy, ix] = 0.0 continue uc = u[iz, iy, ix] ue = u[iz, iy, ix + 1] if ix + 1 <= nx else uc uw = u[iz, iy, ix - 1] un = u[iz, iy + 1, ix] if iy + 1 < ny else uc us = u[iz, iy - 1, ix] if iy - 1 >= 0 else uc ut = u[iz + 1, iy, ix] if iz + 1 < nz else uc ub = u[iz - 1, iy, ix] if iz - 1 >= 0 else uc vn = 0.5 * (v[iz, iy + 1, ix - 1] + v[iz, iy + 1, ix]) if iy + 1 < ny else 0.0 vs = 0.5 * (v[iz, iy, ix - 1] + v[iz, iy, ix]) wt = 0.5 * (w[iz + 1, iy, ix - 1] + w[iz + 1, iy, ix]) if iz + 1 < nz else 0.0 wb = 0.5 * (w[iz, iy, ix - 1] + w[iz, iy, ix]) dux = (uc - uw) / dx if uc > 0 else (ue - uc) / dx duy = (un - uc) / dy if vn > 0 else (uc - us) / dy duz = (ut - uc) / dz if wt > 0 else (uc - ub) / dz conv = uc * dux + 0.5 * (vn + vs) * duy + 0.5 * (wt + wb) * duz lap = (ue - 2.0 * uc + uw) / (dx * dx) + (un - 2.0 * uc + us) / (dy * dy) \ + (ut - 2.0 * uc + ub) / (dz * dz) unew = uc + dt * (-conv + nu * lap) if plate[iz, iy, ix] or plate[iz, iy, ix - 1]: unew /= (1.0 + dt * rporous) u[iz, iy, ix] = unew # v faces: (nz, ny+1, nx) for iz in range(nz): for iy in range(1, ny): for ix in range(nx): if solid[iz, iy, ix] or solid[iz, iy - 1, ix]: v[iz, iy, ix] = 0.0 continue vc = v[iz, iy, ix] ve = v[iz, iy, ix + 1] if ix + 1 < nx else vc vw = v[iz, iy, ix - 1] if ix - 1 >= 0 else vc vn = v[iz, iy + 1, ix] if iy + 1 <= ny else vc vs = v[iz, iy - 1, ix] vt = v[iz + 1, iy, ix] if iz + 1 < nz else vc vb = v[iz - 1, iy, ix] if iz - 1 >= 0 else vc ue = 0.5 * (u[iz, iy, ix + 1] + u[iz, iy - 1, ix + 1]) if ix + 1 <= nx else 0.0 uw = 0.5 * (u[iz, iy, ix] + u[iz, iy - 1, ix]) wt = 0.5 * (w[iz + 1, iy, ix] + w[iz + 1, iy - 1, ix]) if iz + 1 < nz else 0.0 wb = 0.5 * (w[iz, iy, ix] + w[iz, iy - 1, ix]) dvy = (vc - vs) / dy if vc > 0 else (vn - vc) / dy dvx = (ve - vc) / dx if ue > 0 else (vc - vw) / dx dvz = (vt - vc) / dz if wt > 0 else (vc - vb) / dz conv = 0.5 * (ue + uw) * dvx + vc * dvy + 0.5 * (wt + wb) * dvz lap = (ve - 2.0 * vc + vw) / (dx * dx) + (vn - 2.0 * vc + vs) / (dy * dy) \ + (vt - 2.0 * vc + vb) / (dz * dz) vnew = vc + dt * (-conv + nu * lap) if plate[iz, iy, ix] or plate[iz, iy - 1, ix]: vnew /= (1.0 + dt * rporous) v[iz, iy, ix] = vnew # w faces: (nz+1, ny, nx) for iz in range(1, nz): for iy in range(ny): for ix in range(nx): if solid[iz, iy, ix] or solid[iz - 1, iy, ix]: w[iz, iy, ix] = 0.0 continue wc = w[iz, iy, ix] we = w[iz, iy, ix + 1] if ix + 1 < nx else wc ww = w[iz, iy, ix - 1] if ix - 1 >= 0 else wc wn = w[iz, iy + 1, ix] if iy + 1 < ny else wc ws = w[iz, iy - 1, ix] if iy - 1 >= 0 else wc wt = w[iz + 1, iy, ix] if iz + 1 <= nz else wc wb = w[iz - 1, iy, ix] ue = 0.5 * (u[iz, iy, ix + 1] + u[iz - 1, iy, ix + 1]) if ix + 1 <= nx else 0.0 uw = 0.5 * (u[iz, iy, ix] + u[iz - 1, iy, ix]) vn = 0.5 * (v[iz, iy + 1, ix] + v[iz - 1, iy + 1, ix]) if iy + 1 <= ny else 0.0 vs = 0.5 * (v[iz, iy, ix] + v[iz - 1, iy, ix]) dwz = (wc - wb) / dz if wc > 0 else (wt - wc) / dz dwx = (we - wc) / dx if ue > 0 else (wc - ww) / dx dwy = (wn - wc) / dy if vn > 0 else (wc - ws) / dy conv = 0.5 * (ue + uw) * dwx + 0.5 * (vn + vs) * dwy + wc * dwz lap = (we - 2.0 * wc + ww) / (dx * dx) + (wn - 2.0 * wc + ws) / (dy * dy) \ + (wt - 2.0 * wc + wb) / (dz * dz) wnew = wc + dt * (-conv + nu * lap) if plate[iz, iy, ix] or plate[iz - 1, iy, ix]: wnew /= (1.0 + dt * rporous) w[iz, iy, ix] = wnew return @numba.njit(cache=True, fastmath=True) def apply_bc(u, v, w, solid, inlet, outlet, uin): nz, ny, nx = solid.shape # 压力进口: 边界面速度零梯度 (让气流自由进入) for iz in range(nz): for ix in range(nx): if inlet[iz, 0, ix]: v[iz, 0, ix] = v[iz, 1, ix] return @numba.njit(cache=True, fastmath=True) def divergence(u, v, w, solid, dx, dy, dz): nz, ny, nx = solid.shape d = np.zeros((nz, ny, nx)) for iz in range(nz): for iy in range(ny): for ix in range(nx): if solid[iz, iy, ix]: continue d[iz, iy, ix] = (u[iz, iy, ix + 1] - u[iz, iy, ix]) / dx + \ (v[iz, iy + 1, ix] - v[iz, iy, ix]) / dy + \ (w[iz + 1, iy, ix] - w[iz, iy, ix]) / dz return d @numba.njit(cache=True, fastmath=True) def correct(u, v, w, p, solid, dt, rho, dx, dy, dz): nz, ny, nx = solid.shape c = dt / rho for iz in range(nz): for iy in range(ny): for ix in range(1, nx): if solid[iz, iy, ix] or solid[iz, iy, ix - 1]: continue u[iz, iy, ix] -= c * (p[iz, iy, ix] - p[iz, iy, ix - 1]) / dx for iz in range(nz): for iy in range(1, ny): for ix in range(nx): if solid[iz, iy, ix] or solid[iz, iy - 1, ix]: continue v[iz, iy, ix] -= c * (p[iz, iy, ix] - p[iz, iy - 1, ix]) / dy for iz in range(1, nz): for iy in range(ny): for ix in range(nx): if solid[iz, iy, ix] or solid[iz - 1, iy, ix]: continue w[iz, iy, ix] -= c * (p[iz, iy, ix] - p[iz - 1, iy, ix]) / dz return def build_poisson(solid, outlet, inlet): n = NX * NY * NZ idx = np.arange(n).reshape(NZ, NY, NX) fluid = ~solid rows, cols, vals = [], [], [] def add(r, c, v): rows.append(np.asarray(r).ravel()); cols.append(np.asarray(c).ravel()) vals.append(np.asarray(v).ravel()) diag = np.zeros((NZ, NY, NX)) m = fluid[:, :, :-1] & fluid[:, :, 1:] c = 1.0 / DX ** 2 add(idx[:, :, :-1][m], idx[:, :, 1:][m], np.full(m.sum(), -c)) add(idx[:, :, 1:][m], idx[:, :, :-1][m], np.full(m.sum(), -c)) diag[:, :, :-1] += m * c; diag[:, :, 1:] += m * c m = fluid[:, :-1, :] & fluid[:, 1:, :] c = 1.0 / DY ** 2 add(idx[:, :-1, :][m], idx[:, 1:, :][m], np.full(m.sum(), -c)) add(idx[:, 1:, :][m], idx[:, :-1, :][m], np.full(m.sum(), -c)) diag[:, :-1, :] += m * c; diag[:, 1:, :] += m * c m = fluid[:-1, :, :] & fluid[1:, :, :] c = 1.0 / DZ ** 2 add(idx[:-1, :, :][m], idx[1:, :, :][m], np.full(m.sum(), -c)) add(idx[1:, :, :][m], idx[:-1, :, :][m], np.full(m.sum(), -c)) diag[:-1, :, :] += m * c; diag[1:, :, :] += m * c fl = fluid.ravel() df = diag.ravel().copy() df = np.where(fl & (df == 0.0), 1.0, df) df[~fl] = 1.0 add(np.arange(n), np.arange(n), df) A = csr_matrix((np.concatenate(vals), (np.concatenate(rows), np.concatenate(cols))), shape=(n, n)) of = outlet.ravel() inf = inlet.ravel() fix = of | inf if fix.any(): oi = np.where(fix)[0] A = A.tolil(); A[oi, :] = 0.0; A[oi, oi] = 1.0; A = A.tocsr() return A, idx def main(): t_all = time.time() print(f"[rans] staggered grid {NX}x{NY}x{NZ}={NX*NY*NZ} PPM={PPM}") solid, plate, inlet, outlet = build_masks() print(f"[rans] fluid={int((~solid).sum())} plate={int(plate.sum())} R_porous={R_POROUS:.0f}") A, idx = build_poisson(solid, outlet, inlet) print(f"[rans] poisson nnz={A.nnz}, AMG ...") t0 = time.time() ml = pyamg.smoothed_aggregation_solver(A, max_coarse=300) print(f"[rans] AMG {time.time()-t0:.2f}s levels={len(ml.levels)}") u = np.zeros((NZ, NY, NX + 1)) v = np.zeros((NZ, NY + 1, NX)) w = np.zeros((NZ + 1, NY, NX)) p = np.zeros((NZ, NY, NX)) dt = CFL * min(DX, DY, DZ) / (U_IN + 1e-9) print(f"[rans] 恒功率风机: P={P_FAN:.0f} W 额定({Q_RATED*3600:.0f} m3/h, {DP_RATED:.0f} Pa) " f"关死={DP_SHUT:.0f} Pa nu_eff={NU_EFF:.4f} dt={dt:.3e}") p_fan = DP_RATED t0 = time.time() for it in range(MAX_STEPS): predict(u, v, w, solid, plate, inlet, outlet, NU_EFF, dt, DX, DY, DZ, R_POROUS, U_IN) d = divergence(u, v, w, solid, DX, DY, DZ) rhs = (RHO / dt) * d rhs[outlet] = 0.0 rhs[inlet] = p_fan sol = ml.solve(rhs.ravel(), x0=p.ravel(), tol=1e-8, maxiter=200) p = sol.reshape(NZ, NY, NX) correct(u, v, w, p, solid, dt, RHO, DX, DY, DZ) apply_bc(u, v, w, solid, inlet, outlet, U_IN) # 测进口流量 -> 更新风机压力 iz_in = int(round((-0.45 - Y0) / DY)) m_in = ~solid[:, iz_in, :] q_in = abs(v[:, iz_in, :][m_in].sum()) * DX * DZ p_fan = fan_dp(q_in) if it % 100 == 0: uc = 0.5 * (u[:, :, :-1] + u[:, :, 1:]) vc = 0.5 * (v[:, :-1, :] + v[:, 1:, :]) wc = 0.5 * (w[:-1, :, :] + w[1:, :, :]) sp = np.sqrt(uc ** 2 + vc ** 2 + wc ** 2) print(f" it={it:5d} t={time.time()-t0:.1f}s |U|max={sp.max():.2f} " f"Q_in={q_in*3600:.0f} p_fan={p_fan:.0f} p=[{p.min():.0f},{p.max():.0f}]") if time.time() - t0 > TIME_BUDGET: break izp = int(round((PLATE_T - Z0) / DZ)) q_plate = w[izp][~solid[izp]].sum() * DX * DY * 3600.0 iz_in = int(round((-0.45 - Y0) / DY)) m_in = ~solid[:, iz_in, :] q_in = v[:, iz_in, :][m_in].sum() * DX * DZ * 3600.0 iz_out = NZ - 2 q_out = w[iz_out][~solid[iz_out]].sum() * DX * DY * 3600.0 a0, a1 = int(round((-0.30 - Z0) / DZ)), int(round((-0.05 - Z0) / DZ)) p_pl = p[a0:a1][~solid[a0:a1]].mean() b0, b1 = int(round((0.15 - Z0) / DZ)), int(round((0.35 - Z0) / DZ)) p_fb = p[b0:b1][~solid[b0:b1]].mean() print(f"[rans] steps={it+1} solve={time.time()-t0:.1f}s") print(f"[rans] plate Q={q_plate:.0f} inlet Q={q_in:.0f} outlet Q={q_out:.0f} m3/h") print(f"[rans] p_风箱={p_pl:.0f} p_沸腾段={p_fb:.0f} Δp={p_pl-p_fb:+.0f} Pa") out = os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))), "docs", "assets") np.savez_compressed(os.path.join(out, "fbed_rans.npz"), u=u, v=v, w=w, p=p, solid=solid, plate=plate) print(f"[rans] total {time.time()-t_all:.1f}s") if __name__ == "__main__": main()