实际气体状态方程：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/contentcnr/100766.html

七、使用说明

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

八、应用场景

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

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