Skip to content
cfd-lab:~/zh/posts/2026-08-20-franklin-stov…online
NOTE #136DAY THU 유체역학DATE 2026.08.20READ 5 min read#Stack-Effect#Natural-Convection#Buoyancy#Historical#Flow-Phenomena

让烟先下降1 m,抽力就少了7 Pa — 富兰克林壁炉与浮力预算

浮力只能按当地密度逐米买入,热的时候下降的那一米,永远换不回冷端补上的那一米。

1741年,一台故意让烟往下走的壁炉#

本杰明·富兰克林在1741年重新设计了壁炉。目的很简单:把本该直接从烟囱溜走的热量,在屋里再抓一次。 当时欧洲人口激增、木柴短缺,取暖效率几乎是生存问题。

他的办法是把流道折起来。火上升的烟先越过挡板(baffle,横在流道里让流动转向的板),沿铸铁壁 向下走一段,然后才进入烟囱。下降途中铸铁被加热,热量辐射进房间。当时的记载把这个把戏称作 "虹吸现象"。

烟囱顶端的高度没有变,烟最终获得的高差也没有变。可是这台炉子却以往屋里倒烟而出名。本文把这笔 损失出在哪里、有多大,逐条记进浮力账本。结论只有一句 — 浮力只能在气体还热的那个高度上赚到。

抽力是密度差乘以高度

烟囱把烟往上拉的力叫抽力(draft)。它不是泵,而是静压差。外面的空气柱和里面的烟柱重量不同, 差额就是抽力。

对高度为 HH 的烟道(flue),抽力写成

Δpdraft=0H(ρρ(z))gdz\Delta p_{\text{draft}} = \int_0^{H} \left( \rho_\infty - \rho(z) \right) g \, \mathrm{d}z

其中 ρ\rho_\infty 是环境空气密度,ρ(z)\rho(z) 是高度 zz 处的烟气密度,gg 是重力加速度。烟道内外 绝对压力几乎相同,理想气体给出 ρ=p0/(RT)\rho = p_0 / (R T);若烟气温度均匀为 TgT_g,积分立刻可以求出。

Δpdraft=ρgH(1TTg)\Delta p_{\text{draft}} = \rho_\infty g H \left( 1 - \frac{T_\infty}{T_g} \right)

代入 10 °C 的外界空气(ρ=1.247\rho_\infty = 1.247 kg/m³)、烟温 500 K、烟囱高 7 m,得到 37 Pa。 真实壁炉抽力在 10–30 Pa 之间,量级是对的。

关键在于 dz\mathrm{d}z 在积分号里面。浮力不是拿总高度乘一次的数,而是按米、按当地密度一段段 买来的。富兰克林付出的代价就在这里。

在下面的模拟里亲手试一试。

The chimney top stays at 7 m however you drag down-leg, so the elevation gained never changes — but the red bar grows and the net draft falls, because the descent is paid for with the hottest, lightest gas. Now drag fire tempdown: the mouth gauge crosses the amber spill limit far earlier with a deep down-leg than with none.

down-leg 滑块推到底,烟囱顶端仍然钉在 7 m。可是右侧账本的红条在变长,净抽力在下降。接着把 fire temp 往下拉:炉口仪表会越过橙色的倒烟限值,而下降段越深,这一刻来得越早。

浮力只能在气体还热的高度上赚到

把路径切成三段:升到挡板(+0.8+0.8 m)、下降段(hd-h_d)、烟囱(+6.2+hd+6.2 + h_d m)。三段高程变化 相加,不论 hdh_d 取多少都恒等于 7 m。可是每一段的浮力贡献,要用它通过时的烟气密度来算。

烟气沿途向壁面散热而变冷。沿路径坐标 ss 写出来就是一个一阶常微分方程。

dTds=UPm˙cp(TT)\frac{\mathrm{d}T}{\mathrm{d}s} = -\frac{U P}{\dot{m} c_p} \left( T - T_\infty \right)

