基于Matlab的巖石蠕變模型參數(shù)辨識算法設(shè)計

張華賓,閔小東,張頃頃

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

?

基于Matlab的巖石蠕變模型參數(shù)辨識算法設(shè)計

張華賓,閔小東,張頃頃

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

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

蠕變模型；算法設(shè)計；面向?qū)ο?；可視?/p>

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

1　蠕變模型介紹

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

(1)

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

2　算法程序設(shè)計

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

(2)

(3)

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

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

(4)

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

最后，由曲線瞬時蠕變階段的數(shù)據(jù)點，通過

(5)

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

(6)

通過Matlab軟件編譯該算法的M函數(shù)，并編譯完成nnModelMK.dll的COM組件。為了使COM組件脫離開發(fā)環(huán)境，從“Component”菜單中選擇“PackageComponent”選項，將創(chuàng)建包含如表1所示文件的自解壓可執(zhí)行程序。在目標(biāo)計計算機(jī)上運(yùn)行安裝自解壓文件nnModelMK.exe進(jìn)行注冊即可調(diào)用算法程序[18]。

M文件算法源程序關(guān)鍵代碼：

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);

%讀取數(shù)據(jù)賦值到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));

%計算剪切黏性系數(shù)η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);

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

表1　 nnModelMK文件說明

3　實例應(yīng)用

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

以我國云應(yīng)地區(qū)擬建鹽穴儲庫提取鹽樣進(jìn)行常規(guī)三軸蠕變試驗數(shù)據(jù)為例，其中軸壓為σ1=46.2mPa，圍壓為σ2=22.1mPa，泊松比0.34，對試驗數(shù)據(jù)采用M-K模型擬合，并用該軟件計算各蠕變參數(shù)見表2，直接得到彈性體積模量K、彈性剪切模量G0、彈性剪切模量G1、剪切粘性系數(shù)η1、剪切粘性系數(shù)η2。同時給出擬合對比圖如圖3所示，由圖3可知該軟件計算結(jié)果同樣與實驗數(shù)據(jù)很好的吻合。

表2　本構(gòu)方程蠕變參數(shù)表

4　結(jié)論

1)將伯格斯蠕變模型參數(shù)辨識方法結(jié)合Matlab與VB二者混合編程進(jìn)行設(shè)計，以鹽巖蠕變實驗為算例進(jìn)行驗證，結(jié)果表明可以實現(xiàn)巖石蠕變模型的高效快速計算及可視化操作，避免反復(fù)的驗算過程，提高了計算效率。

2)借助Matlab的圖像顯示及數(shù)值計算能力，為相關(guān)巖石力學(xué)參數(shù)計算問題提供了可供借鑒的快速處理方法。

[1]董志宏,丁秀麗.地下洞室軟巖蠕變參數(shù)反分析[J].礦山壓力與頂板管理,2005,3:60-62.

[2]孫 鈞,黃 偉.巖石力學(xué)參數(shù)彈塑性反演問題的優(yōu)化方法[J].巖石力學(xué)與工程學(xué)報,1992，6(3):221-229.

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

[4]林元華,曾德智，施太和,等.巖鹽層蠕變規(guī)律的反演方法研究[J].石油學(xué)報，2005,26(5):111-114.

[5]李青麒.軟巖蠕變參數(shù)的曲線擬合計算方法[[J].巖石力學(xué)與工程學(xué)報，1998,17(5):559-564.

[6]陳炳瑞,馮夏庭,丁秀麗.基于模式搜索的巖石流變模型參數(shù)識別[J].石力學(xué)與工程學(xué)報，2005,24(2):207-211.

[7]袁鴻鵠,李云鵬,唐明明,等.山嶺隧道圍巖流變參數(shù)識別及應(yīng)用研究[J].代隧道技術(shù)，2009,46(5):61-65.

[8]張華賓,王芝銀,趙艷杰，等.鹽巖全過程蠕變試驗及模型參數(shù)辨識[J].石油學(xué)報，2012，33(5):904-908.

[9]沈振中,徐志英.三峽大壩地基花崗巖蠕變試驗研究[J].河海大學(xué)學(xué)報,1997, 25(2):1-7.

[10]何 峰,王來貴,于永江,等.巖石試件非線性蠕變模型及其參數(shù)確定[J].遼寧工程技術(shù)大學(xué)學(xué)報,2005,24(2):181-183.

[11]王芝銀,李云鵬.巖體流變理論及其數(shù)值模擬[M].北京：科學(xué)出版社,2007.

[12]閻 巖,王思敬,王恩志.基于西原模型的變參數(shù)蠕變方程[J].巖土力學(xué),2010, 31(10):3025-3035.

[13]孟力力,楊其長.VB調(diào)用Matlab的COM組件實現(xiàn)二者混合編程[J].電腦開發(fā)與應(yīng)用,2008, 21(6):24-26.

[14]郭慶武, 張 湜.用VB和MATLAB混合開發(fā)軟測量模型生成系統(tǒng)[J].計量與測試技術(shù),2004(1):25-27.

[15]隗燕琳，陳進(jìn)明.基于COM的Matlab混合編程方法[J].計算機(jī)與數(shù)字工程,2013,41(8):1388-1390.

[16]張文軍,萬 宇.基于COM的Matlab混合編程技術(shù)常見問題分析[J].計算機(jī)與現(xiàn)代化,2011(4):153-158.

[17]張 杰，張蕊蕊，胡卜元.煤焦SEM圖像的表面孔洞分形維數(shù)的Matlab實現(xiàn)[J].河北工程大學(xué)學(xué)報:自然科學(xué)版，2007, 24(2):40-44.

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

(責(zé)任編輯李軍)

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

國家自然科學(xué)基金青年項目(51504124)；遼寧省教育廳科學(xué)研究一般項目(L2014142)；國家自然科學(xué)基金青年項目(51404130)

張華賓(1983-)，男，黑龍江呼瑪縣人，講師，博士，主要從事巖石力學(xué)、計算力學(xué)方面的研究工作。

1673-9469(2016)03-0029-04

10.3969/j.issn.1673-9469.2016.03.006

TG333.17

A