三种强度准则的盐穴应力弹塑性解及参数敏感性分析*

张华宾,张顷顷,王来贵

(辽宁工程技术大学 力学与工程学院，辽宁 阜新 123000)

0 引言

盐穴储库为人们在深部盐层中采用水溶造腔方式形成的地下溶洞，可作为石油、天然气等的储存场所。在深埋地下的盐岩层中造腔势必产生人为扰动，在一定范围内引起地应力重新分布。造腔过程会在盐穴周边形成一个弹塑性分区，甚至导致围岩的破坏。

目前，研究人员已对地下空间硐室平面应变问题涉及的弹塑性理论进行了广泛研究。研究人员通过采用不同的岩体本构关系和强度准则，取得了诸多成果[1-11]。例如，倪国荣[1]推导了地下球穴围岩的弹塑性区应力和位移公式,同时给出了喷射混凝土支护抗力的表达公式；张浪等[5]基于M-C准则, 推导了深埋球洞开挖后应力弹性分布、塑性分布和塑性区半径的公式；王武[6]研究了盐穴及含夹层圆柱储气库围岩弹性应力分布的完整解析解并结合H-B破坏准则计算了储库极限内压值范围。由于空间问题的理论求解涉及复杂的数学推演过程，研究人员仅仅采用M-C，D-P等较为简单的准则进行了研究，涉及盐穴储库的弹塑性研究较少见到。此外，盐穴造腔过程采用水溶方式，必定会对硐室围岩造成侵蚀作用，导致盐岩力学性能降低。鉴于此，本文将盐穴储库简化为地下球型硐室，分别采用M-C，D-P准则和考虑岩体结构特征与开采扰动的H-B准则，结合空间平衡方程对盐穴储库围岩进行理想弹塑性分析，并研究塑性半径对准则参数的敏感性。该研究对盐穴储库工程设计和施工具有理论参考作用。

1 盐穴储库简化模型

考虑我国可作为盐穴造腔的地层均在地表1 000 m以下。为简化研究，假设球型围岩处于均质、各向同性的理想弹塑性状态。地层应力趋于均匀分布。设模型边界在无穷远处且受均匀压力p0，盐穴设计的半径为R0，将油气对围岩的作用简化为作用于洞壁表面的均匀压力pi。暂不考虑重力作用与侧压系数影响。以盐穴中心为原点，满足球对称条件，设塑性半径Rp，则盐穴围岩Rp≤r≤∞为在弹性区，在R0≤r≤Rp为塑性区；在半径r=Rp处弹塑性分界面上满足连续性条件。图1为球型盐穴硐室围岩的弹塑性模型。

图1 球型盐穴围岩的弹塑性模型Fig.1 The elastic-plastic model of spherical salt cavern and its surrounding rock

2 强度准则

目前广泛应用的广义H-B准则[9]：

(1)

式中：σci为完整岩石的单轴抗压强度,MPa；mb，s，a为H-B常量,为无量纲。2002年，研究人员引入岩体扰动系数D，并提出了基于地质强度指标GSI的参数mb，s，a取值的新方法。

(2)

式中：GSI为岩体地质强度指标；D作为工程扰动因素的岩体弱化因子，取决于外界因素对原位岩体的扰动程度，D=0时为未受扰动状态；mi反映岩石软硬程度的岩性参数，岩石取值为0～25；s，a为与岩体特征有关的参数，其中0≤s≤1，大多数完整岩石时a=0.5。

3 盐穴球型硐室围岩弹塑性力学分析

3.1 弹性区分析

设主应力与球坐标应力对应关系为σ1=σ2=σθ=σφ,σ3=σr。在Rp≤r≤∝为弹性区，其平衡方程为：

(3)

根据拉梅解，可直接得到盐穴围岩弹性区应力：

(4)

式中：σRP为弹、塑性交界面上弹性区一侧的径向应力。

3.2 H-B准则下盐穴围岩塑性区分析

在塑性内外边界满足条件：

(5)

(6)

为确定塑性半径Rp，考虑在弹塑性交界面上的连续性条件，弹性区一侧的应力满足：

σr+2σθ=3p0

(7)

结合塑性区的应力表达式(6)和式(7)可得：

(8)

(9)

(10)

将式(10)代入式(6)，取r=Rp可确定塑性半径：

(11)

式中：σRp为弹塑性交界面上径向应力。由式(11)知，塑性区半径与H-B准则参数σci,mi,s,a有关，还与盐穴所在地层原始地应力p0及储库运行压力pi有关系。

3.3 M-C，D-P准则下盐穴围岩塑性区分析

