CN103077556B - 油井出砂的三维数值模型设计方法 - Google Patents

Abstract

本发明的目的是提供一种考虑射孔试验特征，反映射孔试验中砂岩的宏细观力学特性的模型设计方法。为实现本发明目的而采用的技术方案是这样的，油井出砂的三维数值模型设计方法，包括以下步骤：1）读入实际油井的几何尺寸，建立柱坐标系下的射孔三维数值模型；2）根据射孔试验过程的实际情况， 将所述射孔三维数值模型划分成若干单元；3）根据所述射孔三维数值模型的岩石颗粒分布，计算出每一个单元的孔隙率；4）根据所述孔隙率计算出流体压力梯度力和拖曳力，所述流体压力梯度力和拖曳力引起所述射孔三维数值模型中岩石颗粒的重新分布；5）重复步骤3～4， 直到计算结果收敛。

Description

油井出砂的三维数值模型设计方法

技术领域

本发明涉及石油工程领域，特别是油井出砂的三维数值模型设计。

背景技术

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

对于实际油藏开采，射孔过程中的流体运动和岩层受力与流体相互作用复杂，油井出砂过程中的砂岩力学特性也较复杂，常规的室内试验方法和等效连续介质分析方法不能反映砂岩颗粒的细观特性，虽然现在一些学者基于土力学理论开始了细观力学分析，但仅仅基于达西渗流的二维离散元模型，不能准确反映油井出砂的发生和发展。特别是常规注水开采中，考虑射孔过程的流体运动主要以径向为主，如果依然在笛卡尔坐标系下研究流体运动，其数值模拟中的流体边界条件很难确定，不能准确地描述射孔过程中的孔隙率变化和流固耦合效应，也无法准确说明出砂的发生和发展过程。

发明内容

本发明的目的是提供一种考虑射孔试验特征，反映射孔试验中砂岩的宏细观力学特性的模型设计方法。

为实现本发明目的而采用的技术方案是这样的，油井出砂的三维数值模型设计方法，包括以下步骤：

1)读入实际油井的几何尺寸，建立柱坐标系下的射孔三维数值模型；

2)根据射孔试验过程的实际情况，将所述射孔三维数值模型划分成若干单元；

3)根据所述射孔三维数值模型的岩石颗粒分布，计算出每一个单元的孔隙率；

4)根据所述孔隙率计算出流体压力梯度力和拖曳力，所述流体压力梯度力和拖曳力引起所述射孔三维数值模型中岩石颗粒的重新分布；

5)重复步骤3～4，直到计算结果收敛。

记录上述过程的所有中间值和最终值，确立射孔出砂模型。

值得说明的是：本发明需要在步骤1)中预先读入油井的相关参数，其来源可以是实际油井中的数据，也可以是实验室模拟的油井数据，其中特别重要的是油井的射孔的数据。步骤1)中，所述射孔包括贯通的内孔，以及内孔周边的砂岩。在油井不断开采过程中，流体不断从内孔中流出，其周边的砂岩会有渗流等情况发生，内孔尺寸发生变化。而射孔出砂的来源主要是砂岩，因此，本发明所要考虑的重点就是内孔周边的砂岩。所述射孔三维数值模型即是将所述内孔和内孔周围的砂岩放在柱坐标系下建立起来的模型。

步骤2)中，划分单元时，主要是将内孔周围的砂岩所在区域划分为若干虚拟的单元，以便于进一步的处理和研究，每一个单元的边界即为柱坐标系下的坐标区间。初始时已经贯通的内孔所在的区域可以不进行划分。

步骤3)中，由于砂岩又是若干岩石颗粒胶结而成，因此，内孔周围的砂岩具有孔隙。本步骤需要计算出每一个单元内砂岩的孔隙率，即单元的孔隙率。

