張華賓,閔小東,張頃頃

(遼寧工程技術大學 力學與工程學院，遼寧 阜新 123000)

?

基于Matlab的巖石蠕變模型參數辨識算法設計

張華賓,閔小東,張頃頃

(遼寧工程技術大學 力學與工程學院，遼寧 阜新 123000)

為減少巖石蠕變模型參數辨識的繁冗計算，探討一種可靠、高效的巖石蠕變模型辨識方法。以M-K(伯格斯)蠕變模型為例，給出蠕變模型參數辨識方法，并編寫Matlab算法函數通過COM組件從而剝離出程序，借助VB環(huán)境下編制面向對象的可視化計算軟件。結果表明，該軟件用于巖石蠕變模型參數辨識的計算過程同樣有較好的效果并提高了計算效率，為實現其他巖石力學參數計算問題提供了可供借鑒的快速處理方法。

蠕變模型；算法設計；面向對象；可視化

巖石蠕變力學性質模型大致可以分為以下4類：經驗模型、元件組合模型、基于損傷力學的蠕變本構模型和基于熱力學理論的蠕變本構模型。其中基本元件模型通過組合可以得到很多模型來描述巖土的流變特性，同時蠕變模型的參數辨識方法也很多，根據蠕變試驗數據選擇合適的蠕變模型確定相應的蠕變參數，識別參數可以獲得較吻合試驗數據的擬合曲線，主要的參數辨識方法有回歸分析，最小二乘，曲線分解等方法。董志宏[1]等基于現場模型洞開挖位移檢測結果，運用均勻設計—神經網羅—遺傳算法對軟巖蠕變參數進行了反演分析，得到該區(qū)軟巖的蠕變參數。孫鈞[2]教授等用最小二乘法建立目標函數，通過全局最優(yōu)化，同時辨識了初始地應力三個分量及變形和強度3個參數，從實踐上得出了非線性逆問題的唯一解。蔣昱州[3]提出 BFGS-LSM 算法對巖石非線性黏彈塑性蠕變模型參數進行辨識，結果和試驗曲線吻合較好。林元華[4]建立了巖鹽本構模型辨識的目標函數，提出了以井徑測量信息為基礎反演確定巖鹽層蠕變特性參數的位移反分析法。李青麒[5]介紹一種根據室內蠕變資料確定流變參數和選擇流變模型的方法。陳炳瑞[6]提出模式搜索與非線性最小二乘法的有機結合識別流變模型參數的方法。袁鴻鵠[7]基于流變反演理論，建立油房灣隧道圍巖流變變形特征的流變模型。但無論何種方法得到準確合理的模型參數，都需要反復提取數據計算[8-12]。 然而數學軟件Matlab進行數據處理和數據顯示方面的優(yōu)勢非常顯著，同時考慮VB編寫的程序執(zhí)行效率較高，可以提供交互性很好的界面，因此二者聯合編程可以充分發(fā)揮出兩者各自的優(yōu)點，大大提高巖石蠕變模型參數辨識的工作效率。采用COM組件處理技術，將Matlab開發(fā)的計算函數文件做成組件，進行VB、VC等開發(fā)工具直接調用，讓算法脫離Matlab環(huán)境，則更具有靈活性和可移植性[13-17]。本文針對巖石流變學蠕變模型參數辨識問題，以最常用的伯格斯模型為實例，來實現巖石蠕變模型參數可視化軟件的設計開發(fā)。為相關巖石蠕變參數辨識問題提供可供參考的快速確定方法。

1　蠕變模型介紹

伯格斯模型(M-K模型)主要由Kelvin體、Maxwell體串聯而成，可以描述經過瞬時彈性變形后進入衰減蠕變階段，隨著時間的增加應變持續(xù)增加，但蠕變速率呈非線性降低趨勢，最終穩(wěn)定于某值即進入穩(wěn)態(tài)蠕變階段，此時的蠕變速率基本趨于不變，應變與時間呈線性遞增關系(如圖1所示)。具體軸向應變公式(1)，其中σ1為軸向應力，σ2為環(huán)向應力，而K、G0、G1、η1、η2為待辨識參數。

(1)

