引言

出砂是弱胶结砂岩油气藏开发中广泛存在的问题，在完井、生产和修井等环节都可能会出现^[1]。特别是对于海上疏松气藏，%由于其一般都是通常采用衰竭式开采，随着生产的进行，地层能量亏空，地质状况越来越差，同时，高速产出气流产生很高的拖拽力，更容易拖曳砂粒进入井筒，使出砂更为严重。出砂会导致一系列的危害，如：冲蚀套管、磨损生产设备、增加冲砂的井下作业次数、影响产量等，严重时甚至导致油管油泵被卡、储层被埋以及油气井的停产^[2-3]。因此，要防止海上疏松砂岩衰竭开采出砂带来的危害，关键在于提前做好出砂预测，适时调整生产压降，进行有效的防和治。

常用的出砂预测主要有现场观察、经验公式计算、室内实验测试和数值模拟等，具体分析如下。

(1) 现场观测法是通过观察岩芯疏松程度或DST测试来判断是否出砂，其分析过程简单，适用于现场对出砂预测精度要求不高的情况^[4]。

(2) 经验公式法则是基于测井数据和岩石力学测试数据，计算声波时差(95 < ∆t_c < 105，轻微出砂；∆t_c ≥105，严重出砂)，组合模量(15 <E_c≤20，轻微出砂；E_c≤ 15，严重出砂)，出砂指数(出砂指数≤ 20 GPa，出砂)，斯伦贝谢比(斯伦贝谢比≤ 59 MPa²，出砂)来预测地层是否出砂，这种方法简单实用，在国内外被广泛应用^[5-6]。

(3) 室内实验测试是以未胶结或弱胶结砂岩为对象，模拟井眼及其生产条件，进行流动实验来判断是否出砂^[7]。目前常用的实验模型是厚壁空芯圆柱模型，可评价砂岩初始破坏，而无法说明孔眼的扩大及破坏后的稳定性，未考虑射孔产生的破碎带、损伤带等情况。

(4) 数值模拟方法是建立岩石力学与流体力学的动态耦合数值模型，研究出砂油层流固耦合渗流情况下，骨架砂发生剥离的临界生产压差或临界流速来预测地层出砂^[8]，模型所需的储层物性参数、地质力学参数及岩石力学参数过多，参数获取困难，且分析过程相当繁琐复杂，所耗时间长。

目前，出砂临界压差预测模型主要有Morita模型、经验模型、剪切破坏模型、Vaziri模型^[9-11]。其中，Vaziri模型、剪切破坏模型、Morita模型预测考虑参数较少，适用性有限，而经验模型是建立在大量油气田开发数据统计上。

以上出砂预测方法针对的是油气井初期生产的情况，未考虑气井生产中后期气藏压力衰竭对出砂的影响，而压力衰竭均会对出砂构成影响，且不可忽略。现场开采经验表明，很多油气藏开采初期地层强度较高，即使生产产量高，初期也不出砂，但进入开发中后期，地层能量得不到有效补充，产层压力衰竭，井壁周围岩石承受的应力增加，地层开始出砂并逐渐严重^[12-13]。因此，对于衰竭开采海上疏松油气藏，开展气井生产中后期储层压力衰竭的出砂预测及临界生产压差研究，对减少出砂危害、有效防砂和提高产量等都尤为重要。

在储层压力衰竭出砂方面，一些学者通过建立目标油田岩石力学模型，结合岩芯强度测试或厚壁圆筒实验，研究储层压力衰竭与井壁坍塌关系，预测储层压力衰竭后临界生产压差。Ispas等以井周切向应力≥岩芯强度为依据判别是否出砂，并计算衰竭后临界生产压降^[14-15]，Elyasi等以临界生产压降取n倍岩芯强度(n=0 ~ 1)计算井底临界生产压差^[16]。