由于流体不断从内孔周围的砂岩中渗出，因此，岩石颗粒必然会受到流体的影响，最终导致孔隙率发生变化。而正是这种变化，带来了本领域实际生产中的诸多难题。现有技术是不能在实验室通过模拟的手段真实反映上述变化的。步骤4)正是需要计算出每一个单元的流体压力梯度力和拖曳力，其依据主要是前一步骤计算出的孔隙率。由于所述流体压力梯度力和拖曳力引起所述射孔三维数值模型中岩石颗粒的重新分布。因此，步骤4)的结果不一定是整个模型的稳定状态。而步骤5)说明了，将步骤4)所获得的岩石颗粒分布作为步骤3)的计算依据，并由此循环若干次，直到计算出的单元的孔隙率不发生变化或其变化小于某个设定的阈值时，模型到达稳定状态，停止循环。

本发明所公开的上述方法，提供了一种能简单、直接和有效的油井出砂三维数值模型设计，该模型通过准确反映每一个流体单元的孔隙率，保证了射孔试验过程中流固耦合作用力计算的有效性，提高了三维细观数值模型的可行性和稳定性。

附图说明

本发明可以通过附图给出的非限定性实施例进一步说明。

图1油井出砂示意图；

图2基于柱坐标系的流体单元划分俯视图；

图3颗粒与环向内边界相交；

图4颗粒与环向外边界相交；

图5与流体单元环向边界相交时颗粒体积计算示意图；

图6颗粒与径向边界相交；

图7颗粒与水平边界相交；

图8颗粒与径向边界和水平边界同时相交；

图9交于环向边界时流体单元内颗粒体积V₀的计算；

图10流体单元孔隙率计算流程图；

图11数值模拟与理论计算的模型切向应力比较；

图12基于PFC3D出砂模型的切向应力分布。

具体实施方式

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

油井出砂的三维数值模型设计方法，包括以下步骤：

1)读入实际油井的几何尺寸，建立柱坐标系下的射孔三维数值模型；

2)根据射孔试验过程的实际情况，将所述射孔三维数值模型划分成若干单元；

3)根据所述射孔三维数值模型中砂岩颗粒分布，计算出每一个单元的孔隙率；

4)根据所述孔隙率计算出流体压力梯度力和拖曳力，所述流体压力梯度力和拖曳力引起所述射孔三维数值模型中岩石颗粒的重新分布；

5)重复步骤3～4，直到计算结果收敛。

参见图1，所述射孔三维数值模型是根据实际油井的射孔的尺寸参数建立的，所述实际油井的射孔包括贯通的内孔和所述内孔周围的砂岩。所述射孔三维数值模型是以所述内孔的中轴线为z轴，来建立的柱坐标系下的三维数值模型。由于射孔被射穿时，内孔周边的砂岩受到一定程度的影响，在出油的过程中，内孔的孔径会不断发生变化。因此，根据岩土力学的基本知识，在模型中出现弹性区域和塑形区域。图1中所示射孔在初始时，内孔的内孔半径为R_i，塑性边界到中心点o的距离为R_c，弹性边界到中心点o的距离为R_o，均记录为柱坐标系下的r值，以便于后续研究。

步骤1)中，读入实际油井的几何尺寸的方式较多，可以选定具体的储层砂岩，得到砂岩颗粒的密度、根据级配曲线得到砂岩的最大颗粒与最小颗粒半径比和最小颗粒半径。根据试验得到平行粘结的法向刚度和切向刚度、颗粒的摩擦系数、颗粒的法向刚度和切向刚度、平行粘结的法向强度和切向强度以及粘结程度的胶结半径比。进一步，计算得到每个流体单元的孔隙率n、沿着径向方向颗粒的平均直径颗粒的平均速度以及流体单元的平均速度u_f,r、考虑流体对颗粒产生的拖曳力和流体压力梯度力，计算每个流体单元的拖曳力系数和流体压力梯度

