实际气体状态方程：Peng-Robinson（P-R）方程计算指南

一、P-R方程基本原理

Peng-Robinson（P-R）方程是立方型状态方程的代表之一，由Peng和Robinson于1976年提出，用于描述实际气体的压力（ $P$ ）、体积（ $V$ ）、温度（ $T$ ）关系，考虑了分子间作用力（通过参数 $a$ ）和分子固有体积（通过参数 $b$ ），比理想气体方程（ $PV=nRT$ ）更准确。

二、P-R方程数学表达式

对于任意物质的量 $n$ 的气体，P-R方程的摩尔体积形式为：
$P = \frac{RT}{V_m - b} - \frac{a(T)}{V_m(V_m + b) + b(V_m - b)}$

其中：

$P$ ：气体压力（Pa）， $V_m = V/n$ ：摩尔体积（m³/mol）， $T$ ：热力学温度（K）；
$R$ ：通用气体常数（ $R=8.314\ \text{J/(mol·K)}$ 或 $0.08206\ \text{L·atm/(mol·K)}$ ）；
$a(T)$ ：温度依赖的分子间吸引力参数（Pa·m⁶/mol²）；
$b$ ：分子固有体积参数（m³/mol）。

三、参数 $a(T)$ 与 $b$ 的计算

参数 $a(T)$ 和 $b$ 由气体的临界性质（ $T_c$ ：临界温度， $P_c$ ：临界压力）和偏心因子（ $\omega$ ）决定：

1. 参数 $b$ （分子固有体积）

$b = 0.07780 \cdot \frac{RT_c}{P_c}$

2. 参数 $a(T)$ （温度依赖的吸引力参数）

$a(T) = a_c \cdot \alpha(T)$
其中：

临界吸引参数 $a_c$ ： $a_c = 0.45724 \cdot \frac{R^2 T_c^2}{P_c}$ ；
温度修正因子 $\alpha(T)$ ： $\alpha(T) = \left[1 + \kappa\left(1 - \sqrt{T_r}\right)\right]^2$ ；
对比温度 $T_r = T/T_c$ （ $T$ 为当前温度， $T_c$ 为临界温度）；
系数 $\kappa$ ： $\kappa = 0.37464 + 1.54226\omega - 0.26992\omega^2$ （ $\omega$ 为偏心因子）。

四、P-R方程求解步骤

实际应用中，常需根据已知条件（ $T,P$ 求 $V_m$ ，或 $T,V_m$ 求 $P$ ）解方程。由于P-R方程是关于 $V_m$ 的三次多项式，需通过代数或数值方法求解实根。

步骤1：整理为三次方程标准形式

将P-R方程展开并整理为关于 $V_m$ 的三次方程：
$V_m^3 - \left(b + \frac{RT}{P}\right)V_m^2 + \left(a(T) - 3b^2 - \frac{2RTb}{P}\right)V_m - \left(ab - b^2 - \frac{RTb^2}{P}\right) = 0$

令：
$A = \frac{a(T)P}{R^2 T^2}, \quad B = \frac{bP}{RT}$
则方程简化为无量纲形式（压缩因子 $Z = PV_m/(RT)$ ）：
$Z^3 - (1 - B)Z^2 + (A - 3B^2 - 2B)Z - (AB - B^2 - B^3) = 0$

步骤2：求解三次方程实根

三次方程可能有1个或3个实根，物理意义如下：

1个实根：超临界流体或单一气相/液相；
3个实根：气液共存区（最大根为气相 $V_m^g$ ，最小根为液相 $V_m^l$ ，中间根不稳定）。

求解方法：

解析法：卡尔达诺公式（复杂，较少用）；
数值法：牛顿迭代法（实用，需给定初值）。

五、计算实例（以甲烷为例）

已知条件：甲烷（ $CH_4$ ）的临界参数 $T_c=190.56\ \text{K}$ ， $P_c=4.5992\ \text{MPa}=4.5992\times10^6\ \text{Pa}$ ，偏心因子 $\omega=0.011$ ；当前状态 $T=300\ \text{K}$ ， $P=1\ \text{MPa}=1\times10^6\ \text{Pa}$ 。求摩尔体积 $V_m$ 。

1. 计算参数 $b$

$b = 0.07780 \cdot \frac{RT_c}{P_c} = 0.07780 \cdot \frac{8.314 \times 190.56}{4.5992\times10^6} \approx 2.69\times10^{-5}\ \text{m}^3/\text{mol}$

2. 计算参数 $a(T)$

