基于MATLAB的巖體結構面三維模擬與表征單元體研究

李鈺潔

摘 要：巖體中的結構面對巖體等效物理力學參數有重要影響，當巖體體積大于表征單元體體積（REV）時，巖體等效參數趨于穩定。以新疆提依爾金礦為工程背景，使用Matlab軟件，采用蒙特卡洛方法進行巖體結構面網絡三維模擬，建立了與實際巖體具有相同統計特征的結構面三維模型；通過定義REV指標，綜合確定巖體表征單元體為13m×13m×13m，約為最大節理跡長的3.7倍，為后續研究提供了基礎數據。

關鍵詞：Matlab 蒙特卡洛方法 結構面三維模擬 表征單元體

中圖分類號：TU452 文獻標識碼：A 文章編號：1672-3791（2015）07（a）-0093-02

研究巖體工程問題前需要合理確定巖體的等效力學參數，然而由于巖體存在復雜的結構面（斷層、節理、裂隙等），其物理力學性質表現出不均勻、不連續性。分析大尺度巖體的力學行為時，應當注意一般情況下巖體的物理力學性質隨巖體體積增大逐漸趨于穩定，即巖體存在表征單元體積（REV）。目前通常可采用數值分析的方法確定巖體的等效力學參數，國內外不少學者對此進行了深入研究[1-2]。使用數值分析方法確定巖體等效參數時，需要先明確巖體的REV尺度，只有當巖體模型體積大于巖體的REV時，才能保證數值分析法得到的力學參數與巖體實際的力學參數具有等效性。

巖體內結構面分布具有隨機性，可用概率分布描述[3]。該文以新疆提依爾金礦露天采場巖質邊坡為工程背景，根據現場結構面測量數據統計分析所得的結構面產狀概率分布及特征參數，利用Matlab軟件，運用蒙特卡洛方法建立了與邊坡巖體具有相同分布特征的結構面三維模型，然后以此為基礎，研究了邊坡巖體的REV大小，為后續邊坡穩定性研究提供了依據。

1 巖體結構面三維模擬

該文的研究對象為新疆提依爾金礦露天采場南部邊坡巖體。該部分邊坡以凝灰巖為主，結構面發育，主要表現為節理，屬于Ⅳ級結構面。該礦最終邊坡高度228 m，為合理確定邊坡巖體的等效參數以便進行下一步研究，需先確定邊坡巖體的表征單元體體積（REV），該文通過建立巖體結構面三維模型來研究巖體的REV大小。

1.1 結構面三維模擬原理

提依爾金礦南部邊坡巖體結構面主要表現為節理，因此，該文針對巖體節理建立結構面三維模型。

該文使用Baecher模型模 擬巖體中的節理，該模型被廣泛用于巖石力學相關分析[4]。該模型假設節理面形狀為薄圓盤，在三維坐標系中，可由6個參數確定：中心點坐標（x，y，z），圓盤直徑D，產狀（傾向、傾角），其中需重點說明的是節理面的圓盤直徑和模型內節理面數目的確定。

（1）節理面圓盤直徑D。

圓盤直徑D的概率分布可用負指數分布或者對數正態分布來描述，Villaesusa（1991）研究得出節理面直徑呈負指數分布時，根據實測節理跡長數據估算節理面直徑概率分布參數的公式，可以在工程中直接使用[3]。

如果節理面直徑的概率分布fD（D）為負指數分布，則有μ=σ=E（D）=πμ′/8；E（D2）=π2μ′2/32。其中，μ是節理面直徑的均值；μ′是節理面跡線長度的均值。

假設節理面出露跡線即為節理面圓盤的弦長，平均跡長μ′即為平均弦長，則有

，其中，r和D表示節理面半徑與直徑。

假設節理面直徑D服從分布fD（D），由μ′=0.25πD可得fD（D）=0.25π·f（0.25πD）。其中，f （x）為節理面跡長的概率分布密度函數。該次研究中，實測節理跡長數據均服從負指數分布，因此假設節理面直徑D也服從負指數分布，然后根據公式μ=σ=E（D）=πμ′/8，求取分布參數。

（2）節理面數目。

伍法權等[5]指出，在相對均質區域內，節理面的體密度相對恒定，與節理面方向無關。因此，模擬空間中節理面數目n可以表示為結構面體積密度λv與模擬空間體積V的乘積，即n=λvV。

節理面體積密度λv=4λL/[πE（D2）]，其中，λL為節理面線密度[3]。

使用節理面間距測量數據的均值μ〃求出λL=1/μ〃后，再通過（1）中的公式求出E（D2），就可得到節理體積密度。

如果巖體中有m組節理面，則總體密度λv總=（λv1+λv2+……+λvm）/2π，其中λv1，λv2，……，λvm分別為第1組至第m組節理的體積密度。

1.2 蒙特卡洛法建立結構面三維模型