UU 是总传热系数,PP 是烟道周长,m˙\dot{m} 是质量流量,cpc_p 是定压比热。解是指数衰减,衰减 长度为 m˙cp/(UP)\dot{m} c_p / (U P),这一点稍后会起作用。

现在损失看得见了。下降段是整条路径上第二热的区段。密度接近最低,(ρρ)(\rho_\infty - \rho) 接近 最大,而它要乘上一个负的 dz\mathrm{d}z。反过来,把这一米补回来的烟囱段,用的是已经把热交给铸铁、 冷下来的烟气。这不是高卖低买,而是高买低卖。

更糟的是,在富兰克林的设计里,下降段就在房间内部,也就是散热最快的区段。为提高取暖效率而 设计的那个特征,恰好先花掉了烟囱要用的浮力预算。

用一个回路闭合动量平衡

流量 m˙\dot{m} 由浮力与损失相等的地方决定。这是绕烟道一圈的动量平衡。

iLi(ρρ)gdz=0LfDm˙22ρA2ds+Km˙22ρeA2\sum_{i} \int_{L_i} (\rho_\infty - \rho) g \, \mathrm{d}z = \int_0^{L} \frac{f}{D} \frac{\dot{m}^2}{2 \rho A^2} \, \mathrm{d}s + K \frac{\dot{m}^2}{2 \rho_e A^2}

左边是浮力水头,右边是摩擦损失加局部损失(弯头与出口)。ff 是 Darcy 摩擦系数,DD 是烟道直径, AA 是截面积,KK 是局部损失系数,ρe\rho_e 是出口密度。管摩擦项的形式与 入口段与 Hagen–Poiseuille 中用的一样。

这个式子没有看上去那么老实。m˙\dot{m} 变大,右边按 m˙2\dot{m}^2 增长,但左边也跟着增长:流量越大 停留时间越短,烟气到烟囱时还更热。由于这个反馈,残差 ΔpΔploss\Delta p - \Delta p_{\text{loss}} 关于 m˙\dot{m} 是上凸的,会出现两个根。下面那个根不稳定, 真正的工作点是上面那个根

用 Python 写出的浮力账本#

把路径细分后推进温度,分段收集浮力水头,再用二分法求上根。把下降高度从 0 扫到 1.2 m。

import math
 
G, R_AIR, CP = 9.81, 287.0, 1100.0
T_AMB, P_ATM = 283.0, 101325.0
D_FLUE = 0.15
A_FLUE = math.pi * D_FLUE ** 2 / 4.0
PERIM = math.pi * D_FLUE
F_DARCY, K_MINOR = 0.030, 3.0
U_ROOM, U_STACK = 40.0, 5.0        # 室内铸铁流道 / 烟囱砌体 [W/m2K]
Z_TOP, Z_BAFFLE = 7.0, 0.80        # 烟囱出口 / 挡板顶端 [m]
A_MOUTH, V_SPILL = 0.12, 0.25      # 炉口开口面积 [m2] / 倒烟临界面速度 [m/s]
 
 
def gas_density(T):
    return P_ATM / (R_AIR * T)
 
 
def march_flue(mdot, h_down, T_fire, n=120):
    # 沿上升 / 下降 / 烟囱三段推进,收集分段浮力水头与总损失
    legs = [(Z_BAFFLE, +1.0, U_ROOM),
            (h_down, -1.0, U_ROOM),
            (Z_TOP - Z_BAFFLE + h_down, +1.0, U_STACK)]
    T, loss, heads = T_fire, 0.0, []
    rho_a = gas_density(T_AMB)
    for L, sgn, U in legs:
        head = 0.0
        if L > 0.0:
            ds = L / n
            for _ in range(n):
                T_new = T_AMB + (T - T_AMB) * math.exp(-U * PERIM * ds / (mdot * CP))
                rho = gas_density(0.5 * (T + T_new))
                head += (rho_a - rho) * G * sgn * ds
                loss += F_DARCY / D_FLUE * 0.5 * mdot ** 2 / (rho * A_FLUE ** 2) * ds
                T = T_new
        heads.append(head)
    loss += K_MINOR * 0.5 * mdot ** 2 / (gas_density(T) * A_FLUE ** 2)
    return heads, loss, T
 
 
