基于差分导数的超临界燃料裂解反应模拟研究

Study on simulation for supercritical fuel cracking reaction flow based on differential derivative

摘要

准确捕捉燃料反应热与动力学特征对先进发动机热管理技术的发展至关重要。现有超临界裂解反应流建模方法通常需要对特定的实际气体状态方程进行复杂解析，且缺乏可扩展性。本文建立一种基于差分导数的超临界反应流模拟方法，实现了多种立方型实际气体状态方程的通用化耦合求解。进一步以正癸烷为代表性燃料，对比研究了理想气体状态方程和4种实际气体状态方程对高压裂解过程的热物性和反应行为的预测差异。研究结果表明，理想气体模型在高压条件下对密度和黏度的预测误差较大，尤其低温区域更为显著；密度和黏度的低估进一步导致反应速率计算的偏差。PRVT和RKPR方程分别引入体积转化和Z_c拟合参数，显著提高了密度、黏度和温度的预测精度。在出口温度的预测中，4种立方型状态方程与实验数据的相对误差均不超过0.3%，而理想气体模型的误差为2.18%；在实际气体模型对裂解转化率的预测中，PR模型与实验数据的相对误差最大，为10.96%；SRK模型的相对误差最小，为6.66%；而理想气体模型的相对误差则高达20.41%。本文为超临界反应流模拟提供了新方法，并可为状态方程的合理选取提供参考。

Abstract

Accurate capture of fuel reaction heat and kinetic characteristics is crucial for the development of thermal management technologies in advanced engines. Existing modeling methods for supercritical pyrolysis reaction flows require complex analytical solutions for specific real-gas equations of state （EOS） and lack scalability. In this study， a differential derivative-based simulation method for supercritical reactive flows is established， enabling universal coupling with multiple cubic real-gas EOS. Using n-decane as a representative fuel， comparative investigations were conducted to evaluate the differences in predicting thermophysical properties and reaction behavior during high-pressure pyrolysis between the ideal gas EOS and four real-gas EOS （PR， SRK， PRVT， and RKPR）. The results show that the ideal gas model exhibits significant prediction errors in density and viscosity under high-pressure conditions， particularly at low temperatures. The underestimation of density and viscosity leads to deviations in reaction rate calculations. The PRVT and RKPR equations significantly improve prediction accuracy for density， viscosity， and temperature through the introduction of volume translation and Z_c fitting parameters， respectively. In outlet temperature predictions， the four cubic EOS show relative errors not exceeding 0.3% compared with experimental data， while the ideal gas model yields a 2.18% error. The results indicate that among real gas models predicting cracking conversion rates， the PR model exhibits the highest relative error of 10.96%， while the SRK model shows the lowest relative error of 6.66%. The ideal gas model， in contrast， demonstrates a significantly higher relative error of 20.41%. This work provides a novel methodology for supercritical reactive flow simulations and offers practical guidance for EOS selection.

Graphical abstract

关键词

碳氢燃料 / 裂解反应 / 超临界 / 反应热力学 / 建模方法

Key words

hydrocarbon fuels / cracking reaction / supercritical / reaction thermodynamics / modeling method

碳氢燃料裂解反应是石油化工及能源领域中的重要过程^［1］。近年来，碳氢燃料裂解所带来的化学热沉使其成为发动机热管理技术中的重要冷却介质^［2-3］。在该技术中，冷却燃料流经发动机壁面毫米级别的通道发生裂解吸热，从而冷却发动机壁面。与传统的石油裂解过程不同，碳氢燃料裂解发生在超临界条件下，这导致燃料物性发生显著变化，理想气体假设不再适用。其次，燃料吸热升温经历多个复杂过程，裂解和流动高度耦合^［4-5］。因此，准确预测燃料在管内流动和反应热力学特性对于发动机热管理技术的发展至关重要。

