一种适用于水力驱动裂缝的自压内聚力本构模型构建方法与流程

1.本发明涉及材料破坏与失效力学描述技术领域，特别涉及一种适用于水力驱动裂缝的自压内聚力本构模型构建方法。


背景技术：

2.对于常规裂缝扩展、材料失效分层、塑性延展的模拟，常采用内聚力模型进行分析，其控制方程简洁，且易于与有限元方法进行结合。并且，依据裂缝扩展和材料损伤的环境特性，内聚力模型存在多种扩展类别，如针对疲劳损伤下的裂缝扩展的循环内聚力模型，针对复合材料渐近损伤的锯齿内聚力模型。
3.裂缝扩展过程中，会生成新的表面，通常新生成的裂缝面为自由表面，不会受到直接的外力作用。但是在液体环境下，如大坝中形成的裂缝、油气储层的水力驱动裂缝、承压设备内表面的裂缝，新生成的裂缝表面会受到随后的水压载荷。但常规内聚力单元失效后结合力为零，若需要考虑裂缝在扩展过程中水力侵入和液体压力的影响，需要进行流体模拟进行流固耦合分析，即孔压内聚力模型。孔压内聚力模型需要在计算过程中增加流体渗流场，当裂缝中流体流动的动态效应可忽略且介质渗透率极低时计算效率低且收敛难度大。


技术实现要素：

