吳佳寧 邢同振 宋義敏

*(北方工業(yè)大學(xué)土木工程學(xué)院，北京100144)

?(北京理工大學(xué)宇航學(xué)院，北京100081)

在人類生產(chǎn)和實(shí)踐過程中，巖石是重要的生產(chǎn)材料，巖石材料的力學(xué)參數(shù)作為采礦、巖土、地下空間等工程在設(shè)計(jì)、施工和運(yùn)行維護(hù)過程中的基礎(chǔ)數(shù)據(jù)，對(duì)于整個(gè)工程和人員的安全保障影響重大。因此，對(duì)于巖石的邊界條件和力學(xué)參數(shù)的研究具有重要的理論和實(shí)踐意義。

對(duì)于巖石的邊界條件和力學(xué)參數(shù)的研究，專家學(xué)者取得了一些有意義的成果。袁風(fēng)波[1]基于地應(yīng)力實(shí)測(cè)數(shù)據(jù)和數(shù)值模擬手段，使用神經(jīng)網(wǎng)絡(luò)算法和遺傳算法進(jìn)行優(yōu)化，反演了巖體地應(yīng)力場(chǎng)；陳亮[2]基于有限元基本原理，建立結(jié)構(gòu)的數(shù)值模型，反演了桁架、梁和方板的力學(xué)邊界條件；邱道宏等[3]基于主應(yīng)力測(cè)量數(shù)據(jù)，建立數(shù)值計(jì)算模型，利用神經(jīng)網(wǎng)絡(luò)反演了初始地應(yīng)力場(chǎng)；黃亞哲[4]將計(jì)算得到的邊界條件和力學(xué)參數(shù)代入到有限元數(shù)值模型中，通過小波神經(jīng)網(wǎng)絡(luò)優(yōu)化求解的方法反演了初始地應(yīng)力場(chǎng)；付玉華等[5]基于實(shí)測(cè)的離散點(diǎn)應(yīng)力值，使用有限元法和有限差分法優(yōu)化目標(biāo)函數(shù)，反演了構(gòu)造應(yīng)力場(chǎng)邊界力。通過以上研究成果可以看出，反演方法是獲取巖石邊界條件和力學(xué)參數(shù)的重要手段。

本文基于有限元方法推導(dǎo)巖石參數(shù)反演方程組，通過數(shù)字散斑相關(guān)方法觀測(cè)巖石單軸壓縮實(shí)驗(yàn)的位移場(chǎng)，將位移場(chǎng)作為已知量，反演巖石單軸壓縮實(shí)驗(yàn)的邊界條件和巖石試件的力學(xué)參數(shù)。

1 反演方程組推導(dǎo)

單元基本方程為[6]

式(1)中，k為單元?jiǎng)偠染仃嚕琿為節(jié)點(diǎn)位移矩陣，F(xiàn)為節(jié)點(diǎn)外載荷矩陣，單元?jiǎng)偠染仃噆展開形式見式(2)。i，j，m為三角形單元的3個(gè)節(jié)點(diǎn)。

以兩個(gè)3節(jié)點(diǎn)三角形單元為例，推導(dǎo)參數(shù)反演方程組。圖1是一個(gè)厚度為t的正方形，劃分為e1和e2兩個(gè)3節(jié)點(diǎn)三角形單元，共4個(gè)節(jié)點(diǎn)，節(jié)點(diǎn)順序按逆時(shí)針排列。

圖1 單元示意圖

單元e1和e2的單元?jiǎng)偠染仃噆e1和ke2表示為

整體剛度矩陣是各單元?jiǎng)偠染仃囍停?/p>

每個(gè)子矩陣都可以由式(6)計(jì)算得到

式中，A為單元面積，t為厚度，μ為泊松比，E為彈性模量，s=i，j，m，r=i，j，m。將式(6)和式(5)代入基本方程，得到參數(shù)反演方程組T為已知參數(shù)矩陣

式中，h代表與i在同一單元的所有節(jié)點(diǎn)號(hào)。C為未知彈性參數(shù)矩陣，F(xiàn)為節(jié)點(diǎn)外載荷矩陣

將兩單元拓展到多單元，假設(shè)節(jié)點(diǎn)數(shù)為n，則未知彈性參數(shù)矩陣C不變，已知參數(shù)矩陣T和節(jié)點(diǎn)外載荷F分別為

式中，i=1，2，···，n。

2 巖石單軸壓縮實(shí)驗(yàn)

2.1 實(shí)驗(yàn)準(zhǔn)備及過程

實(shí)驗(yàn)選取紅砂巖材料，制成長方體試件，長和寬均為50 mm，高為100 mm，通過精細(xì)打磨的方式使試件各個(gè)端面平整，選擇一個(gè)長50 mm，高100 mm的平面制作人工散斑場(chǎng)。