同前一节的推导思路，限于篇幅。参考文献[5],M-C，D-P准则应力解直接给出推导结果，应力解依次为：

(12)

(13)

塑性半径Rp的解依次为：

(14)

(15)

由式(14)、(15)知，M-C，D-P准则的塑性区半径与盐岩黏聚力c、内摩擦角φ均有关，与原始地应力p0的正相关，与运行内压pi的负相关。

4 参数估计与塑性区半径对参数的敏感性

4.1 参数估计

H-B准则考虑了岩体结构、岩块强度及应力状态等多种因素，弥补了M-C，D-P准则的不足，并且主应力的非线性关系更适于处理深埋高应力状态的围岩弹塑性分析。然而应用H-B准则处理岩土工程问题时，参数选择的正确性很关键。其中参数GSI值可以根据岩体结构特征和结构表面条件估计，Cai等[16]提出了基于块体尺寸Vb和结构面条件因子Jc的量化GSI围岩分级方法，该方法使GSI值的确定更具客观性。

参数σci，mi可通过多组室内三轴试验数据采用式(16)、(17)进行估计。

y=miσcix+sσci

(16)

(17)

式中：x=σ3,y=(σ1-σ3)2；较完整岩石可假设s=1，xi，yi为同组试验数据。N为试验总次数。

为了进行3种不同准则下应力解的比较，根据Hoek、Brown的研究，可以将H-B准则与M-C准则进行等效转化，获得盐岩的等效粘聚力c′、内摩擦角φ′。其等效公式为：

(18)

4.2 塑性半径对H-B准则参数的敏感性

塑性半径与各准则参数、盐穴设计指标均有关系。基于单因素分析法，考虑H-B准则参数单一变化时对盐穴围岩塑性区半径的影响情况。

考虑参数GSI，mi，σci及D在2，10 MPa运行压力下塑性半径变化过程。由图2可知，随着GSI，mi，σci这3个参数的增大，Rp/R0减小，均呈衰减的指数形式降低。同等条件下，参数GSI，mi均在较低水平，运行压力由2 MPa增大到10 MPa时，Rp/R0变化较为显著，Rp/R0下降近约一半，但随参数的增大变化幅度趋于平缓。当参数GSI，mi均在较高水平，运行压力由2 MPa增大到10 MPa时，Rp/R0变化较不显著，这种变化随参数值的增大更不明显。根据Rp/R0数值可知参数GSI，mi相对于σci对计算结果的影响更为显著。参数σci的增大使得Rp/R0逐渐变小，变化趋势呈衰减指数形式。运行压力的增大使得Rp/R0减小。同等条件下，低内压运行时，Rp/R0随参数σci增大的降低程度最为显著，如在2 MPa运行压力下，参数σci由8 MPa增大到24.4 MPa时，Rp/R0下降30%。较低σci值时，运行压力由2 MPa增大到10 MPa时，Rp/R0下降28%。参数D的增大使得Rp/R0同时提高，近似为线性递增，但Rp/R0增幅程度相对其他情况并不显著。运行压力的增大使得Rp/R0减小。同等条件下，运行压力由2 MPa增大到10 MPa时，Rp/R0最大下降18%，即为D=0.8时的情况。在较低运行压力时，参数D由0增大至0.8时，Rp/R0

最大增幅6%，即为运行压力为2 MPa时的情况。综上，可以说参数σci，D对计算结果的影响不显著。

若讨论2个因素及以上组合对计算结果的影响，可采用多因素分析法。正交试验法是1种常见的多因素分析方法，其利用数理统计手段，挑选出具有代表性的因素组合合理地安排试验，可以有效减少试验次数，并且能够较全面反映各因素对试验指标的影响程度。现讨论H-B准则参数GSI，mb，D，σci这4个因素对结果的影响关系。将各量化为低、中、高3种水平(见表1)，忽略准则参数之间的交互作用。选择L9(34)正交表安排试验及方差分析结果如表2所示，其中，方差分析中各因子总偏差平方和为5 167.89；总自由度为8。取α=5%，根据表2，查F分布分位数表[17]F0.05(2,8)=4.459，根据表2可知H-B准则参数组合对盐穴围岩塑性区半径的无显著影响。综上，塑性区半径对H-B准则各参数的敏感性排序为：GSI>mi>σci>D。

表1 正交设计水平因素Table 1 Orthogonal design level factor

图2 H-B准则参数与Rp/R0关系曲线Fig.2 H-B criterion parameters and Rp/R0 relationship curve

