一种实时定量化表征纤维封堵岩石孔隙的数值模拟方法

1.本发明涉及一种数值模拟方法，特别是用于封堵岩石孔隙的包含纤维材料的钻井液的数值模拟方法。


背景技术：

2.勘探开发页岩气和深部地质钻探过程中维持井壁稳定离不开钻井与钻井液。维持井内压力大于地层孔隙压力(但不高于破裂压力)是深部地质钻探和海洋钻井较合理的安全开采方式。在此条件下，水基钻井液不可避免地会侵入岩石孔隙。
3.而要解决上述问题，就需要在维持井壁稳定的同时尽可能地降低钻井液侵入程度。加入纤维材料封堵岩石孔隙是增强井壁稳定性的有效方法之一。由于纤维其长径比优势，纤维可以形成网状结构从而形成有效封堵，并阻滞水分的再次侵入。
4.但目前纤维对岩石孔隙的封堵效果多限于物理实验数据，关于纤维封堵岩石孔隙维持井壁稳定的数值模拟研究鲜有报道。钻井液中纤维侵入岩石孔隙后的运移、动态堆积与微观封堵机理并不明确。
5.实验数据问题为物理实验的最大问题，其动态观测难度大、成本高、不易成功，因此完全物理实验来全面掌握纤维封堵微观机理目前不具可行性。
6.此外，孔隙复杂度、流体物理性能和孔隙应力的影响使得纤维封堵岩石孔隙的动态封堵机理研究变的更加复杂和困难。
7.现有技术中对于纤维类在岩石中的应用专利主要集中于纤维的合成和物理实验方面(cn108384520a)。采用的方法为物理封堵实验，通过物理实验测试纤维对于水力裂缝的封堵效果(cn104963673a)。此外，本领域技术人员通过物理实验探究纤维添加模拟纤维对于石油采收率的影响(cn112627788a)，无法实现实时定量分析。此类方法最大的问题为纤维能够改变的参数范围极小，只包含1种纤维参数即纤维长度不同所计算的封堵效率值。
8.同时，现有技术的研究内容不包含流体粘度改变所带来的封堵性能评测和模型验证。以上限制导致不同流体物性特征和多样纤维参数下纤维封堵效率的动态监测和钻井液实时侵入量不易快速确定。


技术实现要素：

