CN103206203B - 油井单一射孔出砂的分析方法 - Google Patents

Abstract

本发明涉及一种油井单一射孔出砂的分析方法，其步骤包括：测定某油井任一个射孔的几何参数；测定上述射孔出砂在数值计算中所需要的参数；在柱坐标系下建立三维数值模型，对射孔的出砂过程进行计算；根据计算得到的颗粒应力和位移等数据进行绘图处理，得到砂岩在不同区域的切向应力σ_θ分布、砂岩的塑性区以及颗粒胶结的破坏范围，并对图形进行分析来确定射孔出砂过程中的砂岩力学响应。本发明在理论和技术两个层面上有效解决了复杂赋存环境下油井单一射孔出砂的数值模拟问题，简化了数值计算渗流边界的确定，准确分析油井出砂过程中的砂岩应力、砂岩胶结状态的变化和颗粒的运动，提高了数值计算的精度和效率，增强了该方法的可操作性和稳定性。

Description

油井单一射孔出砂的分析方法

技术领域

本发明涉及一种油井单一射孔出砂的分析方法，是反映不同结构的油井砂岩在射孔试验过程中的砂岩宏细观力学响应，模拟出砂过程。

背景技术

据国家统计局统计，原油进口依存度自2000年的24.8％飙升到2012年的57.16％，已赫然超过了50％的安全警戒线。而石油开采可持续发展的瓶颈之一便是储层的出砂问题，出砂不仅给采油工艺带来许多麻烦，还影响储层采油速度及油气采收率，严重时造成井壁坍塌、套管损坏，乃至油井报废。因此，储层出砂问题给石油开采带来了巨大的挑战，出砂预测和防治是石油工程和岩土工程领域的研究热点和难点，对我国石油开采的可持续发展意义空前。

但油藏开采过程中，由于孔隙流体运移与储层骨架之间的强流固耦合效应，砂岩颗粒大小、颗粒排列、几何形状、延伸方向、孔隙率以及赋存环境相互作用，原始胶结的地层岩石变成松散砂粒，在流体携带作用下最终形成出砂。目前，对油井出砂预测主要的方法主要为理论分析、室内试验和连续介质数值模拟等方法对油井出砂量和油井稳定性进行研究。但实际油藏开采时的地层受力、剥蚀运移、边界条件等的复杂性，油井出砂的物理过程也十分复杂，基于宏观岩体力学的解析法、经验方法等因其自身局限性忽略了砂岩颗粒的细观特征，只能大致确定砂岩体的宏观破坏，而对砂岩颗粒的运动和赋存状态却不能反映，因此不能解释油井射孔过程中的出砂机理问题，也影响了防砂措施的有效性和针对性。即使有些学者基于离散元方法考虑渗流作用下的颗粒运动，但主要是基于笛卡尔坐标系考虑二维达西渗流效应引起的颗粒运动，不能真实模拟射孔试验中的渗流边界和流固耦合效应。因此，有必要提出一种能反映油井射孔出砂过程中的砂岩颗粒胶结性状和颗粒运动的方法，既能反映出砂过程中砂岩的应力、塑性区和砂岩颗粒胶结的破坏，又能预测出砂速度和出砂量，能为防砂预控提供科学的依据，确保油井正常生产和石油工业的可持续发展。

发明内容

本发明的目的在于提供能简单、直接和高效的分析油井单一射孔出砂的分析方法，该方法能分析不同颗粒级配、孔隙率和胶结程度的油井砂岩在复杂赋存环境中的宏细观力学响应，实现储层砂岩胶结性状的模拟和流固耦合效应的计算，提高三维数值方法的准确性、收敛性和稳定性，为实际储层射孔出砂的预测提供重要的技术手段。

为实现上述发明目的，本发明采用的技术方案为：油井单一射孔出砂的分析方法，包括以下步骤：

1)测定油井任一个射孔的几何参数，包括射孔的内孔直径d(其半径为R_i)、外直径D(其半径为R_o)和高度h，由于砂岩出砂过程中，塑性区是靠近射孔的，远离射孔位置是弹性区域，弹性区域的应力变化是较稳定的，也就是说砂岩模型的塑性区并不是无限大的，为了提高计算效率，对模型的外边界给出一个限值即射孔的外直径D；

