可燃冰沉积物宏细观力学特性真三轴试验离散元模拟

2020-02-27王宏乾周世琛薛世峰

中国石油大学学报（自然科学版） 2020年1期

周博, 王宏乾, 王辉, 周世琛, 薛世峰

(中国石油大学(华东)储运与建筑工程学院,山东青岛 266580)

可燃冰是甲烷气体和水在低温高压下形成的笼型类冰结晶化合物,广泛存在于陆域永冻区、极地大陆架和大陆架边缘的深海沉积区[1-4]。在标准大气压下,1 m3可燃冰分解可产生164 m3甲烷气体和0.87 m3水,是21世纪最具有商业开发前景的战略资源[5]。2017年中国在南海神狐海域实现了可燃冰试采成功,极大推动了中国在可燃冰开发领域的研究[6-7]。可燃冰开采需要以破坏可燃冰的相平衡条件为前提,这容易引起井壁失稳、海底滑坡及地层变形等一系列安全问题[8-12]。研究可燃冰沉积物的力学行为、揭示可燃冰储层的变形规律,对实现可燃冰安全开采具有重要意义。天然可燃冰沉积物试样不易提取和保存,目前主要采用人工合成可燃冰沉积物试样来进行室内试验研究[13-14],但室内试验难以揭示可燃冰沉积物力学行为的细观机制。离散元法(DEM)在计算散粒体和表征沉积物几何特征方面存在优势[15]。贺洁等[16]制备了可燃冰沉积物DEM试样,并对其进行了固结排水三轴压缩DEM模拟试验,研究了可燃冰沉积物的宏微观力学性质。Brugada等[17]利用离散元软件PFC3D对孔隙填充型水合物沉积物进行了常规三轴排水模拟试验,分析了孔隙填充型结构对力学特性的微观影响机制。肖俞等[18]建立一种微观胶结型可燃冰沉积物DEM试样,对其进行了常规三轴排水试验模拟研究。Jung等[19]针对可燃冰在孔隙中的不同生长阶段,对孔隙颗粒型以及块状团簇饱和型的可燃冰沉积物进行了离散元模拟,进一步揭示了可燃冰形成环境对沉积物力学性质的影响。上述常规三轴室内试验和DEM模拟试验均不便于描述复杂应力状态下的可燃冰沉积物力学行为。贺洁等[20]对特定饱和度的孔隙填充型可燃冰沉积物,开展了不同中主应力系数下真三轴压缩DEM模拟试验,分析了中主应力对可燃冰沉积物力学行为的影响。蒋明镜等[21]提出了裹覆型可燃冰沉积物离散元模拟方案,对裹覆型可燃冰沉积物试样进行真三轴压缩试验模拟,探讨了裹覆型可燃冰沉积物在复杂应力状态下的宏微观力学特性。这些研究主要集中在中主应力对可燃冰沉积物的力学特性的影响,而对于可燃冰沉积物的细观和微观力学性质的演化以及其与宏观力学性质的联系,仍有待进一步研究。笔者采用离散元软件PFC3D5.0制备孔隙填充型的可燃冰沉积物DEM试样,对其开展一系列在不同平均应力、不同中主应力系数下的真三轴压缩试验模拟,分析平均应力、中主应力系数对可燃冰沉积物强度、变形特征以及摩擦角的变化规律的影响,并初步探讨中主应力系数对细观和微观力学特性的影响。

1 可燃冰沉积物DEM试样制备

根据文献[22]中所提出方法生成可燃冰沉积物DEM试样。图1是饱和度为10%的可燃冰沉积物DEM试样。试样为边长3 mm的立方体,采用圆球颗粒模拟土体颗粒和可燃冰颗粒,其中浅黄色颗粒表示土颗粒,粒径为0.1～0.4 mm,为更好地与室内试验结果相比较,其颗粒级配曲线与Toyoura砂类似(图2),密度为2.65 g/cm3;深蓝色颗粒表示可燃冰颗粒,密度为0.9 g/cm3,考虑到可燃冰颗粒由于土体孔隙生成及计算效率等因素,将可燃冰颗粒粒径取为0.06 mm。

