CN112883661A - 一种碎软低渗油气储层的压裂模拟方法 - Google Patents
一种碎软低渗油气储层的压裂模拟方法 Download PDFInfo
- Publication number
- CN112883661A CN112883661A CN202110131767.8A CN202110131767A CN112883661A CN 112883661 A CN112883661 A CN 112883661A CN 202110131767 A CN202110131767 A CN 202110131767A CN 112883661 A CN112883661 A CN 112883661A
- Authority
- CN
- China
- Prior art keywords
- fracture
- node
- reservoir
- unit
- fracturing
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 97
- 238000004088 simulation Methods 0.000 title claims abstract description 35
- 239000011159 matrix material Substances 0.000 claims abstract description 47
- 239000011148 porous material Substances 0.000 claims abstract description 43
- 230000008569 process Effects 0.000 claims abstract description 42
- 238000004364 calculation method Methods 0.000 claims abstract description 31
- 239000000463 material Substances 0.000 claims abstract description 17
- 230000000694 effects Effects 0.000 claims abstract description 9
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 7
- 206010017076 Fracture Diseases 0.000 claims description 139
- 239000012530 fluid Substances 0.000 claims description 74
- 230000006378 damage Effects 0.000 claims description 38
- 238000006073 displacement reaction Methods 0.000 claims description 30
- 239000007787 solid Substances 0.000 claims description 29
- 239000003245 coal Substances 0.000 claims description 23
- 239000013598 vector Substances 0.000 claims description 23
- 238000002474 experimental method Methods 0.000 claims description 13
- 230000008878 coupling Effects 0.000 claims description 11
- 238000010168 coupling process Methods 0.000 claims description 11
- 238000005859 coupling reaction Methods 0.000 claims description 11
- 239000011435 rock Substances 0.000 claims description 11
- 238000010008 shearing Methods 0.000 claims description 11
- 230000035699 permeability Effects 0.000 claims description 10
- 239000004576 sand Substances 0.000 claims description 6
- 206010053206 Fracture displacement Diseases 0.000 claims description 4
- 230000006835 compression Effects 0.000 claims description 4
- 238000007906 compression Methods 0.000 claims description 4
- 238000002247 constant time method Methods 0.000 claims description 4
- 238000000926 separation method Methods 0.000 claims description 4
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 4
- 230000001133 acceleration Effects 0.000 claims description 3
- 238000003491 array Methods 0.000 claims description 3
- 238000003776 cleavage reaction Methods 0.000 claims description 3
- 230000005489 elastic deformation Effects 0.000 claims description 3
- 238000001914 filtration Methods 0.000 claims description 3
- 230000005484 gravity Effects 0.000 claims description 3
- 230000005012 migration Effects 0.000 claims description 3
- 238000013508 migration Methods 0.000 claims description 3
- 239000002245 particle Substances 0.000 claims description 3
- 230000000149 penetrating effect Effects 0.000 claims description 3
- 230000007017 scission Effects 0.000 claims description 3
- 239000000758 substrate Substances 0.000 claims description 3
- 238000009864 tensile test Methods 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- 238000009423 ventilation Methods 0.000 claims description 3
- 230000010076 replication Effects 0.000 claims description 2
- 239000004215 Carbon black (E152) Substances 0.000 claims 2
- 229930195733 hydrocarbon Natural products 0.000 claims 2
- 150000002430 hydrocarbons Chemical class 0.000 claims 2
- 208000010392 Bone Fractures Diseases 0.000 description 98
- 239000007789 gas Substances 0.000 description 18
- 238000002347 injection Methods 0.000 description 10
- 239000007924 injection Substances 0.000 description 10
- 239000000243 solution Substances 0.000 description 9
- 238000002591 computed tomography Methods 0.000 description 4
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 description 4
- 230000009471 action Effects 0.000 description 3
- 238000010586 diagram Methods 0.000 description 2
- 230000000977 initiatory effect Effects 0.000 description 2
- 239000007788 liquid Substances 0.000 description 2
- 239000003345 natural gas Substances 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- OKTJSMMVPCPJKN-UHFFFAOYSA-N Carbon Chemical compound [C] OKTJSMMVPCPJKN-UHFFFAOYSA-N 0.000 description 1
- 230000004075 alteration Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 229910052799 carbon Inorganic materials 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 230000008676 import Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000009533 lab test Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000006386 neutralization reaction Methods 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 230000021715 photosynthesis, light harvesting Effects 0.000 description 1
- 230000000644 propagated effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
-
- E—FIXED CONSTRUCTIONS
- E21—EARTH OR ROCK DRILLING; MINING
- E21B—EARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
- E21B43/00—Methods or apparatus for obtaining oil, gas, water, soluble or meltable materials or a slurry of minerals from wells
- E21B43/25—Methods for stimulating production
- E21B43/26—Methods for stimulating production by forming crevices or fractures
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A10/00—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
- Y02A10/40—Controlling or monitoring, e.g. of flood or hurricane; Forecasting, e.g. risk assessment or mapping
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Computer Hardware Design (AREA)
- Geology (AREA)
- Mining & Mineral Resources (AREA)
- Fluid Mechanics (AREA)
- Life Sciences & Earth Sciences (AREA)
- General Engineering & Computer Science (AREA)
- Geometry (AREA)
- Evolutionary Computation (AREA)
- Computing Systems (AREA)
- Pure & Applied Mathematics (AREA)
- Algebra (AREA)
- Mathematical Physics (AREA)
- Environmental & Geological Engineering (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Geochemistry & Mineralogy (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明提供了一种碎软低渗油气储层的压裂模拟方法，包括以下步骤；S1、编制碎软低渗储层结构面三维网络模型程序，建立碎软低渗油气储层压裂的数值计算模型；S2、建立储层韧性破坏‑渗流本构模型，分别针对裂缝单元、实体单元开发基于显示和隐式时间积分的求解算法；S3、数值模型中结构面处以及储层基质内嵌入裂缝单元，合并孔压节点；S4、获取并输入材料参数、设置边界条件和初始条件；S5、求解模型并输出结果，S6、分析模拟效果，优化压裂工艺参数，进而提高碎软低渗储层中的油气产量。本发明可用于分析不同压裂工程参数、不同工程地质参数影响下压裂缝网的形态，优化压裂参数、实现油气资源高效开采，具有较大的应用前景。
Description
技术领域
本发明涉及天然气与石油开发技术领域，具体涉及一种碎软低渗油气储层的压裂模拟方法。
背景技术
我国天然气、石油的对外依存度已经达到43.1％和67.4％，使得能源安全面临巨大风险。而我国2000m以浅的煤层气资源储量高达30万亿立方米，高效开发这一非常规清洁能源将有力促进实现“2060年实现碳中和”的目标。但是，中国大部分地区的煤层具有破碎、软弱、极低渗透率的特征，这些特征会对煤层气高效抽采带来巨大挑战。
在沁水煤田赵庄矿煤层气抽采实践过程中表明：地面直井压裂和水平井分段压裂过程中，水力主裂缝通过煤层中的结构面扩展延伸，次级裂缝主要通过直接压裂煤基质而向前传播；碎软煤显著的韧性破坏特征会导致压裂裂缝网络的面积很小，其很低的力学强度又会引起钻井井壁稳定性很差，导致煤层气产量很低。因而采用合理的数值模拟方法预测储层中的压裂缝网形态，对于优化压裂施工工艺、提高油气产量至关重要。
近年来，国内外学者主要针对页岩、致密砂岩等脆性显著，且无明显结构面的介质开展压裂模拟研究。这些研究大多基于连续介质假设，并采用线弹性断裂力学理论分析在无结构面、或者有极简单结构面分布条件下I型压裂裂缝的传播规律。或者假设储层中缝网的形成是由于压裂裂缝沟通了储层中原有的天然结构面，且天然结构面彼此平行。
但是碎软低渗储层中具有丰富且分布复杂的结构面、以及强烈的韧性破坏特征。前者会导致水力主裂缝沿着结构面传播，而次要裂缝会压穿储层基质，形成复杂的缝网结构，更重要的是，碎软低渗储层的韧性断裂特征会导致大部分注液能量用于消耗与裂缝扩展无关的非弹性能，同时，储层基质一旦其出现塑性损伤会导致压裂液大量漏失，二者共同导致高压压裂液能量耗散，致使储层中压裂裂缝网络面积很小，而现有的压裂工程数值计算模型普遍将碎软低渗储层简化为连续介质、脆性材料，严重影响数值模拟结果的准确性。
发明内容
本发明所要解决的技术问题在于针对上述现有技术的不足，提供一种碎软低渗油气储层的压裂模拟方法，以解决上述背景技术中提出的问题。
为解决上述技术问题，本发明采用的技术方案是：一种碎软低渗油气储层的压裂模拟方法，包括以下步骤；
S1、编制碎软低渗储层结构面三维网络模型程序，确定模型中结构面的几何参数，建立碎软低渗油气储层压裂的数值计算模型；
S2、建立储层韧性破坏-渗流本构模型，分别针对裂缝单元、实体单元开发基于隐式和显式时间积分的求解算法；
S3、数值模型中结构面处以及储层基质内嵌入裂缝单元，合并孔压节点模拟宏细观裂缝的随机分布，具体为在储层结构面处及基质内的单元边界处嵌入零厚度的裂缝单元，合并具有相同节点坐标的裂缝单元上的孔压节点，模拟储层中宏细观裂缝的随机分布；
S4、获取并输入材料参数、设置边界条件和初始条件；
S5、求解模型并输出结果，设定需要输出的参数，包括裂缝张开与剪切位移分布、裂缝传播路径、裂缝内流体压强、裂缝内流体流速、裂缝单元中弹性能与非弹性能、断裂能混合比；完整煤岩体中应力、塑性应变、损伤变量、孔隙流体压强、流体运移速率，以及各节点处孔隙流体压强随时间的演化规律；
S6、分析模拟效果，优化压裂工艺参数，进而提高碎软低渗储层中的油气产量。
优选的，在S1中，是采用CT实验方法确定碎软低渗储层中主要结构面的分布，以及倾角、间距、长度的几何参数，再根据拟建立数值计算模型的大小得出x、y、z三个方向上Voronoi多边形的个数，估算散点在三个方向上的浮动范围，由此在matlab软件中生成一组三维点集，取这组点集中的任意一点，计算距离该点最近的另外一点，以这两点的连线构成Delaunay三角形的一边，然后再找出与这条边距离最近且满足Delaunay三角剖分中的外接圆准则的一点作为第三点，构成一个Delaunay三角形，在此基础上，以Delaunay法的空间外接球准则找出第四点，这四个点构成一个Delaunay四面体；
由上述Delaunay四面体相邻边的垂直平分线所围成的图形即为Voronoi多边形或多面体，通过以上方法，可生成储层结构面的三维网络模型，通过多系修正散点在三个方向上的浮动范围，即可得到符合CT实验结果的3D结构面网络模型。
优选的，在S2中，包括储层基质的塑性损伤及压裂液滤失计算，储层结构面的韧性断裂过程及压裂液流动模拟计算和多物理过程耦合求解计算。
优选的，储层基质的塑性损伤及压裂液滤失计算的具体是：
储层基质在压裂过程中会出现拉、剪两种基本破坏形式，拉、剪有效应力分别为：
其中σi和ni(i＝1，2，3)分别表示主应力大小和方向；
总应力矢量
与总应变矢量εs的关系为：
其中，E0为弹性刚度矩阵；
pw是孔隙流体压强；
I为单位矩阵；
α为有效应力系数；且有α＝1-Kb/Ks，Ks是固体颗粒的体积模量，Kb是多孔介质的排水体积模量；
张量d由拉、压损伤变量d(+)与d(-)表示：
d＝d(+)N(+)+d(-)(I-N(+))
其中，假设储层基质在受拉时出现弹性损伤，以
达到抗拉强度作为张拉损伤起始的判据，且当抗拉强度为0时基质完全产生张拉破坏，通过紧凑拉伸试验获取张拉损伤演化规律；
假设储层在受剪时出现塑性损伤，其判据为
其中，A＝(σb0/σc0-1)/(2σb0/σc0-1)，
C＝3(1-Kc)/(2Kc-1)，peff＝trace(σs,eff)/3,
Seff＝σ-peffI，
是三个方向的最大主应力，σb0/σc0是三轴与单轴压缩条件下屈服应力之比，
和
是等效拉压应力；Kc根据经验取值为0.667。
塑性应变可由塑性势函数通过迭代求解得到，
其中，δ为材料参数，根据经验取值为0.1；σUTS为张拉强度；ψ为剪胀角。
储层基质在发生塑性损伤前后的渗透系数表达式为：
优选的，储层结构面的韧性断裂过程及压裂液流动模拟计算具体为：
地层中压裂裂缝会受到储层结构面的影响而出现弹性变形-混合型断裂的过程，在弹性变形阶段，应力和应变之间的关系为：
σc＝D0,cεc
当达到以下应力条件时，结构面的拉、剪强度开始下降；
其中，σc,l和
中l表示n，s，t，分别为法向和两个切向方向，σc,l和
表示牵引力或者牵引力峰值；
D0,c是弹性刚度矩阵；εc是应变矢量，其与分离位移的关系为εc＝S/T0,T0是裂缝单元的本构厚度，在ABAQUS/CAE的Property/Section模块中设置该值为实体单元特征尺寸的0.01倍；
在峰值应力
之后，储层结构面发生混合型韧性断裂过程，其应力σc,l按照如下规律降低：
其中，Sn，Ss和St是变量，表示法向和两个切向方向的位移；Γn,ΓS是断裂能常数，表达式为：
式中Gn与GS是法向与总的切向断裂能，且有GS＝Gs+Gt，Gs和Gt是两个切向方向的断裂能；β与γ表示纯I型和纯II型断裂实验所得应力-应变曲线的形状参数；m和n与β和γ之间的关系为：
其中，χn＝sn,p/sn,χS＝sS,p/sS，sn,p，sS,p以及sn，sS分别是纯I型和纯II型断裂实验所得应力-应变曲线的峰值位移与断裂位移，分别采用紧凑拉伸、贯穿剪切实验确定材料参数，
引入立方定律反映压裂液的切向流动，考虑煤岩基质塑性损伤引起压裂液的滤失或法向流动，建立压裂液质量守恒方程如下：
其中，Qc是总流量；w是裂缝宽度，可由等式(11)～(13)计算得到；μ是压裂液动力粘度；
是切向流的流体压力梯度；ρw是压裂液密度；g是重力加速度；pn,cen和pn,boun是裂缝中心处和边缘处的流体压力。
优选的，多物理过程耦合求解具体是：压裂模型可分为2个部分：储层基质、储层结构面的流固耦合模型，
对于储层基质，压裂模拟过程中需要全耦合求解的主变量包括：
xT＝[ui,1,ui,2,…,ui,i；ps,1,ps,2,…,ps,i；Qs,1,Qs,2,…,Qs,i]
其中，ui为每个实体单元节点位移，ps,i为每个实体单元孔隙流体压强，Qs,i为流入每个实体单元的流体流量；
固体部分方程采用有限元法离散，流体部分方程采用有限体积法离散，采用Newton-Raphson迭代方法对离散后的方程进行隐式求解，已知第j迭代步求解第(j+1)步未知变量的表达式为：
Jacobi(xj)dxj+1＝-Rj
其中，Jacobi(xj)为第j迭代步时的雅可比矩阵，Rj为第j迭代步时的残差。
迭代收敛的条件：
||R||2＜tolerance并且||x||2＜tolerance
其中||||2为向量的二型范数，当容差tolerance取值为10-8时，实体单元的一个计算步只需要2～5次迭代即可收敛；
对于裂缝单元，其节点力与实体单元的节点力相等，在此基础上显式求解裂缝单元的节点位移sn，ss，st，得到单元张开位移w，结合质量守恒方程以及N-S方程，计算pn,cen值，通过立方定律得到结构面中切向流流速；本发明假设实体单元中的孔隙流体压强ps,i与裂缝边缘处的流体压强pn,boun相等，结合实体单元渗透系数k的表达式，显示求解结构面中法向流的流速。
优选的，在步骤S3中，包括如下内容：确定零厚度裂缝单元的尺寸；零厚度裂缝单元的嵌入和孔压节点合并。
优选的，确定零厚度裂缝单元的尺寸的具体是，裂缝单元的尺寸需小于断裂过程区的长度L，表达式为：
其中，μ为材料泊松比。
优选的，零厚度裂缝单元的嵌入包括更新实体单元节点编号和组装裂缝单元：在步骤S1的基础上建立数值模型，并对其进行全局网格划分，在Job模块中生成对应的.inp脚本文件，在该文件中，记最大节点编号为nmax，最大单元编号为emax；
将第ni(ni≠1)个节点复制为b个，b为共用单元的个数，在坐标相同的前提下，保持第1个节点号与初始值相同，其后的第k个节点的编号修改为：
其中，Ten(nmax)为nmax的十进制位数；
经过上述变换，第k个节点被复制了b次，编号依次为n1,n2,…,nb；而共用第k个节点的单元也有b个，单元编号依次为e1,e2,…,eb，这样在.inp文件中，通过将新的节点编号复制在对应单元编号位置处，即可使两个有序数组Elem＝(e1,e2,…,eb)和Node＝(n1,n2,…,nb)形成一一对应的关系，由此得到更新节点编号的实体单元，
eb→nb,其中eb∈Elem,nb∈Node
遍历原始.inp文件中相邻实体单元之间编号相同的节点，将其复制，并按照下式确定的顺序组成新的数组：
Group＝(ej,ek,nShared01,nShared02,nShared03,nShared04,Node01,Node02,Node03,Node04)
数组Group中的ej和ek为相邻单元的单元编号；nShared01～nShared04为原始.inp文件中的相邻单元之间的相同编号的节点；Node01～Node04为根据等式(2)重新编号后的节点号数组；
在Node01～Node04中，仅保留由nShared01～nShared04更新的节点编号，然后逆序排列，即可得到单元ej和ek之间的裂缝单元的节点编号；
而后按照下式赋予裂缝单元的单元编号：
ecoh＝1000+emax
在.inp文件中写入所有裂缝单元的编号及其节点编号，在该字符前添加关键字“*Element type＝COH3D8”，同时在关键字“*Part”之后“*Assembly”之前添加关键字“*Elset＝cohElem02”，将所有裂缝单元组成一个名称为“cohElem02”的集合；
然后将裂缝单元嵌入步骤S1建立的储层结构面中，具体为将更新后的裂缝单元编号ej，以及对应的节点坐标nj(xj,yj,zj)复制在Excel中，采用内积法求得该单元的方向向量dip,coh＝(d1,d2,d3)。若dip,coh与步骤S1中储层结构面的方向向量dip共线、任一节点在储层结构面(即泰森多边形)上、而且任一节点都落在泰森多边形范围之内，则判定单元ej在储层结构面内，将所有符合上述三个条件的裂缝单元和节点编号复制在.inp文件中，在关键字“*Part”之后“*Assembly”之前添加关键字“*Elset＝cohElem01”，将所有储层结构面上的单元组成一个名称为“cohElem01”的集合，而后将.inp文件导入ABAQUS/CAE，在Mesh/Assign Element type模块下的将所有裂缝单元的类型改为“COH3D8P”。
优选的，孔压节点合并是在更新后的.inp文件中找出裂缝单元最小的孔压节点编号，记为npp,min，将此.inp文件导入ABAQUS/CAE，使用Display group/Node选项并输入npp,min:10000000，显示所有未合并的孔压节点，而后采用Mesh/Node Merge工具，设置合并容差为10-6，即可将相同坐标的孔压节点全部合并，重新输出.inp文件，并在“*Part”之后“*Assembly”之前添加关键字“*Nset＝coh-pp”，将所有裂缝单元的孔压节点组成一个名称为“coh-pp”的集合，完成孔压节点合并。
本发明与现有技术相比具有以下优点：
1、本发明为碎软低渗油气储层压裂提供了有效的模拟工具，可以应用在所有实体单元之间、以及所有储层结构面上嵌入含孔压节点的裂缝单元，实现压裂裂缝在碎软低渗油气储层中的任意扩展，是对简单裂缝模型的重大改进；创建碎软低渗储层的韧性破坏-渗流耦合本构模型，以反映储层结构面及微裂隙的韧性断裂-渗流、储层基质的塑性损伤-渗流、以及二者之间的应力-渗流耦合等特殊过程；通过将上述过程耦合求解，最终获得碎软低渗储层压裂后的缝网形态，为压裂参数的优化提供模拟方法。
2、本发明针对实体单元采用隐式时间积分求解算法，针对裂缝单元采用基于显式时间积分求解算法，可以大大提高数值解的收敛性，同时加快计算效率。
3、本发明可用于分析不同压裂工程参数、不同工程地质参数影响下压裂缝网的形态，优化压裂参数、实现油气资源高效开采，具有较大的应用前景。
附图说明
图1是本发明中方法流程示意图；
图2(a)是实施例2中的CT实验结果；
图2(b)是实施例2中碎软低渗储层中结构面的三维网络几何模型；
图3是裂缝单元的断裂过程区示意图；
图4是实体单元节点复制、编号方法与裂缝单元示意图；
图5是孔压节点合并示意图；
图6(a)是裂缝单元纯I型和纯II型断裂过程中的牵引力-位移曲线；
图6(b)是实体单元损伤演化曲线；
图7是数值计算模型；
图8(a)是在储层结构面处嵌入裂缝单元的结果；
图8(b)是在其余实体单元边界处嵌入裂缝单元的示意图；
图9是水力裂缝穿透煤岩界面时刻顶板岩层的位移矢量；
图10是碎软低渗煤层中水力裂缝的开度云图；
图11是注液点与煤岩界面不同间距条件下水力裂缝穿透界面时塑性功占水力能量的比例以及实体单元中压裂液漏失速率；
图12是储层压裂过程中流体压力-时间曲线的模拟结果。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
实施例1，如图1所示，本发明提供一种技术方案：一种碎软低渗油气储层的压裂模拟方法，包括以下步骤；
S1、编制碎软低渗储层结构面三维网络模型程序，确定模型中结构面的几何参数，建立碎软低渗油气储层压裂的数值计算模型；
具体是采用CT实验方法确定碎软低渗储层中主要结构面的分布，以及倾角、间距、长度的几何参数，再根据拟建立数值计算模型的大小得出x、y、z三个方向上Voronoi多边形的个数，估算散点在三个方向上的浮动范围，由此在matlab软件中生成一组三维点集，取这组点集中的任意一点，计算距离该点最近的另外一点，以这两点的连线构成Delaunay三角形的一边，然后再找出与这条边距离最近且满足Delaunay三角剖分中的外接圆准则的一点作为第三点，构成一个Delaunay三角形，在此基础上，以Delaunay法的空间外接球准则找出第四点，这四个点构成一个Delaunay四面体；
由上述Delaunay四面体相邻边的垂直平分线所围成的图形即为Voronoi多边形或多面体，通过以上方法，可生成储层结构面的三维网络模型，通过多系修正散点在三个方向上的浮动范围，即可得到符合CT实验结果的3D结构面网络模型
S2、建立储层韧性破坏-渗流本构模型，分别针对裂缝单元、实体单元开发基于隐式和显式时间积分的求解算法；
包括储层基质的塑性损伤及压裂液滤失计算，储层结构面的韧性断裂过程及压裂液流动模拟计算和多物理过程耦合求解计算。
储层基质的塑性损伤及压裂液滤失计算的具体是：
储层基质在压裂过程中会出现拉、剪两种基本破坏形式，拉、剪有效应力分别为：
其中σi和ni(i＝1，2，3)分别表示主应力大小和方向；
总应力矢量
与总应变矢量εs的关系为：
其中，E0为弹性刚度矩阵；
pw是孔隙流体压强；
I为单位矩阵；
α为有效应力系数；且有α＝1-Kb/Ks，Ks是固体颗粒的体积模量，Kb是多孔介质的排水体积模量；
张量d由拉、压损伤变量d(+)与d(-)表示：
d＝d(+)N(+)+d(-)(I-N(+))
其中，假设储层基质在受拉时出现弹性损伤，以
达到抗拉强度作为张拉损伤起始的判据，且当抗拉强度为0时基质完全产生张拉破坏，通过紧凑拉伸试验获取张拉损伤演化规律；
储层基质在发生塑性损伤前后的渗透系数表达式为：
储层结构面的韧性断裂过程及压裂液流动模拟计算具体为：
地层中压裂裂缝会受到储层结构面的影响而出现弹性变形-混合型断裂的过程，在弹性变形阶段，应力和应变之间的关系为：
σc＝D0,cεc
当达到以下应力条件时，结构面的拉、剪强度开始下降；
其中，σc,l和
中l表示n,s,t，分别为法向和两个切向方向，σc,l和
表示牵引力或者牵引力峰值；
D0,c是弹性刚度矩阵；εc是应变矢量，其与分离位移的关系为εc＝S/T0,T0是裂缝单元的本构厚度，在ABAQUS/CAE的Property/Section模块中设置该值为实体单元特征尺寸的0.01倍；
在峰值应力
之后，储层结构面发生混合型韧性断裂过程，其应力σc,l按照如下规律降低：
其中，Sn，Ss和St是变量，表示法向和两个切向方向的位移；Γn,ΓS是断裂能常数，表达式为：
式中Gn与GS是法向与总的切向断裂能，且有GS＝Gs+Gt，Gs和Gt是两个切向方向的断裂能；β与γ表示纯I型和纯II型断裂实验所得应力-应变曲线的形状参数；m和n与β和γ之间的关系为：
其中，χn＝sn,p/sn,χS＝sS,p/sS，sn,p，sS,p以及sn，sS分别是纯I型和纯II型断裂实验所得应力-应变曲线的峰值位移与断裂位移，分别采用紧凑拉伸、贯穿剪切实验确定材料参数，
引入立方定律反映压裂液的切向流动，考虑煤岩基质塑性损伤引起压裂液的滤失或法向流动，建立压裂液质量守恒方程如下：
其中，Qc是总流量；w是裂缝宽度，可由等式(11)～(13)计算得到；μ是压裂液动力粘度；
是切向流的流体压力梯度；ρw是压裂液密度；g是重力加速度；pn,cen和pn,boun是裂缝中心处和边缘处的流体压力。
多物理过程耦合求解具体是：压裂模型可分为2个部分：储层基质、储层结构面的流固耦合模型，
对于储层基质，压裂模拟过程中需要全耦合求解的主变量包括：
xT＝[ui,1,ui,2,…,ui,i；ps,1,ps,2,…,ps,i；Qs,1,Qs,2,…,Qs,i]
其中，ui为每个实体单元节点位移，ps,i为每个实体单元孔隙流体压强，Qs,i为流入每个实体单元的流体流量；
固体部分方程采用有限元法离散，流体部分方程采用有限体积法离散，采用Newton-Raphson迭代方法对离散后的方程进行隐式求解，已知第j迭代步求解第(j+1)步未知变量的表达式为：
Jacobi(xj)dxj+1＝-Rj
其中，Jacobi(xj)为第j迭代步时的雅可比矩阵，Rj为第j迭代步时的残差。
迭代收敛的条件：
||R||2＜tolerance并且||x||2＜tolerance
其中||||2为向量的二型范数，当容差tolerance取值为10-8时，实体单元的一个计算步只需要2～5次迭代即可收敛；
对于裂缝单元，其节点力与实体单元的节点力相等，在此基础上显式求解裂缝单元的节点位移sn,ss,st，得到单元张开位移w，结合质量守恒方程以及N-S方程，计算pn,cen值，通过立方定律得到结构面中切向流流速；本发明假设实体单元中的孔隙流体压强ps,i与裂缝边缘处的流体压强pn,boun相等，结合实体单元渗透系数k的表达式，显示求解结构面中法向流的流速。
S3、数值模型中结构面处以及储层基质内嵌入裂缝单元，合并孔压节点模拟宏细观裂缝的随机分布，具体为在储层结构面处及基质内的单元边界处嵌入零厚度的裂缝单元，合并具有相同节点坐标的裂缝单元上的孔压节点，模拟储层中宏细观裂缝的随机分布；
包括如下内容：确定零厚度裂缝单元的尺寸；零厚度裂缝单元的嵌入和孔压节点合并。
确定零厚度裂缝单元的尺寸的具体是，裂缝单元的尺寸需小于断裂过程区的长度L，表达式为：
其中，μ为材料泊松比。
零厚度裂缝单元的嵌入包括更新实体单元节点编号和组装裂缝单元：在步骤S1的基础上建立数值模型，并对其进行全局网格划分，在Job模块中生成对应的.inp脚本文件，在该文件中，记最大节点编号为nmax，最大单元编号为emax；
将第ni(ni≠1)个节点复制为b个，b为共用单元的个数，在坐标相同的前提下，保持第1个节点号与初始值相同，其后的第k个节点的编号修改为：
其中，Ten(nmax)为nmax的十进制位数；
经过上述变换，第k个节点被复制了b次，编号依次为n1,n2,…,nb；而共用第k个节点的单元也有b个，单元编号依次为e1,e2,…,eb，这样在.inp文件中，通过将新的节点编号复制在对应单元编号位置处，即可使两个有序数组Elem＝(e1,e2,…,eb)和Node＝(n1,n2,…,nb)形成一一对应的关系，由此得到更新节点编号的实体单元，
eb→nb,其中eb∈Elem,nb∈Node
遍历原始.inp文件中相邻实体单元之间编号相同的节点，将其复制，并按照下式确定的顺序组成新的数组：
Group＝(ej,ek,nShared01,nShared02,nShared03,nShared04,Node01,Node02,Node03,Node04)
数组Group中的Group中的ej和ek为相邻单元的单元编号；nShared01～nShared04为原始.inp文件中的相邻单元之间的相同编号的节点；Node01～Node04为重新编号后的节点号数组；
在Node01～Node04中，仅保留由nShared01～nShared04更新的节点编号，然后逆序排列，即可得到单元ej和ek之间的裂缝单元的节点编号；
而后按照下式赋予裂缝单元的单元编号：
ecoh＝1000+emax
在.inp文件中写入所有裂缝单元的编号及其节点编号，在该字符前添加关键字“*Element type＝COH3D8”，同时在关键字“*Part”之后“*Assembly”之前添加关键字“*Elset＝cohElem02”，将所有裂缝单元组成一个名称为“cohElem02”的集合；
然后将裂缝单元嵌入步骤S1建立的储层结构面中，具体为将更新后的裂缝单元编号ej，以及对应的节点坐标nj(xj,yj,zj)复制在Excel中，采用内积法求得该单元的方向向量dip,coh＝(d1,d2,d3)。若dip,coh与步骤S1中储层结构面的方向向量dip共线、任一节点在储层结构面(即泰森多边形)上、而且任一节点都落在泰森多边形范围之内，则判定单元ej在储层结构面内，将所有符合上述三个条件的裂缝单元和节点编号复制在.inp文件中，在关键字“*Part”之后“*Assembly”之前添加关键字“*Elset＝cohElem01”，将所有储层结构面上的单元组成一个名称为“cohElem01”的集合，而后将.inp文件导入ABAQUS/CAE，在Mesh/Assign Element type模块下的将所有裂缝单元的类型改为“COH3D8P”。
孔压节点合并是在更新后的.inp文件中找出裂缝单元最小的孔压节点编号，记为npp,min，将此.inp文件导入ABAQUS/CAE，使用Display group/Node选项并输入npp,min:10000000，显示所有未合并的孔压节点，而后采用Mesh/Node Merge工具，设置合并容差为10-6，即可将相同坐标的孔压节点全部合并，重新输出.inp文件，并在“*Part”之后“*Assembly”之前添加关键字“*Nset＝coh-pp”，将所有裂缝单元的孔压节点组成一个名称为“coh-pp”的集合，完成孔压节点合并。
S4、获取并输入材料参数、设置边界条件和初始条件，根据实验室实验数据和测井资料获取计算参数，输入与压裂模拟相关的材料参数：包括储层基质内部和储层结构面上的裂缝单元在法向和两个切向的弹性模量、峰值应力、纯I型与纯II型断裂能、牵引力-分离曲线的形状参数、断裂位移，压裂液的动力粘度，储层基质的弹性模量、泊松比、剪胀角、三轴与单轴屈服应力之比、损伤前后的渗透系数、拉压强度、拉压损伤演化规律。
边界条件参数包括：而后采用位移边界约束数值计算模型的2个(平面模型)或3个(三维模型)方向，同时限定模型边界的孔隙流体压力，在此基础上通过initial condition命令施加地层的初始地应力、初始孔隙比、初始饱和度，在选定的裂缝单元孔压节点处设置注液点，通过cflow命令施加注液速率，在通过Edit keywords中，通过initial gap命令设置初始损伤区的数量、位置、长度、直径，以模拟射孔引起的煤岩层的初始损伤。
S5.求解模型并输出结果。设定需要输出的参数，包括裂缝张开与剪切位移分布、裂缝传播路径、裂缝内流体压强、裂缝内流体流速、裂缝单元中弹性能与非弹性能、断裂能混合比；完整煤岩体中应力、塑性应变、损伤变量、孔隙流体压强、流体运移速率；以及各节点处孔隙流体压强随时间的演化规律。
S6、分析模拟效果，优化压裂工艺参数，进而提高碎软低渗储层中的油气产量，根据数值模拟结果，对碎软低渗储层的压裂效果进行分析；通过研究不同工程、地质参数下储层内压裂裂缝的长度和面积(二维数值模型)，或面积和体积(三维数值模型)，可以得到各工程、地质参数对碎软低渗储层的压裂效果，进而进行油气产能评价，优化压裂方案。
实施例2，S1、通过CT实验确定碎软低渗储层中天然结构面的分布、倾角、间距、方位角、以及各方位角上结构面的长度，具体如图2所示。
S2.在编制好的碎软低渗储层结构面的三维网络几何模型的生成程序中，输入三维散点在x,y,z三个方向上的平均间距分别为20mm,50mm,25mm，在此基础上设置浮动值分别为8mm，8mm，10mm。然后建立x＝10m，y＝14m，z＝20m的碎软低渗储层压裂数值计算模型，具体如图7所示。
S3.确定零厚度裂缝单元的尺寸为0.1m，采用四面体网格对数值计算模型进行全局网格划分，分别在储层结构面位置处和储层基质内任意两个相邻实体单元之间嵌入裂缝单元，如图8所示，将所有裂缝单元上的具有相同坐标的孔压节点合并，并设置为一个名称为coh-pp的集合。通过减选coh-pp集合得到与煤岩界面间距为1m的注液点。
S4.分别约束数值计算模型6个边界的法向位移。然后在实例中输入材料参数和工程地质参数。工程地质参数如表1所示，其中，初始地应力和初始孔隙流体压强分别施加在实体单元和单元节点上；材料参数如表2和图6所示。
表1工程地质参数
表2材料参数
数值计算开始前，根据实际问题设置初始时间步、最小增量步、最大增量步、迭代收敛的孔压条件、初始损伤区域。在该实例中，设置初始时间步0.1s、最小增量步10-12s、最大增量步10s、迭代收敛的孔压条件1015Pa。在注液点下方0.5m长度、0.0009m直径范围内设置初始损伤的裂缝单元。
S5、输入工程参数、材料参数、分析步设定参数，基于储层韧性破坏-渗流本构模型求解算法，进行压裂数值模拟。流体注入时间为60min。
S6、得到碎软低渗储层压裂后的结果文件、数据及图形，包括：
(A)水力裂缝张开位移矢量与云图，如图9和图10所示。从计算结果可以看出，水力压裂顶板引起其张拉位移有助于其下部煤层发生破坏；由于煤层中天然结构面的存在，压裂裂缝形成复杂的缝网结构。
(B)不同注液点位置D条件下，压裂裂缝穿透煤岩界面时塑性功占水力能量的比例以及压裂液漏失速率，如图11所示。从计算结果可以看出，塑性功占水力能量的比例会随着D的增加而快速增加，压裂液漏失速率会随着D的增加而减小。这意味着将水平井布置在距离煤岩界面1m左右处可保证间接压裂效果。
(C)注液点处流体压力随时间的演化规律，如图12所示。从计算结果可以看出，当压裂裂缝在顶板中，起裂压力在45MPa；当压裂裂缝在煤层中，起裂压力在20MPa左右，并随时间呈波动变化。
与现有压裂模型相比，一方面考虑了天然结构面分布对压裂效果的影响，更加符合碎软低渗储层的特点；另一方面，同时考虑了储层结构面的韧性断裂-渗流、储层基质的塑性损伤-渗流，使得压裂模拟结果更好地反映了储层韧性破坏对压裂效果的影响。为碎软低渗储层的压裂模拟提供了有力的工具。
需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。
尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。
Claims (10)
1.一种碎软低渗油气储层的压裂模拟方法，其特征在于：包括以下步骤；
S1、编制碎软低渗储层结构面三维网络模型程序，确定模型中结构面的几何参数，建立碎软低渗油气储层压裂的数值计算模型；
S2、建立储层韧性破坏-渗流本构模型，分别针对裂缝单元、实体单元开发基于隐式和显式时间积分的求解算法；
S3、数值模型中结构面处以及储层基质内嵌入裂缝单元，合并孔压节点模拟宏细观裂缝的随机分布，具体为在储层结构面处及基质内的单元边界处嵌入零厚度的裂缝单元，合并具有相同节点坐标的裂缝单元上的孔压节点，模拟储层中宏细观裂缝的随机分布；
S4、获取并输入材料参数、设置边界条件和初始条件；
S5、求解模型并输出结果，设定需要输出的参数，包括裂缝张开与剪切位移分布、裂缝传播路径、裂缝内流体压强、裂缝内流体流速、裂缝单元中弹性能与非弹性能、断裂能混合比；完整煤岩体中应力、塑性应变、损伤变量、孔隙流体压强、流体运移速率，以及各节点处孔隙流体压强随时间的演化规律；
S6、分析模拟效果，优化压裂工艺参数，进而提高碎软低渗储层中的油气产量。
2.根据权利要求1所述的一种碎软低渗油气储层的压裂模拟方法，其特征在于，在S1中，是采用CT实验方法确定碎软低渗储层中主要结构面的分布，以及倾角、间距、长度的几何参数，再根据拟建立数值计算模型的大小得出x、y、z三个方向上Voronoi多边形的个数，估算散点在三个方向上的浮动范围，由此在matlab软件中生成一组三维点集，取这组点集中的任意一点，计算距离该点最近的另外一点，以这两点的连线构成Delaunay三角形的一边，然后再找出与这条边距离最近且满足Delaunay三角剖分中的外接圆准则的一点作为第三点，构成一个Delaunay三角形，在此基础上，以Delaunay法的空间外接球准则找出第四点，这四个点构成一个Delaunay四面体；
由上述Delaunay四面体相邻边的垂直平分线所围成的图形即为Voronoi多边形或多面体，通过以上方法，可生成储层结构面的三维网络模型，通过多系修正散点在三个方向上的浮动范围，即可得到符合CT实验结果的3D结构面网络模型。
3.根据权利要求1所述的一种碎软低渗油气储层的压裂模拟方法，其特征在于，在S2中，包括储层基质的塑性损伤及压裂液滤失计算，储层结构面的韧性断裂过程及压裂液流动模拟计算和多物理过程耦合求解计算。
4.根据权利要求3所述的一种碎软低渗油气储层的压裂模拟方法，其特征在于，储层基质的塑性损伤及压裂液滤失计算的具体是：
储层基质在压裂过程中会出现拉、剪两种基本破坏形式，拉、剪有效应力分别为：
其中σi和ni(i＝1，2，3)分别表示主应力大小和方向；
总应力矢量
与总应变矢量εs的关系为：
其中，E0为弹性刚度矩阵；
pw是孔隙流体压强；
I为单位矩阵；
α为有效应力系数；且有α＝1-Kb/Ks，Ks是固体颗粒的体积模量，Kb是多孔介质的排水体积模量；
张量d由拉、压损伤变量d(+)与d(-)表示：
d＝d(+)N(+)+d(-)(I-N(+))
其中，假设储层基质在受拉时出现弹性损伤，以
达到抗拉强度作为张拉损伤起始的判据，且当抗拉强度为0时基质完全产生张拉破坏，通过紧凑拉伸试验获取张拉损伤演化规律；
储层基质在发生塑性损伤前后的渗透系数表达式为：
5.根据权利要求3所述的一种碎软低渗油气储层的压裂模拟方法，其特征在于，储层结构面的韧性断裂过程及压裂液流动模拟计算具体为：
地层中压裂裂缝会受到储层结构面的影响而出现弹性变形-混合型断裂的过程，在弹性变形阶段，应力和应变之间的关系为：
σc＝D0,cεc
当达到以下应力条件时，结构面的拉、剪强度开始下降；
其中，σc,l和
中l表示n,s,t，分别为法向和两个切向方向，σc,l和
表示牵引力或者牵引力峰值；
D0,c是弹性刚度矩阵；εc是应变矢量，其与分离位移的关系为εc＝S/T0,T0是裂缝单元的本构厚度，在ABAQUS/CAE的Property/Section模块中设置该值为实体单元特征尺寸的0.01倍；
在峰值应力
之后，储层结构面发生混合型韧性断裂过程，其应力σc,l按照如下规律降低：
其中，Sn，Ss和St是变量，表示法向和两个切向方向的位移；Γn,ΓS是断裂能常数，表达式为：
式中Gn与GS是法向与总的切向断裂能，且有GS＝Gs+Gt，Gs和Gt是两个切向方向的断裂能；β与γ表示纯I型和纯II型断裂实验所得应力-应变曲线的形状参数；m和n与β和γ之间的关系为：
其中，χn＝sn,p/sn,χS＝sS,p/sS，sn,p，sS,p以及sn，sS分别是纯I型和纯II型断裂实验所得应力-应变曲线的峰值位移与断裂位移，分别采用紧凑拉伸、贯穿剪切实验确定材料参数。
引入立方定律反映压裂液的切向流动，考虑煤岩基质塑性损伤引起压裂液的滤失或法向流动，建立压裂液质量守恒方程如下：
其中，Qc是总流量；w是裂缝宽度；μ是压裂液动力粘度；
是切向流的流体压力梯度；ρw是压裂液密度；g是重力加速度；pn,cen和pn,boun是裂缝中心处和边缘处的流体压力。
6.根据权利要求3所述的一种碎软低渗油气储层的压裂模拟方法，其特征在于，多物理过程耦合求解具体是：压裂模型可分为2个部分：储层基质、储层结构面的流固耦合模型，
对于储层基质，压裂模拟过程中需要全耦合求解的主变量包括：
xT＝[ui,1,ui,2,…,ui,i；ps,1,ps,2,…,ps,i；Qs,1,Qs,2,…,Qs,i]
其中，ui为每个实体单元节点位移，ps,i为每个实体单元孔隙流体压强，Qs,i为流入每个实体单元的流体流量；
固体部分方程采用有限元法离散，流体部分方程采用有限体积法离散，采用Newton-Raphson迭代方法对离散后的方程进行隐式求解，已知第j迭代步求解第(j+1)步未知变量的表达式为：
Jacobi(xj)d xj+1＝-Rj
其中，Jacobi(xj)为第j迭代步时的雅可比矩阵，Rj为第j迭代步时的残差。
迭代收敛的条件：
||R||2＜tolerance并且||x||2＜tolerance
其中|| ||2为向量的二型范数，当容差tolerance取值为10-8时，实体单元的一个计算步只需要2～5次迭代即可收敛；
对于裂缝单元，其节点力与实体单元的节点力相等，在此基础上显式求解裂缝单元的节点位移sn,ss,st，得到单元张开位移w，结合质量守恒方程以及N-S方程，计算pn,cen值，通过立方定律得到结构面中切向流流速；本发明假设实体单元中的孔隙流体压强ps,i与裂缝边缘处的流体压强pn,boun相等，结合实体单元渗透系数k的表达式，显示求解结构面中法向流的流速。
7.根据权利要求1所述的一种碎软低渗油气储层的压裂模拟方法，其特征在于，在步骤S3中，包括如下内容：确定零厚度裂缝单元的尺寸；零厚度裂缝单元的嵌入和孔压节点合并。
8.根据权利要求7所述的一种碎软低渗油气储层的压裂模拟方法，其特征在于，确定零厚度裂缝单元的尺寸的具体是，裂缝单元的尺寸需小于断裂过程区的长度L，表达式为：
其中，μ为材料泊松比。
9.根据权利要求7所述的一种碎软低渗油气储层的压裂模拟方法，其特征在于，零厚度裂缝单元的嵌入包括更新实体单元节点编号和组装裂缝单元：在步骤S1的基础上建立数值模型，并对其进行全局网格划分，在Job模块中生成对应的.inp脚本文件，在该文件中，记最大节点编号为nmax，最大单元编号为emax；
将第ni(ni≠1)个节点复制为b个，b为共用单元的个数，在坐标相同的前提下，保持第1个节点号与初始值相同，其后的第k个节点的编号修改为：
其中，Ten(nmax)为nmax的十进制位数；
经过上述变换，第k个节点被复制了b次，编号依次为n1,n2,…,nb；而共用第k个节点的单元也有b个，单元编号依次为e1,e2,…,eb，这样在.inp文件中，通过将新的节点编号复制在对应单元编号位置处，即可使两个有序数组Elem＝(e1,e2,…,eb)和Node＝(n1,n2,…,nb)形成一一对应的关系，由此得到更新节点编号的实体单元，
eb→nb,其中eb∈Elem,nb∈Node
遍历原始.inp文件中相邻实体单元之间编号相同的节点，将其复制，并按照下式确定的顺序组成新的数组：
Group＝(ej,ek,nShared01,nShared02,nShared03,nShared04,Node01,Node02,Node03,Node04)
数组Group中的ej和ek为相邻单元的单元编号；nShared01～nShared04为原始.inp文件中的相邻单元之间的相同编号的节点；Node01～Node04为重新编号后的节点号数组；
在Node01～Node04中，仅保留由nShared01～nShared04更新的节点编号，然后逆序排列，即可得到单元ej和ek之间的裂缝单元的节点编号；
而后按照下式赋予裂缝单元的单元编号：
ecoh＝1000+emax
在.inp文件中写入所有裂缝单元的编号及其节点编号，在该字符前添加关键字“*Element type＝COH3D8”，同时在关键字“*Part”之后“*Assembly”之前添加关键字“*Elset＝cohElem02”，将所有裂缝单元组成一个名称为“cohElem02”的集合；
然后将裂缝单元嵌入步骤S1建立的储层结构面中，具体为将更新后的裂缝单元编号ej，以及对应的节点坐标nj(xj,yj,zj)复制在Excel中，采用内积法求得该单元的方向向量dip,coh＝(d1,d2,d3)。若dip,coh与步骤S1中储层结构面的方向向量dip共线、任一节点在储层结构面(即泰森多边形)上、而且任一节点都落在泰森多边形范围之内，则判定单元ej在储层结构面内，将所有符合上述三个条件的裂缝单元和节点编号复制在.inp文件中，在关键字“*Part”之后“*Assembly”之前添加关键字“*Elset＝cohElem01”，将所有储层结构面上的单元组成一个名称为“cohElem01”的集合，而后将.inp文件导入ABAQUS/CAE，在Mesh/Assign Element type模块下的将所有裂缝单元的类型改为“COH3D8P”。
10.根据权利要求7所述的一种碎软低渗油气储层的压裂模拟方法，其特征在于，孔压节点合并是在更新后的.inp文件中找出裂缝单元最小的孔压节点编号，记为npp,min，将此.inp文件导入ABAQUS/CAE，使用Display group/Node选项并输入npp,min:10000000，显示所有未合并的孔压节点，而后采用Mesh/Node Merge工具，设置合并容差为10-6，即可将相同坐标的孔压节点全部合并，重新输出.inp文件，并在“*Part”之后“*Assembly”之前添加关键字“*Nset＝coh-pp”，将所有裂缝单元的孔压节点组成一个名称为“coh-pp”的集合，完成孔压节点合并。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110131767.8A CN112883661B (zh) | 2021-01-30 | 2021-01-30 | 一种碎软低渗油气储层的压裂模拟方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110131767.8A CN112883661B (zh) | 2021-01-30 | 2021-01-30 | 一种碎软低渗油气储层的压裂模拟方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112883661A true CN112883661A (zh) | 2021-06-01 |
| CN112883661B CN112883661B (zh) | 2023-07-07 |
Family
ID=76052119
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110131767.8A Active CN112883661B (zh) | 2021-01-30 | 2021-01-30 | 一种碎软低渗油气储层的压裂模拟方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112883661B (zh) |
Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101042775A (zh) * | 2007-04-03 | 2007-09-26 | 浙江大学 | 体三维显示中的体素数据生成方法 |
| CN106593393A (zh) * | 2016-12-09 | 2017-04-26 | 太原理工大学 | 一种提高碎软油气储层渗透率的方法 |
| CN108446415A (zh) * | 2017-12-22 | 2018-08-24 | 北京工业大学 | 一种适用于涂层粒子复合材料考虑几何随机分布的鲁棒设计方法 |
| CN109933939A (zh) * | 2019-03-22 | 2019-06-25 | 西南石油大学 | 非常规双重介质储层多裂缝起裂与扩展的数值模拟方法 |
| CN109992864A (zh) * | 2019-03-22 | 2019-07-09 | 成都理工大学 | 非常规双重介质储层体积压裂数值模拟及参数优化方法 |
| CN110147561A (zh) * | 2018-11-05 | 2019-08-20 | 中国石油大学(华东) | 一种含天然裂缝致密油气储层体积压裂缝网预测方法 |
| CN110175419A (zh) * | 2019-05-30 | 2019-08-27 | 新疆大学 | 风机叶片复合材料细观力学损伤演化分析方法 |
| CN110750930A (zh) * | 2019-10-17 | 2020-02-04 | 西南石油大学 | 一种基于裂缝连续体模型预测裂缝性储层应力演化的方法 |
| CN111988985A (zh) * | 2018-02-20 | 2020-11-24 | 欧司朗有限公司 | 受控农业系统和农业的方法 |
-
2021
- 2021-01-30 CN CN202110131767.8A patent/CN112883661B/zh active Active
Patent Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101042775A (zh) * | 2007-04-03 | 2007-09-26 | 浙江大学 | 体三维显示中的体素数据生成方法 |
| CN106593393A (zh) * | 2016-12-09 | 2017-04-26 | 太原理工大学 | 一种提高碎软油气储层渗透率的方法 |
| CN108446415A (zh) * | 2017-12-22 | 2018-08-24 | 北京工业大学 | 一种适用于涂层粒子复合材料考虑几何随机分布的鲁棒设计方法 |
| CN111988985A (zh) * | 2018-02-20 | 2020-11-24 | 欧司朗有限公司 | 受控农业系统和农业的方法 |
| CN110147561A (zh) * | 2018-11-05 | 2019-08-20 | 中国石油大学(华东) | 一种含天然裂缝致密油气储层体积压裂缝网预测方法 |
| CN109933939A (zh) * | 2019-03-22 | 2019-06-25 | 西南石油大学 | 非常规双重介质储层多裂缝起裂与扩展的数值模拟方法 |
| CN109992864A (zh) * | 2019-03-22 | 2019-07-09 | 成都理工大学 | 非常规双重介质储层体积压裂数值模拟及参数优化方法 |
| CN110175419A (zh) * | 2019-05-30 | 2019-08-27 | 新疆大学 | 风机叶片复合材料细观力学损伤演化分析方法 |
| CN110750930A (zh) * | 2019-10-17 | 2020-02-04 | 西南石油大学 | 一种基于裂缝连续体模型预测裂缝性储层应力演化的方法 |
Non-Patent Citations (6)
| Title |
|---|
| LIU Y: "Visualization and simulation of asphalt concrete with randomly generated three-dimensional models", 《JOURNAL OF COMPUTING IN CIVIL ENGINEERING》 * |
| LIU Y: "Visualization and simulation of asphalt concrete with randomly generated three-dimensional models", 《JOURNAL OF COMPUTING IN CIVIL ENGINEERING》, 31 December 2009 (2009-12-31), pages 340 - 347 * |
| 王丹丹: "新场气田储层裂缝特征及其与动态气水分布的关系", 《石油实验地质》 * |
| 王丹丹: "新场气田储层裂缝特征及其与动态气水分布的关系", 《石油实验地质》, vol. 38, no. 06, 28 November 2016 (2016-11-28), pages 748 - 756 * |
| 邓志鑫: "脱轨列车撞击下盾构隧道的破坏力学特性研究", 《中国优秀硕士学位论文全文数据库工程科技Ⅱ辑》 * |
| 邓志鑫: "脱轨列车撞击下盾构隧道的破坏力学特性研究", 《中国优秀硕士学位论文全文数据库工程科技Ⅱ辑》, no. 09, 15 September 2018 (2018-09-15), pages 034 - 234 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112883661B (zh) | 2023-07-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Li et al. | A numerical investigation of the hydraulic fracturing behaviour of conglomerate in Glutenite formation | |
| Alhuthali et al. | Field applications of waterflood optimization via optimal rate control with smart wells | |
| Hattori et al. | Numerical simulation of fracking in shale rocks: current state and future approaches | |
| Jiang et al. | Effects of fractures on the well production in a coalbed methane reservoir | |
| Zhao et al. | Stability evaluation model for high rock slope based on element extension theory | |
| Ren et al. | CDEM-based simulation of the 3D propagation of hydraulic fractures in heterogeneous Coalbed Methane reservoirs | |
| Shahid et al. | A review of numerical simulation strategies for hydraulic fracturing, natural fracture reactivation and induced microseismicity prediction | |
| Vyazmensky | Numerical modelling of surface subsidence associated with block cave mining using a finite element/discrete element approach | |
| CN110705168A (zh) | 构造应力场的模拟方法 | |
| Xiao et al. | Productivity Prediction and Influencing Factors of Low Permeability Reservoirs after Steering Fracturing Stimulation | |
| CN112200423A (zh) | 一种基于bq、数值模拟的围岩稳定性动态评价方法 | |
| Jang et al. | The oil production performance analysis using discrete fracture network model with simulated annealing inverse method | |
| Yang et al. | Numerical study on the use of alternating injection hydraulic fracturing technology to optimize the interaction between hydraulic fracture and natural fracture | |
| Liu et al. | Mining method selection and optimization for hanging-wall ore-body at Yanqianshan Iron Mine, China | |
| CN112883661A (zh) | 一种碎软低渗油气储层的压裂模拟方法 | |
| CN113378410B (zh) | 一种采动覆岩导水通道演化的模拟方法 | |
| Chen et al. | Numerical modeling of fracture process using a new fracture constitutive model with applications to 2D and 3D engineering cases | |
| Zhang et al. | Recapitulation and prospect of research on flow field in coal mine gob | |
| Yin et al. | Numerical simulation study of enhanced geothermal system modification and heat mining performance under cyclic injection | |
| Chang et al. | Research on the impact of pre‐existing geological fractures on hydraulic fracturing in high in situ stress environments | |
| Zhu et al. | Numerical simulation of complex fractures propagation in thin sand-shale interbeded formation: A case study of shale oil reservoir in Shengli Oilfield | |
| Ma et al. | Numerical simulation of hydraulic fracture extension patterns at the interface of coal-measure composite rock mass with Cohesive Zone Model | |
| Shahsavar et al. | Stochastic analysis of rock slope stability considering cracked rock masses | |
| Chen et al. | Application of stochastic joint network simulation to composite strata of shallow-buried long-span metro tunnels | |
| Mohajerani et al. | An efficient computational model for simulating stress-dependent flow in three-dimensional discrete fracture networks |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |