CN116429295A - 一种评估岩体结构面接触应力分布的方法 - Google Patents

一种评估岩体结构面接触应力分布的方法 Download PDF

Info

Publication number
CN116429295A
CN116429295A CN202310466317.3A CN202310466317A CN116429295A CN 116429295 A CN116429295 A CN 116429295A CN 202310466317 A CN202310466317 A CN 202310466317A CN 116429295 A CN116429295 A CN 116429295A
Authority
CN
China
Prior art keywords
expressed
microprotrusions
point cloud
contact
microprotrusion
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202310466317.3A
Other languages
English (en)
Inventor
李育宗
崔跃
甘雨翔
张庆贺
沈树强
高胜
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Anhui University of Science and Technology
Original Assignee
Anhui University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Anhui University of Science and Technology filed Critical Anhui University of Science and Technology
Priority to CN202310466317.3A priority Critical patent/CN116429295A/zh
Publication of CN116429295A publication Critical patent/CN116429295A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01L—MEASURING FORCE, STRESS, TORQUE, WORK, MECHANICAL POWER, MECHANICAL EFFICIENCY, OR FLUID PRESSURE
    • G01L1/00—Measuring force or stress, in general
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B21/00—Measuring arrangements or details thereof, where the measuring technique is not covered by the other groups of this subclass, unspecified or not relevant
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computing Systems (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开一种评估岩体结构面接触应力分布的方法，包括：将结构面的上盘点云和下盘点云进行初步点云配准后，再在坐标系上分区，取均值作为坐标系的横纵竖方向上方格的代表值，对上盘点云和下盘点云进行精简，并对再一次进行精准的点云配准；获取精准的点云配准后的上盘点云和下盘点云的复合形貌；并按照高程矩阵分布确定微凸体的位置及几何特征；根据结构面的力学参数确定微凸体和结构面基体的刚度，以及两者的串联刚度；根据平衡方程确定微凸体的接触变形，以及获取当前应力状态下微凸体引起的接触面积；获取微凸体区域的总接触应力，并遍历所有微凸体，获取结构面的接触应力分布。通过本发明能够估计岩体结构面接触区应力的大小。

Description

一种评估岩体结构面接触应力分布的方法
技术领域
本发明涉及岩体结构面间的应力测量技术领域，具体为一种评估岩体结构面接触应力分布的方法。
背景技术
结构面广泛存在于天然岩体当中，结构面的接触状态对于岩体结构的稳定性、渗流特性和应力波的传播特性具有决定性的影响。结构面上接触应力的大小对于揭示岩体渗流及滑移破坏规律及内在力学机制至关重要。不仅如此，对结构面接触应力的评估对于岩体滑移型灾害，如岩质滑坡和结构面滑移型岩爆的预警预报具有重要意义。
在过往的研究中，学者们对岩体结构面的法向本构做了详细的研究，从经验本构模型到统计本构模型做了大量研究，但并没有针对结构面不同接触区的应力做详细研究和讨论。在法向应力作用下，结构面间的真实接触区域只占结构面名义接触面积的很小一部分，接触区域内的真实接触应力与结构面名义平均接触应力存在较大差异。结构面上接触区域的真实接触应力是导致结构面在法向压缩或剪切滑移过程中微凸体破坏的直接原因，因此对于结构面上真实接触应力的评估是预测结构面压缩或剪切行为的基础。因此计算岩体结构面接触区域的真实接触应力对于预测结构面力学行为、预测预警结构面型灾害具有重要意义。
现有技术，公开号为CN106682347A的发明专利，一种基于岩体结构面随机强度的非连续变形分析方法，通过对岩体中的每一条结构面生成随机强度，基于每一条结构面的随机强度进行非连续变形分析模拟，这种强度参数的取值方法更接近岩体的实际情况，计算及分析结果更为准确、真实。
发明内容
本发明所要解决的技术问题在于：解决是估计岩体结构面接触区应力的大小问题。
为解决上述技术问题，本发明提供如下技术方案：
一种评估岩体结构面接触应力分布的方法，包括：
S100，将结构面的上盘点云和下盘点云进行初步点云配准；
S200，将经过初步点云配准的上盘点云和下盘点云在坐标系上分区，取均值作为坐标系的横纵竖方向上方格的代表值，对上盘点云和下盘点云进行精简，并对精简后的上盘点云和下盘点云再一次进行精准的点云配准；
S300，获取精准的点云配准后的上盘点云和下盘点云的复合形貌；并按照高程矩阵分布确定微凸体的位置及几何特征；
S400，根据所述结构面的力学参数确定所述微凸体和结构面基体的刚度，以及两者的串联刚度；
S500，根据平衡方程确定所述微凸体的接触变形；
S600，根据所述微凸体的接触变形，获取当前应力状态下微凸体引起的接触面积；
S700，获取微凸体区域的总接触应力，并遍历所有微凸体，获取结构面的接触应力分布。
优点：通过结构面形貌上的微凸体特征，将微凸体作为接触的媒介，通过识别微凸体的具体位置，并根据其几何特性计算其接触应力，遍历所有接触的微凸体之后，即可得到此岩体结构面上不同接触区的接触应力，且可以绘制云图，根据微凸体的位置查询此区域的接触应力。
在本发明的一实施例，步骤S200包括：首先将经过初步点云配准的上盘点云和下盘点云在坐标系横纵方向上划分方格，其均值作为横纵方向上方格的代表值；再将每个方格区域内的竖方向坐标取均值，并代表这个方格区域内的高程坐标；最后获取精简后的上盘点云和下盘点云。
在本发明的一实施例，在坐标系横纵方向上划分方格的划分依据通过以下公式获取：
cell(ii,jj)＝{x,y,z|0.5ii≤x≤0.5(ii+1),0.5jj≤y≤0.5(jj+1)}；
式中，x、y和z分别表示为点云在坐标系横纵竖方向上的坐标，cell(ii,jj)表示为一个点云所在方格区域的名称及位置；
所述高程坐标通过以下公式获取：
Figure BDA0004202430890000031
式中，z'(ii,jj)表示为点云在方格区域内的高程坐标，Nc(ii,jj)表示为点云在方格区域点云的数目，z'表示为方格区域内所有点云的竖坐标。
在本发明的一实施例，步骤S300包括：
S310，将精准点云配准后的上盘点云和下盘点云的形貌组合，获取复合形貌；
S320，判断所述高程矩阵内的元素是否满足判断规则，若满足所述判断规则，则对应位置的区域形貌判定为微凸体；
S330，将步骤S320筛选出来的微凸体的顶点坐标以及在所述定点坐标正交方向相邻四个点的坐标，通过最小二乘法拟合成球体，将所述球体的半径作为微凸体的曲率半径。
在本发明的一实施例，所述复合形貌通过以下公式获取：
comb(x,y)＝zu(x,y)+zl(x,y)；
式中，comb(x,y)表示为复合形貌，comb表示为复合后的高程矩阵，zu表示为精准点云配准后的上盘点云，zl表示为精准点云配准后的下盘点云，x、y分别表示为点云在坐标系横纵竖方向上的坐标；
所述判断规则通过以下公式获取：
comb(ii,jj)＞max{comb(ii-1,jj),comb(ii+1,jj),comb(ii,jj+1),comb(ii,jj-1)}；
式中，max表示为最大值；comb表示为复合后的高程矩阵，(ii,jj)表示为高程矩阵的为位置；comb(ii-1,jj),comb(ii+1,jj),comb(ii,jj+1),comb(ii,jj-1)表示为复合后的高程矩阵在正交方向上相邻四个点的高程。
在本发明的一实施例，步骤S400包括以下步骤：
S410，根据Hertz微凸体接触理论，获取单个微凸体发生弹性变形时的刚度；
S420，根据弹性理论，获取当微凸体受到荷载时，在微凸体中心的最大接触压力，以及微凸体和结构面基体作用的部分，在接触区内的沉降；
S430，根据所述沉降获取结构面基体的刚度。
在本发明的一实施例，所述Hertz微凸体接触理论，通过以下公式获取：
Figure BDA0004202430890000041
式中，fan表示为第n个微凸体所受的压力，E表示为岩体的弹性模量，ρn表示为复合形貌中第n个微凸体的曲率半径，δan表示为复合形貌中第n个微凸体的变形量；
所述单个微凸体发生弹性变形时的刚度通过以下公式获取：
Figure BDA0004202430890000042
式中，kan表示为单个微凸体发生弹性变形时的刚度；
所述最大接触压力通过以下公式获取：
Figure BDA0004202430890000043
式中，P0表示为最大接触压力，FN表示为外载荷，π表示为圆周率，rb表示为微凸体与基体的接触半径；
所述接触区内的沉降通过以下公式获取：
Figure BDA0004202430890000044
式中，Uz(r)表示为接触区内的沉降，E*表示为结构面上下盘接触的等效弹性模量，r表示为微凸体与结构面基体接触范围；
所述结构面基体的刚度通过以下公式获取：
Figure BDA0004202430890000045
式中，kbn表示为结构面基体的刚度，δbn表示为结构面基体的变形量；Uz(0)表示为微凸与结构面基体接触范围为0时，其接触区的沉降；
所述微凸体和所述结构面基体的串联刚度，通过以下公式获取：
Figure BDA0004202430890000051
式中，kn表示为串联刚度，kan表示为第n个微凸体的刚度。
在本发明的一实施例，在步骤S500中，所述微凸体的接触变形通过以下公式获取：
Figure BDA0004202430890000052
式中，f(δn)表示为微凸体的接触变形，δn表示为总变形，ρn表示为复合形貌的微凸体的曲率半径，rb表示为微凸体与结构面基体的接触半径。
在本发明的一实施例，步骤S600中，所述获取当前应力状态下微凸体引起的接触面积，通过以下公式获取：
An＝πρnδan＝πρn(δn-δbn)；
式中：An表示为微凸体引起的接触面积，π表示为圆周率，ρn表示为复合形貌的微凸体的曲率半径，δan表示为复合形貌上第n个微凸体的变形，δn表示为总变形，δbn表示为结构面基体的变形量。
在本发明的一实施例，步骤S700中，所述微凸体区域的总接触应力，通过以下公式获取：
Figure BDA0004202430890000053
式中，σn表示为第n个微凸体区域的总接触应力，An表示为微凸体区域的接触面积，δan表示为复合形貌上第n个微凸体的变形，π表示为圆周率，f(δan)表示为复合形貌上第n个微凸体由于变形产生的接触力，ρn表示为复合形貌的微凸体的曲率半径；π表示为圆周率，ρn表示为复合形貌的微凸体的曲率半径，E*表示为结构面上下盘接触的等效弹性模量。
与现有技术相比，本发明的有益效果是：
(1)本发明将上盘点云和下盘点云在坐标系上分区并计算结构面的法向接触应力，用坐标描述结构面上的微凸体形貌更精确，在计算微凸体几何特征时，得到的几何参数更精准。
(2)本发明进行了两次点云配准，且计算过程中了考虑可微凸体与结构面基体的接触半径，使其在计算微凸体形貌收敛效果更好。
(3)本发明有严谨的数学力学推导过程，在接触应力预测方面可更好的与实际工程相结合，预测结构面产生滑移时产生应力较大的区域。
(4)本发明的计算过程较容易与编程程序相结合，产生具体的实施流程，与工程实际相结合，易于形成集收集信息，预测，防护手段为一体的系统。
(5)本发明可根据计算能力的不同设置不同的计算精度，比如不考虑结构面基体的变形或者引入摩擦参数，具体可根据实际工程精度需求设置，提高工作效率。
附图说明
图1为本发明实施例的一种评估岩体结构面接触应力分布的方法的流程图。
图2为本发明实施例结构面的上、下盘点云示意图。
图3为本发明实施例的复合形貌示意图。
图4为本发明实施例的微凸体与其邻域接触示意图。
具体实施方式
为便于本领域技术人员理解本发明技术方案，现结合说明书附图对本发明技术方案做进一步的说明。
术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括一个或者更多个该特征。在本申请的描述中，“多个”的含义是两个或两个以上，除非另有明确具体的限定。
请参图1所示，本发明提供一种评估岩体结构面接触应力分布的方法，包括以下步骤：
S100，将结构面的上盘点云和下盘点云进行初步点云配准。
请参阅图1和图2所示，在本发明的一实施例中，在步骤S100中，用三维扫描仪获取的点云如图2所示，用点云处理软件例如：CloudCompare或Geomagic Studio等软件进行去噪，封装等操作，将结构面上下盘点云配准成上下面吻合的位置，采用ICP(IterativeClosest Point)配准方法进行迭代，迭代直至误差目标小于0.001。具体的，包括以下步骤：
S110，设定参考点云和测试点云，获取沿点云质心刚体平移后的参考点云和测试点云的新点云坐标矩阵。
其中，所述参考点云和测试点云的新点云坐标矩阵，通过以下公式获取：
Figure BDA0004202430890000071
Figure BDA0004202430890000072
式中，P'表示为测试点云平移后新点云坐标矩阵，P表示为测试点云，p”表示为测试点云的三维坐标，A表示为参考点云，A'表示为参考点云平移后新点云坐标矩阵，a'表示为参考点云的三维坐标，N表示为点云的数目，i＝1，2，3......N。
S120，获取参考点云和测试点云的新点云坐标矩阵构成的中间变量，并对所述中间变量按照奇异值分解。
其中，对所述中间变量按照奇异值分解，通过以下公式：
Figure BDA0004202430890000073
式中，W表示为中间变量，T表示为矩阵装置符号，U和V表示为中间变量的奇异向量，S表示为中间变量的奇异值组成的特征向量。
S130，根据所述中间变量的奇异向量获取点云刚体变换的旋转矩阵，并根据所述旋转旋转矩阵，获取点云刚体变换的平移矩阵。
其中，所述点云刚体变换的旋转矩阵和平移矩阵通过以下公式获取：
R＝U*VT；
Figure BDA0004202430890000081
式中，R表示为点云刚体变换的旋转矩阵，T0表示为点云刚体变换的平移矩阵，ai表示为参考点云的坐标向量，pi表示为测试点云的坐标向量
S140，利用变化矩阵，将所述参考点云和所述测试点云在同一坐标系下，进行初步点云配准。
其中，这里所说的变化矩阵为点云刚体变换的旋转矩阵和平移矩阵，以及所述初步点云配准的误差目标函数为：
Figure BDA0004202430890000082
式中，f(R,T0)表示为初步点云配准误差目标函数。
S200，将经过初步点云配准的上盘点云和下盘点云在坐标系上分区，取均值作为坐标系的横纵竖方向上方格的代表值，对上盘点云和下盘点云进行精简，并对精简后的上盘点云和下盘点云再一次进行精准的点云配准。
请参阅图1和图2所示，在本发明的一实施例中，步骤S200包括：首先将经过初步点云配准的上盘点云和下盘点云在坐标系横纵方向上划分方格，其均值作为横纵方向上方格的代表值。再将每个方格区域内的竖方向坐标取均值，并代表这个方格区域内的高程坐标。最后获取精简后的上盘点云和下盘点云。
具体的，在坐标系横纵方向上划分方格的划分依据通过以下公式获取：
cell(ii,jj)＝{x,y,z|0.5ii≤x≤0.5(ii+1),0.5jj≤y≤0.5(jj+1)}；
式中，x、y和z分别表示为点云在坐标系横纵竖方向上的坐标，cell(ii,jj)表示为一个点云所在方格区域的名称及位置。即在本实施例中，经过初步点云配准的上盘点云和下盘点云按照横纵坐标每隔0.5mm的距离划分成方格。
所述高程坐标通过以下公式获取：
Figure BDA0004202430890000091
式中，z'(ii,jj)表示为点云在方格区域内的高程坐标，Nc(ii,jj)表示为点云在方格区域点云的数目，z'表示为方格区域内所有点云的竖坐标，单位为mm。其中，ii，jj∈[1,200]且取整。
通过此步骤获取两个精简之后的点云，再进行一次ICP配准，即可快速精确的达到吻状态。
S300，获取精准的点云配准后的上盘点云和下盘点云的复合形貌；并按照高程矩阵分布确定微凸体的位置及几何特征。
请参阅图1和图3所示，在本发明的一实施例中，步骤S300包括以下步骤：
S310，将精准点云配准后的上盘点云和下盘点云的形貌组合，获取复合形貌。
其中，所述复合形貌通过以下公式获取：
comb(x,y)＝zu(x,y)+zl(x,y)；
式中，comb(x,y)表示为复合形貌，comb表示为复合后的高程矩阵，zu表示为精准点云配准后的上盘点云，zl表示为精准点云配准后的下盘点云，x、y分别表示为点云在坐标系横纵竖方向上的坐标。
S320，判断所述高程矩阵内的元素是否满足判断规则，若满足所述判断规则，则对应位置的区域形貌判定为微凸体。
其中，所述判断规则通过以下公式获取：
comb(ii,jj)＞max{comb(ii-1,jj),comb(ii+1,jj),comb(ii,jj+1),comb(ii,jj-1)}；
式中，max表示为最大值；comb表示为复合后的高程矩阵，(ii,jj)表示为高程矩阵的为位置，comb(ii-1,jj),comb(ii+1,jj),comb(ii,jj+1),comb(ii,jj-1)表示为复合后的高程矩阵在正交方向上相邻四个点的高程。
S330，将步骤S320筛选出来的微凸体的顶点坐标以及在所述定点坐标正交方向相邻四个点的坐标，通过最小二乘法拟合成球体，将所述球体的半径作为微凸体的曲率半径。
S400，根据所述结构面的力学参数确定所述微凸体和结构面基体的刚度，以及两者的串联刚度。
请参阅图1和图4所示，在本发明的一实施例中，步骤S400包括以下步骤：
S410，根据Hertz微凸体接触理论，获取单个微凸体发生弹性变形时的刚度。
其中，所述Hertz微凸体接触理论，通过以下公式获取：
Figure BDA0004202430890000101
式中，fan表示为第n个微凸体所受的压力，E表示为岩体的弹性模量，ρn表示为复合形貌中第n个微凸体的曲率半径，δan表示为复合形貌中第n个微凸体的变形量。
所述单个微凸体发生弹性变形时的刚度通过以下公式获取：
Figure BDA0004202430890000102
式中，kan表示为单个微凸体发生弹性变形时的刚度。
S420，根据弹性理论，获取当微凸体受到荷载时，在微凸体中心的最大接触压力，以及微凸体和结构面基体作用的部分，在接触区内的沉降。
其中，所述最大接触压力通过以下公式获取：
Figure BDA0004202430890000103
式中，P0表示为最大接触压力，FN表示为外载荷，π表示为圆周率，rb表示为微凸体与结构面基体的接触半径。
其中，
Figure BDA0004202430890000104
式中，ρn表示为复合形貌中第n个微凸体的曲率半径，hn表示为复合形貌中第n个微凸体的高度，即hn相当于对应微凸体位置复合形貌的高程comb(ii,jj)。
所述接触区内的沉降通过以下公式获取：
Figure BDA0004202430890000111
式中，Uz(r)表示为接触区内的沉降，E*表示为结构面上下盘接触的等效弹性模量，r表示为微凸体与结构面基体接触范围；rb表示为微凸体与基体的接触半径。其中
Figure BDA0004202430890000112
式中，ν表示为材料的泊松比，为一个常数。
S430，根据所述沉降获取结构面基体的刚度。
其中，所述结构面基体的刚度通过以下公式获取：
Figure BDA0004202430890000113
式中，kbn表示为结构面基体的刚度，δbn表示为结构面基体的变形量，Uz(0)表示为微凸与结构面基体接触范围为0时，其接触区的沉降。
所述微凸体和所述结构面基体的串联刚度，通过以下公式获取：
Figure BDA0004202430890000114
式中，kn表示为串联刚度，kan表示为第n个微凸体的刚度。
S500，根据平衡方程确定所述微凸体的接触变形。
请参阅图1和图4所示，在本发明的一实施例中，在步骤S500中，根据相互作用力大小相等，则有：
kanδan＝kbnδbn＝kbn(δn-δan)；
式中，kan表示为第n个微凸体的刚度，δan表示为复合形貌中第n个微凸体的变形量，kbn表示为结构面基体的刚度，δbn表示为结构面基体的变形量，δn表示为总变形。
对上式进行化简，则有：
Figure BDA0004202430890000121
式中，ρn表示为复合形貌中第n个微凸体的曲率半径，rb表示为微凸体与结构面基体的接触半径。
微凸体的变形与总变形相关，也即δan＝f(δ)，用不动点迭代法可解出，则所述微凸体的接触变形通过以下公式获取：
Figure BDA0004202430890000122
式中，f(δn)表示为微凸体的接触变形，δn表示为总变形，ρn表示为复合形貌的微凸体的曲率半径，rb表示为微凸体与结构面基体的接触半径。
S600，根据所述微凸体的接触变形，获取当前应力状态下微凸体引起的接触面积。
请参阅图1和图4所示，在本发明的一实施例中，步骤S600中，所述获取当前应力状态下微凸体引起的接触面积，通过以下公式获取：
An＝πρnδan＝πρn(δn-δbn)；
式中：An表示为微凸体引起的接触面积，π表示为圆周率，ρn表示为复合形貌的微凸体的曲率半径，δan表示为复合形貌上第n个微凸体的变形，δn表示为总变形，δbn表示为结构面基体的变形量。
S700，获取微凸体区域的总接触应力，并遍历所有微凸体，获取结构面的接触应力分布。
请参阅图1和图4所示，在本发明的一实施例中，步骤S700中，所述微凸体区域的总接触应力，通过以下公式获取：
Figure BDA0004202430890000131
式中，σn表示为第n个微凸体区域的总接触应力，An表示为微凸体区域的接触面积，δan表示为复合形貌上第n个微凸体的变形，π表示为圆周率，f(δan)表示为复合形貌上第n个微凸体由于变形产生的接触力，ρn表示为复合形貌的微凸体的曲率半径；π表示为圆周率，ρn表示为复合形貌的微凸体的曲率半径，E*表示为结构面上下盘接触的等效弹性模量。
在遍历所有微凸体时，计算的微凸体会返回一个具体的应力值，绘制云图即可得到此岩体结构面的接触应力分布，例如可利用matlab的contourf函数，横纵坐标为微凸体的位置坐标，内容即为接触应力大小。
对于本领域技术人员而言，显然本发明不限于上述示范性实施例的细节，而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现本发明。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内，不应将权利要求中的任何附图标记视为限制所涉及的权利要求。
以上所述实施例仅表示发明的实施方式，本发明的保护范围不仅局限于上述实施例，对于本领域的技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明保护范围。

Claims (10)

1.一种评估岩体结构面接触应力分布的方法，其特征在于，包括：
S100，将结构面的上盘点云和下盘点云进行初步点云配准；
S200，将经过初步点云配准的上盘点云和下盘点云在坐标系上分区，取均值作为坐标系的横纵竖方向上方格的代表值，对上盘点云和下盘点云进行精简，并对精简后的上盘点云和下盘点云再一次进行精准的点云配准；
S300，获取精准的点云配准后的上盘点云和下盘点云的复合形貌；并按照高程矩阵分布确定微凸体的位置及几何特征；
S400，根据所述结构面的力学参数确定所述微凸体和结构面基体的刚度，以及两者的串联刚度；
S500，根据平衡方程确定所述微凸体的接触变形；
S600，根据所述微凸体的接触变形，获取当前应力状态下微凸体引起的接触面积；
S700，获取微凸体区域的总接触应力，并遍历所有微凸体，获取结构面的接触应力分布。
2.根据权利要求1所述的评估岩体结构面接触应力分布的方法，其特征在于，步骤S200包括：首先将经过初步点云配准的上盘点云和下盘点云在坐标系横纵方向上划分方格，其均值作为横纵方向上方格的代表值；再将每个方格区域内的竖方向坐标取均值，并代表这个方格区域内的高程坐标；最后获取精简后的上盘点云和下盘点云。
3.根据权利要求2所述的评估岩体结构面接触应力分布的方法，其特征在于，在坐标系横纵方向上划分方格的划分依据通过以下公式获取：
cell(ii,jj)＝{x,y,z|0.5ii≤x≤0.5(ii+1),0.5jj≤y≤0.5(jj+1)}；
式中，x、y和z分别表示为点云在坐标系横纵竖方向上的坐标，cell(ii,jj)表示为一个点云所在方格区域的名称及位置；
所述高程坐标通过以下公式获取：
Figure FDA0004202430880000021
式中，z'(ii,jj)表示为点云在方格区域内的高程坐标，Nc(ii,jj)表示为点云在方格区域点云的数目，z'表示为方格区域内所有点云的竖坐标。
4.根据权利要求3所述的评估岩体结构面接触应力分布的方法，其特征在于，步骤S300包括：
S310，将精准点云配准后的上盘点云和下盘点云的形貌组合，获取复合形貌；
S320，判断所述高程矩阵内的元素是否满足判断规则，若满足所述判断规则，则对应位置的区域形貌判定为微凸体；
S330，将步骤S320筛选出来的微凸体的顶点坐标以及在所述定点坐标正交方向相邻四个点的坐标，通过最小二乘法拟合成球体，将所述球体的半径作为微凸体的曲率半径。
5.根据权利要求4所述的评估岩体结构面接触应力分布的方法，其特征在于，所述复合形貌通过以下公式获取：
comb(x,y)＝zu(x,y)+zl(x,y)；
式中，comb(x,y)表示为复合形貌，comb表示为复合后的高程矩阵，zu表示为精准点云配准后的上盘点云，zl表示为精准点云配准后的下盘点云，x、y分别表示为点云在坐标系横纵竖方向上的坐标；
所述判断规则通过以下公式获取：
comb(ii,jj)＞max{comb(ii-1,jj),comb(ii+1,jj),comb(ii,jj+1),comb(ii,jj-1)}；
式中，max表示为最大值；comb表示为复合后的高程矩阵，(ii,jj)表示为高程矩阵的为位置；comb(ii-1,jj),comb(ii+1,jj),comb(ii,jj+1),comb(ii,jj-1)表示为复合后的高程矩阵在正交方向上相邻四个点的高程。
6.根据权利要求1所述的评估岩体结构面接触应力分布的方法，其特征在于，步骤S400包括以下步骤：
S410，根据Hertz微凸体接触理论，获取单个微凸体发生弹性变形时的刚度；
S420，根据弹性理论，获取当微凸体受到荷载时，在微凸体中心的最大接触压力，以及微凸体和结构面基体作用的部分，在接触区内的沉降；
S430，根据所述沉降获取结构面基体的刚度。
7.根据权利要求6所述的评估岩体结构面接触应力分布的方法，其特征在于，所述Hertz微凸体接触理论，通过以下公式获取：
Figure FDA0004202430880000031
式中，fan表示为第n个微凸体所受的压力，E表示为岩体的弹性模量，ρn表示为复合形貌中第n个微凸体的曲率半径，δan表示为复合形貌中第n个微凸体的变形量；
所述单个微凸体发生弹性变形时的刚度通过以下公式获取：
Figure FDA0004202430880000032
式中，kan表示为单个微凸体发生弹性变形时的刚度；
所述最大接触压力通过以下公式获取：
Figure FDA0004202430880000033
式中，P0表示为最大接触压力，FN表示为外载荷，π表示为圆周率，rb表示为微凸体与基体的接触半径；
所述接触区内的沉降通过以下公式获取：
Figure FDA0004202430880000034
式中，Uz(r)表示为接触区内的沉降，E*表示为结构面上下盘接触的等效弹性模量，r表示为微凸体与结构面基体接触范围；
所述结构面基体的刚度通过以下公式获取：
Figure FDA0004202430880000035
式中，kbn表示为结构面基体的刚度，δbn表示为结构面基体的变形量；Uz(0)表示为微凸与结构面基体接触范围为0时，其接触区的沉降；
所述微凸体和所述结构面基体的串联刚度，通过以下公式获取：
Figure FDA0004202430880000041
式中，kn表示为串联刚度，kan表示为第n个微凸体的刚度。
8.根据权利要求1所述的评估岩体结构面接触应力分布的方法，其特征在于，在步骤S500中，所述微凸体的接触变形通过以下公式获取：
Figure FDA0004202430880000042
式中，f(δn)表示为微凸体的接触变形，δn表示为总变形，ρn表示为复合形貌的微凸体的曲率半径，rb表示为微凸体与结构面基体的接触半径。
9.根据权利要求1所述的评估岩体结构面接触应力分布的方法，其特征在于，步骤S600中，所述获取当前应力状态下微凸体引起的接触面积，通过以下公式获取：
An＝πρnδan＝πρn(δn-δbn)；
式中：An表示为微凸体引起的接触面积，π表示为圆周率，ρn表示为复合形貌的微凸体的曲率半径，δan表示为复合形貌上第n个微凸体的变形，δn表示为总变形，δbn表示为结构面基体的变形量。
10.根据权利要求1所述的评估岩体结构面接触应力分布的方法，其特征在于，步骤S700中，所述微凸体区域的总接触应力，通过以下公式获取：
Figure FDA0004202430880000043
式中，σn表示为第n个微凸体区域的总接触应力，An表示为微凸体区域的接触面积，δan表示为复合形貌上第n个微凸体的变形，π表示为圆周率，f(δan)表示为复合形貌上第n个微凸体由于变形产生的接触力，ρn表示为复合形貌的微凸体的曲率半径；π表示为圆周率，ρn表示为复合形貌的微凸体的曲率半径，E*表示为结构面上下盘接触的等效弹性模量。
CN202310466317.3A 2023-04-26 2023-04-26 一种评估岩体结构面接触应力分布的方法 Pending CN116429295A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310466317.3A CN116429295A (zh) 2023-04-26 2023-04-26 一种评估岩体结构面接触应力分布的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310466317.3A CN116429295A (zh) 2023-04-26 2023-04-26 一种评估岩体结构面接触应力分布的方法

Publications (1)

Publication Number Publication Date
CN116429295A true CN116429295A (zh) 2023-07-14

Family

ID=87085323

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310466317.3A Pending CN116429295A (zh) 2023-04-26 2023-04-26 一种评估岩体结构面接触应力分布的方法

Country Status (1)

Country Link
CN (1) CN116429295A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117237403A (zh) * 2023-08-15 2023-12-15 天津大学 一种岩石的节理面吻合度确定方法及装置

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117237403A (zh) * 2023-08-15 2023-12-15 天津大学 一种岩石的节理面吻合度确定方法及装置
CN117237403B (zh) * 2023-08-15 2024-03-29 天津大学 一种岩石的节理面吻合度确定方法及装置

Similar Documents

Publication Publication Date Title
US7451637B2 (en) Method of evaluating contact characteristics, and computer product
JP5703396B2 (ja) 減数された測定点による公差評価
CN107423462B (zh) 工件考虑三维粗糙表面形貌的疲劳寿命预测方法及系统
CN110532591B (zh) 基于dic-efg联合仿真分析裂纹尖端应变场的方法
CN116429295A (zh) 一种评估岩体结构面接触应力分布的方法
CN108445188B (zh) 基于中智区间函数的岩体结构面粗糙度系数尺寸效应下边坡稳定性表达方法
CN101253385A (zh) 应变评价装置及应变评价方法
CN112818563B (zh) 一种基于摩擦接触面预估的路面抗滑性能评价方法
CN107016166A (zh) 一种基于模态应变能的新型结构刚度损伤定位方法
CN113436223B (zh) 点云数据的分割方法、装置、计算机设备和存储介质
CN115408783A (zh) 一种基于数字孪生的多层次信息装配模型的构建方法
CN106644868A (zh) 一种二维非凸形随机骨料周围界面浓度的测定方法
CN113129275A (zh) 一种基于岩土体材料数字图像三维结构表征方法
CN115935129B (zh) 一种土壤分尺度重金属浓度值确定方法和装置
CN112617817A (zh) 坐姿检测方法、装置、计算机设备及存储介质
CN117173345A (zh) 基于分形理论和蒙特卡罗方法的粗糙表面形貌模拟方法
CN107239629B (zh) 一种岩石结构面实验室合理尺寸确定的分形维数分析方法
CN113343523B (zh) 一种基于空间网格的有限元数值模拟分析方法
CN115035026A (zh) 一种基于三维点云信息的隧道爆破质量评估方法
CN110889246B (zh) 结构面抗剪区域的确定方法
Nakata et al. Experimental data of 3D printed granular material for verification of discrete element modeling simulation
CN112507578B (zh) 一种应用于压印成形仿真的基于改进的接触算法的物质点法
JP2007223504A (ja) タイヤ性能予測方法、装置、及び記憶媒体
CN113792392A (zh) 基于完整叠加微凸体的机械结合面接触刚度建模方法
CN118504283B (zh) 一种基于数字孪生的堆石骨架合理性分析系统

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination