0 引言
近年来，在碳捕集、利用与封存（CCUS）技术进步的推动下，注CO2提高采收率技术在页岩油开发中的应用日益增多[1⇓-3]。页岩储层内通常存在大量的纳米孔隙[4⇓⇓⇓-8]，纳米孔隙中孔壁和流体分子之间的相互作用显著影响流体分子的热运动，从而导致纳米孔隙中流体的相行为与体相流体的相行为存在明显差异，这种差异被称为限域效应[9⇓⇓⇓⇓-14]。限域效应对流体的影响主要体现在毛管压力增大、临界性质参数偏移、饱和压力偏移、密度分布不均匀等方面。准确理解限域效应对于预测页岩油储层中注入流体的行为和设计有效的提高采收率方法十分重要。目前研究纳米孔隙中流体相行为的方法主要有实验研究、分子动力学模拟（MDS）和理论计算[15⇓-17]。
MDS可以精准描述纳米孔隙中不同类型流体的相行为，直观体现流体的密度分布及组分分布，但是不能很好地处理较大孔隙尺度甚至是矿场尺度这些大尺度问题的计算。状态方程（EOS）是一种被广泛应用于流体相行为研究的理论计算方法，其中Peng-Robinson状态方程（简称PR状态方程）[24]近年来被广泛地应用于石油工业。Nojabaei等[25]使用PR状态方程对纳米孔隙内流体相行为进行研究，发现在相平衡计算过程中引入毛管压力后，流体泡点压力随孔径变小而降低；Teklu等[26]发现在PR状态方程中引入毛管压力后，计算过程中较高的露点压力会增加，较低的露点压力会降低，临界点处的偏移程度变化较小，几乎与其体相流体的临界点相同；Song等[27]提出了一种描述纳米孔隙中流体相行为的新方法，修改了PR状态方程中的有效摩尔体积项，使预测的受限流体临界性质变化规律更符合实验结果；Feng等[28]同样对PR状态方程进行了修正，在方程中考虑了孔壁润湿性对碳氢化合物的影响，发现限域效应不仅受孔径大小影响，还受孔壁润湿性影响；Yang等[29]考虑临界性质偏移、毛管压力及分子-孔壁相互作用的影响，对PR状态方程进行修正，发现在10 nm孔隙中流体泡点压力较其体相流体有所下降。
通过对前人研究进行调研总结，发现目前尚没有一种修正PR状态方程能综合考虑临界性质偏移、毛管压力、流体分子-孔壁相互作用以及吸附对纳米孔隙中流体相行为的影响。因此，本文在前人研究的基础上，对PR状态方程进行修正，充分考虑临界性质偏移、毛管压力、流体分子-孔壁相互作用以及吸附对闪蒸计算过程的影响。其次，通过拟合前人的实验数据及MDS数据，验证本文提出的修正PR状态方程。然后，基于此方程结合气相-液相-吸附相三相相平衡的闪蒸计算方法对纳米孔隙中CO2-页岩油相行为进行分析。
1 数学模型
1.1 流体分子与孔壁间的相互作用
使用传统PR状态方程进行计算时通常忽略流体分子与孔壁间的相互作用。然而在纳米孔隙中，流体分子与孔壁间的相互作用力会减弱流体分子间的吸引力，如果忽略将导致计算不够精准。为了解决这一问题，本文将PR状态方程中的压力拆解为引力项和斥力项，在引力项中引入一个新的系数c用来表征流体分子与孔壁间的相互作用，即：
$p={{p}_{\text{A}}}+{{p}_{\text{R}}}$
${{p}_{\text{A}}}=-\frac{a-c}{{{V}_{\text{m}}}\left( {{V}_{\text{m}}}+b \right)+b\left( {{V}_{\text{m}}}-b \right)}$
${{p}_{\text{R}}}=\frac{RT}{{{V}_{\text{m}}}-b}$
其中 $a=0.457\ 24\frac{{{R}^{2}}T_{\text{c}}^{2}}{{{p}_{\text{c}}}}\alpha $ $b=0.077\ 80\frac{R{{T}_{\text{c}}}}{{{p}_{\text{c}}}}$
$\alpha ={{\left[ 1+m\left( 1-\sqrt{{{T}_{\text{r}}}} \right) \right]}^{2}}$
$m=0.374\ 64+1.542\ 26\omega -0.269\ 92{{\omega }^{2}}$
采用范德华流体混合规则[30]和二元相互作用系数来计算分析多组分流体在纳米孔隙中的行为。
1.2 有效摩尔体积修正
在页岩储层纳米孔隙内，吸附在孔壁上的流体分子数占孔隙中流体分子总数的比例较高。吸附在孔壁上的流体分子会造成体相区的有效摩尔体积发生变化，因此需要对传统PR状态方程的摩尔体积进行修正。如图1 所示，孔隙内流体以吸附相和体相的形式存在，吸附相流体存在于孔壁附近。
为辅助计算流体分子的物质的量，引入了Alder等[31]提出的有效分子体积系数：
$\varepsilon =\frac{{{V}_{\text{e}}}}{{{V}_{\text{real}}}}$
纳米孔隙中流体分子总数的计算方法如下[32]：
${{n}_{\text{t}}}=\frac{{{V}_{\text{p}}}}{{{V}_{\text{et}}}}=\frac{\text{ }\!\!\pi\!\!\text{ }{{R}_{\text{p}}}^{2}L}{{{\varepsilon }_{\text{t}}}{{V}_{\text{real}}}}$
吸附在孔壁上的流体分子数可以表示为吸附分子所占用的体积与每个吸附分子的有效体积之比[32]：
${{n}_{\text{a}}}=\frac{{{V}_{\text{a}}}}{{{V}_{\text{ea}}}}=\frac{\text{ }\!\!\pi\!\!\text{ }L\left( {{R}_{\text{p}}}^{2}-{{R}_{\text{f}}}^{2} \right)}{{{\varepsilon }_{\text{a}}}{{V}_{\text{real}}}}=\frac{\text{ }\!\!\pi\!\!\text{ }L\tau \left( 2{{R}_{\text{p}}}-\tau \right)}{{{\varepsilon }_{\text{a}}}{{V}_{\text{real}}}}$
（6）式中吸附层厚度可由经验公式来计算[32]：
$\tau =\frac{P}{\ln \left( {{R}_{\text{p}}}/{{\sigma }_{\text{LJ}}} \right)}+Q\frac{{{\sigma }_{\text{LJ}}}}{{{R}_{\text{p}}}}$
其中 $P=-8.314\ 0\times {{10}^{-14}}{{M}^{2}}+2.047\ 5\times {{10}^{-11}}M+$$3.088\ 6\times {{10}^{-11}}$
$Q=-6.356\ 5\times {{10}^{-14}}{{M}^{2}}+3.155\ 0\times {{10}^{-11}}M+$$5.853\ 8\times {{10}^{-10}}$
${{\sigma }_{\text{LJ}}}={{0.244}^{3}}\sqrt{{{T}_{\text{c}}}/{{p}_{\text{c}}}}$
传统PR状态方程中的摩尔体积可以表示为：
${{V}_{\text{m}}}=\frac{{{V}_{\text{p}}}}{{{n}_{\text{t}}}/{{N}_{\text{A}}}}$
则修正后的PR状态方程中的有效摩尔体积可以表示为：
${{{V}'}_{\text{m}}}=\frac{{{V}_{\text{p}}}}{{\left( {{n}_{\text{t}}}-{{n}_{\text{a}}} \right)}/{{{N}_{\text{A}}}}\;}=\frac{{{n}_{\text{t}}}}{{{n}_{\text{t}}}-{{n}_{\text{a}}}}{{V}_{\text{m}}}=\frac{{{V}_{\text{m}}}}{1-\frac{{{\varepsilon }_{\text{t}}}}{{{\varepsilon }_{\text{a}}}}\frac{\tau }{{{R}_{\text{p}}}^{2}}\left( 2{{R}_{\text{p}}}-\tau \right)}=\lambda {{V}_{\text{m}}}$
其中 $\lambda =\frac{1}{1-\frac{{{\varepsilon }_{\text{t}}}}{{{\varepsilon }_{\text{a}}}}\frac{\tau }{{{R}_{\text{p}}}^{2}}\left( 2{{R}_{\text{p}}}-\tau \right)}$
1.3 多影响因素条件下的PR状态方程修正
通过综合考虑流体分子与孔壁间的相互作用和有效摩尔体积变化对传统PR状态方程进行整合修正，如（10）式所示。修正后的PR状态方程能更好地解释页岩储层纳米孔隙中的独特物理现象。
$p=\frac{RT}{{{{{V}'}}_{\text{m}}}-b}-\frac{a-c}{{{{{V}'}}_{\text{m}}}\left( {{{{V}'}}_{\text{m}}}+b \right)+b\left( {{{{V}'}}_{\text{m}}}-b \right)}$
PR状态方程是一种立方型状态方程。对于立方型状态方程，临界点须满足特定的条件。在保持温度不变的情况下，临界点处压力相对于摩尔体积的一阶偏导数和二阶偏导数等于零，即：
${{\left( \frac{\partial p}{\partial {{{{V}'}}_{\text{m}}}} \right)}_{T={{T}_{\text{c}}}}}={{\left( \frac{{{\partial }^{2}}p}{\partial {{{{V}'}}_{\text{m}}}^{2}} \right)}_{T={{T}_{\text{c}}}}}=0$
将（10）式代入（11）式可得：
$-\frac{R{{T}_{\text{c}}}}{{{\left( {{{{V}'}}_{\text{m}}}-b \right)}^{2}}}+\frac{2\left( a-c \right)\left( {{{{V}'}}_{\text{m}}}+b \right)}{{{\left( {{{{V}'}}_{\text{m}}}^{2}+2{{{{V}'}}_{\text{m}}}b-{{b}^{2}} \right)}^{2}}}=0$
$\frac{2R{{T}_{\text{c}}}}{{{\left( {{{{V}'}}_{\text{m}}}-b \right)}^{3}}}-\frac{2\left( a-c \right)\left( 3{{{{V}'}}_{\text{m}}}^{2}+6{{{{V}'}}_{\text{m}}}b+5{{b}^{2}} \right)}{{{\left( {{{{V}'}}_{\text{m}}}^{2}+2{{{{V}'}}_{\text{m}}}b-{{b}^{2}} \right)}^{3}}}=0$
由（12）式和（13）式可得：
${{{V}'}_{\text{m}}}^{3}-3{{{V}'}_{\text{m}}}^{2}b-3{{{V}'}_{\text{m}}}{{b}^{2}}-3{{b}^{3}}=0$
求解（14）式可得：
$b=0.253\ 1{{{V}'}_{\text{m}}}=0.253\ 1\lambda {{V}_{\text{m}}}$
将（15）式代入（12）式可得：
$a-c=1.487\ 7R{{T}_{\text{c}}}{{{V}'}_{\text{m}}}=1.487\ 7R{{T}_{\text{c}}}\lambda {{V}_{\text{m}}}$
将（15）式和（16）式代入临界点处的（10）式中可得：
$\frac{{{p}_{\text{c}}}{{{{V}'}}_{\text{m}}}}{R{{T}_{\text{c}}}}=0.307\ 4$
结合（10）式—（17）式可得纳米孔隙中流体的临界温度和临界压力：
${{T}_{\text{cb}}}=0.170\ 2\frac{a-c}{\alpha bR}\lambda $
${{p}_{\text{cb}}}=0.013\ 2\frac{a-c}{\alpha {{b}^{2}}}{{\lambda }^{2}}$
对比纳米孔隙中流体临界压力和临界温度与体相流体临界温度和临界压力，得到临界性质变化量：
$\Delta {{T}^{*}}=\frac{{{T}_{\text{c}}}-{{T}_{\text{cb}}}}{{{T}_{\text{c}}}}=\lambda \left( c-1 \right)+1$
$\Delta {{p}^{*}}=\frac{{{p}_{\text{c}}}-{{p}_{\text{cb}}}}{{{p}_{\text{c}}}}={{\lambda }^{2}}\left( c-1 \right)+1$
$\ln \Delta {{T}^{*}}=0.984\ 8\ln \left( {{\sigma }_{\text{LJ}}}/{{R}_{\text{p}}} \right)+0.748\ 0$
$\ln \Delta {{p}^{*}}=1.010\ 9\ln \left( {{\sigma }_{\text{LJ}}}/{{R}_{\text{p}}} \right)+0.782\ 1$
结合（20）式—（23）式可得PR状态方程中两个修正系数与无因次孔径的关系式：
$\lambda =\frac{2.186\ 1{{\left( {{\sigma }_{\text{LJ}}}/{{R}_{\text{p}}} \right)}^{1.010\ 9}}-1}{2.112\ 9{{\left( {{\sigma }_{\text{LJ}}}/{{R}_{\text{p}}} \right)}^{0.984\ 8}}-1}$
$c=\frac{{{\left[ 2.112\ 9{{\left( {{\sigma }_{\text{LJ}}}/{{R}_{\text{p}}} \right)}^{0.984\ 8}}-1 \right]}^{2}}}{2.186\ 1{{\left( {{\sigma }_{\text{LJ}}}/{{R}_{\text{p}}} \right)}^{1.010\ 9}}-1}+1$
则临界性质的变化量可表示为：
$\Delta {{T}^{*}}=\frac{{{T}_{\text{c}}}-{{T}_{\text{cb}}}}{{{T}_{\text{c}}}}=\lambda \left( c-1 \right)+1=2.112\ 9{{\left( \frac{{{\sigma }_{\text{LJ}}}}{{{R}_{\text{p}}}} \right)}^{0.984\ 8}}$
$\Delta {{p}^{*}}=\frac{{{p}_{\text{c}}}-{{p}_{\text{cb}}}}{{{p}_{\text{c}}}}={{\lambda }^{2}}\left( c-1 \right)+1=2.186\ 1{{\left( \frac{{{\sigma }_{\text{LJ}}}}{{{R}_{\text{p}}}} \right)}^{1.010\ 9}}$
1.4 吸附相-气相-液相三相相平衡计算
传统的相平衡计算过程中只考虑了气液平衡，简单地认为吸附相固定附着在孔壁，忽略了吸附相对相平衡过程中组分变化的影响。本文采用多组分Langmuir吸附模型，使吸附相组分的变化以物质守恒的形式参与到相平衡计算中。
1.4.1 多组分Langmuir吸附模型
Langmuir吸附模型[17]可表示为：
${{m}_{\text{abs}}}={{m}_{\max }}\frac{{{k}_{\text{p}}}p}{1+{{k}_{\text{p}}}p}$
吸附平衡常数（kp）是与温度相关的函数，将吸附平衡常数与逸度相结合可以得到（29）式。多组分Langmuir吸附模型的相关参数通过实验数据拟合获得，可参照文献[40]。
${{m}_{\text{abs},i}}={{m}_{\text{max},i}}\left( T \right)\frac{{{g}_{i}}\left( T \right){{f}_{i}}}{1+\sum\limits_{i=1}^{N}{{{g}_{i}}\left( T \right){{f}_{i}}}}$
结合吸附相体积、体相流体密度和吸附量可以计算得到在三相相平衡计算过程中吸附相的相分率：
${{\delta }_{\text{A}}}=\frac{\sum\limits_{i=1}^{N}{{{m}_{\text{abs},i}}}}{{{\rho }_{\text{b}}}\left( {{V}_{\text{p}}}-{{V}_{\text{a}}} \right)+\sum\limits_{i=1}^{N}{{{m}_{\text{abs},i}}}}$
1.4.2 压缩因子和逸度计算
压缩因子和逸度的计算在相平衡计算中是必不可少的。由于本文对PR状态方程进行了修正，压缩因子和逸度的计算方式与传统的计算方式有所不同。压缩因子的计算式如下：
${{Z}^{3}}+\left( B-1 \right){{Z}^{2}}-\left( 3{{B}^{2}}+2B-A+C \right)Z+$$\left( {{B}^{3}}+{{B}^{2}}-AB+BC \right)=0$
其中 $A=\frac{ap}{{{R}^{2}}{{T}^{2}}}$ $B=\frac{bp}{RT}\lambda $ $ C=\frac{cp}{{{R}^{2}}{{T}^{2}}}$
结合理想气体状态方程、热力学关系式和压缩因子计算式可得气相或液相的逸度系数计算式：
$\ln {{\varphi }_{i}}=\frac{{{b}_{i}}}{b}\left( Z-1 \right)-\ln \left( Z-B \right)-\frac{A-C}{2\sqrt{2}B}\times $
$\left[ \frac{2\sum\limits_{j=1}^{N}{{{x}_{j}}\left( 1-{{k}_{ij}} \right)\sqrt{{{a}_{i}}{{a}_{j}}}}}{a}-\frac{{{b}_{i}}}{b} \right]\ln \frac{Z+\left( \sqrt{2}+1 \right)B}{Z-\left( \sqrt{2}-1 \right)B}$
$\left[ \frac{2\sum\limits_{j=1}^{N}{{{x}_{j}}\left( 1-{{k}_{ij}} \right)\sqrt{{{a}_{i}}{{a}_{j}}}}}{a}-\frac{{{b}_{i}}}{b} \right]\ln \frac{Z+\left( \sqrt{2}+1 \right)B}{Z-\left( \sqrt{2}-1 \right)B}$
1.4.3 毛管压力计算
与常规储层不同的是，在页岩储层纳米孔隙中孔隙半径较小，存在着较大的毛管压力，毛管压力引起的气液两相压差不可忽略。因此，本文在相平衡计算中考虑了纳米孔隙下的毛管压力，以提高模型的精度。气相压力、液相压力、毛管压力间存在如下关系[17]：
${{p}_{\text{G}}}={{p}_{\text{L}}}+{{p}_{\text{C}}}$
毛管压力可以使用Young-Laplace方程结合多组分混合物气相和液相之间的界面张力来计算[41]：
${{p}_{\text{C}}}={{10}^{-6}}\times \frac{2\sigma \cos \theta }{{{r}_{\text{e}}}}$
$\sigma ={{\left[ \sum\limits_{i=1}^{N}{{{P}_{\text{ch},i}}\left( {{\rho }_{\text{L}}}{{x}_{i}}-{{\rho }_{\text{G}}}{{y}_{i}} \right)} \right]}^{4}}$
1.4.4 三相相平衡计算
本文基于热力学定律进行纳米孔隙下气相-液相-吸附相三相相平衡计算。假设纳米孔隙中存在气相、液相、吸附相三相。综合考虑纳米孔隙中的毛管压力、限域效应及纳米孔隙中流体吸附对相平衡计算的影响。以平衡常数作为迭代变量，使用简单稳定的逐次替代法进行求解。假设初始状态下体系不存在吸附相，体系内所有流体全部为体相流体。
体系内达到平衡时满足气相、液相及吸附相中各组分逸度相等，即：
${{f}_{\text{G,}i}}={{f}_{\text{L,}i}}={{f}_{\text{A,}i}}$
该体系中，每相的相分率之和等于1，每相中所有组分的物质的量分数之和也等于1，即：
${{\delta }_{\text{G}}}+{{\delta }_{\text{L}}}+{{\delta }_{\text{A}}}=1$
$\sum\limits_{i=1}^{N}{{{y}_{i}}}=\sum\limits_{i=1}^{N}{{{x}_{i}}}=\sum\limits_{i=1}^{N}{{{w}_{i}}}=1$
${{z}_{\text{f},i}}={{\delta }_{\text{G}}}{{y}_{i}}+{{\delta }_{\text{L}}}{{x}_{i}}+{{\delta }_{\text{A}}}{{w}_{i}}$
三相相平衡计算具体计算过程如下。
①利用Wilson经验公式[42]估算气液平衡常数的初值，如（40）式所示。
$\ln {{K}_{i}}=5.37\left( 1+{{\omega }_{i}} \right)\left( 1-\frac{{{T}_{\text{c},i}}}{T} \right)+\ln \left( \frac{p}{{{p}_{\text{c},i}}} \right)$
②吸附相形成前，体相中的气相和液相存在如（41）式所示的关系。随着吸附相在体系中逐渐形成，体系内吸附相与体相存在如（42）式所示的关系。结合（36）式—（39）式及（41）式、（42）式，使用Newton-Raphson迭代法求解Rachford-Rice方程（见（43）式），计算体相中气相、液相的相分率。
${{z}_{\text{b},i}}=\beta {{x}_{i}}+\left( 1-\beta \right){{y}_{i}}$
其中 $\beta =\frac{{{\delta }_{\text{L}}}}{{{\delta }_{\text{L}}}+{{\delta }_{\text{G}}}}$
${{z}_{\text{b},i}}=\frac{{{z}_{\text{f},i}}-{{\delta }_{\text{A}}}{{{m}_{\text{abs},i}}}/{\sum\limits_{i=1}^{N}{{{m}_{\text{abs},i}}}}\;}{1-{{\delta }_{\text{A}}}}$
$\sum\limits_{i=1}^{N}{\frac{{{z}_{\text{b},i}}\left( {{K}_{i}}-1 \right)}{\lambda +{{K}_{i}}\left( 1-\beta \right)}}=0$
③根据质量守恒定律建立体相和吸附相之间的联系，如（44）式所示，然后根据（45）式、（46）式计算气相、液相中各组分的物质的量分数。
${{z}_{\text{b},i}}=\frac{{{z}_{\text{f},i}}\left[ {{\rho }_{\text{b}}}\left( {{V}_{\text{p}}}-{{V}_{\text{a}}} \right)+\sum\limits_{i=1}^{N}{{{m}_{\text{abs},i}}} \right]-{{m}_{\text{abs},i}}}{{{\rho }_{\text{b}}}\left( {{V}_{\text{p}}}-{{V}_{\text{a}}} \right)}$
${{x}_{i}}=\frac{{{z}_{\text{b},i}}}{\beta +\left( 1-\beta \right){{K}_{i}}}$
$ {{y}_{i}}=\frac{{{z}_{\text{b},i}}{{K}_{i}}}{\beta +\left( 1-\beta \right){{K}_{i}}}$
④计算并更新吸附相形成后体系中的毛管压力。
⑤结合修正后的PR状态方程分别计算气相、液相的压缩因子以及各组分的逸度、逸度系数。
⑥根据（47）式、（48）式计算并更新吸附相形成后的气液平衡常数。
${{K}_{i}}=\frac{{{p}_{\text{L}}}{{\varphi }_{\text{L,}i}}}{{{p}_{\text{G}}}{{\varphi }_{\text{G,}i}}}=\frac{{{y}_{i}}}{{{x}_{i}}}$
$K_{i}^{\left( n+1 \right)}=K_{i}^{\left( n \right)}{{\left( \frac{{{f}_{\text{L,}i}}}{{{f}_{\text{G,}i}}} \right)}^{\left( n \right)}}$
⑦进行误差判断，如果满足（49）式所示的条件，则输出气相、液相、吸附相中各组分的物质的量分数以及各相的相分率，否则重复步骤②—⑦直至计算结果满足条件。
${{\sum\limits_{i=1}^{N}{\left( \frac{{{f}_{\text{L,}i}}}{{{f}_{\text{G,}i}}}-1 \right)}}^{2}}<1\times {{10}^{-12}}$
2 模型验证
2.1 临界性质验证
为了测试修正后的PR状态方程的有效性，计算了无因次孔径与临界性质变化量之间的关系，并与前人实验数据[37-38]及MDS数据[21,43]进行对比。如图3 所示，随着无因次孔径的增加，临界性质变化量增大；无因次孔径小于0.2时，与前人实验数据及MDS数据相比，本文模型计算结果误差较小，模型具有一定的准确性。然而，当无因次孔径大于0.2时，模型计算结果误差较大，原因是此时孔隙实际尺寸过小，在极小的孔隙中，分子的尺寸接近或大于孔隙的尺寸，它们将无法自由移动。以CH4分子为例，它的有效直径通常约为0.38 nm，在无因次孔径大于0.2时，孔隙半径约为2 nm，此时CH4分子大部分将以吸附相存在孔隙内，几乎不存在流动性。
2.2 泡点压力验证
为了验证模型模拟纳米孔隙内流体相行为的可靠性，通过实验测试流体在纳米孔隙内的泡点压力，并与理论计算结果进行对比。
2.2.1 实验装置
设计并制作了具有纳米通道的流体芯片，即纳流控芯片，以模拟页岩孔隙结构。图4 展示了实验装置及纳流控芯片中通道的布局。该纳流控芯片尺寸为38.0 mm×21.5 mm，包含两个微米通道和多组平行的纳米通道。微米通道的两端分别为注入端和流出端，而纳米通道一端与微米通道垂直连接，另一端为密闭端，从而模拟地下纳米孔隙的形态。微米通道的设计尺寸为深4 μm、宽20 μm。针对典型的页岩纳米孔隙尺度，设计了2种纳米通道深度，分别为30 nm和100 nm，代表2 种孔径，同时设计了多孔介质和单管这两种类型通道，多孔介质的孔道长2 mm、宽30 μm，单管长2 mm、宽3 μm。为了更好地重现陆相页岩储层的润湿特性，通道表面润湿性被设置为油湿性。
实验装置由芯片夹持器、高分辨率可视化成像系统、温压控制系统及数据采集终端组成。其中，芯片夹持器配备蓝宝石玻璃，具备高光学清晰度和透明度，可同时承受20 MPa压力和150 ℃高温。高分辨率可视化成像系统包括光学显微镜和CMOS（互补金属氧化物半导体）相机，以实现对纳流控芯片中相行为的可视化和实时捕获。压力由泵系统控制，并通过芯片附近的压力传感器进行实时监测。温度控制采用油浴，并在芯片处配有温度传感器，油浴温度可达150 ℃。为确保实验结果的清晰度和可靠性，选择正构烷烃C4H10、C6H14、C8H18和C10H22配制油样，其比例为6∶5∶5∶4，通过活塞容器将其注入芯片内进行实验。选用该油样比例的目的是避免饱和芯片时流体气化造成过大误差。
2.2.2 实验步骤
为准确测量流体在纳米孔隙中的泡点压力，设计了相应的实验步骤，具体操作如下。
①实验开始前，将芯片置于夹持器内，检查实验管线的气密性，然后对装置进行4 h的抽真空，以排除芯片内的空气。
②设置压力为2~3 MPa（大于实验流体泡点压力），随后将实验流体缓慢注入纳流控芯片中。在显微镜下观察芯片的饱和状态。
③在显微镜下观察到芯片中的纳米通道已完全饱和流体后，设置油浴温度，使芯片所处环境升温至150 ℃。稳定20 min后，以0.1 MPa为间隔逐步降低压力。在每个压力下通过相机记录纳米通道中流体相态变化，等待系统稳定后再降低至下一个压力值。
④在显微镜下观察到纳米通道中出现第1个气泡时，压力传感器记录的压力即为泡点压力。随后，保持压力不变，记录纳米通道内流体相态变化。将拍摄的图像导入ImageJ软件进行图像处理。
2.2.3 实验结果与理论计算结果对比
在实验过程中，通过逐步降低压力并实时观察纳米通道内的流体相态变化，记录了不同压力下的气泡形成情况。图5 中可清晰观察到气泡的生成位置及其扩展过程，反映了流体在纳米孔隙中的相行为。当达到特定的压力时，纳米通道内观察到第1个气泡生成时，即标志着流体的相态转变。
利用本文模型分别计算了孔径为30 nm和100 nm时泡点压力随温度的变化曲线，如图6 所示。在30 nm深度的芯片中，单管与多孔介质测得泡点压力均为1.0 MPa，而实验温度（150 ℃）下的泡点压力本文模型计算结果为0.828 MPa，误差为20.8%；在100 nm深度的芯片中，单管与多孔介质测得泡点压力分别为1.4 MPa和1.3 MPa，而实验温度（150 ℃）下的泡点压力本文模型计算结果为1.271 MPa，误差分别为10.1%和2.3%。泡点压力实验测量结果与理论计算结果之间的误差处于合理范围内，验证了本文模型的可靠性。
3 模型计算及分析
3.1 单组分临界性质变化规律
从图7 中还可以看出，孔径的减小对不同大小和复杂性的分子影响各异。对于碳数少的烷烃，随着孔径的减小，临界温度的变化幅度相对较小，这是因为碳数少的分子整体尺寸小，受到孔隙结构复杂性的影响小。相反，对于碳数多的分子，随着孔径的减小，临界温度的变化幅度增大，这是因为相较于碳数少的烷烃分子，这些碳数多的烷烃分子整体尺寸变大，在相同的空间中限域效应更加明显，而限域效应极大地影响了分子间的相互作用以及分子的形状和排列，从而导致临界温度变化更加显著。同时孔径减小对于不同类型分子的临界压力影响也呈现出差异化特点。对于碳数较少的烷烃分子，随着孔径的减小，临界压力显著降低。这是因为这些分子结构简单、尺寸较小，在狭小的纳米孔隙中受到孔壁限域和分子间相互作用变化的影响更为显著，导致其临界压力变化较大。相比之下，对于碳数较多、分子结构更为复杂且尺寸更大的烷烃分子，孔径减小引发的限域效应对其临界压力的影响相对较小。这些较大分子在较大的孔隙中仍能保持较为自由的排布与旋转，分子间及分子-孔壁的相互作用偏离孔外体系状态的程度较低，从而使其临界压力的下降幅度较为有限。
随着孔隙尺寸的减小，纳米孔隙中的限域效应变得更加显著。这是因为流体分子与孔壁间的相互作用距离缩小，导致范德华力等分子间作用力的相对强度大增，这种增强的流体分子-孔壁相互作用显著削弱了流体分子间的吸引力，降低了达到相变所需的能量阈值。同时，随着孔径的减小，被吸附在孔壁上流体分子占孔隙内流体分子总数的百分比增加，使得孔隙内体相流体分子的有效摩尔体积进一步减小，限制了烷烃分子的运动自由度，使得系统内烷烃分子更容易聚集，进而在较低的压力下达到临界状态，降低了临界压力。流体分子运动受到的限制也随着孔径的减小而加剧，减少了分子间的有效碰撞，降低了分子间相互作用的能量，导致临界温度下降。特别是在孔径小于20 nm时，孔径与分子尺寸更为接近，分子无法像在大空间中那样自由移动，对分子的影响变得显著，进而改变分子的能量状态和分布，影响流体分子的相行为。
3.2 CO2-页岩油相行为分析
为了探究CO2-页岩油相行为，建立了A体系（CH4、C4H10、C6H14的比例为2∶1∶1）和B体系（CH4、C8H18、C10H22的比例为2∶1∶1）两种不同轻重的多组分混合体系。在A、B体系中加入CO2，探究CO2对多组分混合物相态特征的影响。各组分参数如表1 所示。
表1 各组分参数统计表 |
| 组分 | 相对分子 质量 | 临界 压力/MPa | 临界 温度/K | 偏心 因子 | 等张比容/ (N1/4•m11/4•mol-1) |
|---|---|---|---|---|---|
| CO2 | 44.010 | 7.382 | 304.19 | 0.198 | 82.00 |
| CH4 | 16.040 | 4.599 | 190.56 | 0.011 | 74.00 |
| C4H10 | 58.123 | 3.648 | 408.14 | 0.177 | 185.32 |
| C6H14 | 86.177 | 3.025 | 507.60 | 0.301 | 276.71 |
| C8H18 | 114.231 | 2.490 | 568.70 | 0.399 | 359.33 |
| C10H22 | 142.285 | 2.110 | 617.70 | 0.492 | 440.69 |
采用本文模型进行计算，并绘制了A体系和B体系添加CO2前后的相图（见图8 ），以及20 nm孔径条件下和体相条件下A体系、B体系添加不同物质的量分数CO2后的相图（见图9 ）。可以看出，随着孔径的减小，多组分混合体系的临界点移动，包络线包围的面积减小；加入CO2后，相同孔径下临界温度降低、临界压力增加，临界点在相图中向左上偏移，包络线包围的面积减小。此外，加入CO2前，A体系体相条件下与10 nm孔径条件下的临界温度之差为120.757 K，临界压力之差为4.193 MPa；加入CO2后，A体系体相条件下与10 nm孔径条件下临界温度之差为112.952 K，临界压力之差为4.859 MPa。加入CO2前，B体系体相条件下与10 nm孔径条件下临界温度之差为197.114 K，临界压力之差为3.720 MPa；加入CO2后，B体系体相条件下与10 nm孔径条件下临界温度之差为184.456 K，临界压力之差为5.305 MPa。这表明在限域效应影响下，CO2对临界温度的降低起到了抑制作用，对临界压力的降低起到了促进作用。此外，随着CO2物质的量分数的增大，临界点左移，包络线包围的面积减小。对比A、B体系，体系内重质组分含量越高，加入CO2后，随着孔径变小临界点偏移越明显，相图两相区变化越明显。
采用本文模型计算了A体系和B体系中加入CO2后临界性质的变化量，如图10 所示。可以看出，与A体系相比，不同孔径条件下B体系加入CO2后临界温度和临界压力的变化量都更大。这表明CO2的加入对重质组分含量更多的体系的临界性质影响更大。产生这一现象的原因是由于CO2是一种非极性分子，它与重质组分的非极性作用更强。重质组分有更长的碳链可以与CO2产生范德华力相互作用，更容易改变混合体系的热力学性质。
4 矿场应用
4.1 古龙A区块页岩储层流体相态特征
大庆古龙A区块页岩储层温度为379 K，原始地层流体压力为36 MPa，目前地层流体压力为13 MPa。流体中各组分相关参数如表2 所示。选取孔隙半径为10~1 000 nm，这一范围涵盖了目前古龙A区块页岩储层相关报道中的孔隙半径。虽然使用恒定的孔隙半径不能模拟真实储层中的实际情况，但是可以通过探究不同孔径下流体相行为的变化趋势对流体相态进行合理的推测。
表2 大庆古龙A区块页岩储层流体组分参数统计表 |
| 组分 | 物质的 量分数/% | 相对分子质量 | 临界 压力/MPa | 临界 温度/K | 偏心 因子 | 等张比容/ (N1/4•m11/4•mol-1) |
|---|---|---|---|---|---|---|
| CO2 | 6.03 | 44.01 | 7.382 | 304.19 | 0.198 | 82.000 |
| CH4 | 43.25 | 16.04 | 4.599 | 190.56 | 0.011 | 74.000 |
| C2—C4 | 23.01 | 42.82 | 4.310 | 363.30 | 0.143 | 145.200 |
| C6—C8 | 19.56 | 98.05 | 2.941 | 531.56 | 0.316 | 291.600 |
| C10+ | 8.15 | 200.00 | 2.180 | 788.74 | 0.687 | 683.300 |
4.2 相平衡状态下流体组分气液平衡常数变化
采用本文模型可以得到限域效应影响下大庆古龙A区块页岩储层流体达到相平衡状态时的气液平衡常数变化。图12 展示了100 nm以下孔径中各组分的气液平衡常数随孔径变化的曲线，由于孔径大于100 nm 时，气液平衡常数受到限域效应的影响较小并且接近于体相的状态，省略了孔径大于100 nm时的曲线。由图12 可知，CO2、CH4等较轻组分的气液平衡常数随着孔径的减小而减小，而C2—C4、C6—C8、C10+等较重组分的气液平衡常数会随着孔径的减小而增大。孔径10 nm时，各组分气液平衡常数均近似变为1，此时流体的气相组成与液相组成几乎相同。此外，随着孔径的减小，各组分气液平衡常数的变化幅度增大，说明气液平衡常数受到的限域效应影响会随着孔径的减小而增大。
为了探究CO2对页岩油相态的影响，进一步对液相和气相中CO2物质的量分数变化进行研究，如图13 所示，同样省略了孔径大于100 nm时的曲线。孔径小于100 nm时，随着孔径的减小，液相中CO2物质的量分数增加，气相中CO2物质的量分数减少。相比于孔径100 nm，孔径10 nm时液相中CO2物质的量分数增加了20.0%，气相中CO2物质的量分数减少了10.8%，这表明限域效应增强了纳米孔隙中CO2的传质。因此可以推测，限域效应有利于CO2的埋存和微观驱油。
5 结论
受到纳米孔隙内独特的限域效应影响，孔隙内流体的临界性质会随着孔径的减小而发生明显的偏移。当体系内加入CO2后，CO2的存在会改变系统的临界性质等相关物理参数的变化规律。
加入CO2后，相同孔径下临界温度减小、临界压力增大，临界点在相图中向左上偏移。在限域效应影响下，CO2对临界温度的降低起抑制作用、对临界压力的降低起促进作用。体系内重质组分含量越高，加入CO2后，随着孔径变小临界温度和临界压力的变化量越大，相图两相区变化越明显。
CO2、CH4等较轻组分的气液平衡常数会随着孔径的减小而减小，而C2—C4、C6—C8、C10+等较重组分的气液平衡常数会随着孔径的减小而增大。气液平衡常数受到的限域效应影响随着孔径的减小而增大。当孔径小于100 nm时，液相中CO2物质的量分数将随着孔径的减小而增加，相应地，气相中CO2物质的量分数将减少，限域效应有利于CO2的埋存和微观驱油。
符号注释：
a——能量系数，MPa·m6/mol2；A，B，C——压缩因子计算式中的系数，无因次；b——体积系数，m3/mol；c——表征流体分子与孔壁间相互作用的系数，MPa·m6/mol2；g（T）——与温度相关的函数；f——逸度，MPa；fA,i，fG,i，fL,i——吸附相、气相、液相中i组分的逸度，MPa；i，j——组分编号；kij——i组分与j组分之间的二元交互系数，无因次；kp——吸附平衡常数，MPa-1；Ki——i组分的气液平衡常数，无因次；L——圆柱形孔隙的长度，m；mabs——吸附量，kg；mmax——最大吸附量，kg；M——相对分子质量，无因次；n——迭代步序号；na——吸附在孔壁上的流体分子数；nt——孔隙中流体分子总数；N——混合体系的总组分数；NA——阿伏伽德罗常数，6.022 140 76×1023 mol-1；p——压力，MPa；pA——压力拆解后的引力项，MPa；pc——临界压力，MPa；pC——毛管压力，MPa；pcb——受限空间（纳米孔隙）中流体临界压力，MPa；Δp*——受限空间中临界压力变化量，MPa；pG——气相压力，MPa；pL——液相压力，MPa；pR——压力拆解后的斥力项，MPa；Pch,i——i组分的等张比容，N1/4·m11/4/mol；re——除去吸附层厚度的孔隙半径，m；R——气体常数，MPa·m3/（mol·K-1）；Rf——孔隙有效半径，m；Rp——孔隙半径，m；T——温度，K；Tc——临界温度，K；Tcb——受限空间（纳米孔隙）中流体临界温度，K；ΔT*——受限空间中临界温度变化量，K；Tr——拟温度，K；Va——孔隙内吸附分子的总体积，m3；Ve——分子的有效体积，m3；Vea——每个吸附分子的有效体积，m3；Vet——孔隙内体相流体分子的有效体积，m3；Vm——体相摩尔体积，m³/mol；${{{V}'}_{\text{m}}}$——有效摩尔体积，m³/mol；Vp——孔隙总体积，m3；Vreal——分子的真实体积，m3；xi，yi，wi——i组分在液相、气相、吸附相中的物质的量分数，%；zb,i——i组分体相物料投量，%；zf,i——i组分总体系物料投量，%；Z——压缩因子，无因次；δA，δG，δL——吸附相、气相、液相的相分率，%；ε——有效分子体积系数，无因次；εa——吸附分子的有效体积系数；εt——体相流体分子的有效体积系数，无因次；θ——气相与液相的接触角，（°）；λ——摩尔体积的修正系数，无因次；ρb——体相流体的密度，kg/m3；ρG，ρL——气相、液相的摩尔密度，mol/m3；σ——界面张力，N/m；σLJ——Lennard-Jones尺寸参数，m；τ——吸附层厚度，m；φ——逸度系数，无因次；φG,i，φL,i——气相、液相中i组分的逸度系数，无因次；ω——偏心因子，无因次。