生成特定饱和度可燃冰沉积物DEM试样的方法如下。

(1)可燃冰颗粒“转化”为土颗粒,二者同时生成。与文献[16]-[17]中试样生成方法不同,首先将特定饱和度的可燃冰颗粒通过计算转化为土颗粒的一部分,然后将两种颗粒统一由新生的颗粒级配曲线生成。这种做法既体现出可燃冰颗粒填充的随机性,又节省了试样生成时间,提高了计算效率。可燃冰饱和度是指可燃冰体积占试样孔隙体积的占比,即

Smh=Vmh/Vv.

(1)

式中,Vmh和Vv分别为可燃冰体积和孔隙体积, mm3。

(2)消除重叠量,控制不平衡力。为了保证模拟的准确性,试样生成后对整体颗粒进行再平衡运动,使颗粒间的不平衡力小于某一特定值,消除生成颗粒产生的重叠量。

(3)试样固结。为模拟可燃冰沉积物初始的应力状态,在试样上施加1.0 MPa的固结压力。

(4)赋予颗粒接触模型。土颗粒间接触赋予线性接触模型,可燃性颗粒间、可燃冰与土颗粒间赋予平行胶结接触模型。线性模型的法向和切向刚度均为1.5×107N/m,摩擦系数为0.5。平行黏结模型的法向和切向胶结刚度为均为1.5×107N/m,黏结强度和拉伸强度均为8.3 MPa,墙体的法向和切向刚度均为100 MN/m。

图2 DEM模拟的颗粒材料和Toyoura砂级配曲线Fig.2 Particle size distributions of granular material in DEM simulations and Toyoura sand

2 可燃冰沉积物真三轴压缩试验DEM模拟分析

2.1 常规三轴压缩试验DEM模拟

为了验证真三轴DEM试验模拟的有效性和可靠性,利用编制的真三轴伺服程序对制备的孔隙填充型可燃冰沉积物试样进行固结排水三轴压缩模拟试验,并将试验模拟结果与室内试验结果进行对比。图3给出的是饱和度为40.9%的可燃冰沉积物试样在有效围压1.0 MPa下的DEM模拟结果与Masui等[23]对可燃冰沉积物在有效围压1.0 MPa条件下开展的固结排水三轴压缩室内试验结果对比。

图3 饱和度为40.9%的水合物沉积物离散元数值模拟试验与室内试验对比Fig.3 Comparison between DEM simulation and laboratory test for hydrate-bearing sediments with saturation of 40.9%

由图3可以看出,DEM数值模拟结果能较好地反映可燃冰沉积物的主要应力-应变特性,数值模拟结果与试验结果吻合较好,从而说明将本文中提出的可燃冰沉积物离散元试样制备方法应用于真三轴压缩试验是可行的。

2.2 真三轴压缩试验DEM模拟

土力学中将平均应力p和中主应力系数b定义为

(2)

(3)

式中,σ1、σ2与σ3分别为第一、第二、第三主应力,MPa;p为平均应力,MPa;b为中主应力系数;当b=0时,σ2=σ3,对应于常规三轴压缩状态;当b=1时,σ2=σ1,对应三轴拉伸状态。

联立式(2)、(3),中主应力σ2和小主应力σ3可以由已知的平均应力p、中主应力系数b和大主应力σ1表示,使得σ2和σ3在真三轴压缩过程中可以监测得到

(4)

(5)

