一种基于水压致裂法确定渗透性岩体最大水平主应力的方法与流程

1.本发明涉及岩石断裂力学领域，尤其涉及一种基于水压致裂法确定渗透性岩体最大水平主应力的方法。


背景技术：

2.传统水压致裂法是一种能够测量深部地应力的有效方法，具有非常突出的优点：如无需知道岩体的力学参数，并且设备简单，操作方便，可在任意深度进行连续或重复测试等，因此在工程实践中应用广泛。但是目前水压致裂法地应力的求解过程是基于试验区域为低渗岩体的背景下开展的，如sl/t 264
‑
2020《水利水电工程岩石试验规程》中规定了水压致裂法适用的岩体渗透性不宜大于1lu。对于广泛存在的渗透性岩体，现有的水压致裂求解方法则表现出了极大的局限性，如何解决非完整岩体的水压致裂法地应力量值正确求解成为进一步提升工程应力测试的关键。
3.目前关于水压致裂法地应力计算的修正研究主要集中在测试系统本身的控制，如系统柔性影响、加载速度、曲线判定方式等对测试结果的影响。针对渗透性岩体的破裂压力计算，相关科研人员通过推导得出了破裂压力的下限值，但是无渗透状态下的理论解和渗透状态下的理论解无法进行统一描述。与此同时，由于水压致裂实际测试过程中的岩芯具有明显的区域变异性，钻孔中构造、变化复杂，因此破裂压力理论解中的岩芯泊松比和biot系数等参数需要结合水压致裂的实际进行确定，现有的计算方法对biot系数的确定方法较为复杂，无法满足水压致裂实际测试的需要。因此，有必要结合水压致裂的实际测试的特点，开展针对性的岩芯泊松比和biot系数测试，并基于渗透性岩体的实际的特点，设计新的渗透性岩体最大水平主应力的方法。


技术实现要素：

