断层对盐岩储库流-固耦合及稳定性的影响王贵君;李常皓;刘朝鹏;李东;蔡晟鑫;张岭【摘要】层状盐岩中天然气储库围岩存在节理裂隙或断层,可能会增大天然气渗漏的风险.为研究断层对盐岩天然气储库流-固耦合特性及稳定性的影响,建立了有断层的盐岩天然气储库数值模型,分析储库围岩的扩容、渗流深度、塑性区扩展范围和流变变形.结果表明,围岩断层附近围岩应力、应变集中和扩容现象明显,渗透率增大;在地应力和渗流压力(洞室天然气内压)共同作用下,断层岩体渗流深度、塑性区深度和变形均比其他部位大;加大天然气内压渗流深度增大,但可提高围岩稳定性.研究成果对遇断层盐岩储气库的设计和安全运行具有一定的指导作用.【期刊名称】《河北工业大学学报》【年(卷),期】2018(047)003【总页数】8页(P93-100)【关键词】断层;天然气储库;流-固耦合;稳定性;流变变形【作者】王贵君;李常皓;刘朝鹏;李东;蔡晟鑫;张岭【作者单位】河北工业大学土木与交通学院,天津300401;河北工业大学土木与交通学院,天津300401;河北工业大学土木与交通学院,天津300401;华夏幸福基业股份有限公司,北京100028;河北工业大学土木与交通学院,天津300401;河北工业大学土木与交通学院,天津300401;河北工业大学土木与交通学院,天津300401【正文语种】中文【中图分类】TU4570 引言与国外海相沉积形成大而厚的盐丘不同，我国盐岩多为湖相沉积而成，多呈层状结构分布，并含有泥岩、石膏等夹层.在水溶造腔建库及储库注采气过程中，储库围岩发生流变和渗流，在泥岩与盐岩接触面上这种流变和渗流作用尤为突出.一旦盐岩天然气、石油储库遭遇倾斜断层，储库围岩的流-固耦合作用和流变变形及稳定性等工程特性更为复杂.20世纪50年代以来，世界上发生了28起较大的盐岩储库工程事故，有23起是由于油气泄漏和腔体失稳造成的.特别是存在穿越盐岩层的断层，由于蠕变速率不同，接触面产生“拖拽”致使剪切破坏，成为高压流体的渗流通道.可见，对存在不良结构面或断层的盐岩储库渗流特性及稳定性的研究有着十分重要的科学与工程意义.Suppe研究了断层与褶皱形态之间的几何关系，建立了断层运动发展模型[1].Stead和Szczepanik通过试验监测盐岩材料中局域源快速释放能量产生的瞬态弹性波来研究盐岩蠕变损伤的过程，并建立了相应的流变损伤强度理论[2].Cristescu基于弹-黏塑性理论研究了盐岩的剪胀性，进而分析盐岩的流变损伤过程[3].杨春和等[4-5]、陈卫忠等[6]、刘江等[7]、梁卫国等[8]、陈剑文等[9]和其他学者对盐岩变形特征进行了一系列的研究，通过分析含夹层盐岩的孔隙率与渗透率之间关系发现，变形前盐岩及夹层结构致密，不会对盐岩腔体的安全性产生影响；但施加荷载后，因夹层和相邻盐岩体的力学特性差异使二者变形不协调，会造成接触面处发生剪切破坏；据此提出了层状盐岩储库极限运行压力的确定原则.本文作者在Carter模型基础上引入“损伤增速界限”，建立了一种新的流变损伤模型，能够很好地描述盐岩流变-损伤的全过程、特别是第三阶段加速流变-损伤的特性[10].近年来，我国开展了大量盐岩储库方面理论与数值模拟研究，也获得了可喜成果.但是，人们对有断层盐岩储库的研究还很少.本文采用有限差分软件FLAC3D建立有断层含泥岩夹层盐岩天然气储库数值模型，从流-固耦合和流变机理出发，分析断层岩体应力、应变集中，研究围岩渗透率变化、渗流深度、塑性区及流变位移.研究成果对盐岩储库设计和安全运行具有一定的指导意义.1 Kozeny-Carman方程的应用描述岩体渗透率K与孔隙率Φ变化之间关系的Kozeny-Carman方程为[11]式中：Kz为常数，一般取Kz=5；Sp为孔隙比表面积.由式（1）知，渗透率的改变与孔隙率的变化之间的关系为式中：K0为初始渗透率；Φ0为初始孔隙率.不考虑介质中固相的变形，结合小变形假设，孔隙率Φ与体积应变εv的关系近似为[12]或取式中，体积膨胀时体积应变εv取正值，体积压缩时取负值.最后，渗透率变化与体积应变的关系为或需要说明的是，式（3a）与式（3b）所得计算结果，式（4a）与式（4b）所得计算结果均十分接近.2 盐岩储库数值模型2.1 储库模型由于含断层盐岩储库模型较为复杂，为提高建模效率，选择在有限元软件ANSYS 中进行三维建模，导入FLAC3D软件中，再进行数值分析.鉴于实际情况，储气库腔室采用鸭梨形.取30°扇形范围作为研究对象，盐岩储库顶板位于地下1 129 m 深处，模型总高1 320 m，扇形半径为600 m，储库腔室高为160 m，最大半径为64 m.腔室范围内由上到下含有厚度分别为4 m和2 m的2个泥岩夹层，距洞室顶部分别为91 m和128 m.一个倾角为50°的断层同时穿过上覆岩层、盐岩和第一层泥岩夹层，如图1所示.建模所用坐标系中，模型宽度方向右向为x正向，竖直向上为z正向.模型边界条件如下：2个竖向扇形对称面为法向位移约束边界，在扇形体弧形侧面上作用法向原岩应力，模型下表面（z=－611m的平面）为竖向位移约束边界，在模型上表面（z=709 m的平面）作用上覆岩土层压力13.55 MPa.根据计算工况不同，储库腔室内边界作用天然气压力不同，分别为6 MPa、10 MPa和15 MPa.2.2 岩体本构模型与渗流模型采用塑－粘弹性模型——pWIPP模型[13]来描述盐岩力学特性.偏应变速率的粘性部分与偏应力张量sij同轴，表达式为图1 数值模型：储库腔体位置（左）和FLAC3D网格图（右）Fig.1 Numerical model式中式中，ε˙p与ε˙s分别为过渡和稳态蠕变速率，按下式计算：式中：D（s－）1，n，A（s－）1，B（s－1）和εs－1）均为材料参数；σ*=1MPa为参考应力；R=8.314 4 J·mol－·1K－1为气体热力学常数；Q （kJ·mol－1）为材料的激活能；T（K）为材料的绝对温度.pWIPP模型中的屈服条件由Drucker－Prager屈服条件（D－P条件）及抗拉屈服条件组成，与渗流作用耦合后，屈服条件为式中：σm为平均应力；p为渗流压力；qφ、qk为材料参数.假设模型中上覆泥岩和储库范围内夹层泥岩为符合D－P条件的弹塑性材料，断层岩体为符合Mohr－Coulomb条件（M－C条件）+不抗拉的弹塑性材料，下卧泥岩层为弹性材料.对于盐岩和泥岩的渗透性，均采用各项同性连续介质渗流模型来研究.因FLAC3D 软件没有直接通过体积应变计算材料孔隙率及渗透率的功能，本文利用软件中FISH语言编制了相应的计算程序，按照式（4）实现Kozeny－Carman方程的应用，随材料的体积应变变化而改变材料的渗透率，从而实现流－固耦合.另外，在数值分析过程中，储库内天然气的密度和模量随储库内压的变化而改变.FLAC3D中的渗透系数与一般渗流力学中的渗透率是有区别的，需要换算后使用[13].2.3 岩体力学参数根据实验成果并参考Asse盐矿（德国）盐岩的流变实验数据[14-15]，盐岩材料流变参数为D=5.79 × 10-36，Q=56 kJ/mol，A=4.56 s-1，B=127 s-1，ε˙*ss=5.39×10-8s-1，n=5.0，T=300 K.盐岩、上覆泥岩、夹层泥岩、下伏泥岩以及断层岩体的弹塑性参数如表1所示. 表1 岩体弹塑性力学参数Tab.1 Elasto-plastic parameters of rock masses岩层强度准则剪切模量G/GPa 体积模量K/GPa qφ qk/MPa 抗拉强度σt/MPa盐岩D-P 5.71 26.67 0.231 4.20 1.0泥岩夹层 D-P 3.85 8.33 0.273 4.68 1.0上覆岩层D-P 35.89 57.05 0.273 3.51 1.0下伏岩层 Elastic 56.20 78.16断层 M-C 35.89 57.05 φ=25° c=2.0 03 部分计算结果根据含断层盐岩天然气储库的特点，本文采用“先加载后开挖”的模拟顺序进行计算，即先设定模型的边界条件和初始条件，计算达到初始平衡状态；然后采用水溶法造腔形成储库腔体，封井后注入天然气排驱腔体内卤水，使井口天然气压力达到设计运行压力.实际计算中，假定造腔瞬时完成.本文仅选择代表性计算成果，从断层附近围岩应力、应变集中、储库围岩渗透率、渗流、塑性区和位移等方面研究含断层盐岩储库天然气内压分别为6 MPa、10 MPa、15 MPa时流-固耦合特性和围岩稳定性及流变变形特性.需要说明的是，由于泥岩夹层与相邻盐岩体之间联结坚固、紧密，计算中未将这些层理界面作特殊处理，各物理场在界面上是连续的.3.1 断层附近应力、应变集中盐岩的强流变性使得盐岩层中的稳定初始应力场呈静水压力状态，稳态流变过程中的储库围岩应力场也具有抽水井周围水压力特征，这些结论与作者先前的研究成果相同[16-17].断层本质上属于软弱不连续体，它的存在对储库围岩应力分布产生较大影响.图2a）为天然气内压为15 MPa时运营30年后的储库围岩最小应力（即σ3，拉为正，压为负）分布局部放大云图，在储库洞室周边断层附近有明显的应力集中；图2b）为同一条件下储库围岩剪应力σxz分布局部放大云图，剪应力集中现象愈加明显.与应力集中相对应，储库围岩断层附近出现应变集中，图3为体积应变和剪切应变分布局部放大云图，其中体积应变为正表示体积膨胀增大，即产生扩容现象.尽管剪应力集中区域较大，达数十米深，剪切应变集中的区域仍局限于几米、十几米深度范围.图2 断层附近围岩应力分布局部放大图Fig.2 Zoom shows of the stress distribution图3 洞室周边围岩应变分布局部放大图Fig.3 Zoom shows of the strain distribution3.2 围岩体积应变与渗透率图4 为盐岩储库天然气内压为Pi=6 MPa时腔体内壁围岩几个典型部位体积应变随时间的变化曲线.储库断层和夹层附近围岩的体积应变为正值，且数值较大，表明这些位置发生了扩容，也与图3a）结果一致；但在储库运营之初的前两年内体积应变就达到最大值，只是断层附近围岩体积应变随时间仍缓慢增大，储库运营25年后体积应变达到2.9‰（图4曲线4）.在天然气内压和流变作用下，储库顶、底板围岩仍处于压缩状态，体积应变为负值（图4曲线1、2）.由式（4）知，孔隙介质渗透率随体积应变增大而增大，将此式应用于数值分析中，就实现了流-固耦合.尽管式（4）形式上为体积应变的非线性函数，围岩渗透率变化（即K/K0）与体积应变的关系也可以近似采用线性函数来描述，这与作者先前的研究成果结论相同 [18].图5是储库内压Pi=6 MPa运行25年后储库周边围岩渗透率局部放大云图，储库周边断层附近岩体的渗透率比其他部位的渗透率大两个数量级，断层渗透率增大的深度（深色区域）约为28.5 m.尽管夹层泥岩附近有轻微应力集中现象（图2），但此处的体积应变和剪应变与相邻盐岩的数值相比差别不大（图3） .3.3 储库渗流在储库腔室内天然气压力和围岩应力共同作用下，流-固耦合作用使天然气向围岩深部渗透流动，围岩中的天然气渗流压力不断变化.图6为储库内压Pi=6 MPa时腔体内壁围岩几个典型部位渗流压力pp随时间的变化曲线.接近腔室顶板处（点1）的渗流压力逐渐增大，5年后达到最大值5.7 MPa（曲线1）.点3位于第二泥岩夹层下部的腔室内壁，这点的渗流压力即储库腔室天然气内压，自始至终保持6 MPa.点2和点4均处于洞壁深部，距洞壁分别为12.7 m和8.3 m，天然气渗流达到此两点的时间分别为6.4年和2.5年，计算到28年时，这两点的渗流压力分别为2.10 MPa和3.56 MPa.当天然气内压Pi=15 MPa时，图6中腔体内壁围岩的几个典型部位渗流压力pp 随时间的变化趋势与Pi=6 MPa时的相同，只是渗流压力值随内压增大而增大.这仍然是储库腔室围岩流-固耦合作用的结果.虽然储库腔室天然气压力增大了很多，但围岩所处应力状态为受压，且使围岩产生扩容的等效应力更小（或者不存在），围岩渗透率不但不增大，反而会有所减小，渗流速率不会增大很多.断层及受其影响围岩孔隙率较大，初始渗透率也较大，随围岩流变变形及渗流共同作用渗透率增大较多（图5），渗流深度也比其他位置大，如图7和图8中较深颜色区域所示.Pi=6 MPa时，断层处渗流深度（渗流影响范围）约为59 m，其他位置渗流深度较均匀，大约为32～39 m（图7）.Pi=15 MPa时，断层处渗流深度约为64 m（2倍腔室最大半径），其他位置渗流深度差别也不大，大约为36～45 m（图8）.Pi=10 MPa时的渗流深度介于二者之间.图9为储库天然气内压为15 MPa时围岩几个部位渗流深度随时间的变化曲线，渗流深度与时间的关系可近似采用线性函数描述.例如，储库运行5年到30年间储库腔室顶板渗流深度Dp （m）与时间t（a）的关系可拟合为拟合系数R=0.988 1，说明吻合很好.因夹层泥岩的渗透率也较低，其数值仅比相邻盐岩大一个数量级，在不同腔室天然气压力作用下的渗流速率及渗流深度与相邻盐岩的数值相比差别也不大（图7、图8）.3.4 围岩塑性区以储库天然气内压Pi=6 MPa和15 MPa两种工况为例，图10为储库运营30年后围岩塑性区分布的局部放大云图，图中彩色阴影部分为储库围岩塑性单元，其中包括前期屈服、而后又变为弹性的单元.需要说明，由于盐岩的特殊性质，前期发生屈服的单元有可能在流变作用下转化为弹性状态；在FLAC3D中，材料的屈服是单元（zone）的整体状态，因而，单元划分精密程度影响塑性区范围的计算精确性.与不考虑渗流作用和断层影响的情况类似，在储库运行初期阶段，洞室围岩中盐岩与夹层泥岩塑性区的深度基本一致，泥岩中塑性区深度甚至略小；随着盐岩的蠕变变形与渗流压力共同作用，盐岩体中的应力势能逐步转化为变形能，洞周及附近的应力集中得以释放，盐岩中塑性区缓慢减小[19].在数值计算过程中没有考虑泥岩的流变特性，泥岩层中水平应力与竖向应力之间的差值较大，即屈服等效应力数值较大，相邻盐岩体较大的流变变形对泥岩产生“拖拽”作用，使泥岩层中塑性区范围不断扩大，内压较小时（例如Pi=6 MPa）甚者可贯穿整个模型，图10a）中第2个泥岩夹层就属于这种情况，塑性区沿夹层方向扩展深度达90 m.图4 储库围岩体积应变变化曲线（Pi=6 MPa）Fig.4 Volume strains of surrounding rock vs.time图6 储库围岩渗流压力变化曲线（Pi=6 MPa）Fig.6 Seepage pressures of surrounding rock vs.time图8 渗流压力分布（Pi=15 MPa，30 a）Fig.8 Distribution of the seepage pressure图5 围岩渗透率（Pi=6 MPa，25 a）Fig.5 Permeability of surrounding rock 图7 渗流压力分布（Pi=6 MPa，30 a）Fig.7 Distribution of the seepage pressure图9 储库围岩渗流深度变化曲线（Pi=15 MPa）Fig.9 Seepage range ofsurrounding rock vs.time断层岩体强度比其他岩体强度低，尤其是抗剪能力更弱，此处的（剪）应力集中使断层岩体率先屈服，塑性区逐渐扩展，30年时塑性区沿断层方向扩展深度达161 m，超出局部放大图形之外，相邻岩体的塑性区深度比其他位置的大.需要注意，因断层很薄，图10a）中阴影不明显.储库天然气内压Pi=15 MPa时，断层屈服深度大大减小，洞周盐岩体的塑性区范围减小，夹层泥岩塑性区深度也大大减小，仅限于洞壁10 m深度范围，如图10b）所示.3.5 围岩位移储库围岩在地应力和腔室天然气内压共同作用下产生流变变形.如果内压较小，围岩向腔室内部的收敛变形会较大；储库长时间保持在低压状态下运行，会使腔室空间迅速减小，降低储库存储能力和使用年限.适当增大储库天然气内压，可减小围岩收敛变形，提高储库存储能力和使用寿命[19].在考虑渗流和断层作用条件下，上述规律依然成立.图11左图为天然气内压Pi=6 MPa时储库洞周几个部位位移d随时间的增长情况，这些曲线充分反映了盐岩的流变特性.图11中、右图分别为Pi=6 MPa和15 MPa 时洞周围岩位移d分布及位移矢量的局部放大图.显然，断层处的位移远大于其他部位的位移，这一点是非常突出的.增大储库天然气内压，洞壁位移可随之减小.当Pi=6 MPa时，储库腔室底板起鼓达到2.67 m，侧壁除断层处外位移值为小于或接近1 m.当Pi=15 MPa时，底板起鼓及侧壁位移均小于0.5 m，顶板下沉量约为0.7 m，断层处的位移仍达到2.03 m. 图10 储库围岩塑性区分布Fig.10 Distribution of plastic zones图11 围岩位移随时间的发展（Pi=6 MPa，左）和30 a后位移分布（中：Pi=6 MPa；右：Pi=15 MPa）Fig.11 Displacements of surrounding rock：histories and distributions4 结论1）在地应力和储库天然气内压联合作用下，在储库洞室周边断层附近有明显的应力、应变集中，流-固耦合作用使此处围岩孔隙率增大，从而使渗透率增大，渗流深度比其他部位大，塑性区深度大，流变位移也大.断层属于储库围岩最薄弱部位，在此处容易发生渗漏、大变形甚者失稳破坏.因此，在选择储库场址时，必须尽量避开断层.2）在较大天然气内压条件下，流-固耦合作用使围岩中渗流深度加大，当Pi=15 MPa时，储库运营30年后渗流深度达两倍洞室最大半径.在多洞室储库群建设时，为避免相邻腔室渗流联通，必须保证足够大的洞室间距，建议大于2倍洞室最大直径.3）较低储库天然气内压时渗流深度较小；但是，围岩的收敛变形和塑性区开展深度均很大，泥岩夹层和断层的变形和塑性区深度尤为突出.当Pi=6 MPa时，储库运营30年后断层塑性区开展深度超过洞室最大半径的5倍.长时间维持储库低压运行，不利于储库经济运行，降低储库存储能力，减小储库使用寿命.较高储库天然气内压时渗流深度较大，但围岩流变变形和塑性区开展深度均较小，可提高储库围岩稳定性，提高储库运行经济性，提高存储能力，延长使用寿命.参考文献：【相关文献】[1] Suppe K R.Geometric relationship between fault pattern and fold shape[J].Rock Mech，1997，34（4）：646-650.[2] SteadD，SzczepanikZ.Time-dependentacousticemissionstudiesonpotash[C]//RoegiesD.Rock Mechanics as a Multidisciplinary Science.1991，471-479.[3] Cristescu N，Hunsche U.Time effects in rock mechanics[M].Chichester：Wiley，1998.[4] 杨春和，陈峰，曾义全.岩盐蠕变损伤关系研究[J].岩石力学与工程学报，2002，21（11）：1602-1604.[5] 杨春和，梁卫国，魏东吼，等.中国盐岩能源地下储存可行性研究[J].岩石力学与工程学报.2005，24（24）：4409-4417.[6] 陈卫忠，伍国军，戴永浩，等.废弃盐穴地下储气库稳定性研究[J].岩石力学与工程学报，2006，25（4）：848-854.[7] 刘江，杨春和，吴文，等.盐岩短期强度和变形特性试验研究[J].岩石力学与工程学报，2006，25（S1）：3104-3109.[8] 梁卫国，杨春和，赵阳升.层状盐岩储气库物理力学特性与极限运行压力[J].岩石力学与工程学报，2008，27（1）：22-27.[9] 陈剑文，杨春和，郭印同.基于盐岩压缩-扩容边界理论的盐岩储气库密闭性分析研究[J].岩石力学与工程学报，2009，28（S2）：3302-3308.[10] 王贵君.一种盐岩流变损伤模型[J].岩土力学，2003，24（S）：81-84.[11]Kleinberg RL，Flaum C，Griffin DD，et al.Deep sea NMR：methane hydrate growth habit in porous media and its relationship to hydraulic permeability，deposit accumulation，and submarine slope stability[J].Journal of Geophysical Research，2003，108（B10）：2508-2522.[12]Chen Yifeng，Hu Ran，Lu Wenbo，et al.Modeling coupled processes of non-steady seepage flow and non-linear deformation for a concrete-faced rockfill dam[J].Computers and Structures，2011，89：1333–1351.[13]Itasca Consulting Group，Inc.FLAC3D，Fast Lagrange Analysis of Continua in 3 Dimension（Version 3.1）User's Manual[Z].Minneapolis，2003.[14]StaudtmeisterK，Rokahr R B.Rock mechanics design of storage caverns for natural gas in rock salt mass[J].Int J Rock Mech Min Sci，1997，34（3/4）：646.[15]Hunsche U，Schulze O.Das Kriechverhalten von Steinsalz[J].Kali und Steinsalz，1994，11（8/9）：238-255.[16]王贵君.盐岩层中天然气存储洞室围岩长期变形特征[J].岩土工程学报，2003，25（4）：431-435.[17]Wang Guijun，Guo Kangmei，Christianson Mark，et al.Deformation characteristics of rock salt with mudstone interbeds surrounding gas and oil storage cavern[J].International Journal of Rock Mechanics&Mining Sciences，2011，48（6）：871-877.[18]王贵君，刘朝鹏，介少龙，等.考虑流-固耦合的天然气盐岩储库渗流特性研究[J].地下空间与工程学报，2016，12（增2）：470-474，480.[19]王贵君.盐岩油气储库与核废料处置库力学问题研究新进展[M].北京：科学出版社，2013.。