一维柱塞流模型（Plug Flow Reactor， PFR）通过将燃料在管内的流动简化为一维特征，能够快速求解并获取燃料沿程的热力学等信息，常用于裂解机理的发展及发动机设计性能评估^［6］。PFR模型需要结合实际气体状态方程（Equation of State， EOS）进行耦合求解，从而描述流体的P⁃V⁃T关系。长期以来，已有多种真实气体状态方程被提出，如范德华状态方程、索阿韦-雷德利希-邝氏（Soave-Redlich-Kwong，SRK）状态方程、彭-罗宾逊（Peng-Robinson， PR）状态方程、彭-罗宾逊体积平移（Peng-Robinson with Volume-Translation method， PRVT）状态方程和雷德利希-邝氏-彭-罗宾逊（Redlich-Kwong-Peng-Robinson， RKPR）状态方程^［7-8］。研究者基于这些实际状态方程的PFR模型，研究了碳氢燃料在超临界条件下裂解反应流动和换热。例如，Luo等^［9］基于PR状态方程研究了真实流体效应在超临界裂解反应过程中的影响；Cheng等^［10］基于RKPR方程，分析了带壁面冷却的冲压喷嘴的热力学特性；Xu等^［11］采用RKPR方程描述流体的P⁃V⁃T关系和热力学特征，利用PFR模型快速评估碳氢燃料在裂解过程中主动冷却通道的对流传热特性；王慧君等^［12］基于PR状态方程，采用PFR模型研究了质量流量对燃料裂解反应的影响；Jiang等^［13-14］基于RKPR状态方程的PFR模型研究了碳氢燃料在多通道流量偏移下的行为。上述研究中的解析式求导耦合PFR方法形式复杂，且缺乏通用性，这限制了PFR模型的进一步拓展和应用。

本文建立了一种基于差分求导的反应流模拟方法，克服了传统PFR模型耦合状态方程时，由复合热力学参数对原始变量解析求导带来的复杂性。基于该方法，本文进一步研究了不同状态方程对超临界燃料裂解反应过程的影响。通过对比不同状态方程在3~10 MPa裂解过程中的密度、黏度、温度和转化率等关键反应参数，为超临界反应流的建模分析提供了建议。

1 方法

1.1 PFR数学模型及求解方法

PFR是理想的一维反应流模型，其基本假设为： 1） 稳态流动； 2） 忽略沿通道的分子扩散和湍流质量扩散； 3） 垂直于流动方向的参数相等，即一维流动。根据这些假设，对于控制体（dx），质量、动量、能量的一维稳态守恒方程如式（1）~式（3）所示。

\frac{\partial \rho u}{\partial x}=\mathrm{0}(1)

\rho u\frac{\mathrm{d}u}{\mathrm{d}x}+\frac{\mathrm{d}p}{\mathrm{d}x}+\frac{\mathrm{1}}{A}\frac{\mathrm{d}F}{\mathrm{d}x}=\mathrm{0}(2)

\rho u\frac{\mathrm{d}\left(e+\frac{u^{\mathrm{2}}}{\mathrm{2}}\right)}{\mathrm{d}x}+\frac{\mathrm{d}(pu)}{\mathrm{d}x}=\frac{a_{\mathrm{e}}{Q}_{\mathrm{e}}}{A}(3)

式中，ρ、u、p、e分别表示流体的密度、速度、压力、内能；A、Q_e、a_e分别表示横截面积、热通量和单位长度的表面积。式（2）左侧第3项代表沿x坐标的黏性阻力的影响，可以表示为局部雷诺数的函数，如式（4）所示。

Re=\frac{\mathrm{4}{q}_{\mathrm{m}}}{\mathrm{\pi }D\mu }(4)

式中，Re、q_m、μ分别表示流体的雷诺数、质量流量和黏度；D为横截面的直径。混合流体的黏度通过Chung等^［15］提出的方法计算。该方法采用经验修正因子来考虑流体高密度对黏度的影响，其计算黏度的主要方程如式（5）所示。

\mu ={\eta }^{\mathrm{*}}\frac{\mathrm{36.344}(M_{\mathrm{m}}{T}_{\mathrm{c}\mathrm{m}}{)}^{\mathrm{1}/\mathrm{2}}}{V_{\mathrm{c}\mathrm{m}}^{\mathrm{2}/\mathrm{3}}}(5)

式中，η^* 为与混合物密度相关的参数，M_m、T_cm、V_cm为与分子量、临界温度、临界比体积及偶极矩相关的参数。对混合物需使用特定的混合规则计算。

