402 lines
16 KiB
Python
402 lines
16 KiB
Python
|
#!/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()
|