理想气体开始说谎的地方 — 压缩因子 Z 与 van der Waals 状态方程

速度相同时，动量通量按密度比的平方走。26%的密度误差，到ρu²里大约变成60%。

100个大气压、200 K的氮气若按理想气体算，密度会偏离实际将近四分之一。压力相同，温度也相同。分子彼此拉扯，还占着自己的体积。

这份偏差就是压缩因子 $Z$ 。van der Waals 用两个常数模仿气体，临界点以下，S形环把液体和气体分开。

把对比温度滑条拉到1以下，并打开 Maxwell 勾选，环会长大，同时出现一条上下面积相等的水平饱和线。

reduced temperature T_r = 0.900Maxwell construction

saturation pressure P_r,sat = 0.647
Z(sat. vapor) = 0.633  ·   Z(sat. liquid) = 0.163

0.85附近，饱和液体的 $Z$ 掉到0.1以下。把滑条贴到1，液体和气体的 $Z$ 收到同一个点。区分消失的地方就在那里。

100 个大气压的氮气比理想气体重#

$Pv = RT$ 把分子当成点，丢掉分子之间的力。稀薄时它是对的。分子离得远。

压力一升高，分子间距变窄。弱吸引力（van der Waals 力）把体积缩小。与此同时，分子自己的体积还在。

两种效应相反。吸引力缩小体积，有限体积放大体积。哪一边赢，由温度和压力决定。所以实际气体比理想气体重还是轻，会随条件翻转。这个密度排序一旦翻转，底部较轻的流体点燃对流的门槛也会一起挪动。

压缩因子 Z — 把偏差收进一个数#

把偏差收成一个数，就是压缩因子。

Z = \frac{Pv}{RT} = \frac{v_\text{actual}}{v_\text{ideal}}

$P$ 是压力， $v$ 是摩尔体积， $R$ 是气体常数， $T$ 是温度。 $Z=1$ 就是理想气体。离1越远，偏得越多。

$Z<1$ 表示吸引力占上风。实际体积比理想值小。 $Z>1$ 表示有限体积占上风。分子把体积顶大。氮气例子里 $Z \approx 0.79$ ，意思是实际密度约为理想值的 $1/0.79 \approx 1.27$ 倍。

对应态原理:所有气体彼此相似的坐标

每种气体的 $Z$ 曲线不一样。把压力和温度各自除以临界值，曲线就几乎叠上。

P_R = \frac{P}{P_c}, \qquad T_R = \frac{T}{T_c}

$P_c$ 、 $T_c$ 是临界压力和临界温度。对应态原理（principle of corresponding states）说的是：对比压力 $P_R$ 和对比温度 $T_R$ 相同， $Z$ 也几乎相同。

不论氮气、甲烷还是二氧化碳，在 $T_R=1.5$ 、 $P_R=2$ 处 $Z$ 都相近。决定行为的不是分子的名字，而是相对临界点的位置。所以一张通用图能盖住很多气体。

代码里的氮气站在另一个格子上。 $T = 200$ ， $T_c = 126.2$ ，所以 $T_R = 200/126.2 \approx 1.585$ 。 $P = 10^7$ ， $P_c = 3.39\times 10^6$ ，所以 $P_R = 10/3.39 \approx 2.95$ 。把打印的 $Z = 0.7905$ 换回对比体积， $v_R = Z\,(T/T_c)\,(P_c/P)\,(8/3) \approx 1.13$ 。把这个点放回对比等温线， $P_R$ 又回到2.95。打印出来的 $Z$ 就在那条超临界等温线上。

van der Waals:用两个常数模拟分子#

1873年，van der Waals 给理想气体加了两项修正。

\left(P + \frac{a}{v^2}\right)(v - b) = RT

$a/v^2$ 是吸引力削掉的压力。体积越小它越大。 $b$ 是分子自己占的体积，自由空间变成 $v-b$ 。

两个常数在临界点上定下来。那里 $P$ – $v$ 等温线是拐点，斜率和曲率同时为0。

\left(\frac{\partial P}{\partial v}\right)_{T_c} = 0, \qquad \left(\frac{\partial^2 P}{\partial v^2}\right)_{T_c} = 0

解出来就是常数和临界压缩因子。

a = \frac{27 R^2 T_c^2}{64 P_c}, \qquad b = \frac{R T_c}{8 P_c}, \qquad Z_c = \frac{P_c v_c}{R T_c} = \frac{3}{8}

$Z_c = 3/8 = 0.375$ 与气体种类无关。实际气体的 $Z_c$ 在0.27到0.29。0.375比0.28附近大约高出三分之一。定量上偏高，环和临界点的图像是对的。

写成对比变量， $a$ 、 $b$ 就消失了。