實(shí)驗(yàn)系統(tǒng)分為加載系統(tǒng)和數(shù)字散斑相關(guān)方法圖像采集系統(tǒng)。加載系統(tǒng)使用RLJW-2000型液壓伺服控制試驗(yàn)機(jī)，對(duì)試件施加單向壓縮載荷，使用位移加載方式，加載速度大小為0.1 mm/min。數(shù)字散斑相關(guān)方法圖像采集系統(tǒng)由CCD相機(jī)、光源和計(jì)算機(jī)三部分組成，圖像采集速率為每秒2幀，圖像分辨率為1600×1200像素，物面分辨率為每像素0.088 mm。

實(shí)驗(yàn)開始前，在試件上下端面和上下壓板上均勻涂抹一層潤滑劑，減小摩擦力的作用效果，調(diào)試試驗(yàn)機(jī)壓頭與試件上端面稍微接觸，調(diào)整CCD相機(jī)的位置、焦距和光圈，使成像效果最佳；實(shí)驗(yàn)開始，試驗(yàn)機(jī)對(duì)試件進(jìn)行單軸壓縮加載，試件上端保持固定，下端產(chǎn)生向上的位移；同時(shí)，圖像采集系統(tǒng)開始采集散斑圖像，直到試件發(fā)生破壞，停止加載和數(shù)據(jù)采集；實(shí)驗(yàn)結(jié)束后，對(duì)采集到的數(shù)據(jù)進(jìn)行計(jì)算分析。

2.2 實(shí)驗(yàn)分析及結(jié)果

使用數(shù)字散斑相關(guān)方法計(jì)算加載過程中試件的變形場(chǎng)，單軸壓縮的應(yīng)力與應(yīng)變曲線如圖2所示，將加載初始時(shí)刻采集的散斑圖像作為參考圖像，在應(yīng)力與應(yīng)變曲線上按照等應(yīng)變間隔選取標(biāo)識(shí)點(diǎn)1～5，計(jì)算紅砂巖加載過程中5個(gè)標(biāo)識(shí)點(diǎn)對(duì)應(yīng)時(shí)刻的位移場(chǎng)，標(biāo)識(shí)點(diǎn)1～5對(duì)應(yīng)的采集時(shí)間、載荷、應(yīng)力和應(yīng)變量值見表1。

圖2 單軸應(yīng)力與應(yīng)變曲線

表1 標(biāo)識(shí)點(diǎn)1～5參數(shù)值

以標(biāo)識(shí)點(diǎn)5時(shí)刻為例，試件位移場(chǎng)云圖如圖3所示。隨著加載進(jìn)行，在水平方向上，試件由中間向左右兩側(cè)移動(dòng)，位移場(chǎng)等值線出現(xiàn)“X狀”分布，這是由試件上下端面與壓機(jī)之間存在摩擦力導(dǎo)致的；在豎直方向上，隨著加載進(jìn)行，等值線呈水平分層分布，在加載過程中，試件上壓板保持固定，下壓板施加向上的位移，豎直方向的等值線量值由上至下逐漸增大。

圖3 標(biāo)識(shí)點(diǎn)1～5位移場(chǎng)云圖

根據(jù)試件尺寸，將試件劃分為數(shù)量一定的3節(jié)點(diǎn)三角形單元，以3×5網(wǎng)格為例，單元?jiǎng)澐志W(wǎng)格示意圖如圖4所示。在水平方向上，將散斑場(chǎng)平均分成3等份，豎直方向上平均分成5等份，共劃分了30個(gè)單元，24個(gè)節(jié)點(diǎn)。以實(shí)驗(yàn)計(jì)算得到的位移場(chǎng)為已知量，將單元內(nèi)的已知點(diǎn)位移值代入到有限元方程中，優(yōu)化求解各個(gè)節(jié)點(diǎn)位移。

圖4 網(wǎng)格劃分示意圖(3×5)

節(jié)點(diǎn)位移計(jì)算方法如下：假設(shè)在任意一個(gè)節(jié)點(diǎn)為i，j，m的三角形單元中，數(shù)字散斑相關(guān)方法計(jì)算已知點(diǎn)的數(shù)量為n，位移為u1，v1，···，un，vn，通過已知點(diǎn)坐標(biāo)計(jì)算的形函數(shù)為Ni1···Nin，Nj1···Njn，Nm1···Nmn，三角形單元3個(gè)節(jié)點(diǎn)的坐標(biāo)在劃分網(wǎng)格時(shí)已經(jīng)確定，將以上數(shù)據(jù)代入到式(15)中，采用最小二乘法進(jìn)行優(yōu)化，計(jì)算得到節(jié)點(diǎn)位移值ui，uj，um，vi，vj，vm。

由于一個(gè)節(jié)點(diǎn)可能存在于多個(gè)單元中，因此在不同單元中，同一個(gè)節(jié)點(diǎn)計(jì)算得到多個(gè)位移值，將多個(gè)位移值的均值作為該節(jié)點(diǎn)的劃分網(wǎng)格節(jié)點(diǎn)位移值。例如，在標(biāo)識(shí)點(diǎn)5中，圖4中的所有網(wǎng)格節(jié)點(diǎn)位移插值結(jié)果見表2。