黏度确定以后，可根据雷诺数分段计算黏性阻力项。对于层流流动（Re<2300），圆管通道的分析结果如式（6）所示。

\frac{\mathrm{d}F}{\mathrm{d}x}=\sqrt[]{A\mathrm{\pi }}\rho u^{\mathrm{2}}\frac{\mathrm{16}}{Re}(6)

而对于湍流流动，可以使用近似光滑通道的Blasius公式，如式（7）所示。

\frac{\mathrm{d}F}{\mathrm{d}x}=\sqrt[]{A\mathrm{\pi }}\rho u^{\mathrm{2}}\frac{\mathrm{0.079}\mathrm{1}}{R{e}^{\mathrm{0.25}}}(7)

化学反应本质上是一个改变物种含量的过程。物种的稳定守恒方程可用于描述这一过程，如式（8）所示。

\rho u\frac{\mathrm{d}{Y}_i}{\mathrm{d}x}=M_i{w}_i(8)

式中，M_i 、Y_i 分别表示物种i的分子量、质量分数；w_i 是物种i的净反应速率，应根据反应机理计算。典型的燃料反应机理由几十到几百种物种和反应组成。对于物种i， w_i 应通过同时求解涉及该物种的所有反应速率给出，如式（9）所示。

w_i={\displaystyle \sum _k^K}(v_{ki,r}-v_{ki,f})q_k(9)

式中，v_{ki,\mathrm{f}}和v_{ki,\mathrm{r}}分别表示第k反应中物种i的正向和反向化学计量系数；q_k 为第k反应的进展速率变量，由正向和反向速率的差值给出。

为使方程可解，式（2）和式（3）热力学参数的复合变量（压力p，内能e）应展开到原始变量（密度ρ、温度T、质量分数Y_i ）。对展开后的导数重新组合，可得到原始变量的可求解的ODEs形式，如式（10）所示。使用常微分方程（Ordinary Differential Equations，ODEs）求解器——变系数常微分方程求解器（Variable-coefficient ODE solver，VODE）即可实现方程组（10）沿x方向推进求解。

