吳佳寧 邢同振 宋義敏
*(北方工業(yè)大學(xué)土木工程學(xué)院,北京100144)
?(北京理工大學(xué)宇航學(xué)院,北京100081)
在人類生產(chǎn)和實踐過程中,巖石是重要的生產(chǎn)材料,巖石材料的力學(xué)參數(shù)作為采礦、巖土、地下空間等工程在設(shè)計、施工和運行維護過程中的基礎(chǔ)數(shù)據(jù),對于整個工程和人員的安全保障影響重大。因此,對于巖石的邊界條件和力學(xué)參數(shù)的研究具有重要的理論和實踐意義。
對于巖石的邊界條件和力學(xué)參數(shù)的研究,專家學(xué)者取得了一些有意義的成果。袁風波[1]基于地應(yīng)力實測數(shù)據(jù)和數(shù)值模擬手段,使用神經(jīng)網(wǎng)絡(luò)算法和遺傳算法進行優(yōu)化,反演了巖體地應(yīng)力場;陳亮[2]基于有限元基本原理,建立結(jié)構(gòu)的數(shù)值模型,反演了桁架、梁和方板的力學(xué)邊界條件;邱道宏等[3]基于主應(yīng)力測量數(shù)據(jù),建立數(shù)值計算模型,利用神經(jīng)網(wǎng)絡(luò)反演了初始地應(yīng)力場;黃亞哲[4]將計算得到的邊界條件和力學(xué)參數(shù)代入到有限元數(shù)值模型中,通過小波神經(jīng)網(wǎng)絡(luò)優(yōu)化求解的方法反演了初始地應(yīng)力場;付玉華等[5]基于實測的離散點應(yīng)力值,使用有限元法和有限差分法優(yōu)化目標函數(shù),反演了構(gòu)造應(yīng)力場邊界力。通過以上研究成果可以看出,反演方法是獲取巖石邊界條件和力學(xué)參數(shù)的重要手段。
本文基于有限元方法推導(dǎo)巖石參數(shù)反演方程組,通過數(shù)字散斑相關(guān)方法觀測巖石單軸壓縮實驗的位移場,將位移場作為已知量,反演巖石單軸壓縮實驗的邊界條件和巖石試件的力學(xué)參數(shù)。
單元基本方程為[6]
式(1)中,k為單元剛度矩陣,q為節(jié)點位移矩陣,F(xiàn)為節(jié)點外載荷矩陣,單元剛度矩陣k展開形式見式(2)。i,j,m為三角形單元的3個節(jié)點。
以兩個3節(jié)點三角形單元為例,推導(dǎo)參數(shù)反演方程組。圖1是一個厚度為t的正方形,劃分為e1和e2兩個3節(jié)點三角形單元,共4個節(jié)點,節(jié)點順序按逆時針排列。
圖1 單元示意圖
單元e1和e2的單元剛度矩陣ke1和ke2表示為
整體剛度矩陣是各單元剛度矩陣之和,即
每個子矩陣都可以由式(6)計算得到
式中,A為單元面積,t為厚度,μ為泊松比,E為彈性模量,s=i,j,m,r=i,j,m。將式(6)和式(5)代入基本方程,得到參數(shù)反演方程組T為已知參數(shù)矩陣
式中,h代表與i在同一單元的所有節(jié)點號。C為未知彈性參數(shù)矩陣,F(xiàn)為節(jié)點外載荷矩陣
將兩單元拓展到多單元,假設(shè)節(jié)點數(shù)為n,則未知彈性參數(shù)矩陣C不變,已知參數(shù)矩陣T和節(jié)點外載荷F分別為
式中,i=1,2,···,n。
實驗選取紅砂巖材料,制成長方體試件,長和寬均為50 mm,高為100 mm,通過精細打磨的方式使試件各個端面平整,選擇一個長50 mm,高100 mm的平面制作人工散斑場。
實驗系統(tǒng)分為加載系統(tǒng)和數(shù)字散斑相關(guān)方法圖像采集系統(tǒng)。加載系統(tǒng)使用RLJW-2000型液壓伺服控制試驗機,對試件施加單向壓縮載荷,使用位移加載方式,加載速度大小為0.1 mm/min。數(shù)字散斑相關(guān)方法圖像采集系統(tǒng)由CCD相機、光源和計算機三部分組成,圖像采集速率為每秒2幀,圖像分辨率為1600×1200像素,物面分辨率為每像素0.088 mm。
實驗開始前,在試件上下端面和上下壓板上均勻涂抹一層潤滑劑,減小摩擦力的作用效果,調(diào)試試驗機壓頭與試件上端面稍微接觸,調(diào)整CCD相機的位置、焦距和光圈,使成像效果最佳;實驗開始,試驗機對試件進行單軸壓縮加載,試件上端保持固定,下端產(chǎn)生向上的位移;同時,圖像采集系統(tǒng)開始采集散斑圖像,直到試件發(fā)生破壞,停止加載和數(shù)據(jù)采集;實驗結(jié)束后,對采集到的數(shù)據(jù)進行計算分析。
使用數(shù)字散斑相關(guān)方法計算加載過程中試件的變形場,單軸壓縮的應(yīng)力與應(yīng)變曲線如圖2所示,將加載初始時刻采集的散斑圖像作為參考圖像,在應(yīng)力與應(yīng)變曲線上按照等應(yīng)變間隔選取標識點1~5,計算紅砂巖加載過程中5個標識點對應(yīng)時刻的位移場,標識點1~5對應(yīng)的采集時間、載荷、應(yīng)力和應(yīng)變量值見表1。
圖2 單軸應(yīng)力與應(yīng)變曲線
表1 標識點1~5參數(shù)值
以標識點5時刻為例,試件位移場云圖如圖3所示。隨著加載進行,在水平方向上,試件由中間向左右兩側(cè)移動,位移場等值線出現(xiàn)“X狀”分布,這是由試件上下端面與壓機之間存在摩擦力導(dǎo)致的;在豎直方向上,隨著加載進行,等值線呈水平分層分布,在加載過程中,試件上壓板保持固定,下壓板施加向上的位移,豎直方向的等值線量值由上至下逐漸增大。
圖3 標識點1~5位移場云圖
根據(jù)試件尺寸,將試件劃分為數(shù)量一定的3節(jié)點三角形單元,以3×5網(wǎng)格為例,單元劃分網(wǎng)格示意圖如圖4所示。在水平方向上,將散斑場平均分成3等份,豎直方向上平均分成5等份,共劃分了30個單元,24個節(jié)點。以實驗計算得到的位移場為已知量,將單元內(nèi)的已知點位移值代入到有限元方程中,優(yōu)化求解各個節(jié)點位移。
圖4 網(wǎng)格劃分示意圖(3×5)
節(jié)點位移計算方法如下:假設(shè)在任意一個節(jié)點為i,j,m的三角形單元中,數(shù)字散斑相關(guān)方法計算已知點的數(shù)量為n,位移為u1,v1,···,un,vn,通過已知點坐標計算的形函數(shù)為Ni1···Nin,Nj1···Njn,Nm1···Nmn,三角形單元3個節(jié)點的坐標在劃分網(wǎng)格時已經(jīng)確定,將以上數(shù)據(jù)代入到式(15)中,采用最小二乘法進行優(yōu)化,計算得到節(jié)點位移值ui,uj,um,vi,vj,vm。
由于一個節(jié)點可能存在于多個單元中,因此在不同單元中,同一個節(jié)點計算得到多個位移值,將多個位移值的均值作為該節(jié)點的劃分網(wǎng)格節(jié)點位移值。例如,在標識點5中,圖4中的所有網(wǎng)格節(jié)點位移插值結(jié)果見表2。
表2 節(jié)點位移計算結(jié)果
試件上端面豎直方向位移值為0,根據(jù)實驗加載中的力和位移關(guān)系,認為外載荷平均分布到試件下端面的各個節(jié)點上。由于單軸壓縮加載位移等值線在水平方向大致呈“X狀”分布,試件上下端面與上下壓板之間存在摩擦力,試驗機和數(shù)字散斑采集系統(tǒng)均不能測量得到邊界條件摩擦力的量值大小,需對試件上下端面水平方向節(jié)點外載荷進行反演。根據(jù)單軸壓縮實驗條件設(shè)置,試件上端面節(jié)點承受水平方向的節(jié)點外力,試件下端面節(jié)點承受水平和豎直方向的節(jié)點外力,其余節(jié)點不承受外力。
在參數(shù)反演方程組式(7)中,假設(shè)巖石的彈性模量E和泊松比μ未知,散斑場劃分3節(jié)點三角形單元網(wǎng)格的節(jié)點數(shù)共為n個,其中,上/下邊界節(jié)點數(shù)為n1個,參數(shù)反演方程組是由2n個方程和2n1+2個未知量組成的超定方程組,使用最小二乘法進行優(yōu)化,計算得出紅砂巖單軸壓縮的彈性模量E、泊松比μ和上下端面的水平方向外載荷值。其中,彈性模量E的初始值為8 GPa,上限為100 GPa,下限為1 GPa,泊松比μ的初始值為0.25,上限為0.5,下限為0,停止迭代條件為未知數(shù)的變化差值小于10-8。優(yōu)化計算部分程序如圖5所示。
圖5 優(yōu)化計算部分程序
將散斑場劃分為5×10的3節(jié)點三角形網(wǎng)格單元,即在水平方向上將散斑場平均分為5等份,共10個單元,在豎直方向上將散斑場平均分為10等份,共10個單元。分別計算5個加載時刻的彈性力學(xué)參數(shù)和巖石單軸壓縮實驗上下端面邊界條件,計算結(jié)果見表3,其中,F(xiàn)X1~FX6為試件上端面水平方向的節(jié)點外力,F(xiàn)X7~FX12為試件下端面水平方向的節(jié)點外力,F(xiàn)Y為已知的試件下端面豎直方向的均布節(jié)點外力。
表3 標識點1~5計算結(jié)果
表3 的計算結(jié)果如圖6和圖7所示。通過圖6可以看出,水平方向的相鄰節(jié)點外載荷方向可能不一致,量值大小也不完全統(tǒng)一,計算結(jié)果與節(jié)點位移量值大小和方向密切相關(guān),由于巖石是非均質(zhì)的材料[7-8],某個位置的微小裂隙或者結(jié)構(gòu)構(gòu)造弱化都會導(dǎo)致其在加載過程中出現(xiàn)較大裂隙,位移量值比該點周圍點的位移量值大,進而導(dǎo)致力的特征差異。從圖7中可以看出,隨著加載進行,彈性模量為18 GPa左右,加載后期趨于穩(wěn)定,泊松比在0.14~0.34之間,隨加載進行呈現(xiàn)逐漸增大趨勢。
圖6 邊界條件計算結(jié)果
本文通過有限單元法和數(shù)字散斑相關(guān)方法對巖石單軸壓縮實驗的邊界條件和彈性參數(shù)進行反演??梢缘玫饺缦陆Y(jié)論:
(1)利用本文所發(fā)展的方法,可以實現(xiàn)材料的邊界條件和彈性參數(shù)反演。
(2)隨著加載進行,反演得到的彈性模量為18 GPa左右,并且加載后期趨于穩(wěn)定;泊松比在0.14~0.34之間,隨加載進行呈現(xiàn)逐漸增大趨勢。
(3)在水平方向上,相鄰節(jié)點外載荷的方向和量值大小具有一定差異,產(chǎn)生此種結(jié)果與巖石材料的非均勻性及各個節(jié)點的位移量值差異較大有關(guān)。