本发明属于油气田增产改造领域，特别是涉及一种基于相场法的弹塑性储层脉冲压裂模拟方法。
背景技术：
1、随着日益增长的油气资源需求和非常规油气资源的枯竭，不得不向非常规油气资源进军。目前，页岩气资源作为全球能源改革的重要组成部分。然而，深层页岩气储层具有低孔、低渗和塑性强的特点，常规开发方式难以满足需求，因此需要采用水力压裂技术在深层页岩储层中建立复杂的裂缝网络，以形成高速流动通道。水力脉冲压裂技术是近年来提出的一种新型压裂技术。目前的研究认为该技术是通过循环动态载荷使得岩石发生疲劳损伤，从而降低破裂压力并产生复杂裂缝网络。然而，水力脉冲压裂技术中水力裂缝起裂和扩展机理不明，因此揭示其机理具有重要意义。
2、在近二十年内，相场法在复杂裂缝扩展建模方面引起了广泛的兴趣，建立了脆性、准脆性、延性、弹塑性、粘弹性等相场模型，并被用于研究岩石、复合材料、混凝土、陶瓷和聚合物等材料的损伤和断裂问题。疲劳断裂的相场模型最近才被开发出来。然而，目前的疲劳相场模型以机械疲劳断裂问题为主，水力压裂领域的岩石疲劳损伤与传统的机械疲劳不同，水力裂缝涉及复杂的流固耦合作用，表现为储层应力和流体压力的相互影响，截止目前为止还未应用到动态水力压裂领域。因此建立一种基于相场法的弹塑性储层脉冲压裂模拟方法对于弹塑性储层压裂优化设计具有重要意义。
技术实现思路
1、本发明基于断裂力学、塑性力学、变分学、固体物理学、疲劳与失效等多学科知识，建立了一种基于相场法的弹塑性储层脉冲压裂模拟方法。
2、一种基于相场法的弹塑性储层脉冲压裂模拟方法，包括以下步骤：
3、s1：收集脉冲压裂模型所需的参数；
4、s2：基于热力学框架推导弹塑性疲劳相场模型；
5、s3：建立控制方程组的弱形式；
6、s4：对控制方程组进行有限元离散；
7、s5：构建控制方程组的线性化及迭代格式；
8、s6：将步骤s1的参数输入步骤s5，完成脉冲压裂的分析。
9、进一步的，所述步骤s1脉冲压裂模型所需的参数包括：最大水平主应力、最小水平主应力、储层埋深、地层压力、岩石的弹性模量和泊松比、岩石的孔隙度、岩石的抗拉强度、岩石和流体的密度、天然裂缝抗拉强度、天然裂缝渗透率、天然裂缝角度、剪胀角、内摩擦角、施工排量、压裂液粘度、压裂液密度。
10、进一步的，所述步骤s2基于热力学框架推导弹塑性疲劳相场模型，包括以下内容：
11、s2.1：平衡定律和热力学第二定律
12、通过热力学第一定律描述内能、动能和外部功率之间相互作用关系，其表示为：
13、
14、式中，ε，和分别表示净内能，动能和外部功率。
15、
16、式中，ρ为弹塑性介质的密度；e和分别为单位质量内能和位移的时间导数；和b为表面牵引力和体力；ω和π分别为裂缝表面的微牵引力和内部微力。
17、将方程(2)代入能量平衡方程(1)得到如下表达式：
18、
19、线性动量平衡方程可以写为：
20、
21、未受损的弹塑性介质总的应力张量σ可以分解为初始应力张量σ0和弹塑性应力张量σep：
22、σ＝σ0+σep (5)
23、将式(5)代入式(4)在域上进行积分再乘以域上位移的时间导数，可以得到：
24、
25、将式(6)代入式(3)得到如下的能量平衡方程：
26、
27、总应变张量ε可以分解为弹性应变张量εe和塑性应变张量εp：
28、
29、假设水力压裂是等温过程，采用克劳修斯-杜昂不等式来描述水力压裂不可逆的热力学过程，即满足单位体积的耗散能量非负，可以写为：
30、
31、将式(8)和(9)带入式(7)，则弹塑性介质的内部能量平衡方程写为：
32、
33、水力压裂过程中弹塑性介质总的亥姆霍兹自由能，包括初始应变能、弹性能、塑性能、断裂能和流体能，可以写为：
34、
35、式中
36、
37、将式(11)和(12)带入式(10)，化简为：
38、
39、上式必须在任意热力学过程成立。因此，该不等式需要满足耗散项的系数非负和非耗散项的系数为零。热力学一致模型成立的需要满足以下条件：
40、
41、s2.2：控制方程和本构方程推导
42、s2.2.1：流体流动方程
43、弹塑性介质的流体连续性方程可以写为：
44、
45、式中，ζ和v分别为流体体积含量增量和流速。
46、增量流体体积含量增量定义如下：
47、
48、式中，tr(·)为二阶张量的迹；tr(ε)为体积应变，包括弹性和塑性应变，也可写成εii；α(c)和m(c)分别为biot的系数和模量，其均是相场的函数。
49、biot的系数可以写为：
50、
51、式中，kc＝g(c)k表示弹塑性介质的退化体积模量，c为未受损的体积模量；g(c)为退化函数；ks为弹塑性介质骨架的体积模量。
52、相场的退化函数可以表示为：
53、g(c)＝(1-k)(1-c)2+k (18)
54、式中，k＝1×10-9。
55、biot的模量可以写为：
56、
57、式中，kf是流体的体积模量，φ是随着相场演化的孔隙度，表达式为：
58、φ＝φ0+tr(ε) (20)
59、式中，φ0为弹塑性介质未受损的基质孔隙度。
60、裂缝宽度可以写为：
61、ω＝hetr(ε) (21)
62、式中，he为网格尺寸。
63、弹塑性介质的流体能量密度可以表示为：
64、
65、采用达西定律来描述流体在基质和裂缝中的流动规律。忽略重力的影响，则流体流动速度v可以表示为：
66、
67、式中，k和μf分别为弹塑性介质的渗透率和流体的粘度。
68、岩石破裂成缝会提高渗透率，我们采用了一个隐式考虑裂缝中的泊泽维尔型流动的各向异性渗透率公式：
69、
70、式中，km为各向同性渗透率；ξ≥1为加权函数，控制渗透率增强强度，确保数值稳定性；nγ沿裂缝面的法向量。
71、将式(16)和(23)代入式(15)流体的连续性方程可以将改写为：
72、
73、s2.2.2：应力平衡方程
74、线性动量平衡方程式(4)忽略体力的情况下可以写为：
75、
76、弹性应变张量分解拉伸部分和压缩部分分解如下：
77、
78、式中，拉伸和压缩应变张量具体形式如下：
79、
80、式中，εea(a＝1，2，3)为主应变；na和nb是主应变方向；算子符号被定义为
81、
82、未受损的拉伸和压缩应变能密度如下：
83、
84、式中，λ和μ为拉梅常数。
85、只有在拉伸情况下岩石才能发生损伤，而压缩不会导致岩石发生损伤，因此受损岩石骨架的弹性应变能密度可以写为：
86、
87、将式(14)1带入式(26)，应力平衡方程改写为：
88、
89、将式(8)带入式(31)，得到应力平衡方程的最终形式，如下：
90、
91、式中，c是受损的四阶弹性张量刚度，表达式如下：
92、
93、式中c为heaviside函数；如果c＞0，hε(c)＝1；如果c≤0，hε(c)＝0；jijkl为四阶张量，jijkl＝δijδkl，其中δij和δkl是kronecker符号，p+和p-为四阶张量。
94、塑性储能密度形式表达如下：
95、
96、式中，χ表示taylor-quinney系数，即塑性功转化为热量的部分，具体形式是一般取0.9；1-χ表示塑性功转化为裂缝驱动力的部分。
97、初始应变能密度表达形式如下：
98、ψ0＝σ0:ε≈g(c)σ0:ε (35)
99、初始应变能密度参照弹性应变能密度分解形式可以分解为：
100、
101、式中
102、
103、s2.2.3：疲劳损伤模型
104、对断裂表面能进行修改，如下所示：
105、
106、式中，和是累计历史变量和退化函数。
107、累计历史变量可以写为：
108、
109、式中，τ为伪时间，由于疲劳损伤是由于岩石的弹性和塑性应变，以及流体能共同驱动的，即其中具体见式(59)。
110、疲劳退化函数是描述材料在循环损伤过程中如何抗断裂退化，表达式可以写为：
111、
112、式中，为疲劳裂纹阈值，当超过该阈值，疲劳损伤才会开始累计。
113、临界能量释放率gc与岩石的抗拉强度的关系可以写为：
114、
115、式中，e是弹塑性介质的弹性模量，l0是长度尺度参数，一般用于描述裂缝的宽度。
116、s2.2.4：相场方程
117、采用微力平衡方程刻画水力裂缝的演化过程，不考虑微惯性和外部微力的影响，则弹塑性岩石的微力平衡方程可以写为：
118、
119、式中，w和π是分别是微观牵引力和内部微观力。
120、将式(22)、(30)、(34)、(38)和带入式(14)3和(14)4，并结合式(42)得到相场方程如下：
121、
122、在流体的影响下，岩石破坏是不可逆转的，因此需要使用历史状态变量来保证损伤的不可恢复性。另外，塑性应变是永久性变形，无法逆转，因此塑性应变可以单调增加，无需额外约束。其表达式如下：
123、
124、将历史状态变量的表达式(44)带入式(43)，相场方程可以改写为：
125、
126、s2.2.5：弹塑性本构方程
127、弹塑性本构关系通过增量的形式来表达：
128、dσ＝c：(dε-dεp)+g(c)σ0-α(c)dp1 (46)
129、考虑孔隙压力对岩石非线性变形的影响，采用biot有效应力进行本构方程的推导，drucker–prager屈服表达式写为：
130、
131、式中，s(σ′)表示偏应力张量；j2和i1分别是偏应力张量的第二不变量和有效应力张量的第一不变量。和κ分别为强化参数和硬化函数，表达式写为：
132、
133、式中，和分别为内摩擦角和等效塑性应变；a、b、c为强化参数，表达式如下：
134、
135、式中，hr为储层埋深。
136、硬化函数可以写为：
137、
138、页岩的力学行为中，屈服条件与实际塑性变形方向之间并不总是一致。因此非关联流动法则可以更准确地描述页岩在加载条件下的应变响应。非关联的流动法则可以写为：
139、
140、式中，θ为岩石剪胀角。
141、根据非关联流动法则，塑性应变增量和等效塑性应变增量表达式为：
142、
143、式中，是偏应力，dγ和分别表示塑性乘数和塑性流动矢量。
144、根据biot有效应力理论，式(46)可以改写为：
145、dσ'＝c:(dε-dεp)+g(c)σ0 (53)
146、当岩石发生屈服时，必须满足一致性条件：
147、
148、将式(52)和(53)代入式(54)得到：
149、
150、整理化简，塑性乘数表达式为：
151、
152、将式(56)带入式(53)可得到有效应力与应变的关系式，表达式如下：
153、
154、式中，ceq为四阶切线弹塑性算子。
155、
156、s2.3：混合模式裂缝的修正驱动力
157、将弹性应变能密度分解为拉伸部分和拉剪部分，初始应变能密度分解为拉伸部分和拉剪部分。流体能密度被视为i型能，塑性耗散能密度视为ii型能：
158、
159、式中
160、
161、进一步的，所述步骤s3建立控制方程组的弱形式，包括以下内容：
162、水力脉冲压裂模型的控制方程强形式：
163、
164、对应的边界条件为：
165、
166、式中，和γ分别为位移场、压力场和裂缝相场的dirichlect边界；和分别是位移场、压力场和裂缝相场的neumann边界，满足以下条件：
167、
168、对控制方程中的应力平衡方程(61)1、流体连续性方程(61)2和相场方程(61)3分别与位移wu、压力wp和相场wc权函数相乘。然后在计算域上积分，采用分部积分法。并利用散度定理以及结合边界条件式(62)，得到控制方程和边界条件的等效积分“弱”形式如下：
169、
170、进一步所述步骤s4对控制方程进行有限元离散，包括以下内容：
171、在离散域上，对位移、压力、相场、权函数和梯度进行以下近似：
172、
173、式中，上标h代表单元节点的值，nu、np和nc分别为位移、压力和裂缝相场的插值形函数；bu和分别为应变矩阵和体积应变矩阵，bp和bc分别为压力和裂缝相场插值形函数的梯度。
174、将式插值函数(65)带入弱形式(64)，可得：
175、
176、进一步所述步骤s5构建控制方程组的线性化及迭代格式，包括以下内容：
177、通过向后隐式差分法，对方程(66)中的时间导数的进行离散。
178、
179、式中，n和n+1分别表示前一个时间步和当前时间步。
180、采用经典的newmark迭代算法求解加速度，即：
181、
182、式中，e1∈[0，0.5]，e2∈[0，1]，取e1＝0.25，e2＝0.5的带隐式的newmark公式。
183、
184、将所有在第n+1个时间步未知量的下标删掉。式(67)和(69)带入式(66)，改写为余量的形式，可得：
185、
186、采用交错迭代格式对水力脉冲压裂数值模型进行求解。其中，采用nr迭代法求解流固耦合方程组。流固耦合方程组在第i个迭代步的nr迭代格式如下：
187、
188、式中
189、
190、将相场进行固定，求解式(71)。得到第i个迭代步的位移增量δuh和压力增量δph后，进而得到第i+1个迭代步位移和压力的试探解如下：
191、
192、然后，固定i+1个迭代步的位移和压力，然后通过求解方程(74)可得到第i+1个迭代步相场的试探解：
193、
194、式中
195、
196、当节点位移、压力和相场都满足以下收敛条件时，该时间步的迭代过程结束，进入下一个时间步：
197、
198、式中，和分别为位移、压力和相场的收敛容差。
199、本发明提出了一个新的热力学一致性的动态流固耦合混合模式疲劳弹塑性相场模型，用于模拟弹塑性岩石在循环流体载荷下的裂缝成核和扩展。我们在相场框架中引入了疲劳历史应变参数来捕捉疲劳效应，结合非关联的drucker–prager本构模型和非线性饱和应变函数来描述页岩的弹塑性行为。相场的驱动力包括塑性耗散能、流体能、岩石的初始应变能和弹性能。采用交错格式求解全耦合问题。模型在有限元框架内，采用newmark积分格式，并通过newton-raphson迭代算法求解导出的非线性方程组。该模型能够准确捕捉复杂裂缝形态。为实际工程应用提供了更为精确和可靠的解决方案，有望在石油和天然气等领域得到广泛应用。
1.一种基于相场法的弹塑性储层脉冲压裂模拟方法，其特征在于，包括以下步骤：
2.根据权利要求1所述的一种基于相场法的弹塑性储层脉冲压裂模拟方法，其特征在于所述步骤s1脉冲压裂模型所需的参数包括：最大水平主应力、最小水平主应力、储层埋深、地层压力、岩石的弹性模量和泊松比、岩石的孔隙度、岩石的抗拉强度、岩石和流体的密度、天然裂缝抗拉强度、天然裂缝渗透率、天然裂缝角度、剪胀角、内摩擦角、施工排量、压裂液粘度、压裂液密度。
3.根据权利要求1所述的一种基于相场法的弹塑性储层脉冲压裂模拟方法，其特征在于所述步骤s2基于热力学框架推导弹塑性疲劳相场模型，包括以下内容：
4.根据权利要求1所述的一种基于相场法的弹塑性储层脉冲压裂模拟方法，其特征在于所述步骤s3建立控制方程组的弱形式，包括以下内容：
5.根据权利要求1所述的一种基于相场法的弹塑性储层脉冲压裂模拟方法，其特征在于所述步骤s4对控制方程进行有限元离散，包括以下内容：
6.根据权利要求1所述的一种基于相场法的弹塑性储层脉冲压裂模拟方法，其特征在于所述步骤s5构建控制方程组的线性化及迭代格式，包括以下内容：