步骤2)中，考虑流体的连续性，在柱坐标系下将整个模型离散成若干个流体单元，每个流体单元由径向、环向和水平边界控制，计算作用在位于流体单元内的岩石颗粒上的流体压力梯度力和拖曳力。

所述岩石颗粒受流体作用发生运动，引起颗粒位置的重分布，每个流体单元孔隙率n的改变，单元孔隙率的变化又引起拖曳力和流体压力梯度的变化，从而引起颗粒细观力学效应和模型宏观力学特性的变化。数值计算中，对于任意给定的较小时间增量步，每个流体单元孔隙率的变化比较小，流固耦合分析时，引入一个更新参数N_up来确定渗流计算的时间。即进行N_up次力学计算后，再进行流固耦合分析，以此来提高计算机的运行速度和计算效率。

选定具体的油井砂岩，确定具体的砂岩参数，数值计算结果与理论计算结果进行对比分析模型的稳定性，并对流体单元进行不同于步骤1)的划分，重新进行数值计算，以检验模型的稳定性和收敛性。具体地，通过分析不同工况下的模型的切向应力σ_θ分布、平行粘结分布以及颗粒转动来确定，这些状态变量通过绘图来反映三维数值模型的稳定性和收敛性。

实施例1：

本实施例公开具体的油井出砂三维数值模拟方法，主要包括下列步骤，即：

第一步：读入实际油井的数据，建立柱坐标系下的油井三维数值模型。其中，最重要的是根据射孔的具体尺寸和位置，建立柱坐标系下的射孔三维数值模型。

所述射孔三维数值模型具有内孔半径和射孔模型的外半径，以及弹塑性边界。实际中的射孔，即是分布于油井周围的砂岩中，包括贯通的内孔和所述内孔周边的砂岩。而靠近内孔的砂岩为塑形区域，塑形区域外围即为弹性区域。参见图1，图中所示即为柱坐标系下的射孔数值模型，所述射孔的内孔轴线为z轴方向，射孔任一轴截面的中心点为o点。

第二步：根据射孔试验过程的实际情况，读入试验数值，考虑流体的连续性，对所述射孔数值模型进行单元划分。本实施例中，通过设置径向、环向和水平边界条件，以实现单元划分。即在柱坐标系下，通过设定r、φ、z各自的区间值，将所述射孔数值模型划分为若干单元。实施例中，仅将所述内孔与外半径之间的区域划分成若干单元(即图1中R_i和R_o之间的区域，实际工程中是内孔周围砂岩所处的区域)，以便于进一步研究。图2即为所述射孔数值模型的一个轴截面。编号i所示即为第i个流体单元，其是一个具有三个边界范围的三维图形。图2中所能见到的是环向边界和径向边界。值得说明的是，划分的每个单元都有径向、环向和水平向边界，划分的单元越多，则边界相对就越多，消耗计算机的内存越大。由于划分过多的单元会导致计算缓慢，因此，应该根据计算机的处理能力划分适宜数量的单元。

第三步：读入时间t时若干砂岩颗粒的坐标信息、颗粒半径。计算出处于流体单元N_i中所有颗粒的体积，进一步计算出t时流体单元N_i的孔隙率n。实施例中，颗粒的信息可以通过实际测量获得，也可以通过经验模拟设定，但由于真实颗粒较小，通常是根据颗粒级配曲线得到其颗粒直径的分布，在计算及模拟中满足这一级配曲线即可。值得说明的是，颗粒级配曲线是根据颗分试验成果绘制的曲线，采用对数坐标表示，横坐标为粒径，纵坐标为小于(或大于)某粒径的砂土重(累计百分)含量，如实际颗粒的粒径为0.1mm～1.5mm，就可以得到小于各粒径的比重，数值计算中按该比重考虑颗粒均匀分布或按实际位置分布，如颗粒小于0.2mm的比重为5％，则数值计算中生成5％该粒径的颗粒，均匀放在射孔模型中。或者通过监测颗粒实际位置，可以按照实际位置将其放置在模型中。

实施例中，由于所划分的单元的体积已知，需要计算出处于每一个单元内所有颗粒的体积，其重点在于处于单元边界上颗粒体积的计算。参见图9～10，通过以下A～D计算孔隙率n：

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

B)判断颗粒i是否与流体单元N_i的环向边界相交：若否，则颗粒i的体积即为V₀；若是，则进一步判断颗粒i中心点的位置，参见图3～图5分情况计算V₀：当颗粒i与环向边界相交时，相交区域的体积计算参见图5，得到 $V_{s1}=\frac{\mathrm{\π}}{3}(r1-a1)(2r_1^2+a_1{r}_1-a_1^2),V_{s2}=\frac{\mathrm{\π}}{3}(r2-a2)(2r_2^2+a_2{r}_2-a_2^2),$ 式中：V_s1表示位于相交区域内流体单元的体积，V_s2表示位于相交区域内颗粒的体积。r1表示流体单元环向边界代表的半径，r2表示颗粒半径，a1表示流体单元的中心点到流体单元与颗粒相交面的距离(由于附图采用了投影方式，图中a1为中心点到相交线的距离)，a2表示的是颗粒中心点到流体单元与颗粒相交面的距离(由于附图采用了投影方式，图中a2为中心点到相交线的距离)。(1)当颗粒与流体单元的内边界交于颗粒中心点之前时，则(2)当颗粒与流体单元的内边界交于颗粒中心点之后时，则V₀＝(V_s2-V_s1)；(3)当颗粒与流体单元的外边界交于颗粒中心点之前时，则V₀＝(V_s1+V_s2)；(4)当颗粒与流体单元的外边界交于颗粒中心点之后时，则；(设颗粒中心点到射孔中心点距离为d_b，流体单元的环向边界到射孔中心点距离为d_f，若d_b＞d_f，则颗粒与流体单元的边界交于颗粒中心点之前，若d_b＜d_f，则为颗粒与流体单元的边界交于颗粒中心点之后)

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

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

第四步，计算流体单元N_i内作用在颗粒上的流体压力梯度力和拖曳力：

流体压力梯度： $\▿p_j=(150\frac{{(1-n)}^2}{n^2{\tilde{d}}_p^2}{\mathrm{\μ}}_f+1.75\frac{(1-n){\mathrm{\ρ}}_f}{n{\stackrel{\‾}{d}}_p}|{\stackrel{\‾}{v}}_j-u_j\left|\right)({\stackrel{\‾}{v}}_j-u_j);$

拖曳力系数为： ${\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}|{\stackrel{\‾}{v}}_j-u_j|;$

式中，μ_f,ρ_f为流体的运动粘滞系数(m²/s)和密度(kg/m³)，为颗粒的平均直径(m)，n为孔隙率，分别为颗粒运动的平均速度(m/s)和流体实际流动速度(m/s)。

完成一次上述步骤后，即完成了一次力学计算。可以根据第四步的结果，进一步进行流固耦合分析。

由于在力的作用下，颗粒的位置不断发生变化。因此，本发明所公开的模型设计需要模拟这种变化，以便为本领域的研究提供可靠的参数。进一步地，根据第四步获得的结果，(例如，对每个颗粒进行编号，如1～20000，得到初始时刻所有颗粒的位置。如初始时刻的颗粒i坐标(x_i,y_i,z_i)，随着数值计算的进行，得到颗粒i的位移增量(Δx_i,Δy_i,Δz_i)，则颗粒i的当前位移即为(x_i+Δx_i,y_i+Δy_i,z_i+Δz_i))。循环执行第三步和第四步，直到第三步计算出的n值不再发生变化或其变化低于一定阈值时结束。

实施例2

本实施例根据实施例1所公开的方法，进行流固耦合分析。