4.为了对披露的实施例的一些方面有一个基本的理解，下面给出了简单的概括。该概括部分不是泛泛评述，也不是要确定关键/重要组成元素或描绘这些实施例的保护范围。其唯一目的是用简单的形式呈现一些概念，以此作为后面的详细说明的序言。
5.本发明的目的在于提供一种水力驱动裂缝起裂破坏的本构模型，解决现有孔压内聚力模型存在的计算量大和收敛难的问题。
6.本发明实施例提供的一种适用于水力驱动裂缝的自压内聚力本构模型构建方法，其改进之处在于，包括：
7.(1)将材料失效分为三个阶段；
8.(2)重构内聚力的应力位移关系；
9.(3)建立自压内聚力模型在数值实现中的刚度矩阵；
10.(4)通过向后时间差分计算刚度，离散化自压内聚力单元的积分方程，并根据虚功原理得到最终离散形式。
11.优选的，所述步骤(1)包括弹性变形阶段、损伤扩展阶段和水力侵入阶段；
12.所述弹性变形阶段为受激励载荷应力弹性上升阶段，所述损伤扩展阶段为超过临界应力/位移逐渐损伤失效阶段，所述水力侵入阶段为完全失效后自由表面阶段。
13.优选的，所述步骤(2)包括重构内聚力的应力位移关系，使得内聚力单元失效后的最初法向牵引力等同于液体压力值，并在第三个阶段，修正损伤参数和法向应力，对于指数内聚力模型的自压特性，采用sigmoid函数对其应力法向分量进行修正。
14.进一步的，所述在第三个阶段表面法向应力为零时：
[0015][0016]
紧凑模式下的内聚力模型中应力位移表达式为：
[0017]
t＝kδ＝(1-d)k0δ；
[0018]
公式变形得：
[0019][0020]
其中，δ是位移向量，包含法相位移δn与切向位移δ
t
，下角标n、t指代参数的法相与切向分量；k为刚度；d为损伤参数，范围从0到1，代表失效状态；当d＝0时，内聚力单元仍在弹性区域，当d开始增加，代表材料开始发生开裂失效，1代表完全失效；代表最大应力下的分离位移，δc代表完全失效时的临界位移；max下角标则用于指代损伤加载过程中已达到的最大值。
[0021]
进一步的，通过流固耦合或渗漏耦合的方式加载在裂缝表面，当裂缝向前扩展时进行流动和形变的迭代计算，将法向应力进行如下修正：
[0022][0023]
损伤参数则修正为：
[0024][0025]
优选的，所述步骤(3)包括将构建好的自压内聚力本构模型与有限元进行结合进
行数值实现，根据本构模型计算得出非线性迭代过程中的刚度矩阵，负值刚度的数值意义代表液压在裂缝形成过程，引入粘度避免收敛问题。
[0026]
进一步的，所述步骤(3)包括自压内聚力模型与有限元进行结合进行数值实现过程中，在非线性迭代过程中的刚度矩阵为：
[0027][0028]
指数自压内聚力本构方程的构建方法如下：
[0029]
指数内聚力模型其应力位移曲线为连续的，应力为势能的导数，势能可表示为：
[0030][0031]
应力通过计算得到：
[0032][0033][0034][0035]
采用sigmoid函数对应力法向分量进行修正：
[0036][0037]
其中，a是与液体压力侵入强度相关的参数，b代表液体侵入新裂缝的程度。
[0038]
优选的，所述步骤(4)包括
[0039]
通过粘性损伤参数dv来代替损伤参数：
[0040][0041]
在数值实现中，通过向后时间差分进行计算：
[0042][0043]
[0044]
刚度计算公式为：
[0045][0046]
根据虚功原理，自压内聚力单元的积分方程如下：
[0047][0048]
通过高斯散度定理进行简化：
[0049][0050]
对自压内聚力单元进行离散化：
[0051][0052]
其中，r为将位移向量从全局坐标系转换为局部坐标系的矩阵：
[0053][0054][0055]
l矩阵则是获得裂缝分离量的变换矩阵：
[0056][0057]
由于方程适用于任意虚位移向量，通过约去δ{ue}
′
，可获得自压内聚力单元的最终离散形式：
[0058][0059]
本发明实施例提供的技术方案可以包括以下有益效果：
[0060]
本发明通过重构自压内聚力模型，简化了受水力驱动裂缝的模拟过程，在不引入流体扩散的物理模型的条件下，仍达到了裂缝扩展新表面受到液体压力作用的效果。
[0061]
本发明中构建的自压内聚力模型，提高了受水力驱动裂缝扩展的计算效率和结果
收敛度。
[0062]
本发明提出的自压内聚力本构模型，适用于受水力驱动因素影响的裂缝扩展模拟。裂缝开裂后受到的液压通过内聚力本构中负值界面粘合力的方式施加，也就是说，开裂的自由表面在不进行流场耦合的条件下自动加载等同于液体压力的法向力。
[0063]
应当理解的是，以上的一般描述和后文的细节描述仅是示例性和解释性的，并不能限制本发明。
附图说明
[0064]
此处的附图被并入说明书中并构成本说明书的一部分，示出了符合本发明的实施例，并与说明书一起用于解释本发明的原理。
[0065]
图1是根据一示例性实施例示出的双线性自压内聚力模型应力位移示意图。
[0066]
图2是根据一示例性实施例示出的指数自压内聚力模型应力位移示意图。
[0067]
图3是根据一示例性实施例示出的自压内聚力模型在不同水压下失效过程示意图。
[0068]
图4是根据一示例性实施例示出的一种适用于水力驱动裂缝的自压内聚力本构模型构建方法流程示意图。
具体实施方式
[0069]
以下描述和附图充分地示出本发明的具体实施方案，以使本领域的技术人员能够实践它们。实施例仅代表可能的变化。除非明确要求，否则单独的部件和功能是可选的，并且操作的顺序可以变化。一些实施方案的部分和特征可以被包括在或替换其他实施方案的部分和特征。本发明的实施方案的范围包括权利要求书的整个范围，以及权利要求书的所有可获得的等同物。在本文中，各实施方案可以被单独地或总地用术语“发明”来表示，这仅仅是为了方便，并且如果事实上公开了超过一个的发明，不是要自动地限制该应用的范围为任何单个发明或发明构思。本文中，诸如第一和第二等之类的关系术语仅仅用于将一个实体或者操作与另一个实体或操作区分开来，而不要求或者暗示这些实体或操作之间存在任何实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素。本文中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的结构、产品等而言，由于其与实施例公开的部分相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
[0070]
下面结合附图及实施例对本发明做进一步描述：
[0071]
本发明从内聚力本构方程出发，开发了一种能够描述液力环境下裂缝扩展的自压内聚力模型。该模型基于双线性及指数内聚力模型构建，用于描述液力环境下的裂缝受压起裂及延展裂缝由于液体侵入而持续受压的过程。
[0072]
如图4所示，本发明提供了一种适用于水力驱动裂缝的自压内聚力本构模型构建方法，具体方案如下：
[0073]
材料开裂之前，潜在的开裂表面存在粘合力。开裂后，裂缝开度超过临界值，液体
进入产产生的裂缝表面，那么，新表面的最终载荷为液体压力。将内聚力的应力位移关系进行重构，使得内聚力单元失效后的最初法向牵引力等同于液体压力值。
[0074]
1、将材料失效分为三个阶段，包括弹性变形阶段、损伤扩展阶段和水力侵入阶段；
[0075]
弹性变形阶段为受激励载荷应力弹性上升阶段，所述损伤扩展阶段为超过临界应力/位移逐渐损伤失效阶段，所述水力侵入阶段为完全失效后自由表面阶段。
[0076]
2、重构内聚力的应力位移关系，使得内聚力单元失效后的最初法向牵引力等同于液体压力值，并在第三个阶段，修正损伤参数和法向应力，对于指数内聚力模型的自压特性，采用sigmoid函数对其应力法向分量进行修正；
[0077]
在第三个阶段表面法向应力为零时：
[0078][0079]
紧凑模式下的内聚力模型中应力位移表达式为：
[0080]
t＝kδ＝(1-d)k0δ；
[0081]
其中，
[0082]
δ是位移向量，包含法相位移δn与切向位移δ
t
，下角标n、t指代参数的法相与切向分量；k为刚度；d为损伤参数，范围从0到1，它代表失效状态。当d＝0时，内聚力单元仍在弹性区域，当d开始增加，代表材料开始发生开裂失效，1代表完全失效。代表最大应力下的分离位移，δc代表完全失效时的临界位移。max下角标则用于指代损伤加载过程中已达到的最大值。
[0083]
对于水力驱动的裂缝，由于液体会侵入裂缝间隙并在扩展过程中持续流动，新产生表面的最终法向应力与液体压力p一致。液体压力一般作为外力载荷，通过流固耦合或渗漏耦合的方式加载在裂缝表面，并且当裂缝向前扩展时进行流动和形变的迭代计算。若忽略流体流动的动态效应，将法向应力进行如下修正：
[0084]
[0085]
损伤参数则修正为：
[0086][0087]
在自压内聚力本构中，可以看到，损伤参数的值可以超过1，这代表液体压力在新的裂缝表面的侵入。同时，裂缝的切向应力不受液压影响。如图1、2所示，从自压内聚力的应力位移图可看出，存在负值刚度，其数值意义代表液压在裂缝形成过程中做的功。可以看到，通过本构方程的构建，计算过程中可省略流体相的物理耦合，新裂缝表面形成后的水力渗入直接通过内聚力中的负值法向力实现。
[0088]
3、将构建好的自压内聚力本构模型与有限元进行结合进行数值实现，根据本构模型计算得出非线性迭代过程中的刚度矩阵，负值刚度的数值意义代表液压在裂缝形成过程，引入粘度避免收敛问题；
[0089]
自压内聚力模型与有限元进行结合进行数值实现过程中，在非线性迭代过程中的刚度矩阵为：
[0090][0091]
指数自压内聚力本构方程的构建方法如下：
[0092]
指数内聚力模型是基于势能概念开发的，其应力位移曲线是连续的，应力为势能的导数。势能可表示为：
[0093][0094]
因此，应力可通过计算得到：
[0095][0096][0097][0098]
构建指数内聚力模型的自压特性，需要考虑保持其应力位移曲线的平滑特性，因
此，采用sigmoid函数对应力法向分量进行修正：
[0099][0100]
其中，a是与液体压力侵入强度相关的参数，b代表液体侵入新裂缝的程度。
[0101]
4、在数值实现中，通过向后时间差分计算刚度，离散化自压内聚力单元的积分方程，并根据虚功原理得到其最终离散形式；
[0102]
从指数自压内聚力模型的应力位移曲线可以看出，裂缝前端损伤开始前的弹性阶段与常规指数模型是相同的。损伤开始后，流体进入新生成的表面，法向应力随后受到影响并逐渐降低为零，表示总体失效，后继续降低为负值，即表示液体全部侵入。其中，剪切应力不受影响。
[0103]
由于在本构模型的构建中出现负刚度，为在隐式求解中获得收敛的解，通过引入一个小的粘度来避免收敛问题。
[0104]
在实际计算中，通过粘性损伤参数dv来代替损伤参数：
[0105][0106]
在数值实现中，通过向后时间差分进行计算：
[0107][0108][0109]
最终的刚度计算公式为：
[0110][0111]
在自压内聚力模型结合有限元方法的数值实现中，自压内聚力单元不含有液体压力自由度，只需要考虑位移分量，大大降低了计算量。且由于液压通过法向表面力的形式参与计算，无需在右侧力向量中额外添加。根据虚功原理，自压内聚力单元的积分方程如下：
[0112][0113]
通过高斯散度定理进行简化：
[0114][0115]
对自压内聚力单元进行离散化：
[0116][0117]
其中，r为将位移向量从全局坐标系转换为局部坐标系的矩阵：
[0118][0119][0120]
l矩阵则是获得裂缝分离量的变换矩阵：
[0121][0122]
由于方程适用于任意虚位移向量，那么通过约去δ{ue}
′
，可获得自压内聚力单元的最终离散形式：
[0123][0124]
实施例1
[0125]
通过自压内聚力模型对i-型水力驱动裂缝进行有限元分析。
[0126]
通过有限元软件中提供的用户自定义单元子程序结构，导入用fortran语言自编的自压内聚力本构模型，从而模拟水力驱动裂缝。
[0127]
在用户自定义单元子程序中，以双线性自压内聚力模型为例，其应力-应变状态依据以下公式定义，即本构模型：
[0128][0129]
单元状态变量用d表示：
[0130][0131]
刚度计算部分代码如下：
[0132][0133]
单元刚度矩阵及右侧向量构建部分代码如下：
[0134][0135]
采用下述方法验证新产生裂缝表面是否会自动加载水力压力：
[0136]
定义裂缝的初始长度，假使初始裂缝长度为10mm，材料杨氏模量为3gpa，泊松比为0.25，裂缝开裂过程应力极值对应临界位移为1mm，液体在位移2mm全部侵入。当初始裂缝受到0.1mpa、0.5mpa和1mpa的水力驱动时，对自压内聚力模型的载荷-位移曲线进行分析。
[0137]
从图3(a)中可以看出随着水力压力的增加，损伤参数d从0开始增加，常规内聚力单元，d的最终值为1，代表没有面应力的自由表面。对于自压内聚力模型，d的最终值高于1，表示表面应力为负值，实际代表新生成裂缝面开始再次受到水力载荷。d的值越高，表示水力载荷越大。图3(b)显示了裂缝前端第一个界面单元的法向应力的变化，最终的法向应力与水压值相符。
[0138]
应当理解的是，本发明并不局限于上面已经描述并在附图中示出的流程及结构，并且可以在不脱离其范围进行各种修改和改变。本发明的范围仅由所附的权利要求来限制。

技术特征：
1.一种适用于水力驱动裂缝的自压内聚力本构模型构建方法，其特征在于，包括：(1)将材料失效分为三个阶段；(2)重构内聚力的应力位移关系；(3)建立自压内聚力模型在数值实现中的刚度矩阵；(4)通过向后时间差分计算刚度，离散化自压内聚力单元的积分方程，并根据虚功原理得到最终离散形式。2.根据权利要求1所述的一种适用于水力驱动裂缝的自压内聚力本构模型构建方法，其特征在于，所述步骤(1)包括弹性变形阶段、损伤扩展阶段和水力侵入阶段；所述弹性变形阶段为受激励载荷应力弹性上升阶段，所述损伤扩展阶段为超过临界应力/位移逐渐损伤失效阶段，所述水力侵入阶段为完全失效后自由表面阶段。3.根据权利要求1所述的一种适用于水力驱动裂缝的自压内聚力本构模型构建方法，其特征在于，所述步骤(2)包括重构内聚力的应力位移关系，使得内聚力单元失效后的最初法向牵引力等同于液体压力值，并在第三个阶段，修正损伤参数和法向应力，对于指数内聚力模型的自压特性，采用sigmoid函数对其应力法向分量进行修正。4.根据权利要求3所述的一种适用于水力驱动裂缝的自压内聚力本构模型构建方法，其特征在于，所述在第三个阶段表面法向应力为零时：紧凑模式下的内聚力模型中应力位移表达式为：t＝kδ＝(1-d)k0δ；公式变形得：其中，δ是位移向量，包含法相位移δ
n
与切向位移δ
t
，下角标n、t指代参数的法相与切向分量；k为刚度；d为损伤参数，范围从0到1，代表失效状态；当d＝0时，内聚力单元仍在弹性区域，当d开始增加，代表材料开始发生开裂失效，1代表完全失效；代表最大应力下的分离位移，δ
c
代表完全失效时的临界位移；max下角标则用于指代损伤加载过程中已达到的最大值。5.根据权利要求4所述的一种适用于水力驱动裂缝的自压内聚力本构模型构建方法，其特征在于，通过流固耦合或渗漏耦合的方式加载在裂缝表面，当裂缝向前扩展时进行流动和形变的迭代计算，将法向应力进行如下修正：
损伤参数则修正为：6.根据权利要求1所述的一种适用于水力驱动裂缝的自压内聚力本构模型构建方法，其特征在于，所述步骤(3)包括将构建好的自压内聚力本构模型与有限元进行结合进行数值实现，根据本构模型计算得出非线性迭代过程中的刚度矩阵，负值刚度的数值意义代表液压在裂缝形成过程，引入粘度避免收敛问题。7.根据权利要求6所述的一种适用于水力驱动裂缝的自压内聚力本构模型构建方法，其特征在于，所述步骤(3)包括自压内聚力模型与有限元进行结合进行数值实现过程中，在非线性迭代过程中的刚度矩阵为：指数自压内聚力本构方程的构建方法如下：指数内聚力模型其应力位移曲线为连续的，应力为势能的导数，势能可表示为：应力通过计算得到：应力通过计算得到：
采用sigmoid函数对应力法向分量进行修正：其中，a是与液体压力侵入强度相关的参数，b代表液体侵入新裂缝的程度。8.根据权利要求1所述的一种适用于水力驱动裂缝的自压内聚力本构模型构建方法，其特征在于，所述步骤(4)包括通过粘性损伤参数dv来代替损伤参数：在数值实现中，通过向后时间差分进行计算：在数值实现中，通过向后时间差分进行计算：刚度计算公式为：根据虚功原理，自压内聚力单元的积分方程如下：通过高斯散度定理进行简化：对自压内聚力单元进行离散化：其中，r为将位移向量从全局坐标系转换为局部坐标系的矩阵：
l矩阵则是获得裂缝分离量的变换矩阵：由于方程适用于任意虚位移向量，通过约去δ{u
e
}
′
，可获得自压内聚力单元的最终离散形式：

技术总结
本发明一种适用于水力驱动裂缝的自压内聚力本构模型构建方法，包括：(1)将材料失效分为三个阶段；(2)重构内聚力的应力位移关系；(3)建立自压内聚力模型在数值实现中的刚度矩阵；(4)通过向后时间差分计算刚度，离散化自压内聚力单元的积分方程，并根据虚功原理得到最终离散形式。本发明通过重构自压内聚力模型，简化了受水力驱动裂缝的模拟过程，在不引入流体扩散的物理模型的条件下，仍达到了裂缝扩展新表面受到液体压力作用的效果。新表面受到液体压力作用的效果。新表面受到液体压力作用的效果。


技术研发人员：宋慧芳 徐明明 张瑾 李风 于洲 王晓君 胥成亮 魏新晨 唐永安 李浩程
受保护的技术使用者：中国石油化工股份有限公司胜利油田分公司
技术研发日：2023.05.31
技术公布日：2023/9/14