双孔双渗介质多场耦合理论
（1）多孔介质模型描述
对于多孔介质可能存在的高渗通道-裂缝，一般视为离散和连续介质方法进行模拟。然而，模拟现场尺度的裂缝岩石似乎从数学和物理上是不可能的（例如，68m3 岩石从图像统计估算有可渗透的5 百万条裂缝）[27-28]。因此，使用连续介质方法模拟裂缝发育地层似乎更为合理。
首先，多孔介质岩石存在不同孔径分布的孔隙结构，如何去模拟他们？由于多孔介质内部的复杂孔隙结构，需要在代表单元体（REV，Representive VolumeElement，存在均匀的孔隙结构）进行研究分析。Biot 指出，REV 与包含的孔隙尺寸相比足够大，但与我们所感兴趣的宏观现象尺度相比又足够小，因此REV 才是均匀多孔介质的代表性单元。当然，研究局限于井眼附近的井壁稳定问题，就不能将大尺度裂缝，如断层等的归入此问题的REV。针对含有裂缝发育的多孔介质的本构模型，根据孔隙类型的不同，已从单重孔隙介质[13-15]扩展至双重孔隙介质[1-3,8,22-26,29]，三重孔隙介质[30]甚至多重（>3）孔隙介质[16-22]。目前多重孔隙介质多场耦合理论已考虑了“流固耦合”现象[19]。
Bai 等人首先基于混合物理论（mixture theory）定义了多重孔隙度-多重渗透率介质用于描述变形-流动耦合过程，并且针对裂缝性地层特性给出了不同重介质模型的适用范围[17]，如下表1.1。之后，Bai 将多孔介质孔隙分为宏观，介观，微观孔，并研究了其溶质运移过程，得到溶质扩散的半解析解[30]。
上表1.1 中，对于典型天然裂缝性地层，当裂缝相渗透率远高于基质相渗透率，才可以使用双孔双渗（双重孔隙度-双重渗透率）模型。但2016 年，Liu 和Abousleiman 对超低渗化学活性页岩的压力传递试验进行数值模拟，表明即使粘土相与非粘土矿物（两相）渗透率仅相差2 倍，双重孔隙介质模型相比单重孔隙介质模型更符合实验结果。因此，为更好模拟多孔介质的孔隙弹性响应，简化双重孔隙理论可能导致误差产生[29]。
Mehrabian 博士论文[18]系统介绍了多重孔隙度-多重渗透率介质理论的构建过程，如图1.2 所示，并将其应用于典型地质力学问题，如CO2 存取，井壁稳定，地层压力衰竭等问题。此后，Mehrabian 和Abousleiman 考虑热动力学势能密度建立多重孔隙度-多重渗透率孔隙介质广义Biot 理论[19]，根据其中单重介质参数确定总体介质孔隙力学参数，并获得其N 重介质的本构模型和Gass mann 方程关系式[20]。
此外，并将其应用于Mandel[18]和Cryer[21]问题的孔隙弹性响应解，即突然施加载荷后多孔介质的孔隙弹性响应，得到其问题的解析解，分析了多重孔隙之间的孔隙压力变化和流体质量交换过程。
Liu 和Abousleiman 等人认为页岩存在不同层次的孔隙介质（见图1.3），进而使用多重孔隙度/渗透率模型分析页岩地层井壁稳定问题，获得了N 重孔隙介质流固耦合的地层孔隙弹性响应解，相比于单重孔隙介质计算结果，双重乃至三重介质模型导致更高的坍塌压力[22]。
多重孔隙介质理论问题在于需要更多参数来描述每重孔隙介质的力学强度，输运（渗透率）参数以及窜流参数等等。Berryman 指出对于N 重孔隙介质（N>2），基于各个孔隙介质相参数叠加无法得到总体介质力学参数，N 重孔隙介质的参数获取仍是一个待解决问题[31]。最终，Mehrabian 和Abousleiman 提出一种基于能量密度函数唯一性原理的方法，解决了总体介质力学参数计算方程的闭合问题[32]。
国内研究方面，陈勉和陈至达[33]、赵颖等[34-35]、陈勉等[36]、蔡新树等[37]使用边界条件叠加方法分别建立各向同性、各向异性单相、多相饱和双重和多重孔隙介质有效应力定律。高彦芳、陈勉等[38]提出了考虑基质吸力的多相非饱和多重孔隙介质的有效应力定律。
（2）双孔双渗介质模型
由于多重孔隙介质模型的复杂性，目前双重孔隙-双重渗透率多孔介质理论应用较为广泛。双重孔隙介质理想模型首先由Barenblatt[39]提出，Warren 和Root[40]给出了模型参数获取方法且多应用于裂缝性油藏渗流、试井等问题。其假设为基质孔隙度大，渗透率低，为主要流体储存空间；而裂缝孔隙度低（一般小于5%）[2]，渗透率大，为主要流动通道。Meng 等人模拟裂缝性地层为双重孔隙-单重渗透率介质，假设基质渗透率为0，得到井周应力场半解析解，并与常用的双孔双渗模型结果对比，发现此模型钻井安全密度窗口更为保守[41]。罗毅等人使用双孔双渗模型分析致密气藏裂缝岩心渗透率，得到其非稳态测试流动特征规律[42]。
对于双重孔隙连续介质的模拟，首先要引入理想裂缝间距的概念[2,43-49]。图1.4中三维方格代表基质，其长度即为其理想裂缝间距。需要指出的是，此为概念模型，为复杂裂缝的等效参数，其与基质与裂缝之间的流体质量交换相关。
Gelet 等人对基于混合物理论的双孔双渗孔隙介质模型的使用范围、理论假设、参数获取方法等做了详细归纳[46]。双孔双渗介质的本构模型框架通过基质和裂缝相的平均化处理得到，不同模型的区别在于平均方法，包括有效介质理论（The effective medium approach），混合物理论（The mixture theory），归一化理论（The homogenization theory）和统计学方法（The statistical approach）,但不管何种平均方法，其本构模型表达式完全一致。Gelet 等人推荐基于混合物理论的模型推导方法，优势之一在于此模型直接可退化至单重孔隙介质。
基于混合物理论，Aifantis[50]，Villiappan 和Khalili[43-45]，Berryman、Wang[51-52]以及Pride[53-54]分别建立双重孔隙介质本构模型。Aifantis 首先构建双重孔隙双重渗透率本构模型，但忽略了基质和裂缝相相互作用，即裂缝流体压力变化不会对基质孔隙压力产生影响[50,55]。但裂缝流体压力、基质孔隙压力、骨架应力通过耦合变形存在相互影响[46,50]。
Berryman 和Wang 首先建立双重孔隙-双重渗透率介质的唯象本构模型，使用试验方法确定总体介质的6 个独立模量参数[51]。唯象理论是根据实验事实得到物理规律（本构模型），但缺乏对参数的相应物理过程的解释。因此，两人[52]使用Biot 单重孔隙介质模型推导方法-热力学方法，得到其相应本构模型系数的物理解释，并将其应用于裂缝性介质声波数据解释。Berryman 和Pride[53]建立裂缝（节理）相为100%流体，或者含有充填材料（裂缝相骨架）的双重孔隙介质本构模型，并通过实验方法和理论推导获得相应参数。之后，Pride 和Berryman 建立各向同性双重孔隙-双重渗透率线性动力学控制方程，基于体积平均方法，用于分析声波在双重孔隙介质的传播特征[54]。
Nguyen 和Abousleiman[55]使用Berryman 等人建立的模型，将总体介质分为两个隔离但重叠的子介质（基质和裂缝相）：各自含有各自的骨架和流体相。Nguyen将各向同性模型扩展至横观各向同性模型，最后得到双重孔隙介质流-固，流-固-化，流-固-热耦合模型的井周应力场解析解[1]。其相应本构模型系数通过基质和裂缝相参数的加权平均得到，最常用的是体积加权平均。模型中裂缝相的（占总介质）体积分数很小（大约1%）且孔隙度接近1（如95%），且裂缝相刚度很小（假设为基质1/50）。
Villiappan 和Khalili 考虑孔隙变形建立裂缝性多孔介质非线性本构方程，相应系数变化可通过实验方法得到。通过数值模拟发现，对于较硬地层，线性本构模型可得到满意结果；对于大变形地层，线性假设会引起显著误差[43]。
Biot 对单重孔隙介质有效应力系数的定义得到学术界的广泛认同，多重孔隙度介质为有效应力定义带来了新问题[48]。今天有两个有效应力定义仍然存在：（1）基于总体介质分为多个子孔隙介质，其对应多个（重）有效应力，如Berryman 和Wang[52]和Abousleiman 等人[1,55];（2）单重孔隙介质的Biot 有效应力理论扩展至双重孔隙介质，对于实体骨架只有一个有效应力，Khalili 和Valliappan[43-45]，Gelet 等人[46-49]，Liu 和Abousleiman 等人[22,56-59]。Cheng[56]指出应使用多重有效应力概念，可分别考虑流体压力分别对基质和裂缝相的影响，基质、裂缝相、总体介质的破坏失效并不一定一致。对于多重有效应力概念，一般使用太沙基（Terzaghi）有效应力系数[55]；单介质有效应力概念需要考虑相应相对总体介质的影响程度，其不等于单重孔隙介质内的有效应力系数。学者Abousleiman 等人针对多重孔隙介质模型的井壁稳定分析使用单有效应力概念较多[22,57-59]。因此，本文使用单有效应力概念进行井壁稳定分析。
为考虑热，化学势，电势等对多孔介质的影响，只考虑多孔介质内骨架应力-应变，流体压力-流体含量变化的流固耦合本构模型已扩展至多场耦合的本构模型。
例如考虑热因素，温度变化引起多孔介质内部的热应力和流体压力发生变化；化学离子浓度变化，同样对固相应力，流体压力产生影响[60]。电势变化通过影响离子浓度分布，进而影响多孔介质孔隙弹性响应[61]。
（3）多场耦合输运理论
耦合输运方程考虑多种驱动力下，如流体压力、温度、化学势和电势梯度[62-64]下多孔介质内部流体流量、溶质浓度，电流量、热流量的实时变化过程。Yeung 和Mitchell[62]基于体积平均方法建立流量与驱动力的线性本构关系（输运方程），如下表1.2。
基于考虑因素种类的区别，结合固体变形，目前已有流固（渗流-应力），流固化（渗流-应力-化学），流固热（渗流-应力-温度），流固热化（渗流-应力-温度-化学），流固化电（渗流-应力-化学-电势）等输运模型。
Yeung 和Mitchell 基于非平衡热动力学理论，建立等温条件下流体流量，离子流量和电流量与水力压力梯度，电势梯度和化学势梯度的输运方程，并给出了输运方程矩阵系数[62]。之后Nguyen 和Abousleiman[64],Tran 和Abousleiman[65]得到横观各向同性含粘土页岩斜井井眼问题的流固化电响应的半解析解，主要考虑由于粘土带负电，导致页岩可移动孔隙流体整体带正电。由于Donnan 平衡，井壁界面阴阳离子浓度与钻井液内浓度并不一致，并对井周应力场和井壁稳定产生影响[64-65]。
此后Liu 和Abousleiman 建立双孔双渗流固化电耦合模型，并与Ewy 等人页岩岩心流体压力传递实验结果进行验证，发现双孔双渗流固化电耦合模型的解相比单重孔隙介质模型解更接近实验结果，同样试验结果可反过来确定页岩孔隙力学参数[29]。
此外，Revil 和Leroy 基于体积平均的局部Nernst-Planck 和Stocks 方程跨尺度建立多孔页岩的离子输运方程[66]，考虑了粘土带负电性和Donnan 平衡导致含粘土多孔介质表面阴阳离子浓度不同的物理过程。此外阴阳离子扩散速度不一致，建立了以阴阳离子浓度，流体压力，温度，电流量为变量的输运方程。
常龙等人基于Revil 和Leory 的模型建立了非饱和页岩膜效率的计算方法，发现随着水饱和度增加，页岩膜效率减小。且CEC（阳离子交换容量）、孔隙流体溶质浓度都会影响非饱和页岩膜效率[67]。
最终，沉积地层存在力学和强度性质的各向异性特征[68-71]，结合单重至多重孔隙介质的本构模型和多场耦合输运方程[60,72]，构建了均匀介质各向异性多重孔隙介质多场耦合模型，应用于不同地层与工况。本文为分析裂缝性地层钻井全周期井壁稳定，分别针对双孔双渗介质流-固耦合动力学，流-固-热耦合和流-固-热-化耦合模型方面展开调研。