在現場勘察中，共得到208條節理測量數據。各組結構面幾何參數的概率分布見表1，其中傾向和傾角數據服從正態分布，跡長數據服從負指數分布。

根據以上數據，該文使用蒙特卡洛法建立巖體結構面（節理面）三維模型的步驟為：（1） 假設生成節理三維模型的空間為立方體，坐標原點為立方體中心。（2）如前文所述，確定節理面數目和節理面直徑的概率分布與參數。（3）由泊松分布假設可知節理面中心點的位置服從均勻分布，利用Matlab軟件的隨機數生成函數，隨機產生節理面中心點坐標（x，y，z）。（4）根據表1中的概率分布參數和步驟（2）中得到的節理面直徑概率分布，利用Matlab軟件隨機產生節理面的產狀參數與直徑D。

通過上述步驟生成節理模型數據后，再編寫程序將這些數據導出為dxf格式，就可在CAD中觀察建立的巖體結構面模型。例如，模型邊長為20 m時，巖體結構面三維模型如圖1所示。

將結構面三維模型與現場巖體照片對比校驗，可以看出現場實際坡面節理分布情況和模型中節理分布情況在產狀和規模上具有相似性，表明該文建立的結構面三維模型是合理的。

2 表征單元體的確定

巖體力學與變形參數等存在明顯的尺寸效應，巖體試樣尺寸較小時，其參數隨著尺寸的改變發生顯著變化；當試樣尺寸增大至某一臨界值時，其參數趨于穩定，該臨界尺寸即為表征單元體體積——REV[6]。可以定義如下三個REV指標[7]。

其中，Rxy、Ryz、Rxz分別為巖體對應于坐標系oxy、oyz、oxz平面的REV指標；Si為第i條節理面的面積；Ni為第i條節理面的單位法向量；Nxy、Nyz、Nxz為oxy、oyz、oxz平面的單位法向向量；V為模型體積；m為巖體中節理面總數。上述三個REV指標的表示巖體所有節理面在oxy、oyz、oxz坐標平面上的投影面積總和與模擬區域體積的比值。

該文編制了Matlab計算程序，自動建立結構面三維模型并求解相應的REV指標值。計算過程中，分別取立方體模型邊長為5 m，7 m，9 m，11 m，13 m，15 m，17 m，20 m，23 m，25 m。對每個尺寸，采用Matlab建立10個巖體節理面立方體模型并求解相應的Rxy、Ryz、Rxz，然后求10個模型Rxy、Ryz、Rxz指標的平均值為該尺寸巖體的最終REV指標。最后，將立方體模型邊長為5～25 m時的REV指標計算結果繪制成曲線圖，見圖2。

從圖2可看出，REV指標在V=13 m時開始趨于穩定，即認為REV尺度為13 m×13 m×13 m。同時從表1中可以看出，應用巖體結構面三維模擬確定的REV大小約為節理最大平均跡長的3.7倍。

提依爾金礦露天采場南部邊坡高度228 m，長度約200 m，遠大于巖體的表征單元體體積，下一步研究邊坡巖體穩定性問題時可將巖體視為等效連續介質；同時使用數值分析法確定巖體等效力學參數時，數值模型尺寸不應小于13m×13m×13m。

3 結語

（1）通過對提依爾金礦南部邊坡巖體進行結構面三維模擬，可看出使用Matlab進行巖體結構面的三維模擬是可行的，所得結構面三維模型與邊坡面實際節理分布的對比校驗結果也說明結構面三維模型是合理的。

（2）借助REV指標，Matlab程序求解出的表征單元體為13m×13m×13m，邊長約為巖體中節理最大平均跡長的3.7倍。

（3）研究區域邊坡的長度和深度遠大于表征單元體尺寸，下一步研究邊坡穩定性問題和數值分析法確定巖體等效力學參數時，所用數值模型尺寸不應小于13m×13m×13m。

參考文獻

[1] 徐光黎，潘別桐，晏同珍.節理巖體變形模量估算新方法[J].地球科學（中國地質大學學報），1991，16（5）：573-580.

[2] 何滿潮，薛延河，彭延飛. 工程巖體力學參數確定方法的研究[J]. 巖石力學與工程學報，2001，20（2）：225-229.

[3] 賈洪彪，等.巖體結構面三維網絡模擬理論與工程應用[M]，北京：科學出版社，2008.

[4] 汪小剛，賈志欣，張發明，等. 巖體結構面網絡模擬原理及其工程應用[M]，北京：中國水利水電出版社，2010.

[5] 伍法權，王思敬，宋勝武等. 巖體力學中的統計方法與理論[J]. 科學通報，1993，38（15）：1345-1354.

[6] 安玉華，王清，基于三維裂隙網絡的裂隙巖體表征單元體研究[J]，巖土力學，2012，33（12）：3775-3780.

[7] 張貴科，徐衛亞. 裂隙網絡模擬與REV 尺度研究[J]，巖土力學，2008，29（6）：1675-1680.