其中K、G0大小與瞬時加卸載相關，即蠕變曲線與變形軸的截距或分級加載時的變形突變值；G1、η1的大小影響蠕變曲線的衰減階段；η2與蠕變穩(wěn)態(tài)階段的曲線斜率相關，確定η2的大小可以得到實際蠕變曲線第二階段的斜率，當η2→∞即為三參量蠕變模型。

2　算法程序設計

圖1所示，常規(guī)三軸蠕變試驗數據擬合時首先取曲線上穩(wěn)態(tài)蠕變階段系列數據擬合漸近線，從而得到斜率及截距。根據公式(1)可得

(2)

(3)

因此可求的η2。其中k1為擬合漸近線斜率，b1為擬合漸近線截距。

其次，取衰減蠕變階段的系列數據，獲得蠕變試驗曲線與漸近線延長線之間的垂直距離，并取自然對數，得到新的擬合直線。根據公式(1)可得則有

(4)

因此可求的G1、η1。其中k2為第二次擬合漸近線斜率，b2為第二次擬合漸近線截距。

最后，由曲線瞬時蠕變階段的數據點，通過

(5)

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

(6)

通過Matlab軟件編譯該算法的M函數，并編譯完成nnModelMK.dll的COM組件。為了使COM組件脫離開發(fā)環(huán)境，從“Component”菜單中選擇“PackageComponent”選項，將創(chuàng)建包含如表1所示文件的自解壓可執(zhí)行程序。在目標計計算機上運行安裝自解壓文件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);

% 線性擬合蠕變穩(wěn)態(tài)階段，參加公式(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個，依次為：軸向應力、圍壓、泊松比、瞬時應變、衰減蠕變數據起始行號、衰減蠕變數據終止行號、穩(wěn)態(tài)蠕變數據起始行號及穩(wěn)態(tài)蠕變數據終止行號。通過提取文件數據，首先擬合圖1所示穩(wěn)態(tài)蠕變階段，即得到公式(2)、公式(3)，然后依據衰減蠕變階段數據得到公式(4)，再選取瞬時蠕變階段數據確定公式(5)、公式(6)中待確定參數。最后函數返回參數依次為：彈性體積模量K、彈性剪切模量G0、彈性剪切模量G1、剪切粘性系數η1、剪切粘性系數η2。并擬合蠕變模型方程具體形式。

表1　 nnModelMK文件說明

3　實例應用

程序在VB6.0、Matlab6.5環(huán)境下開發(fā)，將在Matlab 中生成的算法組件引入VB 中，調用VB 系統，打開“P ro ject”——“Reference”對話框，選中nnModelMK 10 Type Library即可。在開發(fā)軟件功能模塊創(chuàng)建和定義新類即可調用matlab算法完成計算。巖石蠕變模型參數辨識軟件由系統管理(登錄信息加載)，力學參數計算和幫助三部分組成，如圖2(b)所示。操作過程將實測的txt格式數據文件拷貝到D: esult.txt，在圖2所示界面填寫軸向應力、圍壓、泊松比、瞬時應變、衰減蠕變數據起始行號、衰減蠕變數據終止行號、穩(wěn)態(tài)蠕變數據起始行號及穩(wěn)態(tài)蠕變數據終止行號等數據，點擊計算按鈕，即可得到擬合對比圖形及公式(1)各待辨識參數值。點擊保存計算結果按鈕可將擬合圖片默認保存到c盤根目錄下，圖片文件命名規(guī)則采用當前時間命名，同時使用Access數據庫技術存儲數據，便于重復使用。

以我國云應地區(qū)擬建鹽穴儲庫提取鹽樣進行常規(guī)三軸蠕變試驗數據為例，其中軸壓為σ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]孫 鈞,黃 偉.巖石力學參數彈塑性反演問題的優(yōu)化方法[J].巖石力學與工程學報,1992，6(3):221-229.

[3]蔣昱州,張明鳴,李良權.巖石非線性黏彈塑性蠕變模型研究及其參數識別[J].巖石力學與工程學報，2008,27(4):832-839.

[4]林元華,曾德智，施太和,等.巖鹽層蠕變規(guī)律的反演方法研究[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].電腦開發(fā)與應用,2008, 21(6):24-26.

[14]郭慶武, 張 湜.用VB和MATLAB混合開發(fā)軟測量模型生成系統[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