岩土工程材料的数字图像有限元分析

第23卷第6期岩石力学与工程学报23(6)：889～8972004年3月ChineseJournalofRockMechanicsandEngineeringMarch，20042002年8月28日收到初稿，2002年10月19日收到修改稿。*特别行政区政府研究基金委员会(CA99100，EG01，HKU7029102E)和香港赛马会慈善基金资助项目。作者岳中琦简介：男，40岁，博士，1983年毕业于北京大学，现任副教授、博士生导师，主要从事岩土工程、道路工程、滑坡调查、斜坡加固和应用力学等方面的研究工作。E-mail：[e-mail]。岩土工程材料的数字图像有限元分析*岳中琦陈沙郑宏谭国涣(香港大学土木工程系香港)摘要综合数字图像处理理论、几何矢量转换技术与有限元网格自动生成原理，提出了岩土工程材料的数字图像有限元分析方法。以岩土工程材料的图像为研究对象，先采用数字图像算法获得材料的真实细观结构；然后，通过几何矢量转换技术将二元图像的细观结构转化为矢量化的细观结构；在矢量化细观结构的基础上，再利用网格自动生成技术生成材料细观结构的有限元网格；昀后，采用传统的有限元计算理论，对细观结构中的不同材料赋予相应的材料参数，从而实现了岩土工程材料的非均质力学分析，真实地反映了岩土工程材料的力学性能。以花岗岩为例，详细阐明了实现数字图像数值分析方法的整个过程。Brazilian试验的数值模拟表明，材料的细观结构对Brazilian圆盘中轴线上拉应力的分布有着明显的影响，数字图像数值分析方法可真实地实现材料的非均质研究。关键词数值分析，岩土工程材料，数字图像，细观结构，几何矢量化分类号O241文献标识码A文章编号1000-6915(2004)06-0889-09DIGITALIMAGEPROCEEDINGBASEDONFINITEELEMENTMETHODFORGEOMATERIALSYueZQ，ChenS，ZhengH，ThamLG(DepartmentofCivilEngineering，TheUniversityofHongKong，HongKong，China)AbstractThispaperpresentsaninnovativedigitalimageprocessingmethodbasedonamethodtoexaminethemechanicalresponsesofgeomaterialsbytakingintoaccounttheactualmaterialinhomogeneity.Theproposedmethodincorporatesthedigitalimageprocessingprocedures，geometryvectorizationalgorithmsandautomaticfiniteelementmeshgenerationtechniquesintotheconventionalfiniteelementmethods.Usingagraniticrockasanexample,thepaperillustratestheinnovativenumericalmethodformechanicalanalysesofrockmaterials.Numericalresultspresentedinthepapershowthatthismethodcanbeusedtocalculatethemechanicalresponsesofinhomogeneousgeomaterialswithhighefficiencyandaccuracy.Keywordsnumericalanalysis，geomaterial，digitalimage，microstructure，geometryvectorization1概述土、岩石、沥青混凝土和水泥混凝土等岩土工程材料是天然材料或天然材料所组成的非均质材料，它们常常由孔洞、砂砾和石头等不同物质组成。这些物质的物理力学性质各不相同，当岩土工程材料受到外力作用时，它们对外力的响应有着很大的差异，同时，它们之间的相互作用也非常复杂。显然，岩土工程材料的力学性能，如应力分布、裂隙扩展、破坏模式等与材料的非均质性及细观结构有着密切的联系。•890•岩石力学与工程学报2004年文献调查表明，目前大多数岩土工程材料的研究都是在传统表象理论基础上开展的宏观分析，众多力学理论和数值模型忽略了材料的细观结构，在假定材料为均质或分片均质的基础上，岩土工程材料的传统力学性能受到了较多的重视，而其细观力学性质却被大大忽略。昀近，文[1～8]试图开展岩土工程材料的非均质性数值分析，这些研究主要是基于在统计学层次上生成的虚拟细观结构，使生成的细观结构中各种材料的形状和分布与材料的真实情形尽可能一致。虽然以上研究得到了岩土工程材料的一些细观力学性质，然而，似乎很少有文献来说明这类虚拟细观结构能在多大程度上真实地体现岩土工程材料的非均质性。笔者认为，建立在虚拟细观结构上的力学分析是很难得到材料的真实力学性质的。要研究非均质的岩土工程材料，真实地反映材料的细观力学性能，需要开发新的数值计算方法。以往难以对材料实现真实的非均质性研究是由于缺少有效的设备与工具来观测材料的内部结构，以期得到其内部不同性质材料的真实几何分布。由于岩土工程材料在转换为图像时，其内部不同物质可以通过灰色度的变化在图像中得到再现，图像很好地反映了材料的细观结构，保留了材料的非均质性。所以，数字图像处理方法为研究岩土工程材料的力学性能提供了一种新途径。数字图像处理是一门系统地研究数字图像理论、技术和应用的新学科，其主要过程是将研究的物体转化为储存在计算机中的数字图像，并运用计算机对图像信息进行分析和处理从而得到所需要的研究结果。文[9～12]表明，数字图像处理已经在岩土工程领域，如道路状况评估、混凝土细观结构定量观察中得到运用；文[13，14]应用数字图像处理方法成功地对沥青混凝土中骨料的方位、分布和形状进行了定量的分析；文[15，16]应用数字图像处理技术分析了水泥混凝土中粗骨料形状分布，得到了骨料的厚薄率与延伸率；文[5]运用二次复型技术，采用光学显微镜及扫描电镜，研究了香港地区的花岗岩材料的细观破裂模式。目前，在岩土工程领域，数字图像处理主要应用于定量的、几何上的统计分析与观测，对材料的非均质性力学研究尚未得到实现。本文将综合数字图像处理方法、几何矢量转换技术与有限元网格自动生成原理，提出基于数字图像的有限元力学分析方法。本方法可用于研究岩土工程材料内部细观结构对材料力学响应的影响。本文的数字图像数值分析方法是基于二维空间的，但根据立体逻辑变换原理[17]，可将其应用到三维空间的有限元分析。2图像数字化先用圆锯将现场得到的试样或试验室制成的试块锯成垂直或水平的光滑横截面，然后通过数码相机、普通照相机或扫描仪将所研究的试件转换为储存在计算机中的物理图像文件。可将一测量尺放置在试块横截面下，并将试块的实际尺寸拍摄纪录在图像之中。就其自然格式而言，计算机中的物理图像文件是不能直接用来加工与处理的(用已有的图像软件则不属于此类)，它必须先由图形格式转化为数字形式。数字图像在计算机中是由矩形排列的一个个图像元素，亦称像素点组成的。每个像素点是横向和纵向的扫描线组成的交叉区域，这些扫描线均有相同的宽度h。每个像素点有一对应的整数值来代表该像素点的亮度，亦称为该像素点的灰色度。对于常见的256色或二值图像，灰色度值分别为0～255和0～1。实际上，整个图像由有着不同灰色度的像素点阵组成，这个像素点阵的灰色度构成了一个离散函数f(i，j)(i，j代表着笛卡尔系统中的位置，i=1～N，j=1～M)：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=)()2()1()2()22()12()1()21()11()(MNfNfNfMfffMfffjif，，，，，，，，，，LMMMLL(1)图1(a)为岩石试块横截面的图像。岩石断面的实际直径为52mm。它由816条横向扫描线和816条纵向扫描线组成。每个扫描线宽度为0.0637mm，每个像素点的实际面积是0.004mm2。图1(b)为部分数字图像的离散函数f(i，j)示例。图像可以储存为很多种物理格式，如BMP，JPG，GIF等。根据不同的储存格式，编写了程序从物理图像中提取出离散数据函数f(i，j)，这些数据即是进行下一步数字图像处理的基础。3材料的细观结构图1(a)为香港地区典型花岗岩的试样图像，其主要由石英、长石和黑云母组成。3种矿物有不同的力学性能，它们的弹性模量分别为90，70和20第23卷第6期岳中琦等.岩土工程材料的数字图像有限元分析•891•(a)岩石试块横截面图像(b)局部坐标下部分离散函数f(i，j)示例图1香港地区典型花岗岩图像Fig.1TypicaldigitalimageforcylindricalgraniteofHongKongGPa。由于石英和长石的弹性模量相近并比黑云母的弹性模量大很多，因此，可以将黑云母看成一种物质，而将其他的矿物看成另一种物质来进行力学分析。得到花岗岩内部材料的真实分布是实现数字图像有限元计算的第1步，在本节中，将采用2种不同的算法来提取材料的真实细观结构，得到的细观结构实际上是内部材料分界面二值图(其中，灰色度为0的像素点为分界面点)。3.1区域分割算法从图1(b)的离散函数f(i，j)中可以看到，黑云母的灰色度值比较低，而其他矿物的灰色度值比较高。图2为图1(a)的灰色度直方图，其横坐标为0～255，代表着图像的256个灰色度。对于每一灰色度，相应纵坐标上的数字表示在整个图像中有多少个像素点拥有该灰色度值。图2的灰色度直方图曲线上有一明显的波谷，从相应的灰色度位置(灰色度值为100左右)，可以将整个图像分成2个部分。根据下面的算式，得到一个二值图像：⎩⎨⎧=HjifHjifjig)()(01)(，，，(2)式中：H为一整数阈值(对应波谷处的灰色度值)。H值的选定影响着花岗岩的细观结构，如果阈图2图1(a)灰色度直方图Fig.2Histogramofimagegraylevelinfig.1(a)值H设置太低，一部分黑云母将被认定为其他矿物；相反，如果阈值H设定太高，黑云母成分可能会被扩大。因此，可以逐步尝试直到得到一理想的结果。图3为设定H为100后得到的二值图，中间部分灰色度为0的黑色物体代表着黑云母矿物。与图1(a)比较，可以看出，图3正确地检测出了大部分的黑云母矿物。在现阶段，可以只考虑较大面积的黑云母，而将图3中一些小黑色物体转换成其他矿物。如前面所说，图像中每一个像素点的图像面积为1，它相当于h2的实际面积。每个黑色物体的面积可以通过计算它的整个像素点数而得到：≥io6969676464718911313113814114772757065667594119138144144151677269707788107132152159157162737780881001151351601801771741788990103109119143161172187182187194120120136146156173182185194195197199146149164176184195199199203202202200153164175183188195201206210205204200161178181186188190197205207207207202178193192195195192196200196201203199198208203206205199201202193195199195211218210212212206209211200199203196ojii•892•岩石力学与工程学报2004年图3区域分割算法结果(中间黑色区域代表黑云母矿物)Fig.3Resultofregiongrowingmethod(blackareasrepresentbiotite)2nhA=(3)