P_R = \frac{8 T_R}{3 v_R - 1} - \frac{3}{v_R^2}

这一行就是对应态原理的式子。临界点固定在 $(v_R, P_R, T_R) = (1, 1, 1)$ 。分母 $3v_R - 1$ 为零的地方是 $v_R = 1/3$ 。左边比分子体积 $b$ 还小，等温线不往那里画。

临界点之下的 S 形曲线与 Maxwell 作图#

$T_R > 1$ 时等温线单调下降。压力升高，体积减小。

$T_R < 1$ 时出现 S 形环。中间一段 $(\partial P/\partial v)_T > 0$ 。压力升高体积反而变大，力学上不稳定。实际物质会跳过这段。

取而代之的是一条水平线，把液体和气体接起来。高度就是饱和压力。Maxwell 等面积作图（equal-area rule）选取环被切出的上下面积相等的那个高度。这和液体、气体的 Gibbs 自由能相等是一回事。两相能够共存的压力只有这一个。

用 Python 求解 Z#

把 van der Waals 对 $Z$ 整理，得到三次方程。有了临界常数，任意 $P$ 、 $T$ 都能解出 $Z$ 。

import numpy as np
 
R = 8.314  # J/(mol·K)
 
def vdw_constants(Tc, Pc):
    """由临界常数求 van der Waals a, b。"""
    a = 27.0 * R**2 * Tc**2 / (64.0 * Pc)
    b = R * Tc / (8.0 * Pc)
    return a, b
 
def vdw_compressibility(P, T, Tc, Pc):
    """Z^3 - (1+B)Z^2 + A Z - A B = 0 的物理根。"""
    a, b = vdw_constants(Tc, Pc)
    A = a * P / (R * T)**2          # 吸引项的无量纲大小
    B = b * P / (R * T)            # 有限体积项的无量纲大小
    coeffs = [1.0, -(1.0 + B), A, -A * B]
    roots = np.roots(coeffs)
    real = roots[np.abs(roots.imag) < 1e-8].real
    real = real[real > B]          # 只有 v > b  ⇒  Z > B 才是物理的
    return real.max(), real.min(), real.size
 
# 氮气: Tc = 126.2 K, Pc = 3.39 MPa, M = 28 g/mol
Tc, Pc, M = 126.2, 3.39e6, 0.028
P, T = 1.0e7, 200.0                # 100 bar, 200 K
Z_gas, Z_liq, n = vdw_compressibility(P, T, Tc, Pc)
 
rho_ideal = P * M / (R * T)        # 理想气体密度
rho_real = rho_ideal / Z_gas       # 实际密度
print(f"根的个数 = {n},  Z = {Z_gas:.4f}")
print(f"理想气体 密度 = {rho_ideal:6.1f} kg/m^3")
print(f"实际    密度 = {rho_real:6.1f} kg/m^3  ({100*(rho_real/rho_ideal-1):+.0f}%)")

跑出来是三行。

根的个数 = 1,  Z = 0.7905
理想气体 密度 =  168.4 kg/m^3
实际    密度 =  213.0 kg/m^3  (+27%)

$213.0/168.4 = 1.265$ 。增幅是26.5%，输出里的 $+27\%$ 是把这个数进到整数。 $T_R \approx 1.59$ ，所以是超临界。实根只有一个，没有液体和气体之分。即便如此，理想气体密度仍比实际大约轻四分之一。

超临界喷注与 CFD — 状态方程为什么会撼动网格#

液体火箭燃烧室通常跑在推进剂临界压力之上。LOX（液氧）或深冷氮气喷进这个腔室，不会形成熟悉的液滴。

表面张力消失了。过了临界点，液体和气体的边界垮掉。射流不碎成液滴，密度梯度像梳子齿一样散开，和周围掺混。这就是超临界（或跨临界）喷注。

状态方程撼动网格，是因为密度和压力、温度非线性地绑在一起。求解器用理想气体，射流密度会错几十个百分点，动量通量 $\rho u^2$ 和混合长度一起偏掉。实际代码用 cubic EOS（Peng–Robinson、SRK）或 NASG（Noble–Abel Stiffened Gas）。 $Z$ 就是这个模型在每个格子里求解的量。

把密度误差在动量上也写成27%，就会失败#

速度相同， $\rho u^2$ 按密度比的平方走。 $1.265^2 \approx 1.60$ ，所以26.5%的密度误差在动量通量里大约是60%。写成“密度错了27%，所以通量也错27%”，是把平方丢掉了。

若把那个压力原样放进力里，Maxwell 选出的水平线上的界面还会静止吗？把压力直接放进去、静止界面就抖起来的非理想格子玻尔兹曼留下的问题是：这份力要不要改写成自由能，界面才会安静下来。