试验方案GSIDmiσci塑性半径Rp11(30)1(0.2)1(8)1(20)114.294 321(30)2(0.5)2(16)2(50)61.497 1231(30)3(0.8)3(24)3(80)61.436 0942(50)1(0.2)2(16)3(80)35.638 9952(50)2(0.5)3(24)1(20)49.987 8762(50)3(0.8)1(8)2(50)71.267 1573(70)1(0.2)3(24)2(50)32.396 7583(70)2(0.5)1(8)3(80)36.224 8193(70)3(0.8)2(16)1(20)46.967 37K1j79.07660.77773.92970.417K2j52.29849.23748.03455.054K3j38.53059.89047.94144.433Rj40.54611.54025.98825.984偏差平方和2 550.580247.4661 345.9101 023.933自由度2222F值1.9740.1921.0420.793敏感性排序1423

5 算例分析

在国内某拟建地下盐穴储库库区，通过现场钻井取芯分析，通过加载多组三轴压缩试验，实时记录试验全过程应力-应变曲线。试验仪器为长春朝阳仪器厂出产的RLW-2000三轴蠕变试验机(图3所示)。通过式(16)、(17)、(18)对15，16，20，22，25和30 MPa这6种围压下的32组试验数据进行分析，得到参数mi=4.1，σci=24.4 MPa，c=6.83 MPa，φ=34.8°。

为了验证理论解的正确性和有效性。设盐穴储库半径R0=30 m，远场压力为p0=40 MPa。模型涉及岩体力学计算参数见表3。取运行内压分别为0，2，4，6，8和10 MPa的5种工况研究塑性区半径与盐穴半径之比随内压的变化规律，如图4所示。由图4可知，运行内压的增加会使得塑性半径变小，H-B准则下呈指数形式下降，M-C，D-P呈线性降低。其中，D-P(内切圆)结果与M-C相近度最好。H-B准则考虑岩体结构性质的影响，计算塑性半径显著增大。提高H-B准则的GSI参数，塑性半径将会显著变小，这也进一步说明GSI参数对结果影响较为显著。根据当前工程设计条件下的计算结果分析，运行压力8～10 MPa之间，采用M-C，D-P(内切圆)准则时，盐穴围岩将从弹塑性状态过度为完全弹性状态，运行内压达到10 MPa时盐穴围岩内壁不出现塑性区。而H-B准则依然出现塑性区，但塑性区范围显著缩小，10 MPa时，岩体质量GSI=90对应盐穴围岩塑性半径33 m，仅仅在盐穴内腔壁周边3 m范围内出现塑性区。若岩体质量GSI下降到60，对应的围岩塑性区半径增大到44.8 m，即在盐穴内腔壁周边的14.8 m范围内都进入塑性状态。由此可见，在使用H-B准则进行工程设计时，参数的选取更为重要。

图3 试验机及部分标准试样Fig 3 Tester and standard samples

表3 计算参数Table 3 Calculation parameters

图4 Rp/R0随内压的变化曲线Fig 4 Change curves between the ratio of Rp/R0 and the internal pressure

图5 不同屈服准则盐穴应力分布规律比较Fig.5 Comparison of stress distribution law of salt caverns with different yield criteria

图6 H-B准则盐穴围岩的弹塑性分布规律Fig.6 Elastic-plastic distribution law of surrounding rock of salt cavern with H-B criterion

由图6可知，GSI参数的变化导致围岩应力水平有显著差异。 GSI为60时，主应力极值为1.2p0，GSI为90时，主应力极值为1.3p0且围岩应力水平与M-C结果更为接近。GSI值从60变90后，塑性半径降低30%左右。主应力差值σ1-σ3极值提高到原来1.5倍以上。

6 结论

1)将盐穴储库简化为空间球体更符合自然条件，论文给出了M-C，D-P，H-B准则条件，考虑运行内压的盐穴球型硐室的弹塑性应力解析解和塑性半径计算公式，通过算例验证了正确性。

2)采用统计线性回归方法计算盐岩的H-B准则mi，σci值。并通过等效M-C强度理论确定了盐岩的等效粘聚力和等效内摩擦角。

3)塑性半径随H-B准则参数GSI，mi，σci的增大而减小，随着D的增大而增大，影响程度依次为GSI，mi，σci，D。参数组合对盐穴围岩塑性区半径的影响并不显著。塑性半径也与原始地应力p0及运行压力pi有关系。

4)采用H-B准则研究发现，微小内压的波动使得塑性半径的变化更强烈，岩体质量较差时尤为显著。这也更加符合工程的实际情况。M-C，D-P准则相对于H-B准则计算塑性半径偏小，而应力却偏大。