表2 節(jié)點(diǎn)位移計(jì)算結(jié)果

試件上端面豎直方向位移值為0，根據(jù)實(shí)驗(yàn)加載中的力和位移關(guān)系，認(rèn)為外載荷平均分布到試件下端面的各個(gè)節(jié)點(diǎn)上。由于單軸壓縮加載位移等值線在水平方向大致呈“X狀”分布，試件上下端面與上下壓板之間存在摩擦力，試驗(yàn)機(jī)和數(shù)字散斑采集系統(tǒng)均不能測(cè)量得到邊界條件摩擦力的量值大小，需對(duì)試件上下端面水平方向節(jié)點(diǎn)外載荷進(jìn)行反演。根據(jù)單軸壓縮實(shí)驗(yàn)條件設(shè)置，試件上端面節(jié)點(diǎn)承受水平方向的節(jié)點(diǎn)外力，試件下端面節(jié)點(diǎn)承受水平和豎直方向的節(jié)點(diǎn)外力，其余節(jié)點(diǎn)不承受外力。

3 邊界條件和彈性參數(shù)反演結(jié)果

在參數(shù)反演方程組式(7)中，假設(shè)巖石的彈性模量E和泊松比μ未知，散斑場(chǎng)劃分3節(jié)點(diǎn)三角形單元網(wǎng)格的節(jié)點(diǎn)數(shù)共為n個(gè)，其中，上/下邊界節(jié)點(diǎn)數(shù)為n1個(gè)，參數(shù)反演方程組是由2n個(gè)方程和2n1+2個(gè)未知量組成的超定方程組，使用最小二乘法進(jìn)行優(yōu)化，計(jì)算得出紅砂巖單軸壓縮的彈性模量E、泊松比μ和上下端面的水平方向外載荷值。其中，彈性模量E的初始值為8 GPa，上限為100 GPa，下限為1 GPa，泊松比μ的初始值為0.25，上限為0.5，下限為0，停止迭代條件為未知數(shù)的變化差值小于10-8。優(yōu)化計(jì)算部分程序如圖5所示。

圖5 優(yōu)化計(jì)算部分程序

將散斑場(chǎng)劃分為5×10的3節(jié)點(diǎn)三角形網(wǎng)格單元，即在水平方向上將散斑場(chǎng)平均分為5等份，共10個(gè)單元，在豎直方向上將散斑場(chǎng)平均分為10等份，共10個(gè)單元。分別計(jì)算5個(gè)加載時(shí)刻的彈性力學(xué)參數(shù)和巖石單軸壓縮實(shí)驗(yàn)上下端面邊界條件，計(jì)算結(jié)果見表3，其中，F(xiàn)X1～FX6為試件上端面水平方向的節(jié)點(diǎn)外力，F(xiàn)X7～FX12為試件下端面水平方向的節(jié)點(diǎn)外力，F(xiàn)Y為已知的試件下端面豎直方向的均布節(jié)點(diǎn)外力。

表3 標(biāo)識(shí)點(diǎn)1～5計(jì)算結(jié)果

表3 的計(jì)算結(jié)果如圖6和圖7所示。通過圖6可以看出，水平方向的相鄰節(jié)點(diǎn)外載荷方向可能不一致，量值大小也不完全統(tǒng)一，計(jì)算結(jié)果與節(jié)點(diǎn)位移量值大小和方向密切相關(guān)，由于巖石是非均質(zhì)的材料[7-8]，某個(gè)位置的微小裂隙或者結(jié)構(gòu)構(gòu)造弱化都會(huì)導(dǎo)致其在加載過程中出現(xiàn)較大裂隙，位移量值比該點(diǎn)周圍點(diǎn)的位移量值大，進(jìn)而導(dǎo)致力的特征差異。從圖7中可以看出，隨著加載進(jìn)行，彈性模量為18 GPa左右，加載后期趨于穩(wěn)定，泊松比在0.14～0.34之間，隨加載進(jìn)行呈現(xiàn)逐漸增大趨勢(shì)。

圖6 邊界條件計(jì)算結(jié)果

4 結(jié)論

本文通過有限單元法和數(shù)字散斑相關(guān)方法對(duì)巖石單軸壓縮實(shí)驗(yàn)的邊界條件和彈性參數(shù)進(jìn)行反演?？梢缘玫饺缦陆Y(jié)論：

(1)利用本文所發(fā)展的方法，可以實(shí)現(xiàn)材料的邊界條件和彈性參數(shù)反演。

(2)隨著加載進(jìn)行，反演得到的彈性模量為18 GPa左右，并且加載后期趨于穩(wěn)定；泊松比在0.14～0.34之間，隨加載進(jìn)行呈現(xiàn)逐漸增大趨勢(shì)。

(3)在水平方向上，相鄰節(jié)點(diǎn)外載荷的方向和量值大小具有一定差異，產(chǎn)生此種結(jié)果與巖石材料的非均勻性及各個(gè)節(jié)點(diǎn)的位移量值差異較大有關(guān)。