基于Matlab的岩石蠕变模型参数辨识算法设计

张华宾,闵小东,张顷顷

(辽宁工程技术大学 力学与工程学院，辽宁 阜新 123000)



基于Matlab的岩石蠕变模型参数辨识算法设计

张华宾,闵小东,张顷顷

(辽宁工程技术大学 力学与工程学院，辽宁 阜新 123000)

为减少岩石蠕变模型参数辨识的繁冗计算，探讨一种可靠、高效的岩石蠕变模型辨识方法。以M-K(伯格斯)蠕变模型为例，给出蠕变模型参数辨识方法，并编写Matlab算法函数通过COM组件从而剥离出程序，借助VB环境下编制面向对象的可视化计算软件。结果表明，该软件用于岩石蠕变模型参数辨识的计算过程同样有较好的效果并提高了计算效率，为实现其他岩石力学参数计算问题提供了可供借鉴的快速处理方法。

蠕变模型；算法设计；面向对象；可视化

岩石蠕变力学性质模型大致可以分为以下4类：经验模型、元件组合模型、基于损伤力学的蠕变本构模型和基于热力学理论的蠕变本构模型。其中基本元件模型通过组合可以得到很多模型来描述岩土的流变特性，同时蠕变模型的参数辨识方法也很多，根据蠕变试验数据选择合适的蠕变模型确定相应的蠕变参数，识别参数可以获得较吻合试验数据的拟合曲线，主要的参数辨识方法有回归分析，最小二乘，曲线分解等方法。董志宏[1]等基于现场模型洞开挖位移检测结果，运用均匀设计—神经网罗—遗传算法对软岩蠕变参数进行了反演分析，得到该区软岩的蠕变参数。孙钧[2]教授等用最小二乘法建立目标函数，通过全局最优化，同时辨识了初始地应力三个分量及变形和强度3个参数，从实践上得出了非线性逆问题的唯一解。蒋昱州[3]提出 BFGS-LSM 算法对岩石非线性黏弹塑性蠕变模型参数进行辨识，结果和试验曲线吻合较好。林元华[4]建立了岩盐本构模型辨识的目标函数，提出了以井径测量信息为基础反演确定岩盐层蠕变特性参数的位移反分析法。李青麒[5]介绍一种根据室内蠕变资料确定流变参数和选择流变模型的方法。陈炳瑞[6]提出模式搜索与非线性最小二乘法的有机结合识别流变模型参数的方法。袁鸿鹄[7]基于流变反演理论，建立油房湾隧道围岩流变变形特征的流变模型。但无论何种方法得到准确合理的模型参数，都需要反复提取数据计算[8-12]。 然而数学软件Matlab进行数据处理和数据显示方面的优势非常显著，同时考虑VB编写的程序执行效率较高，可以提供交互性很好的界面，因此二者联合编程可以充分发挥出两者各自的优点，大大提高岩石蠕变模型参数辨识的工作效率。采用COM组件处理技术，将Matlab开发的计算函数文件做成组件，进行VB、VC等开发工具直接调用，让算法脱离Matlab环境，则更具有灵活性和可移植性[13-17]。本文针对岩石流变学蠕变模型参数辨识问题，以最常用的伯格斯模型为实例，来实现岩石蠕变模型参数可视化软件的设计开发。为相关岩石蠕变参数辨识问题提供可供参考的快速确定方法。

1　蠕变模型介绍

伯格斯模型(M-K模型)主要由Kelvin体、Maxwell体串联而成，可以描述经过瞬时弹性变形后进入衰减蠕变阶段，随着时间的增加应变持续增加，但蠕变速率呈非线性降低趋势，最终稳定于某值即进入稳态蠕变阶段，此时的蠕变速率基本趋于不变，应变与时间呈线性递增关系(如图1所示)。具体轴向应变公式(1)，其中σ1为轴向应力，σ2为环向应力，而K、G0、G1、η1、η2为待辨识参数。

(1)

其中K、G0大小与瞬时加卸载相关，即蠕变曲线与变形轴的截距或分级加载时的变形突变值；G1、η1的大小影响蠕变曲线的衰减阶段；η2与蠕变稳态阶段的曲线斜率相关，确定η2的大小可以得到实际蠕变曲线第二阶段的斜率，当η2→∞即为三参量蠕变模型。

2　算法程序设计

图1所示，常规三轴蠕变试验数据拟合时首先取曲线上稳态蠕变阶段系列数据拟合渐近线，从而得到斜率及截距。根据公式(1)可得