基于井壁岩石应力分析，引入更适合描述岩石骨架破坏出砂的Mogi-Coulomb准则，建立出砂判别式及临界生产压差力学计算模型，同时考虑地层压力衰竭后地应力变化影响，给出一种适用于海上油气藏压力衰竭开采生产周期出砂预测方法，计算生产井不同生产时期出砂临界生产压差，建立油气藏整个生产周期出砂临界生产压差图版，为生产井配产提供理论依据。

1 岩石骨架破坏出砂临界压差计算

地层砂分为游离砂和骨架砂，游离砂充填在岩石孔喉内，易受地层流体拖拽，随流体进入井筒，但由于孔喉中的游离砂量相对较少，这类出砂不会导致严重的危害。当近井围岩处于高应力状态时，岩石发生屈服塑性变形，达到井壁失稳临界条件，骨架破坏，产生骨架砂粒脱落流出并被带入井筒，会致使地层持续大量产砂^[17]。

骨架破坏出砂是由井壁处逐渐向地层内进行。首先，在井壁应力集中点岩石发生屈服破坏，出砂形成孔穴；然后，孔穴与井筒压力连通，孔穴壁处应力集中且支撑压力最低，此处岩石开始发生屈服破坏，如此周而复始破坏出砂，逐渐向地层内进行形成空穴，产生大量骨架砂，其出砂过程如图 1所示。

1.1 井壁应力状态分析

取井壁处微元段，进行应力分析如图 2所示。假设地层为均质各向同性的弹塑性体，井眼柱坐标系下井周围岩应力如式(1)所示。

$ \left\{ {\begin{array}{*{20}{l}} {{\sigma _r} = \frac{{r_{\rm{w}}^2}}{{{r^2}}}{p_{\rm{w}}} + \frac{1}{2}\left( {{\sigma _x} + {\sigma _y}} \right)\left( {1 - \frac{{r_{\rm{w}}^2}}{{{r^2}}}} \right) + \frac{1}{2}\left( {{\sigma _x} - {\sigma _y}} \right)\left( {1 + \frac{{3r_{\rm{w}}^4}}{{{r^4}}} - \frac{{4r_{\rm{w}}^2}}{{{r^2}}}} \right)\cos 2\theta \\ + {\sigma _{xy}}\left( {1 + \frac{{3r_{\rm{w}}^4}}{{{r^4}}} - \frac{{4r_{\rm{w}}^2}}{{{r^2}}}} \right)\sin 2\theta }\\ {{\sigma _\theta } = - \frac{{r_{\rm{w}}^2}}{{{r^2}}}{p_{\rm{w}}} + \frac{1}{2}\left( {{\sigma _x} + {\sigma _y}} \right)\left( {1 + \frac{{r_{\rm{w}}^2}}{{{r^2}}}} \right) - \frac{1}{2}\left( {{\sigma _x} - {\sigma _y}} \right)\left( {1 + \frac{{3r_{\rm{w}}^4}}{{{r^4}}}} \right)\cos 2\theta \\ - {\sigma _{xy}}\left( {1 + \frac{{3r_{\rm{w}}^4}}{{{r^4}}}} \right)\sin 2\theta }\\ {{\sigma _z} = {\sigma _z} - \mu \left[ {2({\sigma _x} - {\sigma _y})\left( {\frac{{r_{\rm{w}}^2}}{{{r^2}}}} \right)\cos 2\theta + 4{\sigma _{xy}}\left( {\frac{{r_{\rm{w}}^2}}{{{r^2}}}} \right)\sin 2\theta } \right]}\\ {{\tau _{r\theta }} = \left[ {\frac{1}{2}\left( {{\sigma _x} - {\sigma _y}} \right)\sin 2\theta + \frac{1}{2}{\sigma _{xy}}\cos 2\theta } \right]\left( {1 - \frac{{3r_{\rm{w}}^4}}{{{r^4}}} + \frac{{2r_{\rm{w}}^2}}{{{r^2}}}} \right)}\\ {{\tau _{rz}} = \left[ {{\sigma _{xz}}\cos \theta + {\sigma _{yz}}\sin \theta } \right]\left( {1 - \frac{{r_{\rm{w}}^2}}{{{r^2}}}} \right)}\\ {{\tau _{\theta z}} = \left[ {{\sigma _{yz}}\cos \theta - {\sigma _{xz}}\sin \theta } \right]\left( {1 + \frac{{r_{\rm{w}}^2}}{{{r^2}}}} \right)} \end{array}} \right. $

(1)

$ \left( {\begin{array}{*{20}{c}} {{\sigma _x}}&{{\sigma _{xy}}}&{{\sigma _{xz}}}\\ {{\sigma _{xy}}}&{{\sigma _y}}&{{\sigma _{yz}}}\\ {{\sigma _{yx}}}&{{\sigma _y}}&{{\sigma _{yz}}} \end{array}} \right) = H\left( {\begin{array}{*{20}{c}} {{S_{\rm{H}}}}&0&0\\ 0&{{S_{\rm{h}}}}&0\\ 0&0&{{S_{\rm{v}}}} \end{array}} \right){H^{\rm{T}}} $

(2)

$ H = \left( {\begin{array}{*{20}{c}} {\cos \psi \cos \mathit{\Omega }}&{\cos \psi \sin \mathit{\Omega }{\rm{ }}}&{ - \sin \psi }\\ { - \sin \mathit{\Omega }{\rm{ }}}&{\cos \mathit{\Omega }{\rm{ }}}&0\\ {\sin \psi \cos \mathit{\Omega }{\rm{ }}}&{\sin \psi \sin \mathit{\Omega }}&{\cos \psi } \end{array}} \right) $

(3)

以$ \mathit{\boldsymbol{\sigma }} = \left( {\begin{array}{*{20}{c}} {{\sigma _r}}&{{\tau _{r\theta }}}&{{\tau _{rz}}}\\ {{\tau _{r\theta }}}&{{\sigma _\theta }}&{{\tau _{z\theta }}}\\ {{\tau _{rz}}}&{{\tau _{\theta z}}}&{{\tau _z}} \end{array}} \right)$描述井壁围岩一点的应力状态，设I为单位矩阵。当r=r_w时，有

$ \left| {\mathit{\boldsymbol{\sigma }} - \lambda \mathit{\boldsymbol{I}}} \right| = \left| {\begin{array}{*{20}{c}} {{\sigma _r} - \lambda }&{{\tau _{r\theta }}}&0\\ {{\tau _{r\theta }}}&{{\sigma _\theta } - \lambda }&0\\ 0&0&{{\sigma _z} - \lambda } \end{array}} \right| = 0 $

(4)

由式(1)~式(4)，井周岩石微元有效主应力为

$ \left\{ {\begin{array}{*{20}{l}} {{\sigma _r} = {p_{\rm{w}}} - \alpha {p_{\rm{p}}}}\\ {{\sigma _{{\rm{1m}}}} = 0.5({\sigma _\theta } + {\sigma _z}) + }\\ {\;\;\;\;\;\;\;\;\;\;\;\;0.5\sqrt {{{({\sigma _\theta } - {\sigma _z})}^2} + 4\tau _{\theta z}^2} - \alpha {p_{\rm{p}}}}\\ {{\sigma _{{\rm{2m}}}} = 0.5({\sigma _\theta } + {\sigma _z}) - }\\ {\;\;\;\;\;\;\;\;\;\;\;\;0.5\sqrt {{{({\sigma _\theta } - {\sigma _z})}^2} + 4\tau _{\theta z}^2} - \alpha {p_{\rm{p}}}} \end{array}} \right. $

(5)

通过比较$ {\sigma _r}, {\sigma _{{\rm{1m}}}}, {\sigma _{{\rm{2m}}}}$来确定最大主应力、中间主应力和最小主应力

$ \begin{array}{l} {\sigma _1} = \max \left\{ {{\sigma _r}, {\sigma _{{\rm{1m}}}}, {\sigma _{{\rm{2m}}}}} \right\}\\ {\sigma _2} = {\rm{median}}\left\{ {{\sigma _r}, {\sigma _{{\rm{1m}}}}, {\sigma _{{\rm{2m}}}}} \right\}\\ {\sigma _3} = \min \left\{ {{\sigma _r}, {\sigma _{{\rm{1m}}}}, {\sigma _{{\rm{2m}}}}} \right\} \end{array} $

1.2 井壁岩石骨架破坏出砂临界压差

判断岩石是否屈服破坏，目前常用的准则有Mohr-Coulomb准则、Drucker-Prager准则、Modified Lade准则以及Mogi-Coulomb准则(表 1)^[18]。

Mohr-Coulomb准则是分析岩土的经典破坏准则，其假设地层最大原地剪应力由地层的抗剪切强度决定，忽略了中间主应力影响，预测岩石强度偏于保守。Drucker-Prager准则是在Mohr-Coulomb准则和塑性力学中著名Mises准则基础上的扩展和推广，考虑了中间主应力的影响，克服了Mohr-Coulomb准则的主要弱点，但Drucker-Prager准则预测岩石强度偏高。Modified Lade准则也考虑了中间主应力影响，其预测岩石强度破坏结果精确性介于Mohr-Coulomb准则及Drucker-Prager准则之间，但其预测结果仍然偏高；而Mogi-Coulomb准则基于Mohr-Coulomb准则改进而来，充分考虑三向主应力对岩石强度影响，预测结果更符合现场实际效果^[19-23]。

定义F为井壁骨架岩石失稳指数(简称井壁失稳指数)，采用Mogi-Coulomb准则的表达式为

$ F = a + b{\sigma _{\rm{m}}}_{, 2} - {\tau _{{\rm{otc}}}} $

(6)

$ {\tau _{{\rm{otc}}}} = \frac{1}{3}\sqrt {{{\left( {{\sigma _1} - {\sigma _2}} \right)}^2} + {{\left( {{\sigma _2} - {\sigma _3}} \right)}^2} + {{\left( {{\sigma _3} - {\sigma _1}} \right)}^2}} $

(7)

$ {\sigma _{{\rm{m}}, {\rm{2}}}} = \frac{{{\sigma _1} + {\sigma _3}}}{2} $

(8)

$ a = \frac{{2\sqrt 2 }}{3}c\cos \phi $

(9)

$ b = \frac{{2\sqrt 2 }}{3}\sin \phi $

(10)

当岩石微元体所受的应力σ₁，σ₂及σ₃满足式(6)时，F<0时，井壁骨架发生破坏，地层出砂，此时井底流压临界井底流压，可得到井底临界生产压差(Critical Drawdown Pressure，CDP)

$ {p_{{\rm{CD}}}} = {p_{\rm{p}}} - {p_{\rm{w}}} $

(11)

从式(6)~式(11)可见，出砂临界生产压差主要受储层压力、储层岩石所受应力及岩石强度的影响。因此，衰竭式开采气藏，储层孔隙压力衰竭和岩石所受应力改变均会影响地层出砂。

2 孔隙压力衰竭对地应力影响分析

对于构造平缓，面积较大的油气藏，由于压力衰竭引起的水平面内的变形几乎可以忽略，不影响现场工程实际应用效果^[24-28]，有

$ \Delta {\varepsilon _{\rm{H}}} = \Delta {\varepsilon _{\rm{h}}} = 0 $

(12)

根据广义胡克定律，储层开发前的应力-应变关系式为

$ \left\{ {\begin{array}{*{20}{l}} {{\varepsilon _{\rm{H}}} = \frac{1}{E}\left[ {{S_{\rm{H}}}{\rm{ - }}\alpha {p_{\rm{p}}}{\rm{ - }}\mu ({S_{\rm{v}}}{\rm{ - }}\alpha {p_{\rm{p}}} + {S_{\rm{h}}}{\rm{ - }}\alpha {p_{\rm{p}}})} \right]}\\ {{\varepsilon _{\rm{h}}} = \frac{1}{E}\left[ {{S_{\rm{h}}}{\rm{ - }}\alpha {p_{\rm{p}}}{\rm{ - }}\mu ({S_{\rm{v}}}{\rm{ - }}\alpha {p_{\rm{p}}} + {S_{\rm{H}}}{\rm{ - }}\alpha {p_{\rm{p}}})} \right]} \end{array}} \right. $

(13)

孔隙压力改变后，垂直主应力保持不变，而两水平地应力的应力-应变关系为

$ \left\{ {\begin{array}{*{20}{l}} {{\varepsilon _{{\rm{H1}}}} = \frac{1}{E}\left[ {{S_{{\rm{H1}}}} - \alpha {p_{{\rm{p1}}}} - \mu ({S_{\rm{v}}} - \alpha {p_{{\rm{p1}}}} + {S_{{\rm{h1}}}} - \alpha {p_{{\rm{p1}}}})} \right]}\\ {{\varepsilon _{{\rm{h1}}}} = \frac{1}{E}\left[ {{S_{{\rm{h1}}}} - \alpha {p_{{\rm{p1}}}} - \mu ({S_{\rm{v}}} - \alpha {p_{{\rm{p1}}}} + {S_{{\rm{H1}}}} - \alpha {p_{{\rm{p1}}}})} \right]} \end{array}} \right. $

(14)

根据式(12)，有%关系式

$ \left\{ {\begin{array}{*{20}{l}} {{\varepsilon _{{\rm{H1}}}} = {\varepsilon _{\rm{H}}}}\\ {{\varepsilon _{{\rm{h1}}}} = {\varepsilon _{\rm{h}}}} \end{array}} \right. $

(15)

联立式(12)~式(15)，可确定孔隙压力下降后，最大、最小水平主应力变化值，即储层压力衰竭后地层最大、最小水平应力为

$ \left\{ {\begin{array}{*{20}{l}} {{S_{{\rm{H1}}}} = {S_{\rm{H}}} - \Delta {S_{\rm{H}}} = {S_{\rm{H}}} - \frac{{1 - 2\mu }}{{1 - \mu }}\alpha \Delta {p_{\rm{p}}}}\\ {{S_{{\rm{h1}}}} = {S_{\rm{h}}} - \Delta {S_{\rm{h}}} = {S_{\rm{h}}} - \frac{{1 - 2\mu }}{{1 - \mu }}\alpha \Delta {p_{\rm{p}}}} \end{array}} \right. $

(16)

3 实例应用分析 3.1 岩石强度测试

弱胶疏松结气藏F区块内生产井出现不同程度出砂难题，给生产带来很大危害，为保证气井高效生产开发，亟需解决出砂难题。对目标储层井1~440~1~442~m处6块干岩芯进行岩芯强度三轴测试，以获取储层段岩石强度参数，测试结果见表 2。

由表 2计算可得，储层段岩石平均内聚力5.94~MPa，内摩擦角13.97°。

3.2 储层物性及相关参数

结合建立的出砂模型，针对气藏F内的F-1井进行出砂预测，储层及F-1气井参数，见表 3。

3.3 压力衰竭后全生产周期出砂临界生产压差

(1) 孔隙压力衰减后地应力变化趋势

根据压力衰竭地应力变化计算模型，分析最大、最小水平地应力以及有效水平地应力随着地层孔隙压力衰竭变化规律，见图 3。

由图 3可知，储层压力衰竭，最大、最小水平地应力逐渐降低，而有效水平地应力相应增加。

(2) 孔隙压力衰减后井壁失稳指数变化趋势

取井底生产压差2 MPa时，孔隙压力分别衰减至初始孔隙压力的100%，90%，80%及70%时，模拟计算井壁失稳指数，见图 4。

井壁失稳指数越小，井壁越容易发生破坏出砂。图 4中，颜色越红，代表岩石最容易破坏脱落。孔隙压力由原始的100%衰减至70%，井壁失稳指数在逐渐变小，当孔隙压力衰减至70%时，井壁失稳指数小于0，表明该处的岩石已经发生坍塌破坏出砂，图中的红色区域代表井壁骨架在该处最先发生破坏。上述模拟分析了井壁周围骨架的稳定性与孔隙压力大小关系，表明若生产井维持定压差生产，随着孔隙压力衰减骨架容易发生坍塌破坏而导致产砂。

(3) 压力衰竭后临界生产压差

根据CDP计算模型和表 3中F--1井相关参数，计算并绘制出F-1井临界生产压差图版，见图 5。

由图 5可见，(1)暗红色范围内生产压差下气井易出砂，孔隙压力衰竭情况下，淡蓝色范围内生产压差下气井不出砂；(2)储层初期最大临界生产压差为4.90~MPa，随着储层压力衰竭，出砂临界生产压差应逐渐降低以避免产砂；(3)随着压力衰竭，储层临界生产压差应逐渐减小，压力衰竭至临界储层压力11.33~MPa时，井筒内部若存在生产压差，井壁就会发生坍塌破坏，气井需关井不能再生产；(4)从现场实际的配产曲线可以看出，刚开始配产点在淡蓝色区域，现场生产井并没有出砂；随着储层压力衰竭，配产点在暗红色三角形区域内，气井产砂，预测结果与现场实际生产情况一致，验证了模型的有效性}。

4 结论

(1) 基于井壁处岩石应力分析，引进更符合描述井壁处骨架岩石破坏的Mogi-Coulomb准则，建立出砂临界生产压差力学计算模型。

(2) 考虑储层压力衰竭影响，给出了压力衰竭后的地应力计算方式及出砂压力衰竭下临界出砂计算方法，建立生产井全生产周期的临界生产压差模型}。

(3) 应用该模型及相应计算方法，针对南中国海海上衰竭式开采气藏F-1井进行全生产周期出砂预测，计算并绘制出了压力衰竭条件下临界生产压差图版，预测图版与现场实测数据吻合较好，研究结果能为现场生产井配产提供指导意见。

  符号说明

∆t_c——声波时差，µs/ft；

E_c——组合模量，GPa；

S_H——最大水平地应力，MPa；

S_v——垂向地应力，MPa；

p_p——孔隙压力，MPa；

p_w——井筒内压力，MPa；

S_h——最小水平地应力，MPa；

ψ——井斜角，（°）；

Ω——井眼轴线与最大水平地应力夹角，（°）；

σ_r，σ_θ，——r，θ——向的正应力，MPa；

σ_x，σ_y，σ_z——x，y，z向的正应力，MPa；

σ_xy, σ_yz，σ_{xz xy}面，_yz面，_xz面上的应力，MPa；

τ_rθ，τ_rz，τ_θz rθ 面，rz面，θz面上的剪应力，MPa；

r——近井地带一点到井眼中心距离，m；

r_w——井眼半径，m；

θ ——井周角，（°）；

H——变换矩阵；

σ——应力矩阵；

α——Boit系数；

σ_1m，σ_2m井周岩石微元除r向外的另外两个有效主应力，MPa；

σ₁，σ₂，σ₃最大、中间、最小主应力，MPa；

max，median，min——最大、中间、最小值函数；

F——井壁骨架岩石失稳指数，MPa；

a，b——系数；

σ_{m, 2}最大和最小主主应力平均值，MPa；

τ_otc——中间变量，MPa；

c ——内聚力，MPa；

ϕ——内摩擦角，（°）；

∆_εH——最大水平主应变之差，无因次；

∆_εh——最小水平主应变之差，无因次；

ε_H，ε_h——初始最大和最小水平主应变，无因次；

E——杨氏模量，MPa；

µ——泊松比，无因次；

ε_H1, ε_h1 ——孔隙压力衰减后最大和最小水平主应变，无因次；

S_H1，S_h1——孔隙压力衰减后最大和最小水平地应力，MPa；

p_p1——衰减后的孔隙压力，MPa；

ΔS_H, ΔS_h——最大和最小水平主应力变化值，MPa；

Δp_p——孔隙压力改变值，MPa。