加载边界条件的控制是成功进行可燃冰沉积物真三轴压缩DEM模拟试验的关键,通过编写PFC真三轴伺服控制程序对已制备好的可燃冰沉积物试样进行伺服加载并实时监测大主应力σ1,通过计算求出式(4)、(5)中的中主应力σ2、小主应力σ3与大主应力σ1的关系,将计算得到的σ2、σ3重新施加到DEM模型加载板上,保持整个过程中b不变。整个试验过程采用恒应变加载方式,应变率为7%/min。

真三轴压缩试验得到的结果是用偏应力q、体应变εv和偏应变εs的关系来表示的,其中对于偏应力q、偏应变εs和体应变εv的定义如下:

(6)

(7)

εv=ε1+ε2+ε3.

(8)

式中,ε1、ε2、ε3分别为大、中、小主应变,无量纲;q为偏应力,MPa。

真三轴DEM模拟试验主要分为两组:①平均应力p保持定值,中主应力系数b取不同数值0、0.2、0.4、0.5、0.6、0.8、1.0对特定可燃冰饱和度的沉积物试样进行真三轴压缩试验模拟;②中主应力系数b保持定值,平均应力p分别取1.0、1.5、2.0、2.5 MPa,对特定可燃冰饱和度的沉积物试样进行真三轴压缩试验模拟。

3 离散元试验模拟结果

3.1 平均应力取值影响

为了研究平均应力p对DEM模拟结果的影响,在中主应力系数b取值不变的情况下取不同平均应力进行真三轴压缩试验模拟。图4给出了中主应力系数取值相同,平均应力取值不同情况下偏应力与偏应变的关系曲线。图5为偏应力进行归一化后的偏应力比q/p与偏应变εs关系。图6为可燃冰沉积物归一化峰值偏应力比与中主应力系数的关系。由图4可知,当b不变时,可燃冰沉积物的偏应力-偏应变曲线呈弹塑性特性,并且平均应力取值的变化并不影响偏应力-应变曲线的形态,而可燃冰沉积物的偏应力峰值强度以及残余强度会随着平均应力取值增大而增大。由图5、6可知,平均应力的取值对归一化的偏应力比q/p影响较小,b取值相同时不同平均应力下归一化后的偏应力峰值强度几乎相同,说明峰值应力面是过应力空间原点的原锥面。

图4 不同平均应力下可燃冰沉积物的偏应力-偏应变曲线Fig.4 Deviatoric stress-deviatoric strain curve of combustible ice sediments with different mean effective stresses

图5 归一化的偏应力比与偏应变关系Fig.5 Normalized deviatoric stress ratio-deviatoric strain relationship for combustible ice sediments

图6 可燃冰沉积物归一化峰值偏应力比与中主应力系数的关系Fig.6 Relationship of normalized peak stress ratio and intermediate stress ratio for combustible ice sediments

3.2 中主应力系数影响

3.2.1 应力-应变关系

图7 不同中主应力系数对可燃冰沉积物主应力的影响Fig.7 Effect of intermediate stress ratio on principal stress for combustible ice sediments

在主应力系数b的取值对可燃冰沉积物的宏细观力学特性的影响的真三轴压缩试验模拟中,取平均主应力p=1.0 MPa。图7给出了在中主应力系数取值不同时3个主应力σ1、σ2、σ3与偏应变εs的关系曲线。如图7(a)所示,在真三轴压缩过程中大主应力σ1随着偏应变增大而增大,达到峰值强度后会出现轻微的应变软化现象,并且σ1的峰值强度随b增大而减小;由图7(b)可以看出,σ2随着b取值增加表现出截然相反的变化趋势,当b处于0～0.5时,中主应力σ2会随着偏应变εs增大而减小,而其峰值强度则会随着b增大而增大;而当b处于0.5～1.0时,σ2随着εs逐渐增大而增大,且其峰值强度随着b增大而增大;而当b=0.5时,σ2与εs的关系是一条保持在1 MPa的水平直线,属于一条过渡曲线;由图7(c)可以看出,在变形前期小主应力σ3随着εs增大而减小,在εs=0.15时取得峰值,而后σ3会随着εs逐渐增大而略微增大;其峰值强度随着b增大而减小。图8为不同主应力系数时,偏应力q与偏应变εs的关系曲线,图9为偏应力峰值强度与中主应力系数的关系曲线。偏应力-偏应变曲线在变形前期受到b的影响较小,当偏应变εs大于0.1时,中主应力系数对偏应力-偏应变曲线的影响较大。由图9可以看出,b处于0～0.6时,偏应力峰值强度随着b逐渐增大而迅速减小;在b逐渐增大超过0.6之后,偏应力峰值强度趋于一个定值。