(2)

(3)

因此可求的η2。其中k1为拟合渐近线斜率，b1为拟合渐近线截距。

其次，取衰减蠕变阶段的系列数据，获得蠕变试验曲线与渐近线延长线之间的垂直距离，并取自然对数，得到新的拟合直线。根据公式(1)可得则有

(4)

因此可求的G1、η1。其中k2为第二次拟合渐近线斜率，b2为第二次拟合渐近线截距。

最后，由曲线瞬时蠕变阶段的数据点，通过

(5)

求得K，再通过公式(3)式有，因此可求的G0。

(6)

通过Matlab软件编译该算法的M函数，并编译完成nnModelMK.dll的COM组件。为了使COM组件脱离开发环境，从“Component”菜单中选择“PackageComponent”选项，将创建包含如表1所示文件的自解压可执行程序。在目标计计算机上运行安装自解压文件nnModelMK.exe进行注册即可调用算法程序[18]。

M文件算法源程序关键代码：

functionModelMK(q1,q2,vv,e0,e1s,e1e,e2s,e2e)

fid=fopen(′D: esult.txt′, ′r′);

a=fscanf(fid, ′ %f %f′,[2inf]);

a=a′

fclose(fid);

%读取数据赋值到a矩阵(n行2列)

x1=a(e2s:e2e,1);

y1=a(e2s:e2e,2);

yy1=polyfit(x1,y1,1);

% 线性拟合蠕变稳态阶段，参加公式(2)、公式(3)

n2=(q1-q2)/(3*yy1(1));

%计算剪切黏性系数η2

x2=a(e1s:e1e,1);

y2=log(yy1(1)*x2+yy1(2)-a(e1s:e1e,2));

yy2=polyfit(x2,y2,1);

% 线性拟合蠕变衰减阶段，参加公式(4)

g1=(q1-q2)/(3*exp(yy2(2)));

%计算弹性剪切模量G1

n1=-g1/yy2(1);

%计算剪切黏性系数η1

kk=(q1-q2)/(3*(1-2*vv)*e0);

%计算弹性体积模量K

g0=1/((9*kk*yy1(2)-(q1+2*q2))/(3*kk*(q1-q2))-1/g1);

%计算弹性剪切模量G0

ModelMK函数参数有8个，依次为：轴向应力、围压、泊松比、瞬时应变、衰减蠕变数据起始行号、衰减蠕变数据终止行号、稳态蠕变数据起始行号及稳态蠕变数据终止行号。通过提取文件数据，首先拟合图1所示稳态蠕变阶段，即得到公式(2)、公式(3)，然后依据衰减蠕变阶段数据得到公式(4)，再选取瞬时蠕变阶段数据确定公式(5)、公式(6)中待确定参数。最后函数返回参数依次为：弹性体积模量K、弹性剪切模量G0、弹性剪切模量G1、剪切粘性系数η1、剪切粘性系数η2。并拟合蠕变模型方程具体形式。

表1　 nnModelMK文件说明

3　实例应用

程序在VB6.0、Matlab6.5环境下开发，将在Matlab 中生成的算法组件引入VB 中，调用VB 系统，打开“P ro ject”——“Reference”对话框，选中nnModelMK 10 Type Library即可。在开发软件功能模块创建和定义新类即可调用matlab算法完成计算。岩石蠕变模型参数辨识软件由系统管理(登录信息加载)，力学参数计算和帮助三部分组成，如图2(b)所示。操作过程将实测的txt格式数据文件拷贝到D: esult.txt，在图2所示界面填写轴向应力、围压、泊松比、瞬时应变、衰减蠕变数据起始行号、衰减蠕变数据终止行号、稳态蠕变数据起始行号及稳态蠕变数据终止行号等数据，点击计算按钮，即可得到拟合对比图形及公式(1)各待辨识参数值。点击保存计算结果按钮可将拟合图片默认保存到c盘根目录下，图片文件命名规则采用当前时间命名，同时使用Access数据库技术存储数据，便于重复使用。

以我国云应地区拟建盐穴储库提取盐样进行常规三轴蠕变试验数据为例，其中轴压为σ1=46.2mPa，围压为σ2=22.1mPa，泊松比0.34，对试验数据采用M-K模型拟合，并用该软件计算各蠕变参数见表2，直接得到弹性体积模量K、弹性剪切模量G0、弹性剪切模量G1、剪切粘性系数η1、剪切粘性系数η2。同时给出拟合对比图如图3所示，由图3可知该软件计算结果同样与实验数据很好的吻合。