2)测定步骤1)中射孔出砂过程中的砂岩颗粒的几何参数，孔隙率n′，砂岩颗粒的密度ρ_s，砂岩颗粒的力学参数，油井的实际油藏的流体力学参数，射孔的岩层围压σ_confine，划分的流体单元f_element；将所述内孔与外半径之间的区域划分成若干单元(即图2中R_i和R_o之间的区域)，主要可以根据计算机的性能来确定单元的划分，以计算各单元在外荷作用下出现的孔隙率变化。

其中所述砂岩颗粒的几何参数包括最小颗粒半径r_min以及最大颗粒半径r_max与最小颗粒半径r_min的比值r_ratio；所述砂岩颗粒的力学参数包括砂岩颗粒的密度ρ_s，砂岩颗粒的平行粘结法向刚度切向刚度法向强度和切向强度砂岩颗粒的摩擦系数f，砂岩颗粒的法向刚度K_N和切向刚度K_S，砂岩颗粒胶结的程度α(用平行粘结的胶结半径比表示，式中为砂岩颗粒胶结的平均半径，R_A为胶结处颗粒的最小半径)，胶结砂岩颗粒的百分含量m′；所述实际油藏的流体力学参数包括实际油藏流体的密度ρ_f，流体的动力粘滞系数μ_f，根据射孔过程测定的内孔流量q；

3)在柱坐标系下建立三维数值模型，利用步骤2)中的各个参数对步骤1)所述射孔的出砂过程进行计算，得到砂岩颗粒的应力和位移数据；

4)根据步骤3)计算得到的砂岩颗粒的应力和位移数据进行绘图处理，得到砂岩颗粒在不同位置的切向应力σ_θ分布、砂岩颗粒的塑性区以及砂岩颗粒胶结物质的破坏范围；

5)结果分析：若计算结果未显示塑性区，则不会出砂，若出现塑性区，则出砂，且随着流量的增长或流速的增大，塑性区越大，出砂越多。

所述砂岩颗粒胶结的程度α，用平行粘结的胶结半径比表示，式中为砂岩颗粒胶结的平均半径，R_A为胶结处颗粒的最小半径。

具体地，所述步骤3)包括如下步骤：

Step1：利用基于PFC3D程序的Fish语言建立模拟砂岩颗粒的胶结状态的计算模型，输入步骤2)所述的砂岩颗粒的几何参数，对每个砂岩颗粒接触力进行[0,1]区间的均匀设计，每个颗粒接触对应一个随机数，对所有的颗粒接触进行循环，当砂岩中某颗粒接触位置设定的随机数小于实测的砂岩颗粒胶结含量(所有接触的砂岩颗粒上都有胶结物质的话，那么胶结含量为100％，若接触颗粒中只有一半有胶结物质，那么认为胶结含量为50％，若接触的颗粒上没有胶结物质，则胶结含量为0)时，此处的颗粒接触状态确定为胶结，将步骤2)所述的平行粘结法向刚度切向刚度法向强度和切向强度赋给该胶结状态的接触；

Step2：为了提高计算机的运行速度和计算效率，引入一个更新参数N_up来确定渗流计算的时间，初始化力学循环参数N_cycle为1；

Step3：施加岩层围压σ_confine计算初始状态下的颗粒应力和位移，进行力学循环计算，当N_cycle＝N_up时，输入步骤2)所述的油井的实际油藏的流体力学参数，计算每个流体单元f_element中的颗粒位置；力学循环包括颗粒接触以及平行粘结性态判断，流固耦合分析和计算颗粒接触力、位移和转动，并分析砂岩的力学性态等步骤。

Step4：计算每个流体单元f_element的流体单元孔隙率n；

Step5：计算每个流体单元f_element中沿着径向方向颗粒的平均直径式中m为f_element内的颗粒数，d_p,i为第_i颗粒的直径，颗粒的平均速度流体单元的平均速度其中r_i为每个流体单元f_element的中心点到射孔中心的距离；

Step6：根据每个流体单元f_element的流体单元孔隙率n，计算流体对砂岩颗粒产生的压力梯度和拖曳力系数

Step7：根据计算的流体压力梯度力和拖曳力系数计算作用于砂岩颗粒上的作用力 $f_{d_{i_j}}\left(i\right)(i=1,\mathrm{...},n_p,j=r,\mathrm{\θ},z)$ 为

$f_{d_{i_j}}=-\left(\frac{{\mathrm{\β}}_{{\mathrm{int}}_j}}{1-n}\right({\stackrel{\‾}{v}}_j-u_j)+\▿p_j)\frac{\mathrm{\π}}{6}{d}_{pi}^3$

式中，n_p为砂岩颗粒的总数；d_pi为颗粒直径，为流体的运动粘滞系数，u_j为流体流动的表观速度；

Step8：根据和计算作用在颗粒上的作用力赋值N_cycle＝1，对颗粒接触和平行粘结性态进行判断，并将作用力作为外力作用在颗粒上进行力学计算，重新计算颗粒的应力、位移以及每个流体单元孔隙率n，直到N_cycle＝N_up。

具体地，所述步骤4)包括如下步骤：

S1：切向应力σ_θ的分布是根据砂岩颗粒的坐标(x,y,z)确定砂岩颗粒所在的方位角θ_d，再根据每个颗粒的应力张量σ_ij计算颗粒的切向应力σ_s,θ，对每个流体单元f_element计算总体积V以及该单元内颗粒的总切向应力V_m为所计算流体单元f_element内编号为m的颗粒体积，_n1为单元内颗粒总数，计算得到射孔的外边界应力σ_ro，那么每个单元内单位体积上的平均切向应力即可获得；根据切向应力σ_θ的数值，绘图得到应力分布曲线，应力曲线的最大值即为弹塑性边界的分界点，该分界点到油井中心点的距离即为塑性区的半径R_c，根据切向应力σ_θ最大值的位置可确定塑性区的范围；根据数值计算结果就可以得到单颗颗粒的应力σ_ij，根据弹性力学理论就可以得到单颗粒切向方向的主应力值σ_s,θ＝-(σ_xxsin²θ_d+σ_yycos²θ_d-2σ_xysinθ_dcosθ_d)，对应整个模型，为了更好地反映宏观应力变化曲线，需要计算模型中某位置点的应力，这个点的应力可以通过分析某流体单元内的平均应力来获得，对每个流体单元f_element计算总体积V以及该单元内颗粒的总切向应力每个单元内单位体积上的平均切向应力就可以计算出来。计算出不同点的切向应力，就可以得到切向应力的分布曲线，这个曲线有个峰值点，那么这个峰值点对应的位置即为砂岩模型弹塑性区域的分界点。

S2：根据PFC3D确定颗粒胶结的状态，标识无粘结的颗粒为A₁和有粘结的颗粒为A₂、剪切破坏的粘结颗粒数A_1shear和受拉破坏的粘结颗粒数A_1tensile以及对应颗粒中心点的坐标(x,y,z)，对不同的颗粒状态进行分类作图，分别得到A₁和A₂构成的区域，那么A₁构成的区域即代表颗粒粘结破坏的区域，形成砂岩的破坏区域Zone，即获得砂岩颗粒胶结的破坏范围，该破坏区域Zone的半径R_zone代表颗粒胶结的破坏半径，与砂岩的塑性区半径R_c相等，在流体的携带作用下，失去胶结的颗粒即可成为出砂。

本发明采用平行粘结模型、流固耦合理论以及三维颗粒流程序相结合的分析方法，在理论和技术两个层面上有效解决了复杂赋存环境下油井单一射孔出砂的数值模拟问题，简化了数值计算渗流边界的确定，准确分析油井出砂过程中的砂岩应力、砂岩胶结状态的变化和颗粒的运动，提高了数值计算的精度和效率，增强了该方法的可操作性和稳定性。因此，该发明是一种简单、直观和高效的出砂模拟方法，能再现油井出砂的发生和发展过程，并进行出砂速度和出砂量的预测。

附图说明

图1为油井单一射孔出砂的分析流程；

图2为油井出砂物理模型；

图3为柱坐标系下流固耦合模型流体单元划分俯视图；

图4为出砂过程的流固耦合计算流程；

图5为数值模型的可行性对比图；

图6为不同流量对应的切向应力图形；

图7为流量q＝3.14*10^-3m³/s的平行粘结破坏分布图。

具体实施方式

下面结合附图和实施例对本发明作进一步说明，但不应该理解为本发明上述主题范围仅限于下述实施例。在不脱离本发明上述技术思想的情况下，根据本领域普通技术知识和惯用手段，作出各种替换和变更，均应包括在本发明范围内。

本发明提供的油井出砂三维数值模拟方法，本发明中测量的参数都能通过本领域中常规的试验方法和现场测量得到，包括以下步骤：

1)测定某油井任一个射孔的几何参数；

2)测定单一射孔出砂在数值计算中所需要的参数；

3)在柱坐标系下建立三维数值模型，对油井某单一射孔的出砂过程进行计算；

4)根据3)计算得到的颗粒应力和位移等数据进行绘图处理，得到砂岩在不同区域的切向应力分布、砂岩的塑性区以及颗粒胶结的破坏范围，并对图形进行分析来确定射孔出砂过程中的砂岩力学响应。

上述测定某油井任一个射孔的几何参数，确定在数值模拟中某射孔的内孔直径d、外直径D和高度h；

上述测定单一射孔出砂在数值计算中所需要的参数，确定砂岩的最小颗粒半径r_min以及最大颗粒半径r_max与最小颗粒半径r_min的比值r_ratio，r_min取值一般为0.01mm～1.0mm，r_ratio取值1.0～6.0；砂岩颗粒的密度ρ_s，取值为2500kg/m³～2750kg/m³；砂岩颗粒的胶结拟用平行粘结模拟，根据现场试验确定平行粘结的法向刚度和切向刚度取值为1.0*10¹⁰kPa/m～1.0*10¹³kPa/m，平行粘结法向强度和切向强度取值为1.0*10⁵kPa～1.0*10⁸kPa；砂岩颗粒的摩擦系数f，取值为0～0.5；砂岩颗粒的法向刚度K_N和切向刚度K_S，取值为1.0*10⁵kN/m～1.0*10⁷kN/m；砂岩颗粒胶结的程度，拟用平行粘结的胶结半径比表示，式中为砂岩颗粒胶结的平均半径，R_A为胶结处颗粒的较小半径，取值范围为0～1.0；胶结砂岩颗粒的百分含量m′，取值为0～100％；孔隙率n′，取值为0.2～0.8；流体的密度ρ_f，根据实际油藏流体的密度确定，测定为500kg/m³～2000kg/m³；流体的动力粘滞系数μ_f，根据实际油藏的粘性确定，取值为0.5*10^-3Pa.s～50*10^-3Pa.s；内孔流量q，根据射孔过程测定实际流量，一般可取值为1.0*10^-6m/s³～1.0*10^-2m/s³；岩层围压σ_confine，可取值为5MPa～100MPa；流体的单元划分f_element，径向、切向和垂直方向的单元数量需根据计算机内存确定，取值一般可为2～100。

上述在柱坐标系下建立三维数值模型，对油井某单一射孔的出砂过程进行计算，具体步骤如下：

Step1：确定数值模型的力学计算边界条件，输入步骤2)所述的砂岩颗粒几何参数，基于PFC3D程序的Fish语言模拟砂岩颗粒的胶结状态：对每个颗粒接触进行[0,1]区间的均匀设计，每个颗粒接触对应一个随机数，对所有的颗粒接触进行循环，当砂岩中某颗粒接触位置设定的随机数小于实测的砂岩颗粒胶结含量时，此处的颗粒接触状态确定为胶结，将如2)所描述的平行粘结参数赋给该接触。

Step2：为了提高计算机的运行速度和计算效率，引入一个更新参数N_up来确定渗流计算的时间，初始化力学循环N_cycle为1。

Step3：力学计算过程中，如果N_cycle＝N_up，开始渗流计算，输入如2)所描述的流体单元参数、流体力学参数和流量，计算每个流体单元f_element中的颗粒位置。

Step4：计算每个流体单元f_element的流体单元孔隙率n，计算某一具体流体单元(如N_i)的步骤如下：

A)判断砂岩颗粒i的中心点是否处于流体单元N_i内部，是则进入步骤B)，否则舍弃该砂岩颗粒判断下一个砂岩颗粒；

B)判断砂岩颗粒i是否与流体单元N_i的环向边界相交：若否，则砂岩颗粒i的体积即为V₀(根据砂岩颗粒的颗粒半径即可得到)；若是，则进一步判断颗粒i中心点的位置；

C)判断砂岩颗粒i是否与流体单元N_i的水平边界或径向边界相交(或砂岩颗粒i是否与水平边界和径向边界同时相交)。若否，则砂岩颗粒i在流体单元N_i内的体积V₁即为V₀；若是，则进一步计算V₁。当颗粒与径向边界相交时，则V₁＝V₀-(0.25α³-0.75α+0.5)V₀，式中：α为径向边界线的夹角弧度与切向线夹角弧度的比值(切向线表示的是从射孔中心点到颗粒边界作切面，投影后即为切向线)。当颗粒与水平边界相交时，则式中：p为颗粒中心点到水平边界线的距离，r为颗粒半径；当颗粒与径向边界和水平边界同时相交时， $V_1=V_0-\left(0.25{\mathrm{\α}}^3-0.75\mathrm{\α}+0.5\right)V_0-\left(0.25{\left(\frac{p}{r}\right)}^3-0.75\frac{p}{r}+0.5\right)\left(1-\left(0.25{\mathrm{\α}}^3-0.75\mathrm{\α}+0.5\right)\right)V_0;$

D)重复A～C，直到获得处于流体单元N_i内各砂岩颗粒的体积，加和得到流体单元N_i内所有颗粒的总体积V_tot，计算流体单元N_i的孔隙率n＝1-V_tot/V_j，式中：V_j为流体单元N_i的体积。

Step5：计算每个流体单元f_element中沿着径向方向颗粒的平均直径颗粒的平均速度流体单元的平均速度u_f,r。

Step6：根据每个流体单元f_element的孔隙率n，计算流体对砂岩颗粒产生的压力梯度和拖曳力系数其计算步骤如下：

(1)计算流体压力梯度若流体雷诺数R_e较小时，其运动符合达西定律，其值为 $\frac{dp}{{\mathrm{dx}}_j}=-\frac{{\mathrm{\ρ}}_{f}g}{K}{u}_{j,0}\≡-\frac{v_f{\mathrm{\ρ}}_f}{k}{u}_{j,0}.$

式中，ν_f,ρ_f为流体的运动粘滞系数(m²/s)和密度(kg/m³)，g为加速度(m/s²)，u_j,0(j＝r,θ,z)为流体流动的表观速度(m/s)，K为渗透系数(m/s)，k为渗透率(m²)，根据Kozeny-Carmen方程得到k为c为常数，为颗粒的平均直径(m)，n为流体单元孔隙率。

当雷诺数R_e较大时，流体显示出非线性特征而不再符合达西定律，则流体压力梯度为

$\▿p_j(j=r,\mathrm{\θ},z)$ 为 $\frac{dp}{{\mathrm{dx}}_j}=-(150\frac{{\left(1-n\right)}^2}{n^3{\stackrel{\‾}{d}}_p^2}{v}_f{\mathrm{\ρ}}_f{u}_{j,0}+1.75\frac{\left(1-n\right){\mathrm{\ρ}}_f}{n^3{\stackrel{\‾}{d}}_p}{\mathrm{\ρ}}_f{u}_{j,0}^2)$

考虑实际油藏的流体单元孔隙率n的变化，流体压力梯度为

${\mathrm{\Δp}}_j=(150\frac{{\left(1-n\right)}^2}{n^2{\stackrel{\‾}{d}}_p^2}{\mathrm{\μ}}_f+1.75\frac{\left(1-n\right){\mathrm{\ρ}}_f}{n{\stackrel{\‾}{d}}_p}|{\stackrel{\‾}{v}}_j-u_j\left|\right)\left({\stackrel{\‾}{v}}_j-u_j\right)$

式中，分别为颗粒运动的平均速度(m/s)和流体实际流动速度(m/s)。

(2)考虑流体与颗粒运动速度的不同产生拖曳力，需依据三维Navier-Stokes方程 $\frac{\∂n}{\∂t}=-(\▿\·nu)$ 和 $\frac{\∂\left(nu\right)}{\∂t}=-(\▿\·nuu)-\frac{n}{{\mathrm{\ρ}}_f}\▿p-\frac{n}{{\mathrm{\ρ}}_f}\▿\·\mathrm{\τ}+ng+\frac{f_{\mathrm{int}}}{{\mathrm{\ρ}}_f}$ 计算，t为时间(s)，u为流体运动速度(m/s)，p为压强(Pa)，ρ_f为流体密度(kg/m³)，τ为流体粘滞应力(Pa)，g为加速度(m/s²)，f_int为由于颗粒与流体速度差异而产生的单位拖曳力(N/m³)。

那么，不同方向上的拖曳力 $\begin{array}{cc}{f}_{{\mathrm{int}}_j}(j=r,\mathrm{\θ},z)={\mathrm{\β}}_{{\mathrm{int}}_j}({\stackrel{\‾}{v}}_j-u_j),& {\mathrm{\β}}_{{\mathrm{int}}_j}(j=r,\mathrm{\θ},z)\end{array}$ 为拖曳力作用系数(Pa.s/m²)。

考虑油藏实际孔隙率n的变化，即可得到为

${\mathrm{\β}}_{{\mathrm{int}}_j}=150\frac{{(1-n)}^2}{n{\stackrel{\‾}{d}}_p}{\mathrm{\μ}}_f+1.75\frac{(1-n){\mathrm{\ρ}}_f}{{\stackrel{\‾}{d}}_p}\left|{\stackrel{\‾}{v}}_j-u_j\right|$

Step7：根据计算的流体压力梯度力和拖曳力系数计算作用于砂岩颗粒上的作用力为

$f_{d_{i_j}}=-\left(\frac{{\mathrm{\β}}_{{\mathrm{int}}_j}}{1-n}\right({\stackrel{\‾}{v}}_j-u_j)+\▿p_j)\frac{\mathrm{\π}}{6}{d}_{pi}^3$

式中，n_p为砂岩颗粒的总数；d_pi为颗粒直径。

Step8：根据和计算作用在颗粒上的进行力学计算，并设N_cycle＝1，对颗粒接触和平行粘结性态进行判断，并将作为外力作用在颗粒上进行力学计算，重新计算颗粒的应力。位移以及每个流体单元的孔隙率n，直到N_cycle＝N_up。

上述的根据3)计算得到的颗粒应力和位移等数据进行绘图处理，得到砂岩在不同区域的切向应力σ_θ分布、砂岩的塑性区以及颗粒胶结的破坏范围，并对图形进行分析来确定射孔出砂过程中的砂岩力学响应。具体步骤如下：

Step1：切向应力σ_θ的分布是根据颗粒的坐标(x,y,z)确定颗粒所在的方位角θ_d，再根据每个颗粒的应力张量σ_ij计算颗粒的切向应力σ_s,θ，对模型的每个单元f_element计算总体积V以及该单元内颗粒的总切向应力V_m为计算单元f_element内编号为m的颗粒体积，n1为单元内颗粒总数，计算得到模型的外边界应力σ_ro，那么每个单元内单位体积上的平均切向应力σ_θ即可获得。根据切向应力σ_θ的数值，绘图得到应力分布曲线，应力曲线的最大值即为弹塑性边界的分界点，该分界点到油井中心点的距离即为塑性区的半径R_c，根据切向应力σ_θ最大值的位置可确定塑性区的范围。

Step2：颗粒胶结范围的确定是根据PFC3D确定颗粒胶结的状态，标识无粘结的颗粒A₁和有粘结的颗粒A₂、剪切破坏的粘结颗粒数A_1shear和受拉破坏的粘结颗粒数A_1tensile以及对应颗粒中心点的坐标(x,y,z)，通过Matlab软件对不同的颗粒状态进行分类作图，得到A₁和A₂的区域，那么A₁即代表破坏的颗粒，形成砂岩的破坏区域Zone，该Zone的半径R_zone代表颗粒胶结的破坏半径，与砂岩的塑性区半径R_c相等，在流体的携带作用下，失去胶结的颗粒即可成为出砂。

