基于损伤理论的岩石细观尺度微裂隙量化分析

基于损伤理论的岩石细观尺度微裂隙量化分析

1岩石细观结构岩石是由各种矿物结晶、空隙和结构物组成的综合体。经过亿万年的地质演变和多次复杂的构造运动，使岩石含有不同阶次的随机分布的微观孔隙和裂纹。岩石受载后的宏观断裂失稳和破坏与其变形时内部微裂纹的分布以及微裂纹的产生、扩展和贯通密切相关。岩石是一种多尺度的材料，M.E.Kassner等认为只有建立多尺度的力学模型才能真正了解微裂隙起裂的机制，并应在细观尺度上对脆延性材料作进一步研究。因此，本文从岩石的扫描电镜(SEM)图像研究入手，获得了岩石微裂隙的细观结构信息，并从细观角度出发，认为岩石在细观上同样为一结构体；同时，即为岩石细观结构体引入岩石细观结构体量化信息运用于岩体损伤理论中，基于G.Swoboda损伤理论得到岩石细观结构体的损伤本构模型，从而将量化岩石细观结构特征信息运用于损伤理论之中。2岩石开裂程度表征岩石细观结构(rockmesostructure)量化试验研究就是研究岩石体内的微裂隙和微孔洞的三维形状参数。基于微裂隙的长度、方位角、间距和面积等形状参数，得到岩石的初始细观结构损伤张量，从而将细观结构损伤张量引入到损伤理论得到岩石的本构模型。K.I.Kanatani于1984年通过对微结构体的参数——微裂隙密度来表征岩石的开裂程度。D.F.Howarth和J.C.Rowlands于1987年通过对岩石切片图像分析，建立了量化岩石微细颗粒参数(形状、方位角和相对比例等)的方法，并与岩样的力学强度进行了比较，得到两者的相互关系。D.H.Campbell和J.S.Galehouse于1991年讨论了利用光学显微镜进行岩石、砖和水泥等材料微观量化的可行性，并且与X射线方法进行了比较和研究。R.H.CWong等于1996年研究了初始微裂隙密度和岩石微颗粒尺寸对于YuenLong大理石岩样单轴压缩强度的影响。Y.H.Hatzor等于1997年研究了白云石的细观结构与微裂隙起裂的初始应力和试样最终强度之间的关系，发现岩石微结构对其强度极限的影响非常大。X.Y.Wu等于2000年利用光学显微镜和SEM研究了DarleyDale砂岩在压缩破坏中的各向异性损伤的微观力学演化过程，得到了微裂隙密度与应变之间的关系。B.Menendez等于2001年探讨了运用激光共焦点扫描显微镜来研究岩石中微裂隙和微孔洞网络的方法。2.1试样制备和检测岩石试样采自四川锦屏II级水电站工程现场的大理岩，白色、块状、粒状变晶结构，均匀性好。为了保证试样性质的均一性，严格密集采样，岩样严格按照国际岩石力学学会(ISRM)试验规程来加工，加工成3组尺寸为φ50mm×100mm的圆柱体标准试件。其中2组进行单轴压缩试验和常规三轴压缩试验，另外1组按照以下介绍的方法进一步加工成3个切片进行岩石细观结构量化试验。试验切片的制作流程见图1，具体试样加工步骤如下：(1)首先从加工好的大理岩试样中选取1块，另外2块将用来做单轴压缩和常规三轴压缩试验。(2)从试样顶部切割一高大约为10mm左右的圆柱，见图1步骤aㄢ(3)在得到小部分圆柱上，切出一个大致尺寸为30mm×30mm×10mm的切片，其为XY切片(垂直于加载方向)，见图1步骤b,d,eㄢ(4)在剩余的大部分圆柱上，竖向切出一个足够大小的部分，切成大致尺寸为30mm×30mm×10mm，为YZ切片，见图1步骤c,f,h,iㄢ(5)从剩余的部分切出一个垂直于YZ切片的XZ切片，见图1步骤g,jㄢ最后对切片进行相应标号，并在切片的左上角标记，以确保三者之间的空间关系，即观测到恰好是试样内的一个结构立方体的3个面。在进行扫描电镜图像采集前，需要对切片进行清洗和干燥以及镀金，采用的步骤为：(1)在乙醇中超声波处理5min以清除表面杂质。(2)在100˚C下干燥12hㄢ(3)进行镀金。加工好的大理石切片见图2ㄢ2.2图像放大个数采用HitachiS–3000h型扫描电镜对相应的切片进行观测，采集到的图像的分辨率为640像素×480像素。图3为一张放大300倍的大理岩细观结构SEM图像，图中一个像素的距离代表0.6667µm。图像的放大倍数，决定了一张图像上能够显示的微裂隙数量。如果放大倍数太大，则可能出现一张图像上只能显示微裂隙的一部分这种情况；相反，如果放大倍数太小，则会由于微裂隙在图像上显示的太小，在图像处理过程中，丢失图像信息，造成误差过大，所以在参考F.E.Mohamed研究成果的基础上，作者进行相应的尝试，最后决定采用的放大倍数为300倍。图像的采集数量决定了是否能获取到足够多的微裂隙信息，进而形成微裂隙的信息统计分布。在试验中，每个面采集最少200张以上图片。在切片的表面上获取图片，从上而下，避免图像的重复获取。3岩石细观结构的量化分析近年来，随着计算机技术和数字图像处理的发展使得通过分析岩石切片细观图像得到岩石细观结构参数量化信息成为可能。同时数字图像处理技术的发展使得从切片微观图像上获取数据收集更方便、更准确、统计上更有代表性。3.1区域生长——局部分割过程图像分割是将图像细分为构成它的子区域或对象。分割的程度取决于要解决的问题以及如何来识别SEM图像中的微裂隙部分。图像分割算法的一般是基于亮度值的两个基本特性之一：不连续性和相似性。第1个性质的应用途径是基于亮度的不连续变化分割图像，比如图像的边缘；第2个的主要应用途径是依据事先制订的准则将图像分割为相似的区域，如门限处理、区域生长、区域分离和聚合都是这类方法。由于岩石细观SEM图像的特征，决定了采用诸如边缘检测算子如Robison,Roberts,Sobel,Prewitt,Krisch和Log等，无法很好地达到从图像中分割出微裂隙的目的，所以采用后者即采用相似性的原则运用区域生长来实现图像分割的目的。区域生长是一种事先根据定义的准则将像素或子区域聚合成更大区域的过程。基本方法是以一组“种子”点开始将与种子性质相近(诸如灰度级或颜色的特定范围)的相邻像素附加到生长区域的每个种子上。通常根据所解决的问题的性质选择一个或多个起点。当一个先验信息无效时，这一过程将对每个像素计算相同的特性集，最终这个特性集在生长过程中将像素归属某个区域。如果这些计算结果呈现不同簇的值，则那些由于自身性质而处于这些簇中心附近的像素可以作为种子。区域生长的另一个问题是用公式描绘一个终止规则。基本上，在没有像素满足加入某个区域条件时，区域生长就会停止。因此，区域生长的2个基本原则为：(1)任何像素和种子之间的灰度级绝对差必须小于某一个值β，这个数值是根据某一算法得到的，足以区别图上2个不同的区域，如果β超过了其合理大小的值，则分割出来的微裂隙将偏大；反之，则分割出来的微裂隙将偏小。(2)要添加入某一个区域的像素必须与此区域至少一个像素是八连通的。如果某个像素被发现于多于一个区域相联系，就将这些区域合并起来。形象地说，生长的过程中必须有“根”。而β值可以通过多种算法获取，如双峰法、迭代方法(最优方法)、大津法(OTSU法)、灰度拉伸法(改进的OTSU法)以及Krish算子法，具体参见刘爽的研究成果。3.2人工智能以及模糊聚类获取种子点，算法充选取种子点是通过在微裂隙区域人工点击一定种子点，当然在今后可以通过其他一些算法结合人工智能以及模糊聚类也可以实现智能获取种子点，本文不作考虑；图4为设置的区域生长算法的种子点。(2)区域生长运算符使用MATLAB进行编程运算(见图5)，得到微裂隙二值化后的图像(见图5(a))，5(b)为区域生长算法识别结果。(3)边缘有小球质性识别出来的微裂隙图像的图像中微裂隙(图6中白色部分)的边缘存在不光滑以及有一些毛刺，在微裂隙中可能有一些小的孔洞，所以有必要采用形态学的闭操作以及孔洞填充进行处理，处理后的结果见图7ㄢ3.3值化图像的识别在识别出来的二值化微裂隙图像(见图7)上得到微裂隙细观结构参数，MATLAB中有Bwlabel和Regionprops两个函数可以对二值化图像后颗粒或孔隙进行标识和统计各参数值。可以获取的参数值为颗粒或孔隙的总面积、孔隙率、总周长、面积、周长、形状系数、粒径或孔径、扁圆度和方位角。图中一个像素代表的长度为0.6667µm，从计算机识别出来的各量像素单位再转化为微米单位。从图7识别出来的微裂隙细观结构参数如表1所示。微裂隙的平面间距的表达式为式中：L为切片上沿某一方向剖线的长度，NLC为长为L的剖线通过的微裂隙数量。4岩石薄视结构信息分析对所有的3个切片上细观结构图片进行了处理，得到了大量的微裂隙的细观结构参数数据，对这些数据进行相应的统计分析。4.1密度函数的估计图8为XY面微裂隙方位角直方图和理论概率分布函数，从图8中可以看出，微裂隙的方位角服从正态分布，其正态分布的密度函数为采用点估计的方法算出µ和σ的无偏估计量的估计值为得到XY,YZ,XZ面的微裂隙方位角的统计参数的无偏估计量见表2。将表2中相应数值代入式(2)得到理论概率分布函数曲线与实际数据直方图的比较(见图8～10)，其分别为XY,YZ,XZ面的微裂隙方位角数据的直方图和代入统计参数的无偏估计量的理论概率分布函数曲线。4.2密度函数的统计描述图11为XY面微裂隙长度直方图和理论分布密度函数曲线。由图11可看出，微裂隙长度分布基本是服从对数正态分布，其对数正态分布的密度函数为同样得到XY,YZ,XZ面的微裂隙长度的统计参数的无偏估计量见表2。同样将表2中相应数值代入式(4)得到理论概率分布函数曲线与实际数据直方图的比较(见图11～13)，其分别为XY,YZ,XZ面的微裂隙长度数据的直方图和理论概率分布函数曲线。4.3微裂隙间距的统计参数由图14可以看出，XY面的微裂隙间距分布为指数分布，其密度函数为同样得到XY,YZ,XZ面的微裂隙间距的统计参数的无偏估计量见表2。将表2中相应数值代入式(5)得到理论概率分布函数曲线与实际数据直方图的比较(见图14～16)，其分别为XY,YZ,XZ面的微裂隙间距数据的直方图和理论概率分布函数曲线。5岩石细观结构体的模拟微裂隙的方位角服从正态分布，长度服从对数正态分布，间距服从指数分布。模拟岩石细观结构体是通过模拟符合以上3种统计分布的微裂隙来实现的，程序采用MATLAB编制，MATLAB中的随机变量模拟是采用逆变换方法的得到，即所有的分布都是从均匀分布中得到。5.1xy面的微裂隙解析微裂隙的中心坐标是通过模拟微裂隙的间距得到的，将表2中微裂隙间距分布参数输入，编程模拟得到一定量随机微裂隙间距，其为一矩阵S(a,b)。采用式(6)计算得到各模拟微裂隙的中心坐标，即式中：Cx(a,b)，Cy(a,b)为(a,b)处的微裂隙中心。模拟得到XY面的微裂隙中心见图17ㄢ在得到微裂隙中心以后，同样，按照微裂隙方位角的分布参数和长度的分布参数模拟得到相应的微裂隙方位角和长度，微裂隙两端点的坐标分别为式中：为第i条微裂隙的起点坐标；为第i条微裂隙的终点坐标；为第i条微裂隙的中点坐标；iL为第i条微裂隙的长度；θi为第i条微裂隙的方位角。采用式(6)分别计算岩石细观结构的3个表面上微裂隙的起点和终点坐标，运用程序绘制XY,YZ和XZ三个面上的微裂隙，分别见图18～20。最后，模拟细观结构体的三维显示见图21ㄢ5.2微裂隙的尺寸大小由于微裂隙间距的分布情况直接决定了微裂隙的密集程度，也即反映了微裂隙单位面积的数量，所以认为岩石细观REV的大小由微裂隙的间距决定，模拟得到在同一间距分布而不同数量的微裂隙具有的尺寸大小，得到面积与单位面积微裂隙数量之间的关系。图22为以XY面微裂隙间距分布参数得到的以上两者关系图，从图22中可看出，当面积为40×106µm2时，单位面积的微裂隙数趋于稳定，因此得到细观结构体即为REV，即可得到REV各个面上的微裂隙数量、REV的体积以及统计分析得到微裂隙方位角、长度和间距各分布函数的数学期望，见表3ㄢ6岩石细观损伤的本构模型6.1微裂隙损伤张量的计算T.Kyoya等基于Murakami-Ohno蠕变损伤理论给出了节理岩体的损伤基本理论，仿照岩体损伤张量计算方法，计算得到岩石细观结构初始损伤张量。设微裂隙结构体在其体积V内的裂隙平均面积为a，平均微裂隙数为N，则损伤张量可以表示为式中：n为微裂隙的法向矢量，n={n1,n2,n3}T，且有同时，参数λ满足：另外，不同面上的方位角之间满足：式中：Ni(i=1,2)为第in面上平均微裂隙数。对于存在多组的微裂隙，以下为计算存在多组微裂隙的损伤张量的方法，首先依据以上理论求出每组的损伤张量，然后依据能量等效原理，即由于全局损伤张量引起的能量消散等于各个损伤张量引起的能量消散的总合。计算全局损伤张量的表达式为式中：Ωg为全局损伤张量，Ωn为第n组微裂隙的损伤张量，N为微裂隙组别的数量，I为二阶单位张量。将表3中细观结构量化参数数据代入式(8a)和(13)得到初始损伤张量矩阵为6.2g.swoboda损伤理论6.2.1连续介质力学把细观损伤信息引入刚度矩阵是利用连续介质力学方法求解的关键。G.Swoboda和Q.Yang引入了一种计算损伤弹性系数张量的方法。这种方法所得到的损伤弹性系数张量符合对称、正定等连续介质力学条件。其计算公式为式中：λ和µ为Lamé常数；C1,C2和3C均为Ωij的不变量的函数，且1C,C2和C3满足：由于1C,C2和3C为常数，而且Φij为非负性，且有6.2.2主应变损伤张量的确定在节理岩体中，一般认为损伤与应变方向相关联。为了考虑裂缝的开裂与闭合以及由于裂缝中的夹杂材料或表面粗糙不平而引起的不同的裂缝表面接触情况对应力传递的影响，基于主应变的有效损伤张量可定义为式中：pv和εv为第v个的主应力的主矢量和主值；H(·)是修正了的阶跃函数，其表达式为6.2.3．非负材料一阶段的损伤张量的面中,且对于等效状态，基于共轭力的损伤演化法则对于任何应力状态都是适合的，且有其中，的推导在有关研究中有详细介绍，这里不作介绍，可表示为其中，式中：和ωv分别为损伤张量Ω的主方向和主值，且关联的损伤面为式中：β为材料常数，g为非负材料常数。6.2.4本结构方程本构方程可以表示为6.3woda损伤理论由于在G.Swoboda中有一些参数需要事先确定，而其中采用三轴试验可以确定一部分参数，其他的参数可