以实施例1中的流体单元N_i为例，获得t＝t₁时的孔隙率和作用在颗粒上的流体压力梯度力和拖曳力。由于颗粒受力产生运动，引起颗粒位置的重分布，进而流体单元孔隙率发生改变。所述孔隙率的变化又引起拖曳力和流体压力梯度的变化，从而引起颗粒细观力学效应和模型宏观力学特性的变化。作为一种实施方式，时间增量步△t，获得t＝t₁+△t时的孔隙率和作用在颗粒上的流体压力梯度力和拖曳力。根据计算得到的拖曳力系数和流体压力系数得到作用在颗粒上的流固耦合作用力：

作用于砂岩颗粒上的作用力为

$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,$ 式中d_pi为当前颗粒的直径。

在数值计算中，对于任意给定的较小时间增量步，每个流体单元孔隙率的变化比较小。因此，流固耦合分析时，可以引入一个更新参数N_up来确定渗流计算的时间，即是当力学计算N_up次后，才开始渗流计算。如N_up＝10，即进行10次力学计算之后才进行一次渗流分析，以此来提高计算机的运行速度和计算效率。

实施例3

本实施例选定具体的储层砂岩，设置砂岩的最大颗粒与最小颗粒半径比，砂岩颗粒的密度，平行粘结的法向刚度和切向刚度，颗粒的摩擦系数，颗粒的法向刚度和切向刚度、平行粘结的法向强度和切向强度，最小颗粒半径以及粘结程度的胶结半径比，沿着径向方向颗粒的平均直径颗粒的平均速度以及流体单元的平均速度u_f,r。上述参数是根据实际的油井射孔试验得到的，其中的平行粘结刚度和颗粒的刚度是根据工程经验和试验来确定。根据实施例1所公开的方法获得每个流体单元的孔隙率、流体对颗粒产生的拖曳力和压力梯度力，计算每个流体单元的和

流体压力梯度为 $\▿p_j=(150\frac{{(1-n)}^2}{n^2{\stackrel{\‾}{d}}_p^2}{\mathrm{\μ}}_f+1.75\frac{(1-n){\mathrm{\ρ}}_f}{n{\stackrel{\‾}{d}}_p}|{\stackrel{\‾}{v}}_j-u_j\left|\right)({\stackrel{\‾}{v}}_j-u_j)$

拖曳力系数为： ${\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}|{\stackrel{\‾}{v}}_j-u_j|$

式中，μ_f,ρ_f为流体的运动粘滞系数(m²/s)和密度(kg/m³)，为颗粒的平均直径(m)，_n为孔隙率，分别为颗粒运动的平均速度(m/s)和流体实际流动速度(m/s)。

Claims (1)

1.油井出砂的三维数值模型设计方法，其特征在于，包括以下步骤：

1）读入实际油井的几何尺寸，建立柱坐标系下的射孔三维数值模型；

2）根据射孔试验过程的实际情况，将所述射孔三维数值模型划分成N个单元；

3）根据所述射孔三维数值模型的砂岩颗粒分布，计算出每一个单元的孔隙率；

4）根据射孔中流动的流体参数和所述孔隙率，计算出流体压力梯度力和拖曳力：

流体压力梯度：；

拖曳力系数为：；

式中，为流体的运动粘滞系数,为流体的密度，为颗粒的平均直径，n为孔隙率，为颗粒运动的平均速度，为流体实际流动速度；

所述流体压力梯度力和拖曳力引起所述射孔三维数值模型中砂岩颗粒的重新分布；在单元N_k为例，k＝1、2……N，获得时间t＝t₁时的孔隙率和作用在颗粒上的流体压力梯度力和拖曳力，时间增量步△t，获得t＝t₁+△t时的孔隙率和作用在颗粒上的流体压力梯度力和拖曳力;根据计算得到的拖曳力系数和流体压力梯度，得到作用于砂岩颗粒上的作用力为

，式中为当前颗粒的直径；

5）重复步骤3～4，直到计算结果收敛。