图8 不同中主应力系数对可燃冰沉积物的偏应力-偏应变曲线的影响Fig.8 Effect of intermediate stress ratio on deviatoric stress-deviatoric strain curve for combustible ice sediments

图9 可燃冰沉积物偏应力峰值强度随中主应力系数变化Fig.9 Variation of peak strength with intermediate stress ratio for combustible ice sediments

3.2.2 体变关系

图10 不同中主应力系数对可燃冰沉积物的主应变影响Fig.10 Effect of intermediate stress ratio on principal strain for combustible ice sediments

图10为中主应力系数不同时大、中、小主应变ε3与偏应变εs的关系曲线。由图10可以看出,在中主应力系数取值0～1时,大主应力方向始终处于压缩状态,且大主应变ε1随着中主应力系数b增大而减小。当b取值0～0.5时,中主应变ε2处于拉伸状态,且随着中主应力系数增大其拉伸量会逐渐减小;当b=0.5时,中主应变ε2等于0,即这时中主应变方向既没有拉伸也无压缩,处于平面应变状态;当b处于0.5～1.0时ε2处于压缩状态,且随着b不断增大,其压缩量也会增大。当b从0变化到1.0,ε2实现了从拉伸状态到压缩状态的转变。小主应变ε3处于拉伸状态,且随着b逐渐增大,其拉伸量逐渐增大。图11给出了真三轴压缩过程中可燃冰沉积物的体应变εv与偏应变εs的关系。由图11可以看出,在变形初期可燃冰沉积物试样表现为剪缩且随着b增大,其剪缩量逐渐减小;当应变达到某一值时,其剪缩量会达到最大;继续加载,沉积物会转化为剪胀状态,并且随着b增大,其剪胀量也会逐渐增大。

图11 不同中主应力系数对可燃冰沉积物体应变影响Fig.11 Effect of intermediate stress ratio on volume strain for combustible ice sediments

3.2.3 内摩擦角变化规律

由Mohr-Coulomb准则, Mohr-Coulomb的主应力为

(9)

式中,Φ为内摩擦角,(°);c为内聚力,MPa。

根据砂土的强度破坏结果,Matsuoka等[24-25]提出了基于空间滑动面理论的SMP准则,表示为

(10)

其中

I1=σ1+σ2+σ3,I2=σ1σ2+σ2σ3+σ3σ1,I3=σ1σ2σ3.

式中,KSMP为SMP准则下的土体材料的无量纲参数,与土体性质有关,可由试验获取。I1,I2,I3分别为第一、第二和第三主应力不变量。

Lade等[26]根据砂土的真三轴压缩试验提出了Lade-Duncan准则,表达式为

(11)

式中,KLD为Lade-Duncan准则下土体材料的无量纲参数。

假定3种强度准则在b=0时得出的内摩擦角相同,可求得摩尔-库伦准则中内摩擦角为27.92°,KSMP为11.25,KLD为 39.09。

可燃冰沉积物试样受到剪切时的内摩擦角为

(12)