表2　本构方程蠕变参数表

4　结论

1)将伯格斯蠕变模型参数辨识方法结合Matlab与VB二者混合编程进行设计，以盐岩蠕变实验为算例进行验证，结果表明可以实现岩石蠕变模型的高效快速计算及可视化操作，避免反复的验算过程，提高了计算效率。

2)借助Matlab的图像显示及数值计算能力，为相关岩石力学参数计算问题提供了可供借鉴的快速处理方法。

[1]董志宏,丁秀丽.地下洞室软岩蠕变参数反分析[J].矿山压力与顶板管理,2005,3:60-62.

[2]孙 钧,黄 伟.岩石力学参数弹塑性反演问题的优化方法[J].岩石力学与工程学报,1992，6(3):221-229.

[3]蒋昱州,张明鸣,李良权.岩石非线性黏弹塑性蠕变模型研究及其参数识别[J].岩石力学与工程学报，2008,27(4):832-839.

[4]林元华,曾德智，施太和,等.岩盐层蠕变规律的反演方法研究[J].石油学报，2005,26(5):111-114.

[5]李青麒.软岩蠕变参数的曲线拟合计算方法[[J].岩石力学与工程学报，1998,17(5):559-564.

[6]陈炳瑞,冯夏庭,丁秀丽.基于模式搜索的岩石流变模型参数识别[J].石力学与工程学报，2005,24(2):207-211.

[7]袁鸿鹄,李云鹏,唐明明,等.山岭隧道围岩流变参数识别及应用研究[J].代隧道技术，2009,46(5):61-65.

[8]张华宾,王芝银,赵艳杰，等.盐岩全过程蠕变试验及模型参数辨识[J].石油学报，2012，33(5):904-908.

[9]沈振中,徐志英.三峡大坝地基花岗岩蠕变试验研究[J].河海大学学报,1997, 25(2):1-7.

[10]何 峰,王来贵,于永江,等.岩石试件非线性蠕变模型及其参数确定[J].辽宁工程技术大学学报,2005,24(2):181-183.

[11]王芝银,李云鹏.岩体流变理论及其数值模拟[M].北京：科学出版社,2007.

[12]阎 岩,王思敬,王恩志.基于西原模型的变参数蠕变方程[J].岩土力学,2010, 31(10):3025-3035.

[13]孟力力,杨其长.VB调用Matlab的COM组件实现二者混合编程[J].电脑开发与应用,2008, 21(6):24-26.

[14]郭庆武, 张 湜.用VB和MATLAB混合开发软测量模型生成系统[J].计量与测试技术,2004(1):25-27.

[15]隗燕琳，陈进明.基于COM的Matlab混合编程方法[J].计算机与数字工程,2013,41(8):1388-1390.

[16]张文军,万 宇.基于COM的Matlab混合编程技术常见问题分析[J].计算机与现代化,2011(4):153-158.

[17]张 杰，张蕊蕊，胡卜元.煤焦SEM图像的表面孔洞分形维数的Matlab实现[J].河北工程大学学报:自然科学版，2007, 24(2):40-44.

[18]Matlab COM Builder user′s guide[Z].The Mathworks Inc, 2002.

(责任编辑李军)

Parameter identification program design of rock creep model based on Matlab

ZHANG Huabin,MIN Xiaodong,ZHANG Qingqing

(Liaoning Technical University School of Mechanical and Engineering . Liaoning , Fuxin, 123000)

The reliable and efficient identify method of rock creep model was discussed, in order to reduce the onerous calculation of parameter identification of rock creep model . Take M-K creep model for example ,Algorithm was writed use Matlab and then it’s spun off the Matlab by mean of COM.The object-oriented visualization software was writed based on VB. The results suggest that the software used on calculation of parameter identification of rock creep model show up more efficient and worked well , which provide fast process method that can be used for reference in parameter calculations of the others rock mechanics.

M-K creep model; Algorithm design; Object-oriented; Visualization;

2016-05-17

国家自然科学基金青年项目(51504124)；辽宁省教育厅科学研究一般项目(L2014142)；国家自然科学基金青年项目(51404130)

张华宾(1983-)，男，黑龙江呼玛县人，讲师，博士，主要从事岩石力学、计算力学方面的研究工作。

1673-9469(2016)03-0029-04

10.3969/j.issn.1673-9469.2016.03.006

TG333.17

A