9.本发明旨在提供一种实时定量化表征纤维封堵岩石孔隙的数值模拟方法，解决现有针对含有纤维材料的钻井液试验无法实时定量化表征纤维封堵岩石孔隙的问题，同时，能够通过改变多种纤维参数来模拟计算纤维对岩石孔隙的封堵效果。
10.为实现上述目的，本发明采用了下述方案：
11.一种实时定量化表征纤维封堵岩石孔隙的数值模拟方法，以钻井液为液相，以纤维为固相，纤维随钻井液流动，基于cfd
‑
dem耦合方法模拟钻井液中纤维对水平井和垂直井中孔隙的封堵，其中包括，
12.孔隙设定，所述孔隙设定中包括设置至少一条弯折的孔道，孔道的两端分别为入
口和出口，且孔道的孔径在出口处缩小形成孔隙；
13.钻井液设定，所述钻井液设定中包括钻井液粘度设置，钻井液粘度采用钻井液滤液的粘度；
14.纤维设定，所述纤维设定包括纤维几何模型、纤维受力模型和纤维弯折单元，其中：
15.所述纤维几何模型为多个顺序连接的球圆柱体，纤维未发生弯折变形时，多个相同的球圆柱体沿着一条直线依次连接形成一根纤维，且相邻的两个球圆柱体的球面中心重合，相邻的两个球圆柱体的圆柱面轴线同轴；需要说明的是，球圆柱体是指在圆柱体轴线方向的两个端面上各连接一个半球体，且半球体的直径等于圆柱体的直径；
16.所述纤维受力模型包括纤维法向接触力模型、切向接触力模型和法向粘附接触力模型，以及纤维受钻井液粘滞力作用产生的流动阻力模型；
17.所述纤维弯折单元由两个相邻的球圆柱体组成，当纤维与钻井液或另一纤维相互作用发生弯折变形时，纤维弯折单元中两个相邻的球圆柱体发生相对运动，相邻的两个球圆柱体的球面中心产生位移，相邻的两个球圆柱体的圆柱面轴线产生角位移。
18.进一步，所述纤维设定还包括纤维注入方式、纤维运动方向和纤维运动模式，其中，
19.所述纤维注入方式为在孔道入口处随机注入；
20.所述纤维运动方向为在孔道入口处随机方向；
21.所述纤维运动模式包括平移和旋转。
22.进一步，所述钻井液设定中还包括钻井液流动模式的设置，钻井液流动模式为层流。
23.进一步，实时定量化表征纤维封堵岩石孔隙的数值模拟方法还包括模拟计算在不同流体物理性能下，不同纤维材料参数对孔隙封堵效率和钻井液侵入量的影响，其中，
24.所述流体物理性能包括不同数值的钻井液滤液粘度；
25.所述纤维参数包括纤维长径比、纤维直径、纤维浓度、纤维密度、硬度和纤维泊松比。
26.进一步，所述孔隙封堵效率为进入水平井和垂直井孔隙的纤维总量；所述钻井液侵入量为通过水平井和垂直井孔隙出口的钻井液体积。
27.优选的，所述孔隙设定中，孔道与钻井液、纤维接触的壁面边界网格采用楔形网格，其余主体部分为四面体网格，且从孔道入口处到孔道出口处，网格逐渐致密；
28.所述钻井液设定中，钻井液与孔道壁面、纤维接触的边界网格采用楔形网格，其余主体部分为四面体网格；
29.所述纤维设定中，纤维与钻井液、孔道壁面接触的边界网格采用楔形网格，其余主体部分为四面体网格。
30.关于纤维受力模型，所述法向接触力模型中，法向接触力定义为：其中k
n
为法向接触刚度，c
n
为法向阻尼系数，s
n
为重叠区域面积，为重叠区域的时间导数；
31.法向阻尼系数定义为：其中，η为阻尼率，m
*
为纤维碰撞的有效质量，纤维碰撞的有效质量描述为：
[0032][0033]
其中，m1和m2是相互碰撞的两根纤维的质量，单位为kg，m是与孔道边界接触的纤维的质量，单位为kg；
[0034]
所述切向接触力模型中，切向接触力定义为：其中，μ为滑动摩擦系数，f
n
为接触点的法向接触力，为接触重叠区的切向分量；
[0035]
所述法向粘附接触力模型中，法向粘附接触力定义为：其中，f
n，adh
是法向粘附接触力，单位为n，δ
adh
为粘附距离，单位为m，r
adh
是粘附力刚度与接触载荷刚度的比值。
[0036]
关于纤维受钻井液粘滞力产生的流动阻力，所述流动阻力模型定义为施加在纤维上的恒定力矩，
[0037]
其中μ
roll
为滚动阻力系数，d
roll
为纤维滚动直径，单位为m，ω为纤维的角速度，单位为rad/s；
[0038]
滚动直径分为如下类型：
[0039][0040]
其中d
pp1
和d
pp2
为相互接触的两根纤维的滚动直径，单位为m，d
pb
为接触的纤维与孔道壁面的滚动直径，单位为m。
[0041]
关于纤维弯折单元，所述纤维弯折单元中，法向力和切向力分别定义为：
[0042][0043][0044]
其中，和为纤维弯折单元的法向和切向刚度，和为纤维弯折单元的法向和切向阻尼系数，和为纤维弯折单元中两个球圆柱体末端的半球形中心点之间的法向和切向相对位移，和为纤维弯折单元中两个球圆柱体末端的半球形中心点之间的法向和切向相对速度；
[0045]
记和是纤维弯折单元中相连接的两个球圆柱体上圆柱体轴线方向上的单位矢量，则法线方向由单位矢量确定，相对位移矢量d
rel
分解为两个垂直分量，即法向分量和切向分量：
[0046]
[0047][0048][0049]
纤维弯折单元各部分相对速度为：
[0050][0051][0052]
纤维弯折单元的法向刚度和切向刚度定义为：
[0053][0054][0055]
其中，e
vb
为纤维弯折单元杨氏模量，a为纤维弯折单元横截面，单位为m2，l
vb
为纤维弯折单元长度，单位为m，ζ
vb
为纤维弯折单元泊松比；
[0056]
纤维弯折单元杨氏模量设置为：
[0057][0058]
其中，为杨氏模量比，e
p
为纤维的杨氏模量；
[0059]
纤维弯折单元的法向和切向阻尼系数确定为：
[0060][0061][0062]
其中，为纤维弯折单元的有效质量，单位为kg，η
vb
为纤维弯折单元阻尼比；
[0063]
纤维弯折单元的弯矩和扭矩的计算公式为：
[0064][0065][0066]
其中，i和j为惯性矩，θ
ben
和θ
tor
为弯曲和扭转角位移，单位为
°
，和为弯曲和扭曲的相对速度，单位为rad/s。
[0067]
本发明提供了实时定量化表征纤维封堵岩石孔隙的数值模拟方法，与现有技术相比，本发明的数值模拟方法具有如下特点：
[0068]
(1)可实时且定量的表征(定量表征是指计算孔隙出口处的累积纤维总数和孔隙出口处的钻井液侵入量)纤维类材料对于岩石孔隙的封堵过程及封堵效果；
[0069]
(2)纤维参数调节和变化范围广泛且细密，非物理实验所能达到。纤维参数包含：纤维直径、纤维长径比、纤维浓度、纤维密度、纤维旋转、纤维泊松比、纤维硬度和纤维是否可弯折。
[0070]
(3)纤维在流体中(钻井液中)的物理性能可调，主要指流体(钻井液)温度、密度和塑性粘度可调。
[0071]
(4)计算结果准确，且成本消耗低，对于物理实验而言，成本非常高且无法定量实时监测，同时存在非可控因素导致的实验失败。
附图说明
[0072]
图1为本发明中弯折孔道网格划分示意图；
[0073]
图2为本发明中纤维在流体中受力图；
[0074]
图3为本发明中纤维弯折单元设置及方位图；
[0075]
图4为本发明中纤维弯折图；
[0076]
图5为本发明中纤维移动和封堵过程可视化示意图；
[0077]
图6为本发明中不同浓度纤维对于孔隙封堵效果对比示意图。
具体实施方式
[0078]
本发明目的的实现、功能特点及优点将结合实施例以及说明书附图做进一步说明。下面结合附图对本发明的构思做进一步的描述，但本发明的保护范围不局限于以下所述。
[0079]
如图1～图6，为一种实时定量化表征纤维封堵岩石孔隙的数值模拟方法，数值模拟方法包括以下特点：
[0080]
(1)以钻井液为液相，以纤维为固相，纤维随钻井液流动的同时发生弯折和变形，基于cfd
‑
dem耦合方法模拟计算钻井液中纤维对水平井和垂直井孔隙封堵效果(孔隙封堵效果目前没有准确的定义，本发明中采用孔隙出口的纤维封堵数量评价，孔隙出口累积的纤维数越多，孔隙封堵效果越好)，优选不同井况下的流体物理性能参数和纤维材料参数，建立纤维在流体中(钻井液中)的流固耦合模型，以及流体(钻井液)与孔道壁面、纤维与纤维之间的碰撞模型；
[0081]
(2)模拟计算进入水平井和垂直井孔隙的纤维总量和流体总量；
[0082]
(3)由于井壁生成泥饼，进入孔隙中的流体为钻井液滤液，即钻井液经过了过滤，所以基于钻井液滤液粘度进行设定，而非钻井液粘度，基于微量流体粘度仪来确定流体粘度；
[0083]
(4)计算不同流体物理性能(流体粘度变换)下，纤维材料参数(纤维长径比、纤维直径、纤维浓度、纤维密度、硬度和纤维泊松比等)对岩石孔隙封堵效率和钻井液侵入量的影响；
[0084]
(5)模拟计算过程中，发展基于隐式时间积分的四面体网格和楔形网格单元，主要为楔形网格设计，保证流体与纤维接触界面计算准确性，设定流体和纤维边界条件和初始条件；
[0085]
(6)纤维的几何模型为多个顺序连接的球圆柱体，纤维未发生弯折变形时，多个相同的球圆柱体沿着一条直线依次连接形成一根纤维，且相邻的两个球圆柱体的球面中心重合，相邻的两个球圆柱体的圆柱面轴线同轴(如图2)，此外，纤维由于本身的高长径比特征，纤维在流动过程中会发生弯折、拉伸，因此设置弯折单元来建立纤维弯折、移动和拉伸模型；
[0086]
(7)建立释放速度大小、方向以及尺寸各不相同的纤维材料参数模型，计算实时的流体、纤维速度和压力参数；
[0087]
(8)计算不同井况下、不同流体粘度下的模拟结果，通过理论公式推导验证和实验验证，增加模型准确度和适用性，最终对纤维材料封堵垂直井和水平井的孔隙进行评价和
分析。
[0088]
需要说明的是，本发明是基于cfd
‑
dem耦合方法，因此本领域技术人员熟知的流固耦合步骤不再赘述，本发明的核心内容在于纤维的几何模型，受力模型(纤维与纤维的接触碰撞和纤维与孔道壁面的接触碰撞)，纤维在流体中的弯折变形模型(纤维弯折单元)以及不同纤维参数和流体(钻井液)参数下的孔隙封堵效果模拟。
[0089]
其中，孔道的孔隙壁面为无滑移界面的固定壁面，纤维碰撞具有反射性，壁面反射系数分为法向和切向系数。模型中的法向接触力由线性弹力和阻尼力组成。而纤维碰撞之间的切向接触力则使用粘附摩擦系数和滑动摩擦系数计算；模型的网格细度基于纤维碰撞随机度划分；四面体网格是主体网格，主要用于流体流动和纤维迁移，楔形网格被用作边界网格，以更准确地区分边界层的接触和碰撞。
[0090]
纤维的接触力基于cundall和strack模型。模型中的法向接触力由线性弹力和阻尼力组成。法向接触力定义为：
[0091][0092]
其中k
n
为法向接触刚度，cn为法向阻尼系数，s
n
为重叠区域面积，为重叠区域的时间导数。
[0093]
法向阻尼系数可定义为：
[0094][0095]
其中η为阻尼率，m
*
为纤维模型碰撞的有效质量，纤维碰撞的有效质量可以描述为：
[0096][0097]
其中，m1和m2是离散颗粒的质量(kg)，m是与边界接触的粒子的质量(kg)。
[0098][0099]
其中ε是恢复原状的系数。
[0100]
切向力基于摩尔库伦准则，且切向力在接触点可描述为：
[0101][0102]
其中，μ为滑动摩擦系数，f
n
为接触点的法向力，为接触重叠区的切向分量。
[0103]
弹塑性法向接触模型适用于模拟非粘性或干燥的碰撞单元。纤维材料可能具有多种水分含量，导致它们可能与其它表面接触。为了捕获这种行为，斥力法向力由法线吸引力来补充，以准确地预测其流动特性。
[0104]
本质上，线性粘附力模型就像具有引力的线性弹簧。如果接触距离大于最小粘合接触力距离，则粘合力为零。并且，它与该距离和实际纤维接触距离之间的差异成比例增加。比例因子定义为粘合刚度。粘附力模型可以定义为：
[0105][0106]
其中，f
n，adh
是法向粘附接触力(n)，δ
adh
为粘附距离(m)，r
adh
是粘附力刚度与接触载荷刚度的比值。
[0107]
纤维流动阻力设置原理为将恒定力矩施加到离散单元(纤维)上，其公式可表示为：
[0108][0109]
其中μ
roll
为滚动阻力系数，d
roll
为离散单元滚动直径(m)，ω为离散单元的角速度(rad/s)。
[0110]
此外，基于纤维碰撞的类型(纤维与纤维碰撞或纤维与孔道壁面碰撞)，计算过程中滚动纤维(滚动是指纤维平移和旋转的组合)的直径分为如下类型：
[0111][0112]
其中d
pp1
和d
pp2
为接触的两根纤维的滚动直径(m)，d
pb
为接触的纤维与孔道壁面的滚动直径(m)，需要指出的是，此处的纤维是指纤维几何模型整体，同时滚动直径是指纤维几何模型的质心与接触点之间的距离的2倍，即半径的2倍。
[0113]
在流体计算方面，四面体网格优势明显，可快速准确形成模型内部网格，记录和监测流体(钻井液)和离散单元(纤维)的运动信息。但是，对于接触面的稳定性和准确性，楔形网格更加适合。因此，模型网格分为两种类型，分别是楔形网格和四面体网格。对于接触和壁面内容，使用细致的楔形网格，更好的记录离散单元(纤维)运动和碰撞细节。对于其余部分，采用稳定性和准确性兼容的四面体网格。此外，实验观测和记录的重点为离散单元(纤维)堆积和封堵的过程。从孔道入口至孔道出口，网格逐渐细密。因此，在纤维堆积出口处网格设置最为致密，从而保证堆积过程的准确捕捉。经过网格独立性计算，此数量设置的网格可保证计算时间和计算准确度位于最优区间。
[0114]
参数设置主要分为两部分内容，钻井液(流体)参数和纤维(固体)参数。流体设置为粘度不同的水溶液，主要用于分析不同粘度下的流体运移对于离散单元(纤维)的影响。不考虑温度对于离散单元(纤维)运移的影响，因此温度设置为恒温293k。重力方向为y的负方向，取值为
‑
9.81m/s2。固体离散相设置为纤维，默认长径比为1∶10。纤维密度为1300kg/m3，杨氏模量设置为1e6n/m2。纤维注入方式为随机注入，进入孔道时的纤维运移方向是随机的，且运移模式为平移和旋转。流体流动模式为层流，且离散单元追踪模式为非稳态追踪。在流体参数设置一致的情况下，分析纤维长径比(1∶5、1∶10、1∶15、1∶20)、纤维直径(1/20、1/15、1/10、1/5)、纤维浓度(0.5％、1％、1.5％、2％)、纤维密度、硬度和纤维泊松比不同时的离散单元堆积和封堵效果。
[0115]
纤维在流动过程中会自然摆动，为准确计算纤维流动状态，设置纤维弯折单元，从而保证纤维在流动过程中可随流体摆动、弯折和拉伸。纤维弯折单元是具有弹性和粘性性
质的辅助实体，它由两个相邻的球圆柱体元素组成。最初纤维没有发生变形，相邻球圆柱体的半球形末端的中心重合(即图2中两个半球形中心重合形成一个完整的球体)，纤维弯折单元未发生变化。但是，当纤维移动变形时，相邻球形圆柱体之间的相对运动会在纤维弯折单元上产生线性位移和角度位移变形。作为响应，在相邻的球圆柱体上感应并施加力和力矩，以抵抗法向而后切向受力以及弯曲和扭转变形。切向力和法向力计算公式为：
[0116][0117][0118]
其中，和为纤维弯折单元的法向和切向刚度，和为纤维弯折单元的法向和切向阻尼系数，和为纤维弯折单元末端半球形中心点之间的法向和切向相对位移，和为纤维弯折单元中两个球圆柱体末端的半球形中心点之间的法向和切向相对速度。
[0119]
记和是在纤维弯折单元中连接的两个球圆柱体中圆柱体轴线方向上的单位矢量，则法线方向由单位矢量可确定。相对位移矢量d
rel
然后将分解为两个垂直分量，即法向分量和切向分量：
[0120][0121][0122][0123]
纤维弯折单元各部分相对速度为：
[0124][0125][0126]
纤维弯折单元的法向刚度和切向刚度定义为：
[0127][0128][0129]
其中，e
vb
为纤维弯折单元的杨氏模量，a为纤维弯折单元横截面(m2)，l
vb
为纤维弯折单元长度(m)，ζ
vb
为纤维弯折单元泊松比。
[0130]
纤维弯折单元杨氏模量设置为：
[0131][0132]
其中，为杨氏模量比，e
p
为纤维的杨氏模量。
[0133]
纤维弯折单元的法向和切向阻尼系数确定为：
[0134][0135][0136]
其中，为纤维弯折单元的有效质量(kg)，η
vb
为弯折单元阻尼比。
[0137]
因此，纤维弯折单元的弯矩和扭矩的计算公式为：
[0138][0139][0140]
其中，i和j为惯性矩，θ
ben
和θ
tor
为弯曲和扭转角位移(
°
)，和为弯曲和扭曲的相对速度(rad/s)。
[0141]
影响封堵效率的纤维参数主要包含纤维浓度、纤维长径比、纤维直径、纤维密度、软硬度、旋转以及纤维泊松比。可以逐项讨论不同参数下的纤维特征对孔隙封堵效率的影响。通过对比不同参数下的封堵效率，能够量化纤维特征参数对封堵效果起到正相关还是负相关的作用。此外，讨论和计算可以量化不同纤维特征对封堵效率的影响。
[0142]
对于纤维类的离散单元，纤维浓度上升对岩石封堵效果的研究内容并不多见。影响封堵最重要的因素是离散单元的浓度。随着纤维浓度的上升，封堵效果一定是上升的。模型设置的纤维浓度分别为0.5％，1％，1.5％和2％。结果显示，随着纤维浓度的上升，累积的封堵纤维数的确是逐渐上升的。在本发明的模型中，随着纤维浓度从0.5％上升至1％，1.5％和2％，封堵效率分别提升69.40％，139.12％，254.89％。此外，从图6中可以看出，封堵效率的增长并非是线性的。同时，不同浓度下的纤维封堵效率差距并非一成不变。随着计算步长的增加，纤维封堵效率差距逐渐增大。
[0143]
以上仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。