\left\{\begin{array}{l}\frac{\mathrm{d}{Y}_i}{\mathrm{d}x}=\frac{M_i{w}_i}{\rho u},i=\mathrm{1,2},\mathrm{3},\cdots ,k\\ \begin{array}{l}\frac{\mathrm{d}\rho }{\mathrm{d}x}=\left(\frac{n_{\mathrm{1}}\cdot m_{\mathrm{3}}}{m_{\mathrm{1}}}-\frac{n_{\mathrm{1}}\cdot m_{\mathrm{4}}}{m_{\mathrm{1}}}-n_{\mathrm{3}}-n_{\mathrm{4}}\right)/\\ \left(n_{\mathrm{2}}-\frac{n_{\mathrm{1}}\cdot m_{\mathrm{2}}}{m_{\mathrm{1}}}\right)\end{array}\\ \frac{\mathrm{d}u}{\mathrm{d}x}=-\frac{u}{\rho }\cdot \frac{\mathrm{d}\rho }{\mathrm{d}x}\\ \frac{\mathrm{d}T}{\mathrm{d}x}=\frac{\left(m_{\mathrm{4}}-m_{\mathrm{2}}\cdot \frac{\mathrm{d}\rho }{\mathrm{d}x}-m_{\mathrm{3}}\right)}{m_{\mathrm{1}}}\\ \frac{\mathrm{d}p}{\mathrm{d}x}=u^{\mathrm{2}}\cdot \frac{\mathrm{d}\rho }{\mathrm{d}x}-n_{\mathrm{4}}\end{array}\right.(10)

1.2 热力学参数的差分求导方法

在方程组（10）的求解过程中，热力学参数压力p和内能e对原始变量（ρ、T、Y_i ）的偏导数的计算十分重要。这些导数应根据状态方程及对应的热力学关系求解，尤其对于高压反应过程，真实气体状态方程是必需的。立方型状态方程是一种被广泛使用的真实气体状态方程。这些方程通常只需使用较少的性质参数。然而，即便如此，对于多组分的流体混合物，使用立方型状态方程计算压力、内能的偏导数的推导仍十分繁琐^［8］。这在一定程度上限制了模型的可拓展性。因此，为了简化这种计算方法，本文采用一种差分方法来求解压力内能的偏导数，具体步骤如下。对于立方型状态方程，压力、内能均可表示为一组标量的函数，如式（11）所示。

\mathrm{\Psi }=F\left(\phi \right),\phi =\left(T,\rho ,Y_{\mathrm{1}},Y_{\mathrm{2}},\dots ,Y_n\right)(11)

式中，\mathrm{\Psi }代表压力或内能；\phi代表一组标量，由温度、密度及各组分质量分数（Y₁，Y₂，…，Y_n ）组成。根据中心差分方法，因变量\mathrm{\Psi }对任一自变量{\phi }_i\in \phi的偏导数可近似计算为式（12）。

{\left.\frac{\partial \mathrm{\Psi }}{\partial {\phi }_i}\right|}_{{\phi }_{j\ne i}}={\left.\frac{F({\phi }_i+\mathrm{\Delta }{\phi }_i)-F({\phi }_i-\mathrm{\Delta }{\phi }_i)}{\mathrm{2}\mathrm{\Delta }{\phi }_i}\right|}_{{\phi }_{j\ne i}}(12)

特别地，当{\phi }_i=\mathrm{0}时（例如组分的质量分数），偏导数如式（13）所示。

{\left.\frac{\partial \mathrm{\Psi }}{\partial {\phi }_i}\right|}_{{\phi }_{j\ne i}}={\left.\frac{F({\phi }_i+\mathrm{\Delta }{\phi }_i)}{\mathrm{\Delta }{\phi }_i}\right|}_{{\phi }_{j\ne i}}(13)

式中，\mathrm{\Delta }{\phi }_i如式（14）所示。

\mathrm{\Delta }{\phi }_i=\alpha \cdot {\phi }_i(14)

实际计算中可取α=0.001。

通过以上方式，可以方便地利用不同的状态方程计算压力和内能对标量的偏导数。然后，控制式（1）~式（3）及式（8）仍可整理为原始变量温度、密度及各组分质量分数的ODEs形式，并采用VODE求解器沿x方向进行求解。本文对比计算了4种不同的立方型状态方程。

1） SRK EOS

SRK状态方程描述压力与温度、密度及平均分子量之间关系的原始形式如式（15）所示。

p=\frac{\rho RT}{\overline{M}-\rho b}-\frac{a{\rho }^{\mathrm{2}}}{{\overline{M}}^{\mathrm{2}}+b\rho \overline{M}}(15)

相应地，内能计算式如式（16）所示。

e=e_{\mathrm{i}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{l}}+\frac{\mathrm{1}}{b\overline{M}}\left[T\left(\frac{\partial a}{\partial T}\right)-a\right]\mathrm{l}\mathrm{n}\left(\frac{\overline{M}+b\rho }{\overline{M}}\right)(16)

2） PR EOS

PR状态方程是在SRK状态方程基础上修改压缩因子而来，其描述压力与温度、密度及平均分子量关系的原始形式如式（17）所示。

p=\frac{\rho RT}{\overline{M}-\rho b}-\frac{a{\rho }^{\mathrm{2}}}{{\overline{M}}^{\mathrm{2}}+\mathrm{2}b\rho \overline{M}-{\rho }^{\mathrm{2}}{b}^{\mathrm{2}}}(17)

相应地，内能的计算式如式（18）所示。

e=e_{\mathrm{i}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{l}}+\frac{T\frac{\partial a}{\partial T}-a}{\mathrm{2}\sqrt[]{\mathrm{2}}b\overline{M}}\mathrm{l}\mathrm{n}\frac{\overline{M}+\rho b+\sqrt[]{\mathrm{2}}\rho b}{\overline{M}+\rho b-\sqrt[]{\mathrm{2}}\rho b}(18)

式中，e_ideal为理想条件下的内能， \overline{M}为平均分子量。

3） PRVT EOS

原始形式PR方程在预测某些物质的液态密度时往往偏低。为了修正这一偏差，研究者提出了体积平移方法。本文使用Abudour等^［16］提出的体积转换方法对原始PR方程（15）求得的摩尔体积进行直接转换。相应的体积转换方程如式（19）所示。

V_{\mathrm{V}\mathrm{P}\mathrm{T}\mathrm{R}}=V_{\mathrm{P}\mathrm{R}}+c-{\delta }_{\mathrm{c}}\left(\frac{\mathrm{0.35}}{\mathrm{0.35}+d}\right)(19)

式中，V_VTPR是转换后的摩尔体积，V_PR是使用原始PR方程计算得到的摩尔体积。方程中其他几个参数的表达式分别如式（20）~式（23）所示。

{\delta }_{\mathrm{c}}=\left(\frac{R{T}_{\mathrm{c}}}{p_{\mathrm{c}}}\right)\left(z_{\mathrm{c},\mathrm{E}\mathrm{o}\mathrm{S}}-z_{\mathrm{c},\mathrm{E}\mathrm{x}\mathrm{p}}\right)(20)

c=\left(\frac{R{T}_{\mathrm{c}}}{p_{\mathrm{c}}}\right)(c_{\mathrm{1}}-(\mathrm{0.004}+c_{\mathrm{1}})e^{-\mathrm{2}d})(21)

d=\frac{\mathrm{1}}{R{T}_{\mathrm{c}}}{\left(\frac{\partial p_{\mathrm{P}\mathrm{R}}}{\partial \rho }\right)}_T(22)

c_{\mathrm{1}}=\mathrm{0.426}\mathrm{6}{z}_{\mathrm{c}}-\mathrm{0.110}\mathrm{1}(23)

式中，T_c、 p_c和z_c，Exp分别为实验临界温度、临界压力和临界压缩因子。此外，对于PR EOS， z_c，EOS通用值为0.307 4。当利用式（17）及体积转换方程式（19）~式（22）计算压力时，需要通过式（17）和式（23）进行迭代求解，以得到V_PR。

本质上，体积转换方程改变了PR EOS对物质P⁃V⁃T关系的描述，修正后的压力对密度、温度的偏导数已无法显式导出。已有研究者采用迭代方法求解各偏导数，但目前尚无证据表明这些偏导数能够提高热力学计算的精度。因此，在使用体积转换的PR EOS时，本文仍采用式（18）计算内能。

4） RKPR EOS

RKPR EOS是一种三参数立方型状态方程。与PR方程相比，其主要区别在于：RKPR EOS将PR方程中的常量表示为临界压缩因子的函数。因此RKPR EOS被认为在更宽的物性范围内具有更高的精度。其表达式如式（24）所示。

p=\frac{\rho RT}{\overline{M}-b\rho }-\frac{a{\rho }^{\mathrm{2}}}{(\overline{M}+{\delta }_{\mathrm{1}}b\rho )(\overline{M}+{\delta }_{\mathrm{2}}b\rho )}(24)

相应地，内能计算表达式如式（25）所示。

e=e_{\mathrm{i}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{l}}+\frac{\mathrm{1}}{({\delta }_{\mathrm{1}}-{\delta }_{\mathrm{2}})b\overline{M}}\left[T\left(\frac{\partial a}{\partial T}\right)-a\right]\mathrm{l}\mathrm{n}\left(\frac{\overline{M}+{\delta }_{\mathrm{1}}b\rho }{\overline{M}+{\delta }_{\mathrm{2}}b\rho }\right)(25)

2 结果

2.1 不同状态方程的热物性计算结果

密度和黏度作为重要的热物性参数，直接影响反应速率中浓度的换算及雷诺数等关键动力学参数的计算，从而对PFR模型预测的温度、反应转化率及物种分布等结果产生显著影响。为验证不同状态方程在高压条件下的热物性预测能力，本文选取正癸烷（T_c=617.7 K，Z_c=0.256）作为代表性燃料，分别采用4种典型立方型状态方程（SRK、PR、PRVT、RKPR）及理想气体状态方程，对其在压力为3.5 MPa，温度300~800 K范围内的密度和黏度进行了计算，并将结果与NIST数据库中的标准值进行了对比，结果如图1所示。

从图1a可见，在密度预测方面，理想气体状态方程由于完全忽略分子间相互作用，整体低估流体密度。随着温度升高，真实气体效应逐渐减弱，预测误差相对缩小，但在整个温度范围内与实际数据仍存在显著差距。SRK和PR状态方程在低于拟临界温度的区间内，同样存在密度低估问题。这一偏差主要源于它们分别采用固定的临界压缩因子（SRK：0.333，PR：0.307），而未能反映正癸烷较低的实际压缩性（Z_c=0.256），进而高估了体积膨胀效应。相比之下，PRVT状态方程通过引入体积转化方法，有效修正了PR模型在液相区域的压缩偏差，显著提高了密度预测的准确性。在本文考察的温度范围内，PRVT状态方程计算的密度值与NIST数据相对吻合，最大相对误差低于10%。RKPR状态方程则通过显式引入Z_c 拟合参数，自适应调整热力学关系，其密度预测精度虽次于PRVT，但优于SRK和PR，且表现出良好的趋势一致性。图1b为各状态方程对黏度的预测结果。由于本文采用的Chung方法对密度和临界体积高度敏感，状态方程对密度计算中的误差直接传递至黏度计算结果中。由图可见，理想气体模型完全无法反映密度变化趋势；SRK与PR状态方程在此区间内同样低估黏度，这一现象主要源于其密度预测的偏差。PRVT状态方程由于其密度估算精度较高，使得在300~600 K范围内计算出的黏度曲线与NIST数据呈现出高度一致性，整体误差显著低于其他模型。RKPR状态方程在黏度预测上略逊于PRVT状态方程，但远优于未修正的SRK和PR模型。从整体来看，不同状态方程对高压流体的修正主要体现在低温区域，PRVT状态方程在正癸烷热物性预测中表现最优。

2.2 不同状态方程高压下对裂解过程计算的影响

为评估高压条件下不同状态方程对超临界裂解过程非理想效应的捕捉能力，本研究计算了3~10 MPa条件下正癸烷的裂解过程。实验方法见文献［9］。计算采用Jia等^［17］提出的正癸烷裂解机理，此机理共包含842步反应和162种物种，其中自由基的临界性质采用孙彦锦等^［18］发展的基团贡献法估算。模拟条件为：通道直径0.193 mm、通道长度800 mm，进口质量流量为1.2 g/s，热流密度为730 kW/m²。采用不同状态方程计算的出口温度和转化率与实验值的对比结果如图2所示。由图2a可知，理想气体模型无法捕捉高压流体的压缩效应，会过高评估整个流体温度范围内的恒压比热，进而导致计算温度显著高于实验值，其平均出口温度相对误差为2.18%。相反，考虑分子吸引力和体积修正的不同状态方程，能有效捕捉裂解过程的温度变化，其平均出口温度相对误差不超过0.3%。从图2b可知，理想气体模型计算的出口转化率同样高于实验值，平均出口转化率相对误差为20.41%。此误差主要来源于两方面：一是理想气体模型计算得到的过高温度会加快裂解反应速率；二是理想气体模型在超临界裂解条件计算的流体密度会低估流体真实密度，进而导致反应浓度计算偏差。考虑修正的模型中，平均出口转化率相对误差最大的为PR模型计算结果（10.96%），最小的为SRK模型计算结果（6.66%）。在碳氢燃料裂解过程中，随着裂解反应的进行，会逐渐生成多种小分子产物（如H₂、CH₄、C₂H₄等），这些产物的生成导致裂解后混合物的临界温度和临界压力发生显著变化^［19］。而临界温度和临界压力的变化直接影响了流体的相行为，尤其是在高转化率裂解过程中，当流体从压缩液相跨越超临界区再到气相时，温度的突变和转化率的非线性增长变得尤为显著。SRK、PR及其修正形式（PRVT、RKPR）等立方型EOS，能够基于各组分临界常数动态调整压缩因子和体积项，能够反映组分变化对临界点的影响，准确预测相变过程中密度与温度的突变。因此，修正模型可使出口温度和裂解转化率的预测结果更贴近实验值。需要指出的是，实际方程在高压条件下对转化率的预测仍存在误差，这主要是由于其未考虑化学势对净反应速率的影响，且压力越高，偏差越大；但在3~4 MPa压力范围内，模型预测结果与实验值吻合较好，此时化学势影响较小，热物性占据主导。因此实际气体模型预测的转化率显著低于理想气体模型。综上，不同EOS对临界点的描述精度，直接决定了其在高转化率条件下对裂解温度分布和转化率趋势的模拟可靠性。

图3为5 MPa条件下不同状态方程计算的沿管程温度和转化率与实验值的对比结果。从图3a可见，理想气体模型高估流体温度，在拟临界温度附近与实验值误差高达100 K。随着温度升高，裂解反应开始发生，导致温度上升趋势变缓，同时真实流体压缩效应减缓，模型与实验值误差逐渐减少。此外，在770 K以下的非反应阶段，不同立方型状态方程预测的沿程温度无明显差异。该阶段为物理吸热，这说明4种立方型状态方程在比热容预测方面的差距很小。从图3b所知，4种立方型状态方程对裂解反应起始点的预测存在一定影响，这主要是由于物性预测偏差导致反应速率浓度计算偏差。但总体而言，尽管4种立方型状态方程在低温范围内对密度、黏度等物性的预测差异显著，但对最终出口温度的影响却并不明显。具体而言，流体温度预测的差异主要源于实际气体模型对比热容的预测偏差；而转化率的偏差则源于密度和温度预测偏差累积，从而导致反应速率浓度计算出现偏差。在5 MPa下，由于实际模型未考虑化学势对反应速率促进的影响，而理想气体模型又过高预估温度，从而加速了裂解反应，因此理想气体模型高估了转化率，而实际气体模型低估了转化率。

图4为采用中心差分法求解沿程温度对距离的二阶导数的结果。由图可见，理想气体模型的二阶导数曲线仅存在1个极值点，而4种立方型状态方程对应的曲线却有3个极值点。其中，前2个极值点所在区间对应临界温度变化区间，最后1个极值点则标志着燃料进入快速裂解阶段。理想气体模型未能捕捉拟临界温度的转折，过高评估温度，导致裂解提前发生，起始位置为0.46 m。而4种立方型模型因其物性计算差异导致反应速率不同，发生显著裂解的位置约为0.52 m。从图中进一步观察到，在低温阶段，温度变化趋势有所不同，尤其在拟临界温度区域差异更为明显，这是由于各类状态方程自身的特性所导致。此外，从第3个极值点的位置分布可以进一步证实，4种立方型状态方程对裂解反应起始点的预测具有影响。

3 结论

本文建立了一种基于差分求导的超临界燃料裂解反应流模拟方法，此方法有效克服PFR模型与实际状态方程耦合过程中的计算难点，显著简化了内能和压力对原始变量求导的复杂程度，从而实现了不同形式状态方程的通用化求解。基于此方法，本文研究了不同状态方程对超临界燃料裂解反应的影响，主要结论如下：

1） 热力学参数的差分求导法，能够有效实现不同形式的立方型实际气体状态方程与反应流PFR模型的耦合，模拟结果可靠。

2） 与理想气体模型相比，修正后的立方型状态方程（PRVT和RKPR）在高压条件下提供了更准确的热物性预测。PRVT和RKPR模型对密度的预测结果与NIST标准数据相比，最大相对误差低于10%。值得注意的是，不同状态方程对密度的计算误差会直接传递至黏度计算，进而导致模拟结果存在显著差异。

3） 在出口温度预测方面，4种立方型状态方程与实验数据的相对误差不超过0.3%，而理想气体模型的误差为2.18%。在裂解转化率预测方面，实际气体模型的预测精度整体优于理想气体模型（其相对误差为20.41%）。其中，PR模型与实验数据的相对误差最大，为10.96%；SRK模型的相对误差最小，为6.66%。

参考文献

[1] Van Geem K M， Dhuyvetter I， Prokopiev S， et al.Coke formation in the transfer line exchanger during steam cracking of hydrocarbons ［J］.Ind Eng Chem Res， 2009， 48（23）： 10343-10358.

[2] Wang J， Jin H C， Gao H， et al. Cooling capacity optimization of hydrocarbon fuels for regenerative cooling ［J］.Appl Therm Eng， 2022， 200： 117661.

[3] Ward T A， Ervin J S， Striebich R C， et al.Simulations of flowing mildly-cracked normal alkanes incorporating proportional product distributions ［J］.J Propul Power， 2004， 20（3）： 394-402.

[4] Wang C， Du C P， Shang J X， et al.A comprehensive review of the thermal cracking stability of endothermic hydrocarbon fuels ［J］.J Anal Appl Pyrol， 2023， 169： 105867.

[5] Wang Y， Liang C， Zhang Q Y， et al.The study on pyrolysis characteristics of hydrocarbon fuel with fixed heat flux ［J］.Chemical Research and Application， 2017， 29（8）： 1171-1176.

[6] 王雨，梁策，张其翼，等.等热流密度条件下碳氢燃料裂解特性研究［J］.化学研究与应用， 2017， 29（8）： 1171-1176.

[7] Yu W L， Zhou W X， Long L， et al.Effects of different simulation models on pyrolysis performance predictions for regenerative cooling applications ［J］.J Anal Appl Pyrol， 2021， 156： 105136.

[8] Kim J S， Seo J， Han D， et al.Prediction of thermochemical and transport properties of hydrocarbon aviation fuel in supercritical state with thermal decomposition ［J］.Fuel， 2022， 325： 124805.

[9] Kim S-K， Choi H-S， Kim Y.Thermodynamic modeling based on a generalized cubic equation of state for kerosene/LOx rocket combustion ［J］.Combust Flame， 2012， 159（3）： 1351-1365.

[10] Luo L， Chen Y， Li X Y.Modeling and simulation for n-decane cracking flow based on detailed reaction kinetics and complete real gas thermodynamics under supercritical pressure ［J］.Fuel， 2024， 366： 131411.

[11] Cheng K L， Feng Y， Jiang Y G， et al.Thermodynamic analysis for recuperation in a scramjet nozzle with wall cooling ［J］.Appl Therm Eng， 2017， 121： 153-162.

[12] Xu Q， Lin G W， Li H W.Evaluation model for fast convective heat transfer characteristics of thermal cracked hydrocarbon fuel regenerative cooling channel in hydrocarbon-fueled scramjet ［J］.Appl Therm Eng， 2021， 199： 117616.

[13] Wang H J， Shen W H， Chen Y， et al.Effect of mass flow rate on supercritical cracking reaction process ［J］.Chemical Research and Application， 2021， 33（5）： 902-910.

[14] 王慧君，沈文豪，陈玉，等.质量流量对超临界裂解反应过程的影响研究［J］.化学研究与应用， 2021， 33（5）： 902-910.

[15] Jiang Y G， Wang Q， Zhang D， et al.One-dimensional modelling and analysis of flow excursion of cracking hydrocarbon fuel in SCRamjet parallel cooling channels ［J］.Fuel， 2023， 340： 127462.

[16] Jiang Y G， Wang Q， Zhang D， et al.One-dimensional modelling and simulation of multiple solution characteristics of hydrocarbon fuel flow excursion in scramjet parallel cooling channels ［J］.Appl Therm Eng， 2023， 221： 119757.

[17] Chung T H， Ajlan M， Lee L L， et al.Generalized multiparameter correlation for nonpolar and polar fluid transport properties ［J］.Ind Eng Chem Res， 1988， 27（4）： 671-679.

[18] Abudour A M， Mohammad S A， Robinson Jr R L， et al.Volume-translated Peng-Robinson equation of state for saturated and single-phase liquid densities ［J］.Fluid Phase Equilibr， 2012， 335： 74-87.

[19] Jia Z J， Huang H Y， Zhou W X， et al.Experimental and modeling investigation of n-decane pyrolysis at supercritical pressures ［J］.Energy Fuels， 2014， 28（9）： 6019-6028.

[20] Sun Y J， Yao Q， Li Z R， et al.Studies on the calculation of Lennard-Jones parameters for intermediates in hydrocarbon combustion ［J］.Chemical Journal of Chinese Universities， 2016. 37（2）： 328-334.

[21] 孙彦锦，姚倩，李泽荣，等.碳氢化合物燃烧中间体Lennard-Jones系数的计算［J］.高等学校化学学报， 2016， 37（2）： 328-334..

[22] Jiang R P， Liu G Z， You Z Q， et al.On the critical points of thermally cracked hydrocarbon fuels under high pressure ［J］.Ind Eng Chem Res， 2011， 50（15）： 9456-9465.

基金资助

国家自然科学基金(T2441001)