下面结合附图对本发明做进一步的说明，但不限定本发明。

参见图1，一种油井出砂的三维数值模拟方法，本发明包含了油井砂岩计算参数的确定、油井出砂数值模型的建立、流固耦合分析以及砂岩宏细观力学响应的模拟等四个模块。其具体实施步骤如下：

Step1：根据实际砂岩颗粒的物性特征，确定数值计算需要的颗粒的几何参数和颗粒的力学参数，如颗粒摩擦系数、颗粒接触刚度、颗粒的粘结刚度和颗粒的粘结强度。

Step2：根据实际油井的几何参数，建立具有一定内孔径、外半径和高度的物理模型，并根据射孔试验的特征，在柱坐标系下生成砂岩颗粒形成三维颗粒流数值模型，并定义力学边界和渗流边界，如围压、垂直压力和流速等参数。

Step3：计算作用于颗粒的流体压力梯度力和由于颗粒与流体速度差产生的拖曳力，综合为作用于颗粒各方向上的合力f。

Step4：将f作为外荷载进行力学分析，得到颗粒的受力和位移，其力学状态的改变引起颗粒的重分布，最终导致每个流体单元孔隙率n的变化，而n的变化又引起拖曳力和流体压力梯度的变化，从而引起颗粒细观力学效应和模型宏观力学特性的变化，即流固耦合分析。

Step5：根据数值计算结果通过图形软件进行后处理，得到砂岩的宏细观力学响应，如模型的应力、颗粒粘结的应力、颗粒粘结的破坏、颗粒的位移以及颗粒的转动等，并对其特征量进行分析。

参见图2，上述射孔试验包含了数值计算中射孔的内外孔径和高度，参见表1。

参见图3，上述描述的柱坐标系下三维颗粒流数值模型中的流体单元划分，考虑流体的连续性，将三维模型划分为径向、切向和垂直方向的单元，单元数分别为12、12和2，并设置径向、切向和垂向边界。通过设置径向、环向和水平边界条件，以实现单元划分。即在柱坐标系下，通过设定r、θ、z各自的区间值，将所述射孔计算模型划分为若干流体单元。本实施例具体划分中，仅将所述内孔与外半径之间的区域划分成若干单元(即图2中R_i和R_o之间的区域)，以便于进一步研究。图3即为所述射孔数值模型的一个横截面。流体单元是一个具有三个边界范围的三维图形。图3中所能见到的是环向边界和径向边界。值得说明的是，划分的每个流体单元都有径向、环向和水平向边界，划分的单元越多，则边界相对就越多，消耗计算机的内存越大。由于划分过多的流体单元会导致计算缓慢，因此，应该根据计算机的处理能力划分适宜数量的流体单元。

参见图4，上述描述的流固耦合作用是计算作用于颗粒的流体压力梯度力和由于颗粒与流体速度差产生的拖曳力，直到力学平衡为止。

参见图5，上述描述的数值模型的可行性对比图，实线为数值计算得到的切向应力σ_θ曲线，虚线即为某射孔测得的切向应力σ_θ曲线力，两曲线较吻合，该方法可用来分析油井单一射孔出砂过程中的砂岩力学特性。

表1 三维数值模型计算参数

参见图6，上述描述的不同流量对应的切向应力图形，计算参数如表1，考虑流固耦合效应，计算结果显示了对应不同的流量q＝3.14*10^-3m³/s、6.28*10^-3m³/s和7.85*10^-3m³/s，砂岩标准化切向应力σ_θ的最大值为1.8、1.9和2.0，其最大值所在的位置到射孔中心点距离为26mm、31mm和35mm，即为射孔试验中砂岩颗粒被携带出去的区域。

参见图7，上述描述的砂岩胶结的破坏范围，是根据PFC3D计算得到的颗粒胶结状态进行分类绘图而成，图中灰色区域即为颗粒间胶结破坏的范围，计算流量q＝7.85*10^-3m³/s时对应的颗粒胶结破坏的平均半径为35mm，与图6显示的塑性区范围相同，从另一角度说明了出砂过程中的砂岩力学特性。

Claims (1)

1.油井单一射孔出砂的分析方法，其特征是包括以下步骤：

1)测定油井任一个射孔的几何参数，包括射孔的内孔直径d、外直径D和高度h；

2)测定步骤1)中射孔出砂过程中的砂岩颗粒的几何参数，孔隙率n′，砂岩颗粒的密度ρ_s，砂岩颗粒的力学参数，实际油藏的流体力学参数，射孔的岩层围压σ_confine，划分的流体单元f_element；

其中所述砂岩颗粒的几何参数包括最小颗粒半径r_min以及最大颗粒半径r_max与最小颗粒半径r_min的比值r_ratio；所述砂岩颗粒的力学参数包括砂岩颗粒的密度ρ_s，砂岩颗粒的平行粘结法向刚度切向刚度法向强度和切向强度砂岩颗粒的摩擦系数f，砂岩颗粒的法向刚度K_N和切向刚度K_S，砂岩颗粒胶结的程度α，胶结砂岩颗粒的百分含量m′；所述实际油藏的流体力学参数包括实际油藏流体的密度ρ_f，流体的动力粘滞系数μ_f，根据射孔过程测定的内孔流量q；

3)在柱坐标系下建立三维数值模型，利用步骤2)中的各个参数对步骤1)所述射孔的出砂过程进行计算，得到砂岩颗粒的应力和位移数据；

4)根据步骤3)计算得到的砂岩颗粒的应力和位移数据进行绘图处理，得到砂岩颗粒在不同位置的切向应力σ_θ分布、砂岩颗粒的塑性区以及砂岩颗粒胶结物质的破坏范围；

5)结果分析：若计算结果未显示塑性区，则不会出砂，若出现塑性区，则出砂，且随着流量的增长或流速的增大，塑性区越大，出砂越多；

所述砂岩颗粒胶结的程度α，用平行粘结的胶结半径比表示，式中为砂岩颗粒胶结的平均半径，R_A为胶结处颗粒的最小半径；

所述步骤3)包括如下步骤：

Step1：利用基于PFC3D程序的Fish语言建立模拟砂岩颗粒的胶结状态的计算模型，输入步骤2)所述的砂岩颗粒的几何参数，对每个砂岩颗粒的接触进行[0,1]区间的均匀设计，每个颗粒接触对应一个随机数，对所有的颗粒接触进行循环，当砂岩中某颗粒接触位置设定的随机数小于实测的砂岩颗粒胶结含量时，此处的颗粒接触状态确定为胶结，将步骤2)所述的平行粘结法向刚度切向刚度法向强度和切向强度赋给该胶结状态的接触；

Step2：为了提高计算机的运行速度和计算效率，引入一个更新参数N_up来确定渗流计算的时间，初始化力学循环的参数N_cycle为1；

Step3：施加岩层围压σ_c=nfine计算初始状态下的颗粒应力和位移，进行力学循环计算，当N_cycle＝N_up时，输入步骤2)所述的油井的实际油藏的流体力学参数，计算每个流体单元f_element中的颗粒位置；

Step4：计算每个流体单元f_element的流体单元孔隙率n；

Step5：计算每个流体单元f_element中沿着径向方向颗粒的平均直径式中m为f_element内的颗粒数，d_p,i为第i颗粒的直径，颗粒的平均速度流体单元的平均速度其中r_i为每个流体单元f_element的中心点到射孔中心的距离；

Step6：根据每个流体单元f_element的流体单元孔隙率n，计算流体对砂岩颗粒产生的压力梯度和拖曳力系数

Step7：根据计算的流体压力梯度力和拖曳力系数计算作用于砂岩颗粒上的作用力为

式中，n_p为砂岩颗粒的总数；d_pi为颗粒直径，为流体的运动粘滞系数，u_j 为流体流动的表观速度；

Step8：根据和计算作用在颗粒上的作用力赋值对颗粒接触和平行粘结性态进行判断，并将作用力作为外力作用在颗粒上进行力学计算，重新计算颗粒的应力、位移以及每个流体单元孔隙率n，直到 N_cycle=N_up。