本发明公开了一种考虑闭合应力的裂缝岩样酸蚀缝宽数值计算方法，主要步骤：S1、构建未加载闭合应力的天然裂缝初始缝宽；S2、建立碳酸盐岩裂缝受压变形模型，并计算闭合应力下的裂缝受压变形后的缝宽；S3、在步骤S2计算的加载闭合应力的缝宽基础上，建立裂缝酸化的流动反应模型，并计算该时间节点处的酸蚀缝宽；S4、在步骤S3计算得到的酸蚀缝宽的基础上，代入步骤S2，再进行步骤S3重新更新应力平衡后的酸蚀缝宽，得到下一个时间节点的酸蚀缝宽，循环迭代步骤S2和S3直至设定的注酸时间t，最终获取注酸结束时酸蚀缝宽分布。本发明方法能定量评价考虑闭合应力的裂缝岩样酸蚀缝宽，为超深井裂缝型碳酸盐岩储层酸化设计和效果评价提供指导。
一种考虑闭合应力的裂缝岩样酸蚀缝宽数值计算方法
技术领域
[0001] 本发明涉及石油工程技术领域，特别涉及一种碳酸盐岩油气藏酸化增产过程中的考虑闭合应力的裂缝岩样酸蚀缝宽数值计算方法。
背景技术
[0002] 近年来，在深层、超深层海相碳酸盐岩相继发现了多个裂缝型大气藏：如川西观雾山组气藏(平均埋深7500m)、下二叠统气藏(平均埋深7400m)、川东元坝长兴组气藏(平均埋
深7000m)等，展现出巨大的勘探开发前景。这类气藏具有埋藏深(>7000m)、天然裂缝发育、
基质渗透率低(<1mD)、破裂压力高(>90MPa)等特点，受制于现有工程技术条件(井口装置承
压受限)，无法实施大型酸压改造，因此酸化成为确保此类储层获得油气增产的关键技术。
[0003] 深层、超深层裂缝型碳酸盐岩酸化过程中，受制于极高的施工压力，施工排量受限，导致酸液不足以完全开启天然裂缝，从而形成酸液在部分闭合天然裂缝中的流动反应
模式。其酸化机理是：由于基质渗透率极低，天然裂缝将主导酸液流动，酸液进入天然裂缝
增加缝内净压力，同时溶解、破坏原有接触支撑点，破坏缝内应力平衡，导致裂缝内的支撑
应力重新分布；此外，在闭合应力作用下，酸液在天然裂缝内发生流动反应，非均匀刻蚀裂
缝面岩石，最终形成具有一定油气流动能力的通道。详细论述见文献：郭建春,苟波,王坤
杰,等.川西下二叠统超深气井网络裂缝酸化优化设计[J].天然气工业,2017,37(6):34‑
41.)，其中的记载。
[0004] 而现有技术中忽略了“酸蚀裂缝受压变形”与“裂缝内酸流动反应”之间的动态耦合与相互竞争共同作用的影响，容易导致预测的酸蚀缝宽失真，影响深层、超深层海相裂缝
型碳酸盐岩储的层酸化设计和效果评价。
发明内容
[0005] 本发明的目的是提供一种考虑闭合应力的裂缝岩样酸蚀缝宽数值计算方法。
[0006] 本发明提供的考虑闭合应力的裂缝岩样酸蚀缝宽数值计算方法，具体步骤如下：
[0007] S1、构建未加载闭合应力的天然裂缝初始缝宽。具体包括下列步骤：
[0008] S11、将采集自筛选目的层井下岩心或露头，沿天然裂缝延伸方向加工成标准API小岩柱或岩板。
[0009] S12、将岩柱或岩板放置于岩心剖缝器上，采用巴西劈裂方式沿着天然裂缝剖缝，获取两个天然裂缝面。
[0010] S13、获取两个裂缝面高程特征数据，再将获取的裂缝面高程数据的最小值作为基底位置，去掉裂缝面高程数据的基底，仅保留粗糙部分的特征，最后利用数据网格化软件将
高程数据进行标准网格化。
[0011] S14、沿着裂缝面长度、宽度、高度方向建立三维直角坐标系，构成裂缝面形貌高程矩阵，接着，将建立的两个裂缝面形貌高程矩阵在三维空间内重新错位组合，建立起天然裂
缝空间，得到每个网格单元内的初始缝宽b。
[0012] S2、考虑碳酸盐岩弹性‑弹塑性变形特征，建立碳酸盐岩裂缝受压变形模型，并计算闭合应力下的裂缝受压变形后的缝宽。包括以下步骤：
[0013] S21、将缝宽分布数据存储为ASCII文件，利用程序读取并存储为缝宽分布矩阵bi,j。
[0014] S22、输入包括裂缝面岩石的杨氏模量E、泊松比v、最大抗压强度σmax在内的岩石力学参数以及法向加载的闭合应力σc。
[0015] S23、假设一个裂缝网格的形变量di,j，并通过裂缝的形变量di,j与每个网格对应位置的受力大小fi,j的关系，计算每处受力情况fi,j，判断各处的变形特征；然后建立裂缝受压
变形平衡方程如下：
[0016]
[0017] 式中，D0—初始时刻相应的半空间距离，Wi,j—微凸体半空间处的形变， —微凸体长度，ΔLi,j—微凸体受力时自身的形变量。
[0018] S24、计算裂缝表面所有微凸体的支撑应力σ0，并判断与加载闭合应力σc是否满足以下平衡条件：
[0019] |σ0‑σc|≤ε1
[0020] 式中，ε1是计算应力差值，MPa。
[0021] S25、如果不满足平衡条件，则重新修正形变量di,j，重复步骤S23和S24；如果满足条件，则将形变量di,j带入缝宽分布计算公式中，计算得到受压变形后的缝宽，公式如下：
[0022] (bstress)i,j＝bi,j‑di,j。
[0023] S3、在步骤S2计算的加载闭合应力的缝宽基础上，建立裂缝酸化的流动反应模型，并计算该时间节点处的酸蚀缝宽。具体步骤如下：
[0024] S31、首先输入酸液黏度μ、岩石密度ρrock、岩石孔隙度 酸液浓度 酸溶解能力β、酸岩反应速率kg，设定注酸排量Qset、流出端压力Pout、注入端压力Pinj_guess；设定注酸时间
t，并将时间离散为步长为Δt＝1s的微小时间步；然后根据步骤S25计算的缝宽(bstress)i,j
以及裂缝酸化的流动物理模型计算各个网格节点位置处对应的压力场P。
[0025] S32、假设一个初始注入压力Pinj_guess，求出在该注入压力下的压力场，随后再通过局部立方定律式计算注入流量Qcal，接着比较计算注入流量Qcal与设定注入流量Qset大小，最
后通过牛顿循环迭代不断修正Pinj，直到满足条件|Qcal‑Qset|≤ε2，ε2是计算注入流量差值，
获得对应速度场分布。若不满足条件|Qcal‑Qset|≤ε2，则进行以下修正具体修正公式如下：
[0026] 如果是第1次修正注入压力，并且计算注入流量Qcal＜0，那么：
[0027]
[0028] 如果是第1次修正注入压力，并且计算注入流量Qcal≥0，那么：
[0029]
[0030] 如果不是第1次修正注入压力，那么采用如下公式修正注入压力：
[0031] Pinj＝C1Pold‑C2Pold_before
[0032] C1＝(Qold‑Qset)/(Qold‑Qcal)
[0033] C2＝(Qcal‑Qset)/(Qold‑Qcal)
[0034] 式中：Qset—模型设定注入流量，ml/min；
[0035] Qcal—模型计算注入流量，ml/min；
[0036] Pinj—模型带入计算的注入压力，MPa；
[0037] Pold—上一次通过修正带入计算的注入压力，MPa；
[0038] Pold_before—Pold之前一次通过修正带入计算的注入压力，MPa；
[0039] Qold—上一次通过修正计算出的注入流量，ml/min；
[0040] ω—权重，无因次；
[0041] ε2—计算注入流量差值，ml/min。
[0042] S33、通过速度场分布求解各个网格位置处的酸液浓度场分布；浓度场计算公式如下：
[0043]
[0044] 式中： —沿缝宽方向上的平均酸液质量浓度，kg/m3；μ—酸液黏度，kg—酸岩反应速率。
[0045] S34、在求解计算得到裂缝内的酸液浓度场分布之后，再结合缝宽分布、酸液速度场分布以及酸岩反应速率等参数，计算裂缝面的溶解并计算该时间节点处的酸蚀缝宽
(bs)i,j；
[0046]
[0047] 式中：β—酸的溶解能力，kg/kg；
[0048] ρrock—碳酸盐岩基质的密度，kg/m3；
[0049] —碳酸盐岩的基质孔隙度，％；
[0050] Δt—时间离散步长，s；
[0051] kg—酸岩反应速率，kg/(m2·s)；3
[0052] —沿缝宽方向上的平均酸液质量浓度，kg/m ；
[0053] Δx—网格的长度，m；
[0054] Δy—网格的宽度，m；
[0055] An为单元控制体四周与酸液流动方向垂直的四个面的面积：
[0056]
[0057] 为单元控制体四周酸液沿缝宽方向上的平均酸液质量浓度：
[0058]
[0059] S4、在步骤S3计算得到的酸蚀缝宽的基础上，将酸蚀缝宽(bs)i,j代入步骤S2，计算更新受压变形后的缝宽
[0060]
[0061] 再基于缝宽 进行步骤S3重新更新应力平衡后的酸蚀缝宽 即得到下一个时间节点的酸蚀缝宽，计算公式如下：
[0062]
[0063] 式中：β—酸的溶解能力，kg/kg；
[0064] ρrock—碳酸盐岩基质的密度，kg/m3；
[0065] —碳酸盐岩的基质孔隙度，％；
[0066] Δt—时间离散步长，s；
[0067] kg—酸岩反应速率，kg/(m2·s)；3
[0068] —沿缝宽方向上的平均酸液质量浓度，kg/m ；
[0069] Δx—网格的长度，m；
[0070] Δy—网格的宽度，m；
[0071] 循环迭代步骤S2和S3直至设定的注酸时间t，最终获取注酸结束时酸蚀缝宽分布。
[0072] 与现有技术相比，本发明的有益之处在于：
[0073] 本发明的方法考虑了酸化过程中天然裂缝受压变形与酸液在天然裂缝中流动反应耦合的作用，采用数值方法计算考虑酸液的流动反应和裂缝面应力下的酸蚀缝宽。该方
法原理可靠，能定量评价考虑闭合应力的裂缝岩样酸蚀缝宽，为超深井裂缝型碳酸盐岩储
层酸化设计和效果评价提供指导。
[0074] 本发明的其它优点、目标和特征将部分通过下面的说明体现，部分还将通过对本发明的研究和实践而为本领域的技术人员所理解。
附图说明
[0075] 图1是本发明中裂缝半空间和微凸体的获取。
[0076] 图2是本发明中裂缝受压变形计算物理模型。
[0077] 图3是本发明中裂缝中酸液的流动模型。
[0078] 图4是本发明中单元控制体四周与酸液流动方向垂直的四个面。
[0079] 图5为考虑闭合应力的裂缝岩样酸蚀缝宽数值计算方法的流程框图。
[0080] 图6是本发明中天然裂缝初始缝宽分布。
[0081] 图7是本发明中注酸结束后酸蚀缝宽的模拟计算结果。
具体实施方式
[0082] 以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。
[0083] 在本发明的考虑闭合应力的裂缝岩样酸蚀缝宽数值计算方法，包括的四个步骤如下：
[0084] 步骤S1，构建未加载闭合应力的天然裂缝初始宽度，包括以下子步骤：
[0085] 1)将采集自筛选目的层井下岩心或露头，沿天然裂缝延伸方向加工成标准API小岩柱(25mm×50mm)或岩板(175mm×35mm×50mm)；
[0086] 2)将岩柱或岩板放置于岩心剖缝器上，采用巴西劈裂方式沿着天然裂缝剖缝，获取两个天然裂缝面；
[0087] 3)采用专利CN201510319382.9方法获取两个裂缝面高程特征数据，再将获取的裂缝面高程数据的最小值作为基底位置，去掉裂缝面高程数据的基底，仅保留粗糙部分的特
征，最后利用数据网格化软件将高程数据进行标准网格化；
[0088] 4)沿着裂缝面长度、宽度、高度方向建立三维直角坐标系，构成裂缝面形貌高程矩阵，接着，将建立的两个裂缝面形貌高程矩阵在三维空间内重新错位组合，建立起天然裂缝
空间，得到每个网格单元内的初始缝宽b。
[0089] 步骤S2、数值计算闭合应力下的天然裂缝受压变形后的缝宽，包括下裂步骤：
[0090] 1)将缝宽分布数据存储为ASCII文件，利用程序读取并存储为缝宽分布矩阵bi,j，输入相应的网格参数，并计算获取初始时刻相应的半空间距离D0以及微凸体长度 (见图
1)：
[0091] D0＝max(bi,j) (1)
[0092]
[0093] 2)输入裂缝面岩石的杨氏模量E、泊松比v、最大抗压强度σmax等岩石力学参数以及法向加载的闭合应力σc。
[0094] 3)假设一个裂缝网格的形变量di,j，并通过裂缝的形变量di,j与每个网格对应位置的受力大小fi,j的关系，计算每处受力情况fi,j，判断各处的变形特征。
[0095] 定义裂缝形变量d为法向应力由零增加至σc后，两个半空间之间缩短的距离Di,j：
[0096] Di,j＝D0‑di,j (3)
[0097] 则裂缝空间受压与变形平衡控制方程可以成如下形式：
[0098]
[0099] 其中，为Wi,j微凸体半空间处的形变：
[0100]
[0101] 其中， 项代表半空间的形变与距离接触微凸体中心径向距离r的函数形式(见图2)：
[0102]
[0103] 利用辛普森(Simpson)公式求解式(5)的积分，计算对应半空间的形变量：
[0104]
[0105] 其中，ΔLi,j为微凸体受力时自身的形变量：当受力fi,j小于微凸体抗压强度σmax时为弹性变形阶段，当受力fi,j大于微凸体抗压强度σmax时在塑性变形阶段。
[0106] 在弹性变形阶段，微凸体的形变量ΔLi,j可以采用胡克定律进行计算：
[0107]
[0108] 式(3)‑(8)中： ——微凸体(i，j)的初始高度，mm；
[0109] E—杨氏模量，MPa；
[0110] v—泊松比，无因次；
[0111] ΔLi,j—微凸体(i，j)自身的形变量，mm；
[0112] a—微凸体的半径，mm；
[0113] fi,j—微凸体(i，j)所受的力，N。
[0114] 在塑性变形阶段，微凸体及对应半空间位置处的最大应力为岩石抗压强度σmax，且随着形变量增加应力不再变化。
[0115] 将式(7)、(8)代入式(4)整理得到法向应力下裂缝受压与变形平衡控制方程：
[0116]
[0117]
[0118] 在模型限定条件下，当裂缝微凸体形变量为di，j时，比较微凸体(i，j)处所受的应力σi,j与岩石抗压强度σmax的大小，确定微凸体处存在三种接触状态：塑性接触(σi,j≥σmax)、
弹性接触(σi,j＜σmax)和未接触(σi,j＝0)。
[0119]
[0120] 4)计算微凸体的支撑应力σ0并判断与加载闭合应力σc是否满足式(13)的平衡条件，
[0121]
[0122] |σ0‑σc|≤ε1 (13)式中：A—裂缝表面的投影面积，m2；
[0123] σc—所施加的法向应力，MPa；
[0124] σ0—裂缝表面所有微凸体的支撑应力，MPa；
[0125] ε1—计算应力差值(推荐取值0.0001)，MPa。
[0126] 如果不满足上述条件，则重新修正形变量di,j，重复步骤(3)、(4)；如果满足条件，则将形变量di,j带入缝宽分布计算公式中，重新计算并更新缝宽。
[0127] (bstress)i,j＝bi,j‑di,j (14)
[0128] 步骤S3、数值计算闭合应力下各个注酸时刻的酸蚀缝宽。包括下裂步骤：
[0129] 1)输入酸液黏度μ、岩石密度ρrock、岩石孔隙度 酸液浓度 酸溶解能力β、酸岩反应速率kg等参数，设定注酸排量Qset、流出端压力Pout、注入端压力Pinj_guess等；设定注酸时
间t，并将时间离散为步长为Δt＝1s的微小时间步；
[0130] 2)通过步骤S2所得的加载应力下的缝宽矩阵(bstress)i,j以及流动物理模型(见图3)计算各个网格节点位置处对应的压力场P：
[0131]
[0132] 再结合局部立方定律，计算裂缝内部各个网格节点位置处的速度；
[0133]
[0134] 式中：vx—网格节点处沿缝长方向酸液流速，m/s；
[0135] vy—网格节点处沿缝高方向酸液流速，m/s。
[0136] 3)假设一个初始注入压力Pinj_guess，利用式(17)求出在该注入压力下的压力场，随后再通过局部立方定律式(18)计算注入流量Qcal，接着比较计算注入流量Qcal与设定注入流
量Qset大小，最后通过牛顿循环迭代不断修正Pinj，直到满足条件|Qcal‑Qset|≤ε2，获得对应
速度场分布。若不满足条件|Qcal‑Qset|≤ε2，则进行以下修正具体修正公式如下：
[0137] 如果是第1次修正注入压力，并且计算注入流量Qcal＜0，那么：
[0138]
[0139] 如果是第1次修正注入压力，并且计算注入流量Qcal≥0，那么：
[0140]
[0141] 如果不是第1次修正注入压力，那么采用如下公式修正注入压力：
[0142] Pinj＝C1Pold‑C2Pold_before (19)
[0143] C1＝(Qold‑Qset)/(Qold‑Qcal) (20)
[0144] C2＝(Qcal‑Qset)/(Qold‑Qcal) (21)
[0145] 式中：Qset—模型设定注入流量，ml/min；
[0146] Qcal—模型计算注入流量，ml/min；
[0147] Pinj—模型带入计算的注入压力，MPa；
[0148] Pold—上一次通过修正带入计算的注入压力，MPa；
[0149] Pold_before—Pold之前一次通过修正带入计算的注入压力，MPa；
[0150] Qold—上一次通过修正计算出的注入流量，ml/min；
[0151] ω—权重，无因次(推荐取值0.5)；
[0152] ε2—计算注入流量差值(推荐取值0.1)，ml/min。
[0153] 4)通过上述计算得到的速度场分布，求解各个网格位置处的酸液浓度场分布，浓度场计算公式如下：
[0154]
[0155] 式中： —沿缝宽方向上的平均酸液质量浓度，kg/m3；
[0156] 5)在求解计算得到裂缝内的酸液浓度场分布之后，再结合缝宽分布、酸液速度场分布以及酸岩反应速率等参数，计算裂缝面的溶解并计算该时间节点处的酸蚀缝宽
(bs)i,j；
[0157]
[0158] 式中：β—酸的溶解能力，kg/kg；
[0159] ρrock—碳酸盐岩基质的密度，kg/m3；
[0160] —碳酸盐岩的基质孔隙度，％；
[0161] Δt—时间离散步长，s；
[0162] kg—酸岩反应速率，kg/(m2·s)；
[0163] —沿缝宽方向上的平均酸液质量浓度，kg/m3；
[0164] Δx—网格的长度，m；
[0165] Δy—网格的宽度，m；
[0166] An为单元控制体四周与酸液流动方向垂直的四个面的面积(见图4)：
[0167]
[0168] 为单元控制体四周酸液沿缝宽方向上的平均酸液质量浓度：
[0169]
[0170] 步骤S4、在步骤S3计算得到的酸蚀缝宽的基础上，将酸蚀缝宽(bs)i,j代入步骤S2，计算更新受压变形后的缝宽
[0171]
[0172] 再基于缝宽 进行步骤S3重新更新应力平衡后的酸蚀缝宽 即得到下一个时间节点的酸蚀缝宽，计算公式如下：
[0173]
[0174] 式中：β—酸的溶解能力，kg/kg；
[0175] ρrock—碳酸盐岩基质的密度，kg/m3；
[0176] —碳酸盐岩的基质孔隙度，％；
[0177] Δt—时间离散步长，s；
[0178] kg—酸岩反应速率，kg/(m2·s)；
[0179] —沿缝宽方向上的平均酸液质量浓度，kg/m3；
[0180] Δx—网格的长度，m；
[0181] Δy—网格的宽度，m；
[0182] 循环迭代步骤S2和S3直至设定的注酸时间t，最终获取注酸结束时酸蚀缝宽分布。图5是上述的考虑闭合应力的裂缝岩样酸蚀缝宽数值计算方法的流程框图。
[0183] 实施案例：
[0184] 1.基于目标工区某探井X井酸化设计基础参数(表1)，采用步骤S1建立X井初始天然裂缝缝宽分布(见图6)，采用步骤S2建立的裂缝受压变形物理模型，得到加载闭合应力下
的裂缝缝宽。
[0185] 2.以加载闭合应力的裂缝缝宽为初始条件，采用步骤S3中的酸液裂缝内流动反应模型计算裂缝流动场、浓度场，进一步计算得到酸蚀缝宽分布。
[0186] 3.基于上述计算得到的酸蚀缝宽分布，依据步骤S2和S3重新进行应力平衡修正，得到下一个时间节点的缝宽分布；循环迭代直至设定的注酸时间t，最终获取注酸结束时酸
蚀缝宽分布(见图7)，即为考虑闭合应力的裂缝岩样酸蚀缝宽。
[0187] 表1X井酸压设计部分基础参数
[0188]
[0189] 以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人
员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容作出些许更动或修饰
为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质
对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。