def solve_mdot(h_down, T_fire, lo=1e-3, hi=0.4, n=100):
    # 浮力 = 损失的上(稳定)根;无解时返回 0.0
    grid = [lo * (hi / lo) ** (i / (n - 1.0)) for i in range(n)]
    res = []
    for m in grid:
        heads, loss, _ = march_flue(m, h_down, T_fire)
        res.append(sum(heads) - loss)
    k = max(range(n), key=lambda i: res[i])
    if res[k] <= 0.0:
        return 0.0
    a, b = grid[k], grid[-1]
    for _ in range(40):
        m = 0.5 * (a + b)
        heads, loss, _ = march_flue(m, h_down, T_fire)
        a, b = (m, b) if sum(heads) - loss > 0.0 else (a, m)
    return 0.5 * (a + b)
 
 
print("h_down[m]  mdot[kg/s]  V_exit[m/s]  draft[Pa]  T_exit[K]")
for i in range(0, 7):
    h = 0.2 * i
    m = solve_mdot(h, 700.0)
    heads, loss, T_e = march_flue(m, h, 700.0)
    print(f"  {h:4.2f}     {m:8.4f}     {m / (gas_density(T_e) * A_FLUE):6.2f}"
          f"     {sum(heads):7.2f}    {T_e:7.1f}")
 
for h in (0.0, 1.0):
    m = solve_mdot(h, 700.0)
    heads, loss, T_e = march_flue(m, h, 700.0)
    print(f"\nbuoyancy ledger  h_down = {h:.1f} m, mdot = {m:.4f} kg/s")
    for name, val in zip(["rise to baffle", "descent      ", "chimney      "], heads):
        print(f"  {name}  {val:+7.2f} Pa")
    print(f"  net             {sum(heads):+7.2f} Pa   (loss {loss:.2f} Pa)")
 
mdot_spill = gas_density(T_AMB) * A_MOUTH * V_SPILL
print(f"\nspillage threshold: mdot < {mdot_spill:.4f} kg/s")
for h in (0.0, 1.0):
    lo, hi = 290.0, 700.0
    for _ in range(30):
        mid = 0.5 * (lo + hi)
        lo, hi = (mid, hi) if solve_mdot(h, mid) < mdot_spill else (lo, mid)
    print(f"  h_down = {h:.1f} m -> spills below T_fire = {0.5 * (lo + hi):6.1f} K")
h_down[m]  mdot[kg/s]  V_exit[m/s]  draft[Pa]  T_exit[K]
  0.00       0.0629       5.59       44.72      554.5
  0.20       0.0623       5.37       43.48      537.2
  0.40       0.0617       5.15       42.16      520.6
  0.60       0.0610       4.93       40.76      504.6
  0.80       0.0602       4.72       39.29      489.2
  1.00       0.0593       4.51       37.72      474.2
  1.20       0.0583       4.30       36.05      459.6
 
buoyancy ledger  h_down = 0.0 m, mdot = 0.0629 kg/s
  rise to baffle    +5.57 Pa
  descent          +0.00 Pa
  chimney         +39.15 Pa
  net              +44.72 Pa   (loss 44.72 Pa)
 
buoyancy ledger  h_down = 1.0 m, mdot = 0.0593 kg/s
  rise to baffle    +5.56 Pa
  descent          -6.16 Pa
  chimney         +38.32 Pa
  net              +37.72 Pa   (loss 37.72 Pa)
 
spillage threshold: mdot < 0.0374 kg/s
  h_down = 0.0 m -> spills below T_fire =  335.2 K
  h_down = 1.0 m -> spills below T_fire =  377.0 K

只要比较账本的两行就够了。1 m 的下降段直接拿走 6.16-6.16 Pa。而烟囱虽然长了 1 m,贡献却从 +39.15+39.15 降到 +38.32+38.32 Pa,因为进入它的烟气已经更冷了。两笔损失合起来是 6.99 Pa,与净抽力 下降的 7.00 Pa 精确吻合。

先撑不住的是炉口 — 42 K 的余量没了#

抽力从 44.7 降到 37.7 Pa,本身还不至于让烟进屋。真正的限值出在炉口开口。一旦穿过火口的面速度 低于约 0.25 m/s,热烟就开始从开口上沿漏出来。把这个判据换算成流量,是 0.0374 kg/s。

输出的最后两行把这个阈值换算成了温度。没有下降段时,烟温要掉到 335 K 才失守。有 1 m 下降段时, 377 K 就已经撑不住。42 K 的余量凭空消失了。

火烧得旺的时候,两种设计都抽得很好。麻烦出在点火那一刻和火将熄的那一刻。富兰克林壁炉被记成 "会冒烟的炉子",原因就在这里。最终在 1780 年代,戴维·里滕豪斯把下降段去掉,让烟道以 L 形直接 上到烟囱,这个形式后来成了标准。

中和面 — 同一笔账在一个房间里重演

同样的静压逻辑,在没有烟囱的房间里照样成立。内外密度不同,两条静压直线的斜率就不同,而斜率不同 的两条直线恰好只在一个高度相交。那个高度就是中和面(neutral plane)znz_n

Δp(z)=poutpin=(ρρi)g(znz)\Delta p(z) = p_{\text{out}} - p_{\text{in}} = (\rho_\infty - \rho_i) g \left( z_n - z \right)

若上下各有一个开口,面积为 A1,A2A_1, A_2,高度为 z1,z2z_1, z_2,由两个开口质量流量相等的条件,znz_n 有 闭式解。

zn=ρA12z1+ρiA22z2ρA12+ρiA22z_n = \frac{\rho_\infty A_1^2 z_1 + \rho_i A_2^2 z_2}{\rho_\infty A_1^2 + \rho_i A_2^2}

中和面以下外面的空气进来,以上室内空气出去。在下图里改变开口高度、面积比和室内温度。

Enlarge the upper opening and the amber neutral plane climbs toward it — the big hole is the one that sits near zero pressure difference. Then drag inside T below 283 K: the blue line tips the other way, both arrows reverse, and the room that was drawing at the floor starts drawing at the ceiling.

把上开口放大,中和面会被拉向它,因为大孔用更小的压差就能通过同样的流量。然后把室内温度降到 283 K 以下:蓝线朝另一边倾斜,两个箭头同时反向。

在自然对流计算里,这笔预算记在哪一栏

把这套计算搬进 CFD 时,有三处容易出错。

第一是压力边界条件的参考高度。在开口上直接压一个 totalPressure = 0,等于把中和面钉死在 边界面上。实际上要减去 ρgz\rho_\infty g z、用动压来提法,外侧气柱的静压梯度才能保留下来。忽略 这一点,解会照样收敛,但流量能差好几倍。

第二是Boussinesq 近似的适用范围。上面算例的温差超过 200 K,密度变化超过 50 %,此时按 ρρ0(1βΔT)\rho \approx \rho_0 (1 - \beta \Delta T) 线性化会把浮力算得很离谱。这条界线在哪里, Boussinesq 在 30 K 就开始撒谎的地方 里已经梳理过。烟囱问题应该走低马赫可压缩或变密度的提法。

第三是定常解不止一个。回路平衡有两个根这件事,会原样出现在迭代里。把初始流量给到接近零, 求解器就会滑向下根,收敛到"没有抽力"。自然对流问题里普遍存在的这种多解性质,与 自然对流与 Rayleigh 数 中讨论的分岔 结构同源。把初值放在上根附近,或者先固定流量再慢慢松开,都更稳妥。

富兰克林确实把热量留在了房间里。只不过那份热量,本来是烟囱要花的预算。在浮力驱动的系统里决定 换热器放在哪一段之前,先看清楚钱是从账本的哪一栏支出去的,会省很多事。

如果对您有帮助,请分享。