对比温度 $T_r = T/T_c = 300/190.56 \approx 1.574$ ；
系数 $\kappa = 0.37464 + 1.54226\times0.011 - 0.26992\times0.011^2 \approx 0.391$ ；
温度修正因子 $\alpha(T) = \left[1 + 0.391(1 - \sqrt{1.574})\right]^2 \approx 0.839$ ；
临界吸引参数 $a_c = 0.45724 \cdot \frac{(8.314)^2 \times (190.56)^2}{4.5992\times10^6} \approx 2.283\times10^{-1}\ \text{Pa·m}^6/\text{mol}^2$ ；
$a(T) = a_c \cdot \alpha(T) \approx 2.283\times10^{-1} \times 0.839 \approx 0.191\ \text{Pa·m}^6/\text{mol}^2$ （注：此处单位需注意， $a_c$ 实际计算应为 $2.283\times10^5\ \text{Pa·m}^6/\text{mol}^2$ ，修正后 $a(T) \approx 1.91\times10^5\ \text{Pa·m}^6/\text{mol}^2$ ）。

3. 代入P-R方程求 $V_m$

将已知参数代入P-R方程：
$1\times10^6 = \frac{8.314\times300}{V_m - 2.69\times10^{-5}} - \frac{1.91\times10^5}{V_m(V_m + 2.69\times10^{-5}) + 2.69\times10^{-5}(V_m - 2.69\times10^{-5})}$

通过牛顿迭代法求解（初值取理想气体体积 $V_m^0 = RT/P \approx 2.494\times10^{-3}\ \text{m}^3/\text{mol}$ ），最终得 $V_m \approx 2.52\times10^{-3}\ \text{m}^3/\text{mol}$ （气相）。

六、Matlab代码实现

以下代码实现P-R方程计算，支持已知 $T,P$ 求 $V_m$ 或已知 $T,V_m$ 求 $P$ ：

function result = peng_robinson(T, P, Tc, Pc, omega, mode, Vm)
    % P-R方程计算函数
    % 输入:
    %   T: 温度(K), P: 压力(Pa), Tc: 临界温度(K), Pc: 临界压力(Pa), omega: 偏心因子
    %   mode: 'V' (已知T,P求Vm) 或 'P' (已知T,Vm求P)
    %   Vm: 摩尔体积(m³/mol) (mode='P'时需输入)
    % 输出:
    %   result: 摩尔体积(m³/mol) 或 压力(Pa)
    
    R = 8.314; % 气体常数 J/(mol·K)
    
    % 1. 计算参数a(T)和b
    Tr = T / Tc; % 对比温度
    kappa = 0.37464 + 1.54226*omega - 0.26992*omega^2; % 系数κ
    alpha = (1 + kappa*(1 - sqrt(Tr)))^2; % 温度修正因子
    ac = 0.45724 * (R^2 * Tc^2) / Pc; % 临界吸引参数a_c
    a = ac * alpha; % 温度依赖参数a(T)
    b = 0.07780 * (R * Tc) / Pc; % 分子体积参数b
    
    if strcmp(mode, 'V') % 已知T,P求Vm
        % 定义三次方程 f(Vm) = Vm³ + c2*Vm² + c1*Vm + c0 = 0
        c2 = -(b + R*T/P);
        c1 = a - 3*b^2 - 2*R*T*b/P;
        c0 = -(a*b - b^2 - R*T*b^2/P);
        
        % 牛顿迭代法求解（初值取理想气体体积）
        Vm0 = R*T/P; % 理想气体体积
        tol = 1e-6; max_iter = 100;
        for iter = 1:max_iter
            f = Vm0^3 + c2*Vm0^2 + c1*Vm0 + c0;
            df = 3*Vm0^2 + 2*c2*Vm0 + c1;
            Vm1 = Vm0 - f/df;
            if abs(Vm1 - Vm0) < tol, break; end
            Vm0 = Vm1;
        end
        result = Vm1;
        
    elseif strcmp(mode, 'P') % 已知T,Vm求P
        term1 = R*T / (Vm - b);
        denominator = Vm*(Vm + b) + b*(Vm - b);
        term2 = a / denominator;
        result = term1 - term2;
    else
        error('模式错误！请选择''V''(求Vm)或''P''(求P)');
    end
end

参考代码用于计算实际气体状态方程P-R方程 www.youwenfan.com/contentcsr/100766.html

七、使用说明

单位统一：压力 $P$ 用Pa，温度 $T$ 用K，体积 $V_m$ 用m³/mol，临界参数 $T_c$ 、 $P_c$ 单位对应；
常见气体参数：可通过NIST数据库或化工手册查询（如甲烷： $T_c=190.56\ \text{K}$ ， $P_c=4.5992\ \text{MPa}$ ， $\omega=0.011$ ）；
多根选择：气液共存区需根据实际状态（气相/液相）选择对应实根（气相取最大根，液相取最小根）。

八、应用场景

P-R方程广泛用于化工热力学计算，如：

气体压缩因子 $Z$ 计算；
气液平衡（VLE）预测；
压缩机/换热器设计中实际气体物性估算。

实际气体状态方程：Peng-Robinson（P-R）方程计算指南