考虑多重流动机制各向异性的页岩气运移数值模拟

赵静 古斌

摘要：为研究多重流动机制各向异性对页岩气运移的影响，基于Biot线弹性孔隙介质模型，考虑黏性流动、气体流动动态效应和表面吸附扩散等流动机制的各向异性，采用各向同性应力依赖模型，建立流-固耦合的视渗透率模型，对应力约束边界条件下的页岩气运移过程进行数值模拟，定量分析流动各向异性对页岩气运移的影响。结果表明各种流动机制的各向异性均对页岩气运移有显著影响：固有渗透率各向异性可明显改变页岩气运移开始时间和流通量大小，其影响随页岩储层渗透率增加而放大；当固有渗透率较小时，动态效应各向异性会削弱气体流动的动态效应，对页岩气运移开始时间和流通量的影响十分显著，但在固有渗透率较大时其影响可以忽略；表面吸附扩散各向异性的影响与固有渗透率各向异性相似。准确评估页岩气的运移能力和产能应考虑运移过程中的流动各向异性。

关键词：页岩气运移多重流动机制流动各向异性视渗透率模型数值模拟

中图分类号：TE3文献标志码：A文章编号：1671-8755（2023）01-0040-07

Abstract：BasedonBiotsmodeloflinearelasticporousmedium，takingintoaccountfortheanisotropyofflowmechanismssuchasviscousflow，dynamiceffectsofgasflowandsurfaceadsorptivediffusion，andadoptingtheisotropicstressdependencemodel，afluid-solidcouplingapparentpermeabilitymodelwasestablishedtoinvestigatetheeffectofanisotropyofmultipleflowmechanismsonshalegastransport.Numericalsimulationontheshalegasmigrationprocessunderthestressconstraintboundaryconditionswasthenconductedtoquantitativelyanalyzetheinfluenceofflowanisotropyonshalegastransfer.Theresultsshowthattheanisotropyofeachflowmechanismhassignificantimpactonshalegastransport.Theanisotropyofintrinsicpermeabilitycangiverisetoremarkablechangesonthebeginningofshalegastransferandthegasflux，andsuchinfluenceisenhancedbyincreasingtheshalepermeability.Theanisotropyofdynamiceffectscandiminishthedynamiceffectsofgasflowandhasasignificantinfluenceonthebeginningofshalegastransferandthegasfluxwhentheintrinsicpermeabilityislow.Whiletheinfluencecanbeneglectedinthecaseoflargeintrinsicpermeability.Theeffectofanisotropyofsurfaceadsorptivediffusionissimilartothatofintrinsicpermeabilityanisotropy.Therefore，itisnecessarytoconsidertheflowanisotropyinthemigrationprocesswhenaccuratelyevaluatingthemigrationabilityandproductivityofshalegas.

Keywords：Shalegasmigration;Multipleflowmechanisms;Flowanisotropy;Apparentpermeabilitymodel;Numericalsimulation

页岩储层构造复杂、孔隙类型和大小多样、承受的地应力高、呈现超低孔隙率和渗透率、各向异性和非均质等特性。页岩气是自生自储在页岩中的非常规天然气，其在储层中的传输是一个多尺度和多场耦合的复杂过程［1-2］。目前，利用连续介质力学理论建立视渗透率模型是研究页岩气运移规律的主要手段之一。现有的视渗透率模型主要通过不同程度地考虑黏性流动、气体分子与孔壁碰撞引起的气体流动动态效应、表面吸附扩散等多种流动机制和固体变形（即应力依赖）的综合影响，对页岩气运移的机制进行阐述和分析［3-8］。研究发现页岩气运移能力的改变取决于3个因素的相互竞争：表征黏性流动的固有渗透率、动态效应和表面吸附扩散。本课题组的前期研究［7］建立了一个流-固耦合的视渗透率模型，全面考虑固体变形和上述多种流动机制，分析了常应力假设、初始变形状态以及自由气与吸附气权重对页岩气运移能力的影响。在此基础上，进一步讨论了吸附效应引起的孔隙体积模量改变以及该变化对页岩气运移的影响［8］。但该模型与其他大多数视渗透率模型类似，假设与页岩运移相关的力学性质均为各向同性。

众多的实验研究表明，页岩的力学性质和页岩气流动，如渗透率、弹性模量、气体扩散等，都存在明显的各向异性。对不同的页岩，平行和垂直于层理方向的渗透率比值可以达到几至几百［9-10］；对近似横观各向同性的页岩，平行和垂直于各向同性面方向的弹性模量比值最大可达2.7［11］；气体扩散系数在平行和垂直于层理方向上表现出很大的差别［12-13］。数值模拟结果也揭示气体渗透率和扩散系数具有显著的各向异性［14-15］。现有的视渗透率模型中只考虑了黏性流动的各向异性，即采用各向异性固有渗透率。大多数研究关注由固体变形引起的固有渗透率各向异性。基于“火柴棒”等理想化的离散物理介质模型，通过建立孔隙率和固有渗透率与应变的关系，多种适用于裂缝型孔隙的各向异性应力依赖模型已被提出和应用［16-21］。在这些应力依赖模型中，储层、基体和孔隙可处理为具有各向同性或各向异性弹性参数的介质，但这些研究极少关注页岩气运移中的多重流动机制。也有学者［22-23］利用不同的渗透率主值和孔隙取向来描述固有渗透率各向异性，结合经验性的各向同性应力依赖模型［22］或忽略应力依赖性［23］，对页岩气的运移过程进行分析。尽管现有的视渗透率模型存在一定的不足，相应的研究结果均发现流动各向异性对非常规天然气的产气速率和产量有不可忽略的影响。鉴于固有渗透率各向异性的研究存在不足，頁岩气运移中动态效应和表面吸附扩散各向异性影响的研究至今尚未见报道，本工作在前期研究的基础上，利用固有渗透率、Knudsen数和表面吸附扩散系数的各向异性分别表征黏性流动、动态效应和表面吸附扩散等流动机制的各向异性特性，建立流-固耦合的视渗透率模型，对应力约束边界条件下的页岩气运移过程进行数值模拟，定量分析不同流动各向异性对产能的影响。

1视渗透率模型和气体运移控制方程

假设平面应变状态下的页岩气储层流动相关的性质为各向异性，固体变形相关的性质为各向同性。考慮黏性流动、动态效应和表面吸附扩散3种流动机制和固体变形的影响，前期建立的各向同性流-固耦合视渗透率模型［7］可扩展为以下各向异性形式：

数值计算中输入参数的取值见表1。为了表征固有渗透率和表面扩散系数主值的各向异性，引入参数tk=kx∞0/ky∞0，tad=Dxad0/Dyad0。当tk=tad=1时，退化为各向同性情况。为了讨论取向角的影响，本文考虑所有取向角相同的状态（固有渗透率和表面扩散系数的取向角可不一样），即Pk（θ）=δ（θ-θ0），Pad（θ）=δ（θ-θ′0）。

2.1固有渗透率各向异性的影响

为了研究固有渗透率各向异性对流通量的独立影响，假设吸附气流动相关的参数保持不变，相应的表面扩散系数取值为：Dxad0=10-6m2/s，tad=2，Pad（θ）=δ（θ-π/4）。同时设动态效应为各向异性，为了降低动态效应的影响，固有渗透率主值取较大值。图2和图3分别表示当kx∞0=5×10-19m2，tk=0.5时，页岩气运移过程中x方向和y方向的流通量在不同固有渗透率取向角情况下随时间的变化规律。

此外，图中还给出了当取向角均匀分布（即Pk（θ）=1/π）时流通量的变化。由图可知，页岩气的运移主要发生在压力变化的某段时间内，流通量的大小随

时间先增大至峰值，然后逐渐减小。运移开始时间随取向角变化，而且流通量的峰值大小也与取向角明显相关。在所研究的例子中，x方向的流通量在取向角为0时最大，取向角为π/2时最小。y方向的流通量在取向角为0时最小，取向角为-π/4时最大。同时注意到，随着取向角的变化，y方向的流通量发生流入和流出孔隙的转变。由于模型中假设右边界为井筒压力作用边界，因此x方向为页岩气的主要运移方向，y方向流通量的大小预期要远小于x方向流通量，这与数值结果相符。在后续讨论中，我们重点讨论主要运移方向上页岩气的产能，即x方向的流通量。图4表示当kx∞0=5×10-19m2，Pk（θ）=δ（θ-π/6）时，不同初始固有渗透率主值比对x方向流通量的影响。随着主值比的增加，页岩气运移开始越早，x方向流通量峰值越高。而且与各向同性情况（tk=1）相比，页岩气运移开始时间和x方向流通量大小的差异都随着初始固有渗透率主值各向异性程度增加而增加。综合图2和图4可知，固有渗透率各向异性对页岩气运移的产能有着重要的影响，而且各向异性的影响随着整体固有渗透率的增大而放大，如图5所示。图5给出了4种不同的固有渗透率各向异性程度当x方向初始固有渗透率主值由kx∞0=10-14m2增加为kx∞0=10-18m2时流通量的变化。结果显示实心标记曲线间的差异（kx∞0=10-19m2）远小于空心标记曲线间（kx∞0=10-18m2）的差异。

2.2表面吸附扩散各向异性的影响

为了研究表面吸附扩散各向异性对流通量的影响，假设动态效应为各向异性且自由气流动相关的参数保持不变，相应的固有渗透率等参数取值为

