3-D analysis method for pile-supported geosynthetic reinforcements based on analog equation method
-
摘要: 桩承式加筋体的载荷–变形是复杂的三维力学问题,目前三维理论分析方法尚不完善。鉴于此,以三维情况下方形布桩的加筋体为研究对象,将其视为受上部荷载和桩间土反力共同作用下的空间薄膜,并基于最小势能原理推导了其三维平衡偏微分方程组。鉴于该方程组具有形式复杂及高度非线性的特点,针对荷载均匀分布及无地基反力的情况,应用薄膜分析中的比拟方程法进行求解,据此计算加筋体的变形与拉力,初步建立了一种桩承式加筋体的三维分析方法。最后,通过一个现场试验算例计算和结果对比分析初步验证了三维分析方法的合理性。Abstract: The load-deformation behavior of pile-supported geosynthetic reinforcements (GRs) is a complicated three-dimensional mechanical problem, and so far the 3-D analytical methods have not been well-established. Therefore the GRs with piles in square-type layout under 3-D conditions are taken as the research objects, which are regarded as the spacial membranes under the combined effect of the upper load and the subsoil reaction among piles. The partial differential equilibrium equations are derived based on the minimum potential energy principle. For the conditions of uniform load and without subsoil reaction, the analog equation method in the membrane analysis is used to solve the equations in view of their complex forms and highly nonlinear characteristics. Then the deformation and tensile force of GRs are accordingly calculated, and the 3-D analytical method for the pile-supported GRs is established. Finally, the rationality of the method is validated primitively by comparing the calculated and observed results of a field test example.
-
0. 引言
水力压裂技术被广泛应用于解决非常规油气资源开发中的增透提渗、矿山治理中的注浆加固及地应力测量等问题[1-4]。水力压裂破裂压力是指井内液体压力不断增大使地层产生水力裂缝或张开原有裂缝时所达到的峰值压力,与围岩物理力学特性、地应力状态、流体压力及天然裂缝发育情况密切相关,是控制水力压裂设计和施工的重要参数[5-7]。现有水力压裂破裂压力计算模型多基于应用广泛的Hubbert-Willi(H-W)模型和Haimson-Fairhurst(H-F)模型[8-9],其中H-W模型未考虑围岩的渗透效应,可得到破裂压力的上限值;H-F模型考虑了围岩的高渗透性,得到破裂压力的下限值[10]。王六鹏等[11]基于最优化理论扩展了H-W模型,利用井周围岩应力分布和张性破裂准则,分别建立了裸眼和射孔完井条件下破裂压力非线性约束优化计算模型;郭天魁等[12]在H-F模型基础上考虑了渗流引起的应力,提出了新的水平井射孔孔壁的应力分布,建立了考虑页岩射孔水平井水力压裂在围岩内破裂、沿天然裂缝剪切破裂和沿天然裂缝张性破裂3种破裂压力计算模型。李传亮等[13-14]基于双重有效应力概念,集合了H-W模型和H-F模型,建立了裸眼井完井条件下围岩破裂压力计算模型,由此提出了射孔完井条件下产生垂向裂缝和水平裂缝情况下围岩破裂压力计算模型。杨兆中等[15]在井筒围岩应力场计算中,考虑套管和水泥环的影响,建立了煤层直井井周应力场解析模型。上述模型虽考虑了诸多影响因素,但水力压裂过程中注入流体渗流和围岩应力的耦合作用尚未涉及。
水力压裂是井内流体压力、围岩内流体渗流和围岩应力相互耦合的过程[16-17];水力压裂过程中,随着流体注入,井内流体压力逐渐升高;压力驱动作用下流体持续注入井壁围岩内,导致围岩流体压力改变,围岩有效应力场亦改变,影响围岩内部结构的演化,引起流体应力场再次改变;当围岩应力状态达到破裂准则,射孔围岩破裂,此时井内流体压力即为水力压裂破裂压力。故水力压裂过程中,射孔围岩应力场和流体场处于相互联系、相互影响和相互制约的水力耦合状态。渗透率和孔隙度的应力敏感性是研究水力耦合问题的关键。黄远智等[18]基于渗流–应力耦合试验,认为渗透率演化在于有效围压作用下岩石孔隙变形,且渗透率会随有效围压的增加呈指数下降。卢家亭等[19]研究了低渗砂岩微观结构与渗透率应力敏感性的关系,认为乘幂函数可表述渗透率和有效应力的关系。薛永超等[20]探讨了不同渗透率岩石的应力敏感性规律,认为中高渗储层岩石应力敏感性符合“平缓型”模式,而低渗透储层岩石应力敏感性呈“先陡后缓”的规律,且低渗透岩石具有更强的应力敏感性。上述研究表明,岩石渗透率和孔隙度的应力敏感性研究取得了较多的成果,但将其考虑至水力压裂耦合的研究中鲜有成果。
综上所述,现存水力压裂射孔围岩破裂模拟研究多考虑单纯的应力作用,与流体渗流耦合研究成果较少。射孔围岩流体渗流具有非稳态性,导致与围岩应力的水力耦合作用变得异常复杂,且相关围岩应力状态亦产生较大改变。本文基于渗透率和孔隙度的应力敏感性,构建了考虑流体渗流和围岩应力耦合的射孔围岩水力压裂破裂力学模型,并基于有限容积法提出了水力耦合作用下射孔围岩水力压裂破裂数值模拟方法,进而验证了方法的正确性。
1. 射孔围岩应力分布
1.1 物理模型及假设
储层岩石处于垂向地应力σV、最大水平应力σH和最小水平应力σh状态下。假设射孔井是由井筒和射孔组成的正交系统,在远场地应力与井内流体压力作用下,近井区域围岩应力重分布,引起双重应力集中。
为了准确分析射孔围岩应力状态,提出以下假设:①储层岩石为均匀各向同性、线弹性多孔材料,且井壁围岩处于平面应变状态;②井筒与井壁胶结良好时,水力压裂过程中光滑接触;井筒与井壁胶结不好时,不考虑井筒与井壁的相互作用,认为井壁光滑且与井筒间存在微环隙,流体较容易流入微环隙;③在渗流计算时不考虑流体的压缩性;④井筒内和射孔内流体压力相等,储层内均产生渗流。
1.2 井筒围岩应力分布
考虑围岩在远场地应力和井筒内流体压力共同作用下(图 1),结合井壁渗流状态,利用叠加原理得到极坐标系下围岩的应力分布[21-22]:
σr=(σH+σh)2(1−r2wr2)+(σH−σh)2(1−4r2wr2+3r4wr4)cos2θ+δ1⋅Pw2(1−ν′)nr2br2[1+(1−2ν′)n]r2br2a−(1−n)+(1−δ1)Pwrw2r2+δ2[α(1−2ν)(1−ν)1r2∫rrw(P(r)−P0)rdr−ϕ(P(r)−P0)], (1) σθ=(σH+σh)2(1+r2wr2)−(σH−σh)2(1+3r4wr4)cos2θ− δ1⋅Pw2(1−ν′)nr2br2[1+(1−2ν′)n]r2br2a−(1−n)−(1−δ1)Pwr2wr2−δ2[α(1−2ν)(1−ν)1r2∫rrw(P(r)−P0)rdr−(α1−2ν1−ν−ϕ)(P(r)−P0)], (2) σz=σv−2ν(σH−σh)cos2θ+(α1−2ν1−ν−ϕ)(P(r)−P0), (3) τrθ=−(σH−σh)2(1+2r2wr2−3r4wr4)sin2θ, (4) τrθ=τzθ=0。 (5) 式中,σr,σθ,σz为井筒围岩的三向正应力,τrθ,τrz,τzθ为剪切应力,rw为井筒半径,θ是与σH的极坐标角,ν是泊松比,n=E(1+ν′)/(E′(1+ν))为围岩和井筒剪切模量之比,E和E′分别为围岩和井筒的弹性模量,ν和ν′分别为围岩和井筒的泊松比,ra和rb为井筒内径和外径,P0为初始流体压力,P(r)为实际流体压力,ϕ为有效孔隙度,α为有效应力系数。δ1为井筒与井壁胶结状态系数,胶结不好时取0,胶结良好时取1;δ2为井壁的渗透性系数,不渗透时取0,可渗透时取1。当井筒与井壁胶结良好时,井筒围岩不考虑渗流应力。
1.3 射孔围岩应力分布
假设射孔孔眼和井眼正交,力学模型及坐标转换如图 2所示,基于井筒围岩应力可得到射孔围岩应力:
σs=(σz+σθ)2(1−s2ws2)+ (σθ−σz)2(1−4s2ws2+3s4ws4)cos2φ+Pws2ws2+α(1−2ν)(1−ν)1r2∫rrw(P(r)−P0)rdr−ϕ(P(r)−P0), (6) σφ=(σz+σθ)2(1+s2ws2)−(σθ−σz)2(1+3s4ws4)cos2φ− Pws2ws2−α(1−2ν)(1−ν)1r2∫rrw(P(r)−P0)rdr+ (α1−2ν1−ν−ϕ)(P(r)−P0), (7) σzz=σr−2ν(σθ−σz)s2ws2cos2θ+(α1−2ν1−ν−ϕ)(P(r)−P0), (8) τzzφ=2τrθsinφ, (9) τszz=τsφ=0。 (10) 式中,σs,σφ分别为射孔横截面上的径向应力和周向应力,σzz为轴向应力,τzzφ,τszz,τsφ为剪切应力;sw是射孔半径;φ是σs与σθ方向的夹角。
2. 射孔围岩流体压力
2.1 非Darcy渗流方程
射孔围岩不是完全刚性的,流体压力改变引起有效应力变化,导致围岩渗透率K和孔隙度ϕ呈现敏感性。流体在多孔介质中的渗流可用Darcy定律描述;对于非Darcy渗流(雷诺数Re > 10[23]),Forchheimer [24]公式被用于代替Darcy定律:
−∂P∂x=μKxvx+βρvx2 ,−∂P∂y=μKyvy+βρvy2 ,−∂P∂z=μKzvz+βρvz2 ,} (11) 式中,vx,vy,vz分别为流体的三向渗流速度,Kx,Ky,Kz分别为介质的三向渗透率,ρ为流体密度,μ为流体黏度。该方程比Darcy定律多一个渗流速度的二次方项,代表惯性力影响,β被称为惯性系数。本文采用适用于致密和松散的砂岩惯性系数公式[25]:
β = 1.59×103√Kϕ5.5, (12) 式中,渗透率K的单位为mD。
为了与Darcy定律对比分析,引入非Darcy渗透率KN:
KN = K1+βρKvμ。 (13) 式(11)可简化为
vx=−KNxμ∂P∂x ,vy=−KNyμ∂P∂y ,vz=−KNzμ∂P∂z 。} (14) 2.2 渗流偏微分方程
储层围岩中的流体压力可通过渗流偏微分方程得到;将围岩中的流体视为单相微可压缩牛顿流体,在等温无源非稳态渗流下的连续性方程可表示为
∂(ρϕ)∂t+∇⋅(ρ→v)=0。 (15) 考虑Darcy定律和非Darcy定律,认为流体黏度μ与空间位置无关,得到二维非稳态渗流偏微分方程:
ctμϕ∂P∂t=∂∂x(Kx∂P∂x)+∂∂y(Ky∂P∂y), (16) ctμϕ∂P∂t=∂∂x(KNx∂P∂x)+∂∂y(KNy∂P∂y)。 (17) 式中,ct=cf+cϕ为综合压缩系数,cf为流体压缩系数,cϕ为围岩孔隙压缩系数。渗流开始时,认为储层围岩内任意位置的流体压力均为初始孔隙压力,即初始条件为P(r,0)=P0;其边界条件如表 1所示。
表 1 渗流分析的边界条件Table 1. Boundary conditions of fluid flow analysis外边界 自由边界 P(re,t)=0 内边界 射孔壁和井壁压力 P(rw,t)=Pw(t) , P(sw,t)=Pw(t) 2.3 渗透率和孔隙度的应力敏感性
渗透率和孔隙度的应力敏感性直接影响围岩内流体压力。现有研究结果表明围岩渗透率随有效应力增加呈指数降低[26-27]:
K=K0e−Mσ′=K0e−M(σ−αP), (18) 式中,K0为初始渗透率,M为渗透率应力敏感性系数,σ′为Terzaghi有效正应力,P为孔隙压力,α为Biot’s系数。基于多孔介质双重有效应力概念,孔隙度与α紧密相关[28]:
α=(φφc)23φ≤φc ,α=1 φ>φc ,} (19) 式中,ϕc为多孔介质材料的临界孔隙度。Nur[29]基于大量室内试验研究,给出了ϕc的建议值为0.4。基于Kozeny-Carman方程[30],得到渗透率和孔隙度的关系:
K=ϕk2SP2。 (20) 式中SP为固体颗粒表面积与孔隙体积的比值;水力压裂过程中,假设注入流体后围岩变形为小变形,SP认为不变。k2为常数,由此得到渗透率K与孔隙度ϕ成正比:
ϕ=ϕ0e−Mσ′=ϕ0e−M(σ−αP)。 (21) 2.4 破裂力学准则
Hossain等[31]认为最大拉应力准则对于预测破裂压力是准确的。由于井筒和射孔相交的双重应力集中,最大拉应力发生在近交点围岩区域内;当水平方向存在主应力差,射孔方位与最大水平应力方向垂直时,破裂发生在井筒围岩最大水平应力处。针对射孔围岩应力分布,需采用复合应力理论先求出最大拉应力,即第三主应力σ3对于角度φ的极值σ3,min,然后利用最大拉应力准则(22)判定围岩破裂和位置[32]。
σ3,min−αP=−σt, (22) 随着射孔方位角θ的增大,σ3,min逐渐增大,破裂压力亦逐渐增大;当井筒与井壁胶结不好时,最大水平应力方向井壁处围岩应力不会随射孔方位角θ变化;当射孔方位角θ增大时,最大水平应力方向的井壁处围岩会首先达到破裂条件,破裂准则也应该包括如下判据:
(σH+σh)−2(σH−σh)−Pw−δ2(α1−2ν1−ν−ϕ)P−αP=−σt。 (23) 3. 水力压裂破裂数值模拟方法
3.1 有限容积法空间区域及控制方程离散化
流体在围岩内渗流属非稳态过程,采用有限容积法(FVM)开展数值模拟计算[33]。首先,将计算区域离散为互不重叠的控制容积,且使每个网格节点周围有一个控制容积(U,W,E,S,N点)。采用单元中心法对整个计算区域开展离散,即单元位于子区域中心,此时子区域就是控制容积(N1,N2,N3,N4点组成的阴影区域),且网格线与界面线(线段w,e,s,n)重合,故采用单元代表整个控制容积上的力学属性(图 3)。其次,将守恒型的控制方程(微分方程)对每个控制容积积分,离散化后转化成各节点上的离散方程组,离散方程中的未知数就是节点上因变量的数值。
控制方程积分和离散后,需要假定被求函数的型线,求出控制容积积分。本模型选用分段线性分布型线,即假定界面线上被求函数本身及其一阶导数(包括对时间和空间)均在节点间的线性分布上。
3.2 数值计算模型及参数
取垂直井,考虑对称布置的两个射孔,射孔方位角为θ,选取井筒和射孔围岩的四分之一(φ=90°)构建模型(图 4),压裂液体为水,注入流量为15 mL/min,泵容量为200 mL;模型计算参数见表 2。
表 2 FVM模拟计算模型参数Table 2. Model parameters of numerical simulation in FVM井筒外径/mm 井筒内径/mm 射孔直径/mm 射孔长度
/mm模型尺寸/(mm×mm) 网格尺寸/(mm×mm) 6 4 1 4 50×50 0.5×0.5 最大水平主应
力/MPa最小水平主应力/MPa 垂向主应力
/MPa围岩泊松比 井筒弹性模量
/GPa井筒泊松比 8.53 6.57 9.85 0.25 200 0.15 弹性模量/MPa 抗拉强度/MPa 初始孔隙度 应力敏感性模数M 初始渗透率
/10-17m2368.79 0.715 0.124 0.3 27.194 求解过程中,该数值方法亦存在一定的假设;①模型为二维渗流问题,控制容积界面线上的渗流是均匀的,即渗透速度在界面线上处处相同;②模型左边界和下边界为对称边界,对称边界不发生垂直于边界的渗流;射孔井壁存在渗流;③应力边界在模型无穷远处。
3.3 二维非稳态渗流方程离散
渗流偏微分方程是由控制容积的质量守恒推导的连续性方程,对渗流控制方程离散,即把渗流偏微分方程转换为控制容积节点上的离散方程组(图 5)。
时间间隔[t,t+Δt]内,在控制容积U上(图 6),对渗流偏微分方程(16)积分:
∫ns∫ew∫t+Δttctμϕ∂P∂tdxdydt=∫t+Δtt∫ns∫ew∂∂x(Kx∂P∂x)dxdydt+∫t+Δtt∫ns∫ew∂∂y(Ky∂P∂y)dxdydt, (24) 左侧非稳态项可直接离散为
∫ns∫ew∫t+Δttctμϕ∂P∂tdxdydt=(μϕtct)U×(Pt+ΔtU−PtU)ΔxΔy, (25) 右侧扩散项可离散为
∫t+Δtt∫ns∫ew∂∂x(Kx∂P∂x)dxdydt+∫t+Δtt∫ns∫ew∂∂y(Ky∂P∂y)dxdydt=[KxePE−PU(δx)e−KxwPU−PW(δx)w]ΔyΔt+[KynPN−PU(δy)n−KysPU−PS(δy)s]ΔxΔt, (26) 由此得到离散方程:
aUPU=aEPE+aWPW+aNPN+aSPS+atUPtU ,aU=aE+aW+aN+aS+atU ,} (27) 式(27)中离散方程系数计算方法如下:
aE=Δy(δx)e/Kxe,aW=Δy(δx)w/Kxw,aN=Δx(δy)n/KynaS=Δx(δy)s/Kys,atU=(μφtct)UΔxΔyΔt。 (28) 3.4 应力方程的迭代
对于每个控制容积U(图 6),其应力状态σt+ΔtU是关于井内流体压力Pw(t)、孔隙度ϕt和流体压力PtU的函数,可直接参与每一时间步的迭代计算:
σt+ΔtU=f(Pw(t),φt,PtU)。 (29) 为耦合围岩应力分布与流体渗流状态,需修正渗透率和孔隙度与有效应力的关系:
Kt+ΔtU=K0e−M(σtU−αtUPtU), (30) ϕt+ΔtU=ϕ0e−M(σtU−αtUPtU), (31) αt+Δt=(ϕt+Δtϕc)23 ϕt+Δt≤ϕc ,α=1 ϕt+Δt>ϕc 。} (32) 且考虑非Darcy渗透时,应对渗透率予以非Darcy修正:
vtU=KtNUμ(PtU−1−PtU), (33) vtU>vL ,Re>10。 (34) βt+ΔtU=1.59×103√KtU(KtU)5.5 ,Kt+ΔtNU=KtU1+βt+ΔtUρKtUvtUμ ,} (35) 式中,vtU为控制单元流速,vL为临界流速,可由vL=Reϕcμ/(ρd)得到,其中d为岩石孔隙直径。
3.5 边界条件及初始条件
针对上述离散方程,须考虑相应的边界条件予以求解;包括对称边界条件、内边界条件和外边界条件。其中对称边界处不发生垂直于边界的渗流,故边界节点系数满足:asym=0;针对内边界,由于井筒和射孔内节点的整个控制容积内流体压力均为Pw(t),且直接作用在界面线上(图 7)。
与内边界相邻节点的(δx)w或(δy)s为原来的一半,其系数变为
aIn,W=2Δy(δx)In,w/KIn,xw,aIn,S=2Δx(δy)In,s/KIn,ys。 (36) 外边界节点实际为宏观内部节点,故可补充边界节点代数方程对其处理:
PM=PM+1+PM−12,aM+1=aM。 (37) 针对初始条件,当t=0,ϕt=ϕ0,那么对于确定的控制容积尺寸,atU=a0U, PtU=P0U=P0;对于井内和射孔内流体压力始终满足:P(rw,t)=Pw(t),P(sw,t)= Pw(t),且Pw(0)=0。
3.6 数值模拟迭代流程
根据上述离散方程、边界条件和初始条件,结合破裂准则,可迭代计算出井内流体压力、围岩流体压力和应力,由此得到破裂时间和破裂位置。每一时间步的具体迭代流程如下:
(1) 输入围岩力学参数及流体参数,计算临界流速;判断流体是否符合Darcy定律,进而修正非Darcy渗透率KN。
(2) 考虑参数应力敏感性,修正渗透率Kt和孔隙度ϕt(若不考虑应力敏感性,可跳过参数修正);基于井筒内流体压力Pw,计算围岩流体压力Pt,得到渗流场分布。
(3) 基于流体压力Pt,计算围岩应力分布,判断第三有效主应力和最小有效周向应力是否满足破裂条件。
(4) 若满足破裂条件,则输出破裂压力和破裂时间等结果,结束计算;若不满足则重复迭代计算,直至满足破裂条件。
4. 数值模拟结果
4.1 模型验证
考虑达西渗流的情况下,基于广泛使用的Hubbert-Willis(H-W)和Haimson-Fairhurst(H-F)解析模型对本文数值模拟方法予以验证。H-W模型不考虑岩石的渗透性,其结果可作为验证的上限,计算公式如下[8]:
pb=3σh−σH+σt−P0。 (38) 式中,pb为破裂压力。H-W模型计算得到的破裂压力为11.90 MPa。
H-F模型考虑岩石的高渗透性,结果可作为验证的下限,计算公式如下[9]:
pb=3σh−σH+σt−2ηP02(1−η), (39) η=ϕ0(1−2ν)2(1−ν)。 (40) 式中,ϕ0为储层初始孔隙度,η为模型参数。相同条件下,H-F模型计算得到的破裂压力为6.20 MPa。
在射孔方位角为0°,考虑井筒与井壁胶结良好、井筒与井壁胶结不好且不考虑井壁渗透、井筒与井壁胶结不好且考虑井壁渗透条件下,本文数值模拟得到的破裂压力分别为11.56,6.25,6.74 MPa。对比上述结果,当井筒与井壁胶结良好,流体经射孔直接作用于围岩上,破裂压力计算结果接近于H-W模型;当井筒与井壁胶结不好时,井内流体易进入井壁裂隙中,流体直接作用在井壁围岩上,导致破裂压力的降低,计算结果接近于H-F模型。且井筒与井壁胶结不好时,考虑井壁渗透条件下计算得到的破裂压力高于不考虑井壁渗透;主要原因在于井壁渗流条件下考虑了渗透作用引起的附加压应力。3种条件下得到的破裂压力计算结果均处于H-W模型和H-F模型破裂压力计算结果之间,验证了模型的正确性。
4.2 井筒与井壁胶结状态的影响
不考虑渗透率和孔隙度应力敏感性,对井筒与井壁不同胶结状态条件下的射孔围岩水力压裂破裂压力、破裂时间和破裂位置开展模拟研究;主要考虑井与壁胶结良好、井与壁胶结不好但井壁不可渗、井与壁胶结不好且井壁可渗3种工况。3种工况下,射孔围岩破裂压力和破裂时间随射孔方位角增大而增加(图 8);破裂压力越大,破裂需要的流体压力越高,破裂时间越长。当射孔方位角小于45°,破裂压力和破裂时间增长幅度较大,当射孔方位角大于45°,破裂压力和破裂时间增长速度趋于稳定,由此认为水力压裂宜采用略低的射孔方位角。
井与壁胶结良好时的破裂压力最大,井筒内流体压力对围岩的影响可以忽略,破裂位置位于射孔顶端。井筒与井壁胶结不好时,不考虑井壁渗透作用,流体通过围岩裂隙直接作用于井壁上,与射孔内的流体在井筒与射孔相交处形成双重应力集中,故在低射孔方位角下,破裂位置会在射孔末端;随着射孔方位角的增大,最大水平应力方向近井处的孔隙流体压力会逐渐减小,破裂位置开始转向最大水平应力方向的井壁处,故高射孔方位角时破裂压力呈现增大趋势。井筒与井壁胶结不好时,考虑井壁渗透作用,最大水平应力方向井壁处的破裂条件与射孔方向角无关,破裂均发生在最大水平主应力方向井壁处;且井壁渗透导致井内流体压力增长缓慢,渗透作用有效延迟了围岩破裂时间。
此外,射孔方位角0°条件下,井筒与井壁胶结良好时,流体通过射孔进入地层,流体压力呈关于射孔段的椭圆形分布(图 9);随着压裂时间的增长,流体压力分布沿椭圆形扩展,渗流影响范围增加。当井筒与井壁胶结不好且不考虑井壁渗透时,虽然最大流体压力不同,但渗流规律与井壁胶结良好时基本一致。当井筒与井壁胶结不好且考虑井壁渗透时,流体渗透面积增加,不仅在射孔处产生渗流,井壁处亦产生渗流,故流体压力分布是沿井眼向外扩散的,渗流影响范围大幅度增加(图 9)。
监测不同压裂时刻典型截面(图 4,截面1—1,坐标:y=0.5 mm)流体压力的时空演化(图 10)。流体进入射孔围岩后,流体压力大幅度降低;当达到破裂时间t0,射孔围岩流体压力显著增加,直接影响射孔围岩的破裂。当井筒与井壁胶结不好且考虑井壁渗透时,流体压力时空演化曲线较井筒与井壁胶结良好时平缓,印证了井筒与井壁胶结不好且考虑井壁渗透时渗流影响范围增加显著。
4.3 应力敏感性的影响
渗流影响射孔围岩水力压裂的破裂压力和破裂时间,考虑渗透率和孔隙度应力敏感性的条件下,井筒与井壁胶结良好时,破裂压力和破裂时间均随射孔方位角的增大而增加(图 11)。但考虑应力敏感性时,围岩破裂压力和破裂时间均小于不考虑应力敏感性的状态;且考虑应力敏感性时,随着射孔方位角增大,破裂压力增长幅度降低显著,射孔方位角越大,考虑应力敏感性引起的围岩破裂压力和破裂时间降低越明显。其原因在于围岩流体压力增大,有效应力系数α增加,井壁围岩有效周向应力减小,产生应力集中,直至达到抗拉强度准则开始破裂,导致破裂压力和破裂时间的降低。井筒与井壁胶结不好且不考虑井壁渗透作用下,应力敏感性对不同射孔方位角下围岩破裂压力和破裂时间基本不存在影响。井筒与井壁胶结不好且考虑井壁渗透作用下,考虑应力敏感性,围岩破裂压力和破裂时间均降低;且在高射孔方位角下,降低明显(图 12)。
射孔方位角0°条件下,考虑渗透率和孔隙度应力敏感性,井筒与井壁胶结良好时,围岩内流体压力亦呈关于射孔段的椭圆形分布(图 13)。但相对于不考虑应力敏感性,渗流影响范围大幅度增加,原因在于围岩流体压力增大,导致有效应力降低,由于渗透率应力敏感性的存在,其渗透率不断增大,流体渗流范围随之增大,且渗透率演化规律和流体压力分布规律基本一致。井筒与井壁胶结不好且考虑井壁渗透时,流体压力分布规律与不考虑应力敏感性时基本一致,但渗流影响范围扩大,渗透率的演化规律亦与流体压力分布规律基本一致,均沿水平方向变化(图 14)。
监测不同压裂时刻典型截面(图 4,截面1—1,坐标:y=0.5 mm)流体压力的时空演化。考虑渗透率和孔隙度应力敏感性,井筒与井壁胶结良好和井筒与井壁胶结不好时(图 15),由于近井围岩渗透率增大,围岩流体压力时空演化曲线均趋于平缓,流体压力分布更加均匀,导致流体压力的降低速率变缓。
5. 结论
考虑流体渗流和围岩渗透率、孔隙度的应力敏感性,构建考虑水力耦合的下射孔围岩水力压裂破裂力学模型的基础上,采用有限容积法提出了相关数值模拟方法;主要结论如下:
(1) 考虑初始地应力和流体渗流对射孔围岩的影响,应用力学叠加原理得到了射孔围岩的应力分布;且提出了射孔围岩水力压裂破裂准则不仅应包括满足围岩破裂的最大拉应力准则,亦应考虑最大水平应力方向井壁处围岩的破裂判据。
(2) 考虑围岩渗透率和孔隙度的应力敏感性,确定了射孔围岩的流体压力;由此构建了水力耦合作用下射孔围岩水力压裂破裂力学模型;并基于有限容积法(FVM)提出了水力耦合作用下射孔围岩水力压裂破裂数值模拟方法。
(3) 基于上述方法开展数值模拟计算,射孔围岩破裂压力和破裂时间随射孔方位角的增大而增加;且井筒与井壁胶结良好时破裂压力大于井筒与井壁胶结不好时;考虑渗透率和孔隙度的应力敏感性,围岩破裂压力和破裂时间均降低;考虑渗透率和孔隙度的应力敏感性导致流体压力分布更加均匀,近井区域流体压力增大,围岩渗透率增加,引起渗流影响范围扩大。
-
表 1 算例计算方案
Table 1 Calculation schemes of example
计算方案编号 填土内摩擦角φ/(°) 边界单元总数N 内部配点总数M 加筋体荷载B/(kN·pile-1) 加筋体均布荷载p0/kPa 水平位移分量u/m 水平位移分量v/m 1.1 42 228 65 59.849 10.668 0.05 0.05 1.2 50 228 65 50.905 9.074 0.05 0.05 1.3 60 228 65 43.430 7.742 0.05 0.05 1.4 68 228 65 40.445 7.209 0.05 0.05 表 2 算例计算结果
Table 2 Calculated results of example
计算方案编号 填土内摩擦角φ/(°) 竖向位移/m 应变值/% 拉力值/(kN·m-1) wA wC εyA εyB εxC εyC NyA NxB NyB NxC NyC 1.1 42 0.142 0.224 1.64 1.78 1.01 1.01 26.55 24.43 36.10 23.27 23.37 1.2 50 0.134 0.212 1.47 1.59 0.90 0.90 23.82 22.07 32.27 20.84 20.84 1.3 60 0.127 0.201 1.32 1.42 0.81 0.81 21.42 20.02 28.90 18.70 18.70 1.4 68 0.124 0.196 1.26 1.34 0.77 0.77 20.42 19.09 27.42 17.82 17.82 表 3 加筋体竖向位移值及应变值对比
Table 3 Comparison between vertical displacement and strain values of GRs
对比项目 测量仪器 实测值 模拟值 计算值 竖向位移/m SP03 0.36 0.215~0.245 0.196~0.224 SP07 0.17 0.139~0.159 0.124~0.142 SP08 0.17 应变/% ε1 2.05 0.97~1.27 1.34~1.78 ε2 1.73 ε3 1.50 ε6 1.50 1.13~1.47 1.26~1.64 ε10 1.36 ε7 1.14 0.75~0.98 0.77~1.01 ε8 0.97 表 4 加筋体拉力值对比
Table 4 Comparison of tensile force values of GRs
位置 方向 实测值/(kN·m-1) 模拟值/(kN·m-1) 计算值/(kN·m-1) 桩间条带中心 y 29.34 18.27~23.81 20.42~26.55 y 25.84 18.27~23.81 20.42~26.55 四桩中心点 x 25.40 17.45~22.77 17.82~23.27 y 23.28 17.45~22.77 17.82~23.27 -
[1] JONES B M, PLAUT R H, FILZ G. M. Analysis of geosynthetic reinforcement in pile-supported embankments. Part I: 3D plate model[J]. Geosynthetics International, 2010, 17(2): 59-67. doi: 10.1680/gein.2010.17.2.59
[2] HALVORDSON K A, PLAUT R H, FILZ G M. Analysis of geosynthetic reinforcement in pile-supported embankments: Part II 3D cable-net model[J]. Geosynthetics International, 2010, 17(2): 68-76. doi: 10.1680/gein.2010.17.2.68
[3] PLAUT R H, FILZ G M. Analysis of geosynthetic reinforcement in pile-supported embankments. Part Ⅲ: Axisymmetric model[J]. Geosynthetics International, 2010, 17(2): 77-85. doi: 10.1680/gein.2010.17.2.77
[4] 饶为国, 江辉煌, 侯庆华. 桩–网复合地基工后沉降的薄板理论解[J]. 水利学报, 2002(4): 23-27. doi: 10.3321/j.issn:0559-9350.2002.04.005 RAO Wei-guo, JIANG Hui-huang, HOU Qing-hua. Deformation of sheet plate due to residual settlement of pile-net composite foundation[J]. Journal of Hydraulic Engineering, 2002(4): 23-27. (in Chinese) doi: 10.3321/j.issn:0559-9350.2002.04.005
[5] 郑俊杰, 张军, 马强, 等. 双向增强体复合地基桩土应力比三维分析[J]. 华中科技大学学报(自然科学版), 2010, 38(2): 83-86. doi: 10.13245/j.hust.2010.02.012 ZHENG Jun-jie, ZHANG Jun, MA Qiang, et al. Three dimensional analysis of pile-earth stress ratio of biaxial reinforcement composite foundation[J]. Journal of Huazhong University of Science and Technology(Natural Science Edition), 2010, 38(2): 83-86. (in Chinese) doi: 10.13245/j.hust.2010.02.012
[6] 张军, 郑俊杰, 马强. 路堤荷载下双向增强体复合地基受力机理分析[J]. 岩土工程学报, 2010, 32(9): 1392-1398. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201009016.htm ZHANG Jun, ZHENG Jun-jie, MA Qiang. Mechanical performance of biaxial reinforcement composite foundation under embankment loads[J]. Chinese Journal of Geotechnical Engineering, 2010, 32(9): 1392-1398. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201009016.htm
[7] 赵明华, 刘猛, 张锐, 等. 路堤荷载下双向增强复合地基荷载分担比及沉降计算[J]. 岩土工程学报, 2014, 36(12): 2161-2169. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201412003.htm ZHAO Ming-hua, LIU Meng, ZHANG Rui, et al. Calculation of load sharing ratio and settlement of bidirectional reinforced composite foundation under embankment loads[J]. Chinese Journal of Geotechnical Engineering, 2014, 36(12): 2161-2169. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201412003.htm
[8] 闫澍旺, 郎瑞卿, 孙立强, 等. 基于薄板理论的刚性桩网复合地基桩土应力比计算[J]. 岩石力学与工程学报, 2017, 36(8): 2051-2060. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201708025.htm YAN Shu-wang, LANG Rui-qing, SUN Li-qiang, et al. Calculation of pile-soil stress ratio in composite foundation with rigid pile-net based on plate theory[J]. Chinese Journal of Rock Mechanics and Engineering, 2017, 36(8): 2051-2060. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201708025.htm
[9] 李波, 黄茂松, 叶观宝. 加筋桩承式路堤的三维土拱效应分析与试验验证[J]. 中国公路学报, 2012, 25(1): 13-20. https://www.cnki.com.cn/Article/CJFDTOTAL-ZGGL201201002.htm LI Bo, HUANG Mao-song, YE Guan-bao. Analysis of three –dimensional soil arching effect of pile-supported embankment with geosynthetics and its test verification[J]. China Journal of Highway and Transport, 2012, 25(1): 13-20. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-ZGGL201201002.htm
[10] 徐超, 林潇, 沈盼盼. 桩承式加筋路堤张力膜效应模型试验研究[J]. 岩土力学, 2016, 37(7): 1825-1831. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201607001.htm XU Chao, LIN Xiao, SHEN Pan-pan. Model tests of tensile membrane effect of geosynthetic-reinforced piled embankments[J]. Rock and Soil Mechanics, 2016, 37(7): 1825-1831. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201607001.htm
[11] 庄妍, 崔晓艳, 王康宇. 考虑加筋体三维变形的桩承式路堤沉降计算方法[J]. 天津大学学报(自然科学与工程技术版), 2019, 52(12): 1227-1234. https://www.cnki.com.cn/Article/CJFDTOTAL-TJDX201912002.htm ZHUANG Yan, CUI Xiao-yan, WANG Kang-yu. A method to analyze the settlement of reinforced piled embankment considering the three-dimensional deformation of geosynthetic reinforced structure[J]. Journal of Tianjin University(Science and Technology), 2019, 52(12): 1227-1234. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-TJDX201912002.htm
[12] KANG Fei. A simplified method for analysis of geosynthetic reinforcement used in pile supported embankments[J]. The Scientific World Journal, 2014: 1-9, doi: 10.1155/2014/273253.
[13] VAN EEKELEN S J M, BEZUIJEN A, VAN TOL A F. An analytical model for arching in piled embankments[J]. Geotextiles and Geomembranes, 2013, 39: 78-102.
[14] TSIATAS G C, KATSIKADELIS J T. Large deflection analysis of elastic space membranes[J]. International Journal for Numerical Methods in Engineering, 2006, 65: 264-294.
[15] KATSIKADELIS J T. The analog equation method. A powerful BEM-based solution technique for solving linear and non-linear engineering problems[J]. Boundary Element Method Ⅷ, 1994: 167-182.
[16] KATSIKADELIS J T, NERANTZAKI M S, TSIATAS G C. The analog equation method for large deflection analysis of membranes: a boundary-only solution[J]. Computational Mechanics, 2001, 27: 513-523.
[17] KATSIKADELIS J T, TSIATAS G C. The analog equation method for large deflection analysis of heterogeneous orthotropic membranes: a boundary-only solution[J]. Engineering Analysis with Boundary Elements, 2001, 25: 655-667.
[18] GOLDBERG M A, CHEN C S, KAPUR S P. Improved multiquadric approximation for partial differential equations[J]. Engineering Analysis with Boundary elements, 1996, 18: 9-17.
[19] ALMEIDA M S S, EHRLICH M, SPOTTI A P, et al. Embankment supported on piles with biaxial geogrids[C]//Proceedings of the Institution of Civil Engineers, Geotechnical Engineering 160, 2007, Issue GE4: 185-192.
[20] ALMEIDA M S S, MARQUES M ES, ALMEIDA M C F, et al. Performance of two “low”piled embankments with geogrids at rio de janeiro[C]//Proceedings of the First Pan American Geosynthetics Conference and Exhibition, 2008, Cancun: 1285-1295.
[21] VAN EEKELEN S J M, ALMEIDA M S S, BEZUIJEN A. European analytical calculations compared with a full-scale Brazilian piled embankment[C]//Proceedings of IGS10, 2014, Berlin.
[22] VAN EEKELEN S J M, BEZUIJEN A, VAN TOL A F. Validation of analytical models for the design of basal reinforced piled embankments[J]. Geotextiles and Geomembranes, 2015, 43: 56-81.