基于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.