kx∞0=10-19m2，tk=0.5，Pk（θ）=δ（θ-π/6）。图6和图7分别给出了表面扩散系数取向角（Dxad0=6×10-6m2/s，tad=2保持不变）和表面扩散系数主值比（Dxad0=6×10-6m2/s，Pad（θ）=δ（θ-π/4）保持不变）对x方向流通量的影响规律。由图可知，页岩气运移开始时间和流通量的大小均受表面吸附扩散各向异性的影响。由于此时y方向的表面扩散系数主值大于x方向的主值，表面扩散系数取向角对流通量的影响规律与固有渗透率各向异性情况相反，如流通量峰值在取向角为0时最小，取向角为π/2时最大。当tad<1时，两个参数取向角的影响规律相似（相应图形省略）。表面扩散系数主值比的影响规律与固有渗透率各向异性情况相似，如主值比越高，流通量越大。比较图6与图2或图7与图4可以发现，表面吸附扩散各向异性引起的流通

量变化要小于固有渗透率各向异性，这是因为在计算中参数取值条件下，吸附气流动的贡献要小于自由气流动。此外，页岩储层表面扩散系数的增大也将放大表面吸附扩散各向异性的影响（相关图形略）。

2.3动态效应各向异性的影响

通过比较动态效应为各向异性和各向同性时x方向流通量的差别可分析动态效应各向异性的影响。由于表征动态效应的Knudsen数与固有渗透率大小相关，固有渗透率越小，动态效应越强。本文保持x方向初始固有渗透率主值和取向角不变（kx∞0=10-19m2，Pk（θ）=δ（θ-π/6），调整y方向初始固有渗透率主值的大小，考虑固有渗透率主值比tk=0.01，0.5，5等3种情况。此外，在数值计算中以下参数取值不变：Dxad0=6×10-6m2/s，tad=2，Pad（θ）=δ（θ-π/4）。图8显示了对不同程度的固有渗透率各向异性、动态效应分别处理为各向异性和各向同性时x方向流通量在页岩气运移过程中的演化。从图8可以看出，当固有渗透率主值比较大（tk=5）时，动态效应各向异性的结果与各向同性时几乎一致，表明此时动态效应各向异性的影响可以忽略。但是当固有渗透率主值比较小（tk=0.01）时，动态效应各向异性时的流通量远小于动态效应各向同性的情况，而且页岩气运移也更晚开始。同时注意到，当动态效应为各向同性时，流通量峰值大小随着tk取值的增加先减小后增加，这是由于动态效应与固有渗透率的影响相互竞争的结果。当动态效应为各向异性时，流通量峰值大小随着tk取值的增加而单调增加，说明固有渗透率较小时的动态效应受到抑制，从而改变了固有渗透率主值比的影响规律。以上结果表明动态效应各向异性对页岩气运移产能有明显影响。

頁岩气气藏产能的现场数据由于多种原因很难获取，而文献中也只提供储层的各向异性力学性质参数［9-13］，同时基于各向同性视渗透率模型的研究未能提供相应的储层各向异性参数，因此本工作未通过实例验证提出的考虑多重流动机制各向异性的视渗透率模型的有效性。但本文建立的视渗透率模型可退化为实例验证的各向同性模型［7］，而且分析中采用的各向异性参数均依据页岩储层的实验数据，所以本工作中的理论模型和定量分析结果具有合理性和有效性。

3结论

（1）基于Biot线弹性孔隙连续介质模型，利用张量形式的固有渗透率、Knudsen数和表面扩散系数分别表征黏性流动、动态效应和表面吸附扩散3种流动机制的各向异性，同时考虑固有渗透率的应力依赖性，发展了流-固耦合的视渗透率模型。利用数值方法，该新模型可对页岩气运移过程中的运移能力和产能进行分析。

（2）固有渗透率各向异性可明显改变页岩气运移开始时间和流通量大小，其影响随页岩储层渗透率增加而放大。表面吸附扩散各向异性的影响与固有渗透率各向异性相似。

（3）当固有渗透率较小时，动态效应各向异性会削弱气体流动的动态效应，对页岩气运移开始时间和流通量的影响十分显著，但在固有渗透率较大时，其影响可以忽略。

参考文献

［1］CURTISJB.Fracturedshalegassystems［J］.AAPGBulletin，2002，86（11）：1921-1938.

［2］GENSTERBLUMY，GHANIZADEHA，ROBERTCJ，etal.Gastransportandstoragecapacityinshalegasreservoirs—Areview.PartA：Transportprocesses［J］.JournalofUnconventionalOilandGasResources，2015，12：87-122.

［3］BESKOKA，KARNIADAKISGE.Amodelforflowsinchannels，pipes，andductsatmicroandnanoscales［J］.MicroscaleThermophysicalEngineering，1999，3（1）：43-77.

［4］JAVADPOURF.Nanoporesandapparentpermeabilityofgasflowinmudrocks（shalesandsiltstone）［J］.JournalofCanadianPetroleumTechnology，2009，48（8）：16-21.

［5］WUKL，LIXF，GUOCH，etal.Aunifiedmodelforgastransferinnanoporesofshalegasreservoirs：couplingporediffusionandsurfacediffusion［J］.SPEJournal，2016，21（5）：1583-1611.

［6］CAOP，LIUJS，LEONGY-K.Combinedimpactofflowregimesandeffectivestressontheevolutionofshaleapparentpermeability［J］.JournalofUnconventionalOilandGasResources，2016，14：32-43.

［7］YANGX，ZHANGHB，WUWB，etal.Gasmigrationinthereservoirsofultralowporosityandpermeabilitybasedonanimprovedapparentpermeabilitymodel［J］.JournalofPetroleumScienceandEngineering，2020，185：106614.

［8］YINGQ，ZHANGHB，ZHAOJ，etal.InfluenceofconsideringthesorptioneffectintheBetti-Maxwellreciprocaltheoremongastransportcapacityinunconventionalreservoirs［J］.TransportinPorousMedia，2021，137（3）：451-469.

［9］BHANDARIAR，FLEMINGSPB，POLITOPJ，etal.AnisotropyandstressdependenceofpermeabilityintheBarnettshale［J］.TransportinPorousMedia，2015，108（2）：393-411.

［10］PANZJ，MAY，CONNELLLD，etal.Measuringanisotropicpermeabilityusingacubicshalesampleinatriaxialcell［J］.JournalofNaturalGasScience&Engineering，2015，26：336-344.

［11］CHOJW，KIMH，JEONS，etal.DeformationandstrengthanisotropyofAsangneiss，BoryeongshaleandYeoncheonschist［J］.InternationalJournalofRockMechanics&MiningSciences，2012，50：158-169.

［12］WANY，PANZJ，TANGSH，etal.AnexperimentalinvestigationofdiffusivityandporosityanisotropyofaChinesegasshale［J］.JournalofNaturalGasScienceandEngineering，2015，23：70-79.

［13］李赟.C7页岩气藏吸附扩散微观机理研究［D］.成都：西南石油大学，2017.

［14］YINY，QUZG，ZHANGT，etal.Threedimensionalporescalestudyofmethanegasmassdiffusioninshalewithspatiallyheterogeneousandanisotropicfeatures［J］.Fuel，2020，273：117750.

［15］李滔，李闽，荆雪琪，等.孔隙尺度各向异性与孔隙分布非均质性对多孔介质渗透率的影响机理［J］.石油勘探与开发，2019，46（3）：569-579.

［16］WANGGX，MASSAROTTOP，RUDOLPHV.AnimprovedpermeabilitymodelofcoalforcoalbedmethanerecoveryandCO2geosequestration［J］.InternationalJournalofCoalGeology，2009，77（1/2）：127-136.

［17］PANZJ，CONNELLLD.Modellingofanisotropiccoalswellinganditsimpactonpermeabilitybehaviourforprimaryandenhancedcoalbedmethanerecovery［J］.InternationalJournalofCoalGeology，2011，85（3/4）：257-267.

［18］GUFG，CHALATURNYKR.Permeabilityandporositymodelsconsideringanisotropyanddiscontinuityofcoalbedsandapplicationincoupledsimulation［J］.JournalofPetroleumScienceandEngineering，2010，74：113-131.

［19］WUY，LIUJS，DEREKE，etal.Developmentofanisotropicpermeabilityduringcoalbedmethaneproduction［J］.JournalofNaturalGasScienceandEngineering，2010，2（4）：197-210.

［20］WANGJG，LIUJS，KABIRA.Combinedeffectsofdirectionalcompaction，nonDarcyflowandanisotropicswellingoncoalseamgasextraction［J］.InternationalJournalofCoalGeology，2013，109：1-14.

［21］WANGK，ZANGJ，WANGGD，etal.Anisotropicpermeabilityevolutionofcoalwitheffectivestressvariationandgassorption：Modeldevelopmentandanalysis［J］.InternationalJournalofCoalGeology，2014，130：53-65.

［22］HUANGT，TAOZW，LIE，etal.Effectofpermeabilityanisotropyontheproductionofmultiscaleshalegasreservoirs［J］.Energies，2017，10（10）：1549.

［23］單娴，姚军.基于渗透率张量的各向异性油藏两相渗流数值模拟［J］.中国石油大学学报，2011，35（2）：101-106.

［24］AL-YOUSEFHY.Permeabilityanisotropymeasurementonwholecores—analyticalsolutionandapplication［C］∥14thSPEMiddleEastOil&GasShowandConference.Bahrain，2005.