4.本发明为了解决以上问题，提供了一种基于水压致裂法确定渗透性岩体最大水平主应力的方法。
5.为实现上述目的，本发明所采用的技术方案如下：
6.一种基于水压致裂法确定渗透性岩体最大水平主应力的方法，包括以下步骤：
7.步骤1，通过水压致裂法测得岩体测试段的破裂压力p
b
，岩体的瞬时关闭压力p
s
，岩体重张压力p
r
，岩体中孔隙压力p0，并计算未经修正的最大水平主应力s
h0
和最小水平主应力s
h0
；
8.步骤2，根据水压致裂法所需的测试段深度h，取对应测试段的岩心制成标准测试试样；
9.步骤3，获得标准测试试样的体积v、试样截面积a、岩石孔隙率np以及岩石饱和密度ρ；
10.步骤4，将标准测试试样放置到三轴压缩试验装置中进行等向压缩试验，获得加载应力σ和轴向应变ε1、环向应变ε2、体积应变ε
v
的关系曲线；
11.步骤5，根据加载应力σ和体积应变ε
v
的关系曲线，并基于遗传算法获得biot系数α；
12.步骤6，获得泊松比υ；
13.步骤7，将得到的biot系数α和泊松比υ带入如下公式确定修正后的最大水平主应力的修正s
h
14.s
h
＝3s
h0
‑
k1p
r
‑
k2p0[0015][0016]
k2＝1 α
[0017]
可选的，步骤4中，获得加载应力σ和轴向应变ε1、环向应变ε2、体积应变ε
v
的关系曲线包括：
[0018]
步骤401，在加载应力之前对标准测试试样进行饱水；
[0019]
步骤402，将饱水后的标准测试试样用塑料薄膜包裹一层，并使用热缩管加热箍紧，保证标准测试试样表面无明显褶皱；
[0020]
步骤403，在标准测试试样上安装轴向应变测试仪和环向应变测试仪，并放置到三轴压缩试验装置中；
[0021]
步骤404，以0.05mpa/s的加载速率同步施加轴向压力σ1和围压σ3，试验过程中保持σ1＝σ3＝σ至预定压力值σ
p
，预定压力值其中s
v
是根据覆盖层厚度估算的自重应力，其值大小为ρgh；
[0022]
步骤405，根据加载过程中测得的轴向应变ε1、环向应变ε2以及计算体积应变ε
v
＝ε1 2ε2，获得加载应力σ和轴向应变ε1、环向应变ε2、体积应变ε
v
的关系曲线；
[0023]
步骤406，当加载到预定压力值σ
p
后，以0.05mpa/s的速率进行围压卸载，将围压σ3卸载至0.8σ
p
，卸载过程中保持轴向压力σ1不变，围压卸载完成后，计算卸载后和卸载前的轴向应变的变化量δε1和环向应变的变化量δε3。
[0024]
可选的，步骤5包括：
[0025]
步骤501，取水的压缩系数为2.19gpa，孔隙水应变占总应变的比例系数nw＝np；
[0026]
步骤502，生成包含50个个体的初始种群，每个个体包含一组k、n
c
、k
c
，设定k、n
c
、k
c
的取值范围分别为[k
min k
max
]、[n
cmin n
cmax
]、[k
cmin k
cmax
]，每个个体中的初始数值在取值范围内随机生成；
[0027]
步骤503，应力
‑
体积应变曲线中，包含的σ～ε
v
数据数量为j组，通过公式
[0028][0029]
其中，k为岩石骨架的体积模量，k
w
为孔隙水的体积模量，k
c
为孔隙结构的体积模量，n
w
为孔隙水应变占总应变的比例系数，n
c
为孔隙结构应变占总应变的比例系数；
[0030]
计算第i个个体的对应的加载应力
‑
体积应变中j组应力σ
i
对应的ε
vi
，通过公式计算个体i的适应度
[0031]
[0032]
步骤504，选择操作，第i个个体被选中的概率共选择50次，形成选择后种群；
[0033]
步骤505，交叉操作，从选择形成的种群中随机选择10对个体，进行交叉操作，其中每一对个体中，1个是保留个体a，其对应的参数数值为k
a
、n
ca
、k
ca
；另一个是替换个体b，其对应的参数值为k
b
、n
cb
、k
cb
；个体a和个体b交叉形成后的新个体为c，个体c作为新的个体放到种群中，替换掉个体b，个体c对应的参数值为k
c
、n
cc
、k
cc
，其计算公式如下
[0034]
k
c
＝mk
a
(1
‑
m)k
b
[0035]
n
cc
＝nn
ca
(1
‑
n)n
cb
[0036]
k
cc
＝rk
ca
(1
‑
r)k
cb
[0037]
其中，m、n、r为(0 1)之间随机生成的随机数，重复上述操作，完成10对个体的交叉操作，形成新的种群；
[0038]
步骤506，变异操作，从交叉操作形成的种群中，随机选择5个个体，进行变异操作，对于具体的个体s，其对应的参数数值为k
s
、n
cs
、k
cs
，变异后形成的新的个体为t，其对应的参数数值为k
t
、n
ct
、k
ct
，变异后的参数计算公式如下
[0039]
k
t
＝k
s
m(k
max
‑
k
s
)或k
t
＝k
s
‑
m(k
s
‑
k
min
)
[0040]
n
ct
＝n
cs
n(n
cmax
‑
n
cs
)或n
ct
＝n
cs
‑
n(n
cs
‑
n
cmin
)
[0041]
k
ct
＝k
cs
r(k
cmax
‑
k
cs
)或k
ct
＝k
cs
‑
r(k
cs
‑
k
cmin
)
[0042]
其中，每个参数都对应两个计算公式，具体计算时，每个参数对应的计算公式按50％概率随机选择计算，其中m、n、r为[0(1
‑
n/n
max
)]之间生成的随机数，其中n为当前遗传算法迭代次数，n
max
为遗传算法总的迭代次数；
[0043]
步骤507，重复步骤503～505一次为遗传算法迭代一次，直至最大的迭代次数n
max
，遗传算法结束，计算最后一代种群中适应度最大的个体，并输出对应的参数k、n
c
、k
c
；
[0044]
步骤508，确定k、n
c
、k
c
参数，根据公式(1)计算σ
p
和0.95σ
p
时对应的体积应变ε
vp
和ε
vp0.95
，按照下式计算k
σ
[0045]
k
σ
＝0.05σ
p
/(ε
vp
‑
ε
vp0.95
)
[0046]
步骤509，根据如下公式得到biot系数参数α
[0047]
α＝1
‑
k
σ
/k
[0048]
可选的，步骤6包括，由轴向应变的变化量δε1和环向应变的变化量δε3根据以下公式得到泊松比υ
[0049][0050]
可选的，还包括：步骤8，有多个标准测试试样时，最大水平主应力s
h
为多个试样计算结果剔除异常值后的平均值。
[0051]
可选的，标准测试试样尺寸为φ50
×
100cm。
[0052]
本发明与现有技术相比，所取得的技术进步在于：
[0053]
发明基于渗透岩体的破坏模式，通过推导建立水压致裂最大水平主应力的理论计算公式，并针对理论公式中的泊松比和biot系数等参数提出考虑水压致裂实际特点的室内试验方法，通过室内岩石试验和特定的应力路径获得中间参数，中间参数的求解过程采用
遗传算法理论，可有效避免了误差的叠加，同时具有计算快速，结果精确的优点。本发明根据弹性力学公式严格推导，理论严密，因此方法有严格的理论依据；综合考虑了岩体孔隙内气体和液体的影响程度，计算结果更加真实准确；现场试验与室内试验相结合，数据可靠性更高；运用遗传算法理论确定参数，避免了误差累积，计算更加便捷准确。
附图说明
[0054]
附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。
[0055]
在附图中：
[0056]
图1为本发明总体的实施流程图。
[0057]
图2为本发明实施例的加载应力σ和轴向应变ε1、环向应变ε2、体积应变ε
v
的关系曲线。
[0058]
图3为本发明实施例的遗传算法最优个体适应度随迭代次数的变化曲线。
具体实施方式
[0059]
下面这几个具体的实施例可以相互结合，对于相同或相似的概念或过程可能在某些实施例中不再赘述。下面将结合附图，对本发明的实施例进行描述。
[0060]
水压致裂的测试过程、压力时间曲线特征值参数与sl/t 264
‑
2020《水利水电工程岩石试验规程》相关条款保持一致。对于特定深度的钻孔而言，由于岩体抗拉强度不易确定，因此理论推导主要考虑重张阶段的应力过程。重张阶段，最大水平主应力方向上的环向应力可以分3个部分，原始应力的重分布、水压的环向加压和水压的渗透影响，根据弹性力学理论，三个环向应力的分量如下：
[0061]
s
θθ1
＝3s
h
‑
s
h
‑
p0[0062]
s
θθ2
＝
‑
p
w
[0063][0064]
式中sh、sh为钻孔横断面上的大、小平面主应力，mpa；p0为孔隙水压力；p
w
为水压致裂加压压力；s
θθ1
为原始应力作用下钻孔横截面在最大主应力方向上的环向应力分量；s
θθ2
为水压的环向加压形成的环向应力分量；s
θθ3
为水压的渗透影响形成的环向应力分量；α为钻孔对应深度岩石的biot系数；υ为钻孔对应深度岩石的泊松比。
[0065]
环向应力s
θθ
是总的环向应力，因此作用在钻孔孔壁上的环向有效拉应力σ
θθ
如下式所示
[0066][0067]
考虑重张阶段，岩体抗拉强度等于0，即当有效应力σ
θθ
＝0时，水压致裂裂缝重新张开，重张压力pr带入上式可得渗透岩体最大水平主应力表达式如下式所示
[0068]
s
h
＝3s
h
‑
k1p
r
‑
k2p0ꢀꢀ
(1)
[0069][0070]
k2＝1 α
[0071]
当钻孔岩壁无法渗透时，α＝0，上式退化为sl/t 264
‑
2020《水利水电工程岩石试验规程》中的计算方程，如下所示
[0072]
s
h
＝3s
h
‑
p
r
‑
p0[0073]
通过特定的室内岩石试验，求得参数α和υ，即可利用上式对传统水压致裂计算的最大水平主应力量值进行修正。
[0074]
其中参数α为biot系数，其的数学表达式如下式所示
[0075]
α＝1
‑
kσ/k
ꢀꢀ
(2)
[0076]
式中k
σ
—为球应力状态为σ时的岩石体积模量；k为岩石骨架的体积模量。其中岩石体积模量k
σ
受3部分影响，岩石骨架、孔隙水和孔隙结构本身，当岩石致密无孔隙时，k
σ
＝k，参数α＝0。
[0077]
根据岩石的组成，可将岩石变形分为3个部分的变形—岩石骨架、孔隙水和孔隙结构，其中孔隙水和孔隙结构的体积模量分别为k
w
和k
c
。考虑自然应变概念，则在球应力σ作用下，岩石本身的体积应变ε
v
可用下式求得
[0078][0079]
式中，n
w
为孔隙水应变占总应变的比例系数，n
c
为孔隙结构应变占总应变的比例系数。
[0080]
参数υ为试样泊松比，可通过三轴压缩试验获得。
[0081]
实施例一
[0082]
结合图1，步骤1，依据sl/t 264
‑
2020《水利水电工程岩石试验规程》，通过水压致裂法测得岩体测试段的破裂压力p
b
，岩体的瞬时关闭压力p
s
，岩体重张压力p
r
，岩体中孔隙压力p0，并计算未经修正的最大水平主应力s
h0
和最小水平主应力s
h0
。
[0083]
步骤2，按照水压致裂的测试段深度h，取对应测试段的岩心，按照sl/t 264
‑
2020《水利水电工程岩石试验规程》制成φ50
×
100cm的标准试样。
[0084]
步骤3，对于具体的试样，按照sl/t 264
‑
2020《水利水电工程岩石试验规程》相关章节，测试试样的体积v、试样截面积a、岩石孔隙率np、岩石饱和密度ρ。
[0085]
步骤4，将岩石放置到三轴压缩试验装置中进行等向压缩试验，获得等向压缩应力和体积应变ε
v
的关系曲线。
[0086]
步骤4中的等向压缩应力和体积应变关系曲线的可分为以下几个步骤：
[0087]
步骤401，对圆柱形试样加载之前，按照sl/t 264
‑
2020《水利水电工程岩石试验规程》进行饱水。
[0088]
步骤402，饱水后的试样，用塑料薄膜包裹一层，然后使用热缩管加热箍紧，保证表面无明显褶皱，安装轴向应变和环向应变测试仪器，放置到三轴压缩试验装置中。
[0089]
步骤403，以0.05mpa/s的加载速率同步施加轴向压力σ1和围压σ3，试验过程中保持σ1＝σ3＝σ至预定压力值σ
p
。预定压力值其中s
v
是根据覆盖层厚度估算的自重应力，其值大小等于ρgh。
[0090]
步骤404，根据加载过程中测得的轴向应变ε1和环向应变ε2计算体积应变ε
v
＝ε1 2ε2。最终可获得加载应力σ和轴向应变ε1、环向应变ε2、体积应变ε
v
的关系曲线。
[0091]
步骤405，当加载到预定压力值σ
p
后，以0.05mpa/s的速率进行围压卸载，将围压σ3卸载至0.8σ
p
，卸载过程中保持轴压σ1不变。围压卸载完成后，计算卸载后和卸载前的轴向应变的变化量δε1和环向应变的变化量δε3。
[0092]
步骤5，根据应力σ和体积应变ε
v
关系曲线，基于遗传算法理论计算参数α。
[0093]
所述步骤5中的计算参数α可分为以下几个步骤：
[0094]
步骤501,取水的压缩系数为2.19gpa，孔隙水应变占总应变的比例系数n
w
＝n
p
，则公式(2)中未知参数主要有3个：岩石骨架体积模量k、孔隙结构应变占总应变的比例系数n
c
、孔隙结构的体积模量分别为k
c
[0095]
步骤502，生成包含50个个体的初始种群，每个个体包含一组k、n
c
、k
c
，设定k、n
c
、k
c
的取值范围分别为[k
min k
max
]、[n
cmin n
cmax
]、[k
cmin k
cmax
]，每个体中的初始数值在取值范围内随机生成。
[0096]
步骤503，应力
‑
体积应变曲线中，包含的σ～ε
v
数据数量为j组，通过公式(3)计算第i个个体的对应的应力
‑
体积应变中j组应力σ
i
对应的ε
vi
，根据下式计算个体i的适应度
[0097][0098]
步骤504，选择操作，第i个个体被选中的概率共选择50次，形成选择后种群。
[0099]
交叉操作，从选择形成的种群中随机选择10对个体，进行交叉操作，其中每一对个体中，1个是保留个体a，其对应的参数数值为k
a
、n
ca
、k
ca
；另一个是替换个体b，其对应的参数值为k
b
、n
cb
、k
cb
；个体a和个体b交叉形成后的新个体为c，个体c作为新的个体放到种群中，替换掉个体b，个体c对应的参数值为k
c
、n
cc
、k
cc
，其计算公式如下
[0100]
k
c
＝mk
a
(1
‑
m)k
b
[0101]
n
cc
＝nn
ca
(1
‑
n)n
cb
[0102]
k
cc
＝rk
ca
(1
‑
r)k
cb
[0103]
其中，m、n、r为(0 1)之间随机生成的随机数，重复上述操作，完成10对个体的交叉操作，形成新的种群；
[0104]
步骤506，变异操作，从交叉操作形成的种群中，随机选择5个个体，进行变异操作，对于具体的个体s，其对应的参数数值为k
s
、n
cs
、k
cs
，变异后形成的新的个体为t，其对应的参数数值为k
t
、n
ct
、k
ct
，变异后的参数计算公式如下
[0105]
k
t
＝k
s
m(k
max
‑
k
s
)或k
t
＝k
s
‑
m(k
s
‑
k
min
)
[0106]
n
ct
＝n
cs
n(n
cmax
‑
n
cs
)或n
ct
＝n
cs
‑
n(n
cs
‑
n
cmin
)
[0107]
k
ct
＝k
cs
r(k
cmax
‑
k
cs
)或k
ct
＝k
cs
‑
r(k
cs
‑
k
cmin
)
[0108]
其中，每个参数都对应两个计算公式，具体计算时，每个参数对应的计算公式按50％概率随机选择计算，其中m、n、r为[0(1
‑
n/n
max
)]之间生成的随机数，其中n为当前遗传算法迭代次数，n
max
为遗传算法总的迭代次数；
[0109]
步骤507，重复步骤503～505一次为遗传算法迭代一次，直至最大的迭代次数n
max
，遗传算法结束，计算最后一代种群中适应度最大的个体，并输出对应的参数k、n
c
、k
c
；
[0110]
步骤508，确定k、n
c
、k
c
参数，根据公式(1)计算σ
p
和0.95σ
p
时对应的体积应变ε
vp
和ε
vp0.95
，按照下式计算k
σ
[0111]
k
σ
＝0.05σ
p
/(ε
vp
‑
ε
vp0.95
)
[0112]
步骤509，根据如下公式得到biot系数参数α
[0113]
α＝1
‑
k
σ
/k
[0114]
步骤6，根据弹性力学公式，围压卸载完成后，轴向应变的变化量δε1和环向应变的变化量δε3可用来计算泊松比υ，在此次试验的条件下，泊松比υ可用下式计算得到
[0115][0116]
步骤7，将计算得到的参数α和泊松比υ带入公式1完成水压致裂最大水平主应力的修正，修正后的最大主应力为s
h
。
[0117]
步骤8，有多个试样时，最大主应力s
h
为多个试样计算结果剔除异常值后的平均值。
[0118]
实施例二
[0119]
(1)本实施例的钻孔水压致裂测试在一花岗岩钻孔中进行，其中测试深度为h＝424m，测试及计算方法按照sl/t 264
‑
2020《水利水电工程岩石试验规程》进行，通过水压致裂法测得岩体测试段的破裂压力p
b
为20.47mpa，岩体的瞬时关闭压力p
s
为12.38mpa，岩体重张压力p
r
为14.66mpa，岩体中孔隙压力p0为3.37mpa，计算未经修正的最大水平主应力s
h0
和最小主水平应力s
h0
分别为19.10mpa和12.38mpa。
[0120]
(2)在水压致裂的测试段深度h＝424m，取对应测试段的岩心，按照sl/t 264
‑
2020《水利水电工程岩石试验规程》制成φ50
×
100cm的标准试样。
[0121]
(3)对于具体的试样，按照sl/t 264
‑
2020《水利水电工程岩石试验规程》相关章节，测试试样的基本参数，其中体积v＝196.35cm3、试样截面积a＝19.6cm2、岩石孔隙率np＝1.3％、岩石饱和密度ρ＝2.66g/cm3。
[0122]
(3)检查三轴压缩试验装置：装置是否密闭、管路是否畅通正常。将饱水后的岩石试样用塑料薄膜包裹一层，然后使用热缩管加热箍紧，保证表面无明显褶皱，安装轴向应变和环向应变测试仪器，放置到三轴压缩试验装置中。
[0123]
(4)试验：以0.05mpa/s的加载速率同步施加轴向压力σ1和围压σ3，试验过程中保持σ1＝σ3＝σ至预定压力值σ
p
。根据覆盖层厚度估算的自重应力s
v
，其值大小等于ρgh，其中g取9.8m/s2，计算得到s
v
＝11.01mpa。根据计算得14.16mpa。
[0124]
(5)根据加载过程中测得的轴向应变ε1和环向应变ε2计算体积应变ε
v
＝ε1 2ε2。最终可获得加载应力σ和轴向应变ε1、环向应变ε2、体积应变ε
v
的关系曲线，具体曲线如图2所示。
[0125]
(6)加载到预定压力值σ
p
＝14.16mpa后，以0.05mpa/s的速率进行围压卸载，将围压σ3卸载至0.8σ
p
，即卸载后的围压σ3＝11.33mpa，卸载过程中保持轴压σ1不变。围压卸载完成后，计算卸载后和卸载前的轴向应变的变化量δε1和环向应变的变化量δε3。实施例中轴向应变的变化量δε1＝5.44
×
10
‑
5，环向应变的变化量δε3＝8.16
×
10
‑
5。
[0126]
(7)取水的压缩系数为2.19gpa，孔隙水应变占总应变的比例系数n
w
＝n
p
，则公式
(2)中未知参数主要有3个：岩石骨架体积模量k、孔隙结构应变占总应变的比例系数n
c
、孔隙结构的体积模量分别为k
c
[0127]
(8)生成包含50个个体的初始种群，每个个体包含一组k、n
c
、k
c
，设定k、n
c
、k
c
的取值范围分别为[10gpa 40gpa]、[0 10％]、[0gpa 1gpa]，每个体中的初始数值在取值范围内随机生成。
[0128]
(9)本实施例中，加载应力
‑
体积应变曲线中，包含的σ～ε
v
数据数量为83组，通过公式(3)计算第i个个体的对应的应力
‑
体积应变中j组应力σ
i
对应的ε
vi
，根据下式计算个体i的适应度
[0129][0130]
(10)进行选择操作，第i个个体被选中的概率共选择50次，形成选择后种群。
[0131]
(11)进行交叉操作，从选择形成的种群中随机选择10对个体，进行交叉操作。其中每一对个体中，1个是保留个体a，其对应的参数数值为k
a
、n
ca
、k
ca
；另一个是替换个体b，其对应的参数值为k
b
、n
cb
、k
cb
；个体a和个体b交叉形成后的新个体为c，个体c作为新的个体放到种群中，替换掉个体b，个体c对应的参数值为k
c
、n
cc
、k
cc
，其计算公式如下
[0132]
k
c
＝mk
a
(1
‑
m)k
b
[0133]
n
cc
＝n n
ca
(1
‑
n)n
cb
[0134]
k
cc
＝rk
ca
(1
‑
r)k
cb
[0135]
其中，m、n、r为(0 1)之间随机生成的随机数。重复上述操作，完成10对个体的交叉操作，形成新的种群。
[0136]
(12)进行变异操作，从交叉操作形成的种群中，随机选择5个个体，进行变异操作，对于具体的个体s，其对应的参数数值为k
s
、n
cs
、k
cs
；变异后形成的新的个体为t，其对应的参数数值为k
t
、n
ct
、k
ct
，变异后的参数计算公式如下所示
[0137]
k
t
＝k
s
m(k
max
‑
k
s
)或k
t
＝k
s
‑
m(k
s
‑
k
min
)
[0138]
n
ct
＝n
cs
n(n
cmax
‑
n
cs
)或n
ct
＝n
cs
‑
n(n
cs
‑
n
cmin
)
[0139]
k
ct
＝k
cs
r(k
cmax
‑
k
cs
)或k
ct
＝k
cs
‑
r(k
cs
‑
k
cmin
)
[0140]
其中，每个参数都对应两个计算公式，具体计算时，每个参数对应的计算公式按50％概率随机选择计算。其中m、n、r为[0(1
‑
n/n
max
)]之间生成的随机数，其中n为当前遗传算法迭代次数，nmax为遗传算法总的迭代次数，本次实施例中n
max
＝20000。
[0141]
(13)重复步骤(9)～(12)一次为遗传算法迭代一次，直至最大的迭代次数n
max
＝20000，遗传算法结束，计算最后一代种群中适应度最大的个体，并输出对应的参数k、n
c
、k
c
。本次实施例计算得到的参数k、n
c
、k
c
分别为19.709gpa、3.347
×
10
‑4和3.5mpa，实施例的遗传算法最优个体适应度随迭代次数的变化曲线如图3所示。
[0142]
(14)根据确定的k、n
c
、k
c
等参数，根据公式3计算σ
p
和0.95σ
p
时对应的体积应变ε
vp
和ε
vp0.95
，按照k
σ
＝0.05σ
p
/(ε
vp
‑
ε
vp0.95
)计算k
σ
。实施例中计算得到的(ε
vp
‑
ε
vp0.95
)＝4.08
×
10
‑
5，k
σ
＝17.279gpa。
[0143]
(15)计算参数α＝1
‑
k
σ
/k，其中k
σ
＝17.279gpa，k＝19.709gpa，计算得到的参数α＝0.1233。
[0144]
(16)根据轴向应变的变化量δε1和环向应变的变化量δε3用来计算泊松比υ，计算得到的泊松比υ＝0.25。
[0145]
(17)将计算得到的参数α和泊松比υ带入s
h
＝3s
h
‑
k1p
r
‑
k2p0中计算修正后的最大水平主应力，其中k2＝1 α，实施例中计算得到的k1和k2分别为1.041和1.123，计算修正后的s
h
＝18.09mpa。
[0146]
最后应说明的是：以上所述仅为本发明的优选实施例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明权利要求保护的范围之内。