图12为可燃冰沉积物试样内摩擦角与中主应力系数的关系曲线。由图12可得,在b处于0～0.5时,可燃冰沉积物试样的内摩擦角随着中主应力系数增大而增大,在b=0.5处Φ达到最大值34°,之后内摩擦角开始减小,直到b=1.0时,内摩擦角减小为32°,仍然要比b=0时要大。通过DEM模拟结果与3种强度准则的对比可知,Mohr-Coulomb准则由于未考虑中主应力的影响,其计算结果为一条水平直线,并不受中主应力系数的影响。SMP准则与Lade-Duncan准则均考虑了σ2对土体强度的影响,当b较小时,DEM数值结果与SMP准则、L-D准则的预测结果比较接近。当b大于0.1时,SMP准则的预测结果与DEM数值结果差距较大;而L-D准则预测的内摩擦角随中主应力系数的变化曲线在b取0～0.8内与本文模拟结果较为接近;在b取0.8～1.0时略小于本文模拟结果。综合可得,L-D准则更能描述本文的模拟结果。

图12 内摩擦角与中主应力系数关系Fig.12 Relationship of internal friction angle and intermediate stress ratio

3.3 细观机制

可燃冰沉积物内部接触的分布特征可由接触组构张量描述。Satake[27]提出可以通过接触法向分布定义组构张量,表达式为

(13)

式中,N为试样接触总数;ni、nj为接触点处的单位法向向量在i,j方向的分量(i,j=x,y,z)。

接触组构张量为对称的二阶张量,并且接触组构张量在某方向的分量值越大,则证明该方向的接触分布越集中。对于真三轴压缩试验的三维问题,可燃冰沉积物的组构张量φij的矩阵形式表示为

(14)

从而能够较容易地求出沉积物中的颗粒接触点方向的主值φ1,φ2,φ3。可燃冰沉积物内部的颗粒平均接触法向力定义为

(15)

对于可燃冰沉积物来说,根据颗粒接触处各法向接触力和平均法向接触力的关系,可将颗粒间的接触分为强接触与弱接触。当接触点处法向接触力大于试样的平均接触力时,将此类接触定义为强接触,反之,则是弱接触。

Thornton等[28]和张伏光等[29]通过研究净砂与胶结砂土等岩土材料后得出:由颗粒构成的材料宏观力学特性主要通过强接触来体现,然而可燃冰沉积物在这方面的研究比较匮乏,因此,针对可燃冰沉积物DEM试样的宏观力学特性与颗粒间强接触之间的关系作出初步讨论。

图13为DEM试样内部颗粒的强接触点组构的主值与偏应变的关系曲线。由图13可以看出,大主值φ1随着偏应变增大而逐渐增大,达到峰值后出现轻微的减小现象,且大主值φ1受中主应力系数b的影响明显,其值随着b增大而减小,且减小的幅度会越来越大。当0

在可燃冰沉积物试样颗粒体系内部,力链连接形成的力链网络承担着颗粒体系中复杂的动力学响应,从而决定了整体试样宏观力学性质。土体的破坏主要是由强力链的屈服所引起,而正交于强力链的侧向弱力链对强力链的稳定具有辅助作用。在三向不等应力状态下,σ1方向形成的力链承担着大部分荷载,而σ2与σ3方向的力链提供侧向支持。

图13 强接触组构主值与偏应变关系Fig.13 Fabric principal value of strong contacts-deviatoric strain relationship

图14为颗粒间法向接触力投影到XY(σ2-σ3)平面上的力链分布。从图14中可以看出,b=0时的XY平面上的法向接触力分布均匀,随着b增加,X方向的法向接触力分布开始变得不均匀,接触力分布集中在试样中部,而两侧集中力分布较稀疏。对法向接触力的大小和方向进行统计分析,图15给出了p=1 MPa时的法向接触力在XY(σ2-σ3)平面的投影玫瑰分布。从图15中可以看出,b影响了法向接触力的分布,b较小时,法向接触力分布较均匀;随着b增加,中主应力方向的接触力增大,小主应力方向的接触力减小,这表示力链的侧向支持在不断减小,从而使材料的强度降低,稳定性减弱。