一种循缝找洞酸压数值模拟方法及模拟系统与流程

0 0]
t
；pm为岩石基质孔隙流体压力；
18.基质岩体中流体渗流连续性方程为：
[0019][0020]
式中，αm为biot系数；φm岩石基质孔隙度；ks固体岩石体积模量；k
l
流体体积模量；pm为岩石基质孔隙流体压力；t时间；im为单位张量，二维情况下为[1 1 0]
t
，三维情况下为[1 1 1 0 0 0]
t
；ε应变张量；ρ
l
为流体密度；为散度算子；km基质渗透率张量；μ流体粘度；g为重力加速度向量；q汇源项；
[0021]
在裂缝壁面上建立局部坐标系，用x表示沿裂缝长度方向，用y表示沿裂缝高度方向，裂缝宽度较小，忽略沿垂直于裂缝壁面方向上的流动，裂缝内流体流动的连续性方程表示为：
[0022][0023]
式中，v
x
、vy表示酸液流速；w为裂缝宽度；t时间；
[0024][0025]
式中，v
x
、vy裂缝壁面局部坐标x、y的方向流体速度；w裂缝宽度；μ流体粘度；pf裂缝内流体压力；
[0026]
将(4)带入(3)，得到：
[0027][0028]
酸蚀裂缝内酸液浓度方程为：
[0029][0030]
裂缝壁面上局部反应方程可用溶蚀速度rr(cs)表示为：
[0031]
kg(c
a-cs)＝rr(cs)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0032]
其中，
[0033][0034]
式中，ca裂缝内酸液浓度；cs裂缝壁面酸液浓度；v
x
、vy表示酸液流速；w 裂缝宽度；t时间；kg传质系数；μ流体粘度；r气体常数；t温度；对于不同的岩石矿物，x1、x2、x3取不同的值；pf裂缝内流体压力；
[0035]
酸蚀裂缝宽度变化方程为：
[0036]
[0037]
式中，βi为酸液对白云岩/石灰岩的溶解能力，下标i分别表示白云岩/石灰岩；ρi表示白云岩/石灰岩的密度；φ表示孔隙度；rr(cs)表示酸液对白云岩/石灰岩的溶蚀速度；w为裂缝宽度；t时间；
[0038]
初始时刻位移、孔隙流体压力通过初始条件给出：
[0039][0040]
式中，u为位移，u0为初始时刻位移，pm为孔隙流体压力，为初始时刻孔隙流体压力；
[0041]
边界上的位移、流体压力、应力、流量、酸液浓度通过强制边界条件给出：
[0042][0043]
式中，u为位移，为位移边界；pm为孔隙流体压力，为流体压力边界；σ为应力，n为裂缝壁面向量，为应力边界上施加的已知力，t为时间；km为基质渗透率张量，μ为流体粘度，为流体流量边界上施加的已知流量，ρ
l
为流体流量边界上已知的流体密度；c为酸液浓度。γd为位移边界，γ
p
为流体压力边界，γ
t
为应力边界，γq为流体流量边界；
[0044]
裂缝表面受到裂缝内的压力为：
[0045][0046]
式中，σ为应力，为裂缝两个壁面垂向单位向量，pf为裂缝内流体压力，为裂缝两个壁面；
[0047]
裂缝边界处的孔隙流体压力连续，在裂缝边界与基质之间存在流体交换：
[0048][0049]
式中，pm为岩石基质孔隙流体压力；pf为为裂缝内流体压力；km基质渗透率张量；μ流体粘度；为裂缝两个壁面垂向单位向量；q为流量；ρ为流体密度；为裂缝两个壁面；
[0050]
采用最大应变能释放率准则，对裂纹扩展进行判断；
[0051]
采用加权余量法将控制方程、边界条件结合，得到控制方程的等效积分弱形式：
[0052]
基质岩体中流体渗流连续性方程(2)的等效积分“弱”形式为：
[0053][0054]
式中，w为裂缝宽度；αm为biot系数；φm岩石基质孔隙度；ks固体岩石体积模量；kw流体体积模量；pm为岩石基质孔隙流体压力；t时间；im为单位张量，二维情况下为[1 1 0]
t
，三维情况下为[1 1 1 0 0 0]
t
；ε应变张量，km基质渗透率张量；μ流体粘度；q汇源项；ρ为流体密度；为流体流量边界上施加的已知流量；
[0055]
裂缝内流体流动连续性方程(5)的等效积分“弱”形式表示为：
[0056][0057]
式中，w裂缝宽度；t时间；μ流体粘度；pf为裂缝内孔隙流体压力；ρ为流体密度；为流体流量边界上施加的已知流量；
[0058]
将基质岩体中流体渗流连续性方程的等效积分弱形式(14)和裂缝内流体流动连续性方程的等效积分“弱”形式(15)相加，以消除基质和裂缝之间的流体交换项，得到：
[0059][0060]
采用扩展有限元法对控制方程进行离散，得到扩展有限元离散方程的步骤如下：
[0061]
采用扩展有限单元法近似逼近，其中位移可表示为：
[0062][0063]
式中，u为位移；ni、nj、nk、nh插值函数；ua常规位移节点自由度，ub裂纹面贯穿单元节点的额外位移自由度，uc裂纹前缘单元节点的额外位移自由度，ud为溶洞加强节点的额外位移自由度；
[0064]
h(x)为裂纹面函数：
[0065][0066]
线增函数r
α
(x)定义为：
[0067][0068]
溶洞加强函数r(x)定义为：
[0069][0070]
位移表达式(17)中的函数ψ(x)定义为：
[0071][0072]
将扩展有限元位移场逼近模式带入虚功方程，离散得到应力平衡扩展有限元控制方程：
[0073]
ku＝f
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(22)
[0074]
其中，k为整体劲度矩阵；f为等效载荷列阵；u为位移；
[0075]
同理可得到孔隙压力的加强模式：
[0076][0077]
式中，p为孔隙压力；ni、nj、nk插值函数；pa常规压力节点自由度，pb裂纹面贯穿单元节点的额外压力自由度，pc裂纹前缘单元节点的额外压力自由度； h(x)为裂纹面函数；r
α
(x)为线增函数；
[0078]
将扩展有限元压力场模式带入基质-裂缝渗流的等效积分“弱”形式，则流体渗流的扩展有限元控制方程为：
[0079][0080]
其中，hm和hf分别为基质流体整体渗流矩阵和裂缝流体整体渗流矩阵；dm和df分别为基质整体压缩性矩阵和裂缝整体压缩性矩阵；p为压力矩阵；f为等效载荷列阵；
[0081]
水力裂缝宽度w的表达式为：
[0082]
w＝tncd
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(25)
[0083]
其中，w为水力裂缝宽度；t为局部坐标和整体坐标的转换矩阵；nc为节点的形函数；d为单元节点参数矩阵；
[0084]
建立含天然裂缝、溶洞的物理模型：
[0085]
孔洞内弹性模量、泊松比等材料参数设置为零，同时设置位移逼近为零；
[0086]
基于扩展有限元离散方程，对所述物理模型进行计算。
[0087]
本发明提供的循缝找洞酸压数值模拟方法，具有以下有益效果：
[0088]
循缝找洞酸压数值模拟方法，建立了一套包括应力-流动-酸岩反应三个模块的裂缝孔洞型储层酸蚀裂缝扩展模型，采用扩展有限元方法对控制方程进行了离散，利用最大应变能释放率准则对裂纹扩展进行判断，研究酸蚀裂缝扩展路径。
[0089]
本发明还提供了循缝找洞酸压数值模拟系统，包括应力-流动-酸岩反应耦合数学模型建立模块、离散方程获取模块、物理模型建立模块和计算模块，所述应力-流动-酸岩反应耦合数学模型建立模块、离散方程获取模块、物理模型建立模块和计算模块依次相连。
[0090]
作为优选，所述应力-流动-酸岩反应耦合数学模型建立模块执行如下操作：建立应力-流动-酸岩反应耦合控制方程，设置初始条件、边界条件；采用最大应变能释放率准则，对裂纹扩展进行判断；采用加权余量法将控制方程、边界条件结合，得到控制方程的等效积分弱形式，如下：
[0091]
建立应力-流动-酸岩反应耦合控制方程：
[0092]
岩石骨架受力方程，对其采用应力平衡方程进行描述：
[0093][0094]
式中，为散度算子；σ为柯西应力张量；ρ
l
为流体密度；g为重力加速度向量；σ'为有效应力；；αm为biot系数；im为单位张量，二维情况下为[1 1 0]
t
，三维情况下为[1 1 1 0 0 0]
t
；pm为岩石基质孔隙流体压力；
[0095]
基质岩体中流体渗流连续性方程为：
[0096][0097]
式中，αm为biot系数；φm岩石基质孔隙度；ks固体岩石体积模量；k
l
流体体积模量；pm为岩石基质孔隙流体压力；t时间；im为单位张量，二维情况下为[1 1 0]
t
，三维情况下为[1 1 1 0 0 0]
t
；ε应变张量；ρ
l
为流体密度；为散度算子；km基质渗透率张量；μ流体粘度；g为重力加速度向量；q汇源项；
[0098]
在裂缝壁面上建立局部坐标系，用x表示沿裂缝长度方向，用y表示沿裂缝高度方向，裂缝宽度较小，忽略沿垂直于裂缝壁面方向上的流动，裂缝内流体流动的连续性方程表示为：
[0099][0100]
式中，v
x
、vy表示酸液流速；w为裂缝宽度；t时间；
[0101][0102]
式中，v
x
、vy裂缝壁面局部坐标x、y的方向流体速度；w裂缝宽度；μ流体粘度；pf裂缝内流体压力；
[0103]
将(4)带入(3)，得到：
[0104][0105]
酸蚀裂缝内酸液浓度方程为：
[0106][0107]
裂缝壁面上局部反应方程可用溶蚀速度rr(cs)表示为：
[0108]
kg(c
a-cs)＝rr(cs)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0109]
其中，
[0110][0111]
式中，ca裂缝内酸液浓度；cs裂缝壁面酸液浓度；v
x
、vy表示酸液流速；w 裂缝宽度；t时间；kg传质系数；μ流体粘度；r气体常数；t温度；对于不同的岩石矿物，x1、x2、x3取不同的值；pf裂缝内流体压力；
[0112]
酸蚀裂缝宽度变化方程为：
[0113][0114]
式中，βi为酸液对白云岩/石灰岩的溶解能力，下标i分别表示白云岩/石灰岩；ρi表示白云岩/石灰岩的密度；φ表示孔隙度；rr(cs)表示酸液对白云岩/石灰岩的溶蚀速度；w为裂缝宽度；t时间；
[0115]
初始时刻位移、孔隙流体压力通过初始条件给出：
[0116][0117]
式中，u为位移，u0为初始时刻位移，pm为孔隙流体压力，为初始时刻孔隙流体压力；
[0118]
边界上的位移、流体压力、应力、流量、酸液浓度通过强制边界条件给出：
[0119][0120]
式中，u为位移，为位移边界；pm为孔隙流体压力，为流体压力边界；σ为应力，n为裂缝壁面向量，为应力边界上施加的已知力，t为时间；km为基质渗透率张量，μ为流体粘度，为流体流量边界上施加的已知流量，ρ
l
为流体流量边界上已知的流体密度；c为酸液浓度。γd为位移边界，γ
p
为流体压力边界，γ
t
为应力边界，γq为流体流量边界；
[0121]
裂缝表面受到裂缝内的压力为：
[0122][0123]
式中，σ为应力，为裂缝两个壁面垂向单位向量，pf为裂缝内流体压力，为裂缝两个壁面；
[0124]
裂缝边界处的孔隙流体压力连续，在裂缝边界与基质之间存在流体交换：
[0125][0126]
式中，pm为岩石基质孔隙流体压力；pf为为裂缝内流体压力；km基质渗透率张量；μ流体粘度；为裂缝两个壁面垂向单位向量；q为流量；ρ为流体密度；为裂缝两个壁面；
[0127]
采用最大应变能释放率准则，对裂纹扩展进行判断；
[0128]
采用加权余量法将控制方程、边界条件结合，得到控制方程的等效积分弱形式：
[0129]
基质岩体中流体渗流连续性方程(2)的等效积分“弱”形式为：
[0130][0131]
式中，w为裂缝宽度；αm为biot系数；φm岩石基质孔隙度；ks固体岩石体积模量；kw流体体积模量；pm为岩石基质孔隙流体压力；t时间；im为单位张量，二维情况下为[1 1 0]
t
，三维情况下为[1 1 1 0 0 0]
t
；ε应变张量，km基质渗透率张量；μ流体粘度；q汇源项；ρ为流体密度；为流体流量边界上施加的已知流量；
[0132]
裂缝内流体流动连续性方程(5)的等效积分“弱”形式表示为：
[0133][0134]
式中，w裂缝宽度；t时间；μ流体粘度；pf为裂缝内孔隙流体压力；ρ为流体密度；为流体流量边界上施加的已知流量；
[0135]
将基质岩体中流体渗流连续性方程的等效积分弱形式(14)和裂缝内流体流动连续性方程的等效积分“弱”形式(15)相加，以消除基质和裂缝之间的流体交换项，得到：
[0136][0137]
作为优选，所述离散方程获取模块执行如下操作：采用扩展有限元法对控制方程进行离散，得到扩展有限元离散方程的步骤如下：
[0138]
采用扩展有限单元法近似逼近，其中位移可表示为：
[0139][0140]
式中，u为位移；ni、nj、nk、nh插值函数；ua常规位移节点自由度，ub裂纹面贯穿单元节点的额外位移自由度，uc裂纹前缘单元节点的额外位移自由度， ud为溶洞加强节点的额外位移自由度；
[0141]
h(x)为裂纹面函数：
[0142][0143]
线增函数r
α
(x)定义为：
[0144][0145]
溶洞加强函数r(x)定义为：
[0146][0147]
位移表达式(17)中的函数ψ(x)定义为：
[0148][0149]
将扩展有限元位移场逼近模式带入虚功方程，离散得到应力平衡扩展有限元控制方程：
[0150]
ku＝f
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(22)
[0151]
其中，k为整体劲度矩阵；f为等效载荷列阵；u为位移；
[0152]
同理可得到孔隙压力的加强模式：
[0153][0154]
式中，p为孔隙压力；ni、nj、nk插值函数；pa常规压力节点自由度，pb裂纹面贯穿单元节点的额外压力自由度，pc裂纹前缘单元节点的额外压力自由度； h(x)为裂纹面函数；r
α
(x)为线增函数；
[0155]
将扩展有限元压力场模式带入基质-裂缝渗流的等效积分“弱”形式，则流体渗流的扩展有限元控制方程为：
[0156][0157]
其中，hm和hf分别为基质流体整体渗流矩阵和裂缝流体整体渗流矩阵；dm和df分别为基质整体压缩性矩阵和裂缝整体压缩性矩阵；p为压力矩阵；f为等效载荷列阵；
[0158]
水力裂缝宽度w的表达式为：
[0159]
w＝tncd
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(25)
[0160]
其中，w为水力裂缝宽度；t为局部坐标和整体坐标的转换矩阵；nc为节点的形函数；d为单元节点参数矩阵。；
[0161]
作为优选，所述物理模型建立模块执行如下操作：建立含天然裂缝、溶洞的物理模型，如下：孔洞内弹性模量、泊松比等材料参数设置为零，同时设置位移逼近为零。
[0162]
作为优选，所述计算模块执行如下操作：基于扩展有限元离散方程，对所述物理模型进行计算。
[0163]
该模拟系统采用上述提到的模拟方法，所以具有同样的技术效果。
附图说明
[0164]
图1是本发明提供的一种循缝找洞酸压数值模拟方法的流程图；
[0165]
图2是本发明实施例建立的物理模型；
[0166]
图3是本发明实施例裂缝扩展延伸计算流程图；
[0167]
图4是本发明人工裂缝与溶洞相遇后扩展规律俯视图；
[0168]
图5为本发明提供的一种循缝找洞酸压数值模拟系统的结构示意图。
[0169]
图1-5中附图标记如下：
[0170]
1天然裂缝，2溶洞，3井筒，4人工裂缝，5循缝找洞酸压数值模拟系统， 6应力-流动-酸岩反应耦合数学模型建立模块，7离散方程获取模块，8物理模型建立模块，9计算模块。
具体实施方式
[0171]
如图1-5所示，图1是本发明提供的一种循缝找洞酸压数值模拟方法的流程图；图2是本发明实施例建立的物理模型；图3是本发明实施例裂缝扩展延伸计算流程图；图4是本发明人工裂缝与溶洞相遇后扩展规律俯视图；图5为本发明提供的一种循缝找洞酸压数值模拟系统的结构示意图。
[0172]
结合图1，一种循缝找洞酸压数值模拟方法，包括如下步骤：
[0173]
建立应力-流动-酸岩反应耦合数学模型：建立应力-流动-酸岩反应耦合控制方程，设置初始条件、边界条件；采用最大应变能释放率准则，对裂纹扩展进行判断；采用加权余量法将控制方程、边界条件结合，得到控制方程的等效积分弱形式；
[0174]
采用扩展有限元法对控制方程进行离散，得到扩展有限元离散方程；
[0175]
建立含天然裂缝1、溶洞2的物理模型；
[0176]
基于扩展有限元离散方程，对物理模型进行计算。
[0177]
下面结合实施例对本发明做进一步详细的说明：
[0178]
建立应力-流动-酸岩反应耦合数学模型的步骤如下：
[0179]
建立应力-流动-酸岩反应耦合控制方程：
[0180]
岩石骨架受力方程，对其采用应力平衡方程进行描述：
[0181][0182]
式中，为散度算子；σ为柯西应力张量；ρ
l
为流体密度；g为重力加速度向量；σ'为有效应力；；αm为biot系数；im为单位张量，二维情况下为[1 1 0]
t
，三维情况下为[1 1 1 0 0 0]
t
；pm为岩石基质孔隙流体压力。
[0183]
基质岩体中流体渗流连续性方程为：
[0184][0185]
式中，αm为biot系数；φm岩石基质孔隙度；ks固体岩石体积模量；k
l
流体体积模量；pm为岩石基质孔隙流体压力；t时间；im为单位张量，二维情况下为[1 1 0]
t
，三维情况下为[1 1 1 0 0 0]
t
；ε应变张量；ρ
l
为流体密度；为散度算子；km基质渗透率张量；μ流体粘度；g为重力加速度向量；q汇源项。
[0186]
在裂缝壁面上建立局部坐标系，用x表示沿裂缝长度方向，用y表示沿裂缝高度方向，裂缝宽度较小，忽略沿垂直于裂缝壁面方向上的流动，裂缝内流体流动的连续性方程表示为：
[0187]
[0188]
式中，v
x
、vy表示酸液流速；w为裂缝宽度；t时间。
[0189][0190]
式中，v
x
、vy裂缝壁面局部坐标x、y的方向流体速度；w裂缝宽度；μ流体粘度；pf裂缝内流体压力。
[0191]
将(4)带入(3)，得到：
[0192][0193]
酸蚀裂缝内酸液浓度方程为：
[0194][0195]
裂缝壁面上局部反应方程可用溶蚀速度rr(cs)表示为：
[0196]
kg(c
a-cs)＝rr(cs)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0197]
其中，
[0198][0199]
式中，ca裂缝内酸液浓度；cs裂缝壁面酸液浓度；v
x
、vy表示酸液流速；w 裂缝宽度；t时间；kg传质系数；μ流体粘度；r气体常数；t温度；对于不同的岩石矿物，x1、x2、x3取不同的值；pf裂缝内流体压力。
[0200]
酸蚀裂缝宽度变化方程为：
[0201][0202]
式中，βi为酸液对白云岩/石灰岩的溶解能力，下标i分别表示白云岩/石灰岩；ρi表示白云岩/石灰岩的密度；φ表示孔隙度；rr(cs)表示酸液对白云岩/石灰岩的溶蚀速度；w为裂缝宽度；t时间。
[0203]
初始时刻位移、孔隙流体压力通过初始条件给出：
[0204][0205]
式中，u为位移，u0为初始时刻位移，pm为孔隙流体压力，为初始时刻孔隙流体压力。
[0206]
边界上的位移、流体压力、应力、流量、酸液浓度通过强制边界条件给出：
[0207][0208]
式中，u为位移，为位移边界；pm为孔隙流体压力，为流体压力边界；σ为应力，n为裂缝壁面向量，为应力边界上施加的已知力，t为时间；km为基质渗透率张量，μ为流体粘度，为流体流量边界上施加的已知流量，ρ
l
为流体流量边界上已知的流体密度；c为酸液浓度。γd为位移边界，γ
p
为流体压力边界，γ
t
为应力边界，γq为流体流量边界。
[0209]
裂缝表面受到裂缝内的压力为：
[0210][0211]
式中，σ为应力，为裂缝两个壁面垂向单位向量，pf为裂缝内流体压力，为裂缝两个壁面。
[0212]
裂缝边界处的孔隙流体压力连续，在裂缝边界与基质之间存在流体交换：
[0213][0214]
式中，pm为岩石基质孔隙流体压力；pf为为裂缝内流体压力；km基质渗透率张量；μ流体粘度；为裂缝两个壁面垂向单位向量；q为流量；ρ为流体密度；为裂缝两个壁面。
[0215]
采用最大应变能释放率准则，对裂纹扩展进行判断；
[0216]
采用加权余量法将控制方程、边界条件结合，得到控制方程的等效积分弱形式：
[0217]
基质岩体中流体渗流连续性方程(2)的等效积分“弱”形式为：
[0218][0219]
式中，w为裂缝宽度；αm为biot系数；φm岩石基质孔隙度；ks固体岩石体积模量；kw流体体积模量；pm为岩石基质孔隙流体压力；t时间；im为单位张量，二维情况下为[1 1 0]
t
，三维情况下为[1 1 1 0 0 0]
t
；ε应变张量，km基质渗透率张量；μ流体粘度；q汇源项；ρ为流体密度；为流体流量边界上施加的已知流量。
[0220]
裂缝内流体流动连续性方程(5)的等效积分“弱”形式表示为：
[0221][0222]
式中，w裂缝宽度；t时间；μ流体粘度；pf为裂缝内孔隙流体压力；ρ为流体密度；为流体流量边界上施加的已知流量。
[0223]
将基质岩体中流体渗流连续性方程的等效积分弱形式(14)和裂缝内流体流动连续性方程的等效积分“弱”形式(15)相加，以消除基质和裂缝之间的流体交换项，得到：
[0224][0225]
采用扩展有限元法对控制方程进行离散，得到扩展有限元离散方程的步骤如下：
[0226]
采用扩展有限单元法近似逼近，其中位移可表示为：
[0227][0228]
式中，u为位移；ni、nj、nk、nh插值函数；ua常规位移节点自由度，ub裂纹面贯穿单元节点的额外位移自由度，uc裂纹前缘单元节点的额外位移自由度， ud为溶洞加强节点的额外位移自由度；
[0229]
h(x)为裂纹面函数：
[0230][0231]
线增函数r
α
(x)定义为：
[0232][0233]
溶洞加强函数r(x)定义为：
[0234]
[0235]
位移表达式(17)中的函数ψ(x)定义为：
[0236][0237]
将扩展有限元位移场逼近模式带入虚功方程，离散得到应力平衡扩展有限元控制方程：
[0238]
ku＝f
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(22)
[0239]
其中，k为整体劲度矩阵；f为等效载荷列阵；u为位移。
[0240]
同理可得到孔隙压力的加强模式：
[0241][0242]
式中，p为孔隙压力；ni、nj、nk插值函数；pa常规压力节点自由度，pb裂纹面贯穿单元节点的额外压力自由度，pc裂纹前缘单元节点的额外压力自由度； h(x)为裂纹面函数；r
α
(x)为线增函数。
[0243]
将扩展有限元压力场模式带入基质-裂缝渗流的等效积分“弱”形式，则流体渗流的扩展有限元控制方程为：
[0244][0245]
其中，hm和hf分别为基质流体整体渗流矩阵和裂缝流体整体渗流矩阵；dm和df分别为基质整体压缩性矩阵和裂缝整体压缩性矩阵；p为压力矩阵；f为等效载荷列阵。
[0246]
水力裂缝宽度w的表达式为：
[0247]
w＝tncd
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(25)
[0248]
其中，w为水力裂缝宽度；t为局部坐标和整体坐标的转换矩阵；nc为节点的形函数；d为单元节点参数矩阵。
[0249]
建立含天然裂缝1、溶洞2的物理模型：孔洞内弹性模量、泊松比等材料参数设置为零，同时设置位移逼近为零，建立的物理模型如图2所示，图中共7 条天然裂缝，2个溶洞，溶洞半径为2，溶洞中心分别设置为(36,18,20)、(5,7,20)。
[0250]
基于扩展有限元离散方程，对所述物理模型进行计算：裂缝网络扩展延伸计算流程如图3所示，从计算结果来看，如图4所示，溶洞2对人工裂缝4具有一定的吸引作用，两个溶洞2都在井筒3附近的两条天然裂缝1延长线上，天然裂缝1沿最大主应力方向扩展，沟通了溶洞2。
[0251]
结合图5，本发明还提供了循缝找洞酸压数值模拟系统5，包括应力-流动
‑ꢀ
酸岩反应耦合数学模型建立模块6、离散方程获取模块7、物理模型建立模块8 和计算模块9；
[0252]
所述应力-流动-酸岩反应耦合数学模型建立模6执行如下操作：建立应力
‑ꢀ
流动-酸岩反应耦合控制方程，设置初始条件、边界条件；采用最大应变能释放率准则，对裂纹扩展进行判断；采用加权余量法将控制方程、边界条件结合，得到控制方程的等效积分弱形式，如下：
[0253]
建立应力-流动-酸岩反应耦合控制方程：
[0254]
岩石骨架受力方程，对其采用应力平衡方程进行描述：
[0255][0256]
式中，为散度算子；σ为柯西应力张量；ρ
l
为流体密度；g为重力加速度向量；σ'为有效应力；；αm为biot系数；im为单位张量，二维情况下为[1 1 0]
t
，三维情况下为[1 1 1 0 0 0]
t
；pm为岩石基质孔隙流体压力。
[0257]
基质岩体中流体渗流连续性方程为：
[0258][0259]
式中，αm为biot系数；φm岩石基质孔隙度；ks固体岩石体积模量；k
l
流体体积模量；pm为岩石基质孔隙流体压力；t时间；im为单位张量，二维情况下为[1 1 0]
t
，三维情况下为[1 1 1 0 0 0]
t
；ε应变张量；ρ
l
为流体密度；为散度算子；km基质渗透率张量；μ流体粘度；g为重力加速度向量；q汇源项。
[0260]
在裂缝壁面上建立局部坐标系，用x表示沿裂缝长度方向，用y表示沿裂缝高度方向，裂缝宽度较小，忽略沿垂直于裂缝壁面方向上的流动，裂缝内流体流动的连续性方程表示为：
[0261][0262]
式中，v
x
、vy表示酸液流速；w为裂缝宽度；t时间。
[0263][0264]
式中，v
x
、vy裂缝壁面局部坐标x、y的方向流体速度；w裂缝宽度；μ流体粘度；pf裂缝内流体压力。
[0265]
将(4)带入(3)，得到：
[0266][0267]
酸蚀裂缝内酸液浓度方程为：
[0268][0269]
裂缝壁面上局部反应方程可用溶蚀速度rr(cs)表示为：
[0270]
kg(c
a-cs)＝rr(cs)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0271]
其中，
[0272][0273]
式中，ca裂缝内酸液浓度；cs裂缝壁面酸液浓度；v
x
、vy表示酸液流速；w 裂缝宽度；t时间；kg传质系数；μ流体粘度；r气体常数；t温度；对于不同的岩石矿物，x1、x2、x3取不同的
值；pf裂缝内流体压力。
[0274]
酸蚀裂缝宽度变化方程为：
[0275][0276]
式中，βi为酸液对白云岩/石灰岩的溶解能力，下标i分别表示白云岩/石灰岩；ρi表示白云岩/石灰岩的密度；φ表示孔隙度；rr(cs)表示酸液对白云岩/石灰岩的溶蚀速度；w为裂缝宽度；t时间。
[0277]
初始时刻位移、孔隙流体压力通过初始条件给出：
[0278][0279]
式中，u为位移，u0为初始时刻位移，pm为孔隙流体压力，为初始时刻孔隙流体压力。
[0280]
边界上的位移、流体压力、应力、流量、酸液浓度通过强制边界条件给出：
[0281][0282]
式中，u为位移，为位移边界；pm为孔隙流体压力，为流体压力边界；σ为应力，n为裂缝壁面向量，为应力边界上施加的已知力，t为时间；km为基质渗透率张量，μ为流体粘度，为流体流量边界上施加的已知流量，ρ
l
为流体流量边界上已知的流体密度；c为酸液浓度。γd为位移边界，γ
p
为流体压力边界，γ
t
为应力边界，γq为流体流量边界。
[0283]
裂缝表面受到裂缝内的压力为：
[0284][0285]
式中，σ为应力，为裂缝两个壁面垂向单位向量，pf为裂缝内流体压力，为裂缝两个壁面。
[0286]
裂缝边界处的孔隙流体压力连续，在裂缝边界与基质之间存在流体交换：
[0287][0288]
式中，pm为岩石基质孔隙流体压力；pf为为裂缝内流体压力；km基质渗透率张量；μ流体粘度；为裂缝两个壁面垂向单位向量；q为流量；ρ为流体密度；为裂缝两个壁面。
[0289]
采用最大应变能释放率准则，对裂纹扩展进行判断；
[0290]
采用加权余量法将控制方程、边界条件结合，得到控制方程的等效积分弱形式：
[0291]
基质岩体中流体渗流连续性方程(2)的等效积分“弱”形式为：
[0292][0293]
式中，w为裂缝宽度；αm为biot系数；φm岩石基质孔隙度；ks固体岩石体积模量；kw流体体积模量；pm为岩石基质孔隙流体压力；t时间；im为单位张量，二维情况下为[1 1 0]
t
，三维情况下为[1 1 1 0 0 0]
t
；ε应变张量，km基质渗透率张量；μ流体粘度；q汇源项；ρ为流体密度；为流体流量边界上施加的已知流量。
[0294]
裂缝内流体流动连续性方程(5)的等效积分“弱”形式表示为：
[0295][0296]
式中，w裂缝宽度；t时间；μ流体粘度；pf为裂缝内孔隙流体压力；ρ为流体密度；为流体流量边界上施加的已知流量。
[0297]
将基质岩体中流体渗流连续性方程的等效积分弱形式(14)和裂缝内流体流动连续性方程的等效积分“弱”形式(15)相加，以消除基质和裂缝之间的流体交换项，得到：
[0298][0299]
采用扩展有限元法对控制方程进行离散，得到扩展有限元离散方程的步骤如下：
[0300]
采用扩展有限单元法近似逼近，其中位移可表示为：
[0301][0302]
式中，u为位移；ni、nj、nk、nh插值函数；ua常规位移节点自由度，ub裂纹面贯穿单元节点的额外位移自由度，uc裂纹前缘单元节点的额外位移自由度， ud为溶洞加强节点的额外位移自由度；
[0303]
h(x)为裂纹面函数：
[0304][0305]
线增函数r
α
(x)定义为：
[0306][0307]
溶洞加强函数r(x)定义为：
[0308][0309]
位移表达式(17)中的函数ψ(x)定义为：
[0310][0311]
将扩展有限元位移场逼近模式带入虚功方程，离散得到应力平衡扩展有限元控制方程：
[0312]
ku＝f
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(22)
[0313]
其中，k为整体劲度矩阵；f为等效载荷列阵；u为位移。
[0314]
同理可得到孔隙压力的加强模式：
[0315][0316]
式中，p为孔隙压力；ni、nj、nk插值函数；pa常规压力节点自由度，pb裂纹面贯穿单元节点的额外压力自由度，pc裂纹前缘单元节点的额外压力自由度； h(x)为裂纹面函数；r
α
(x)为线增函数。
[0317]
将扩展有限元压力场模式带入基质-裂缝渗流的等效积分“弱”形式，则流体渗流的扩展有限元控制方程为：
[0318][0319]
其中，hm和hf分别为基质流体整体渗流矩阵和裂缝流体整体渗流矩阵；dm和df分别为基质整体压缩性矩阵和裂缝整体压缩性矩阵；p为压力矩阵；f为等效载荷列阵。
[0320]
水力裂缝宽度w的表达式为：
[0321]
w＝tncd
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(25)
[0322]
其中，w为水力裂缝宽度；t为局部坐标和整体坐标的转换矩阵；nc为节点的形函数；d为单元节点参数矩阵。
[0323]
所述物理模型建立模块8执行如下操作：建立含天然裂缝、溶洞的物理模型，如下：孔洞内弹性模量、泊松比等材料参数设置为零，同时设置位移逼近为零；
[0324]
所述计算模块9执行如下操作：基于扩展有限元离散方程，对所述物理模型进行计算。
[0325]
以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容作出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。