CN112258643A - 岩土边坡任意形状落石运动轨迹三维分析方法 - Google Patents

岩土边坡任意形状落石运动轨迹三维分析方法 Download PDF

Info

Publication number
CN112258643A
CN112258643A CN202010976543.2A CN202010976543A CN112258643A CN 112258643 A CN112258643 A CN 112258643A CN 202010976543 A CN202010976543 A CN 202010976543A CN 112258643 A CN112258643 A CN 112258643A
Authority
CN
China
Prior art keywords
slope
rockfall
model
falling rocks
falling
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
CN202010976543.2A
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.)
Army Engineering University of PLA
Original Assignee
Army Engineering University of PLA
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 Army Engineering University of PLA filed Critical Army Engineering University of PLA
Priority to CN202010976543.2A priority Critical patent/CN112258643A/zh
Publication of CN112258643A publication Critical patent/CN112258643A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/05—Geographic models

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Business, Economics & Management (AREA)
  • Geometry (AREA)
  • Software Systems (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Economics (AREA)
  • Human Resources & Organizations (AREA)
  • Strategic Management (AREA)
  • Computer Graphics (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Game Theory and Decision Science (AREA)
  • Development Economics (AREA)
  • Remote Sensing (AREA)
  • Marketing (AREA)
  • Operations Research (AREA)
  • Quality & Reliability (AREA)
  • Tourism & Hospitality (AREA)
  • General Business, Economics & Management (AREA)
  • Devices Affording Protection Of Roads Or Walls For Sound Insulation (AREA)

Abstract

本发明公开了一种岩土边坡任意形状落石运动轨迹三维分析方法，包括以下步骤：在系统坐标系下建立任意形状落石和坡面的三维模型；通过接触搜索算法判断每一时刻落石是否与坡面发生接触；通过接触碰撞模型计算落石与坡面发生接触时的接触力；根据接触力和重力作用建立落石的运动方程；采用时间步长法求解运动方程以获得落石的运动轨迹。本发明通过建立任意形状的落石和坡面模型，根据落石与坡面间的接触搜索算法和接触碰撞模型，建立了落石运动方程，模拟了落石沿坡面的三维运动轨迹，缩小了落石运动轨迹预测结果与实际情况的差距。

Description

岩土边坡任意形状落石运动轨迹三维分析方法
技术领域
本发明属于山区落石运动轨迹模拟分析领域，具体涉及一种岩土边坡任意形状落石 运动轨迹三维分析方法。
背景技术
落石运动轨迹的预测方法是指，根据危岩体的位置以及复杂的山区地形，采用现场 监测或计算机模拟的方法，分析落石沿坡面运动的过程，预测落石的运动轨迹和运动速度。落石运动轨迹的预测结果，决定了落石灾害的威胁范围，对山区线路工程(公路、 铁路、输油输气管道等)建设前的选线绕线设计，以及重要基础设施(房屋、桥梁、隧 道洞口等)的防护范围设计极具参考价值；落石冲击速度的预测结果，直接影响落石冲 击力的大小，可以为防护结构的极限承载力设计提供参考依据。
而当前落石运动轨迹的分析方法研究仍存在以下几点不足，主要体现在：目前模型 仍以二维为主，考虑三维6自由度的分析模型较少；以球形和特定形状的落石模型为主，未能反映落石形状的任意性和随机性的影响，而落石的数字化建模方法虽然可以详细地表征落石的几何形状，但是通过扫描的方法得到大量落石模型却是非常艰巨的工作；落 石与坡面间接触模型仍以回弹系数法为主，十分依赖于经验系数而且需要根据大量试验 确定，而计算接触力的分析模型尚无法模拟落石与坡面间的面-面接触状态；大多数模 型只能模拟落石的飞落和反弹，只有少数模型还可以模拟落石的滚动和滑动，落石多种 运动模式之间的转换尚难以准确地预测，可能导致落石运动轨迹预测结果与实际情况有 较大的差距。
发明内容
针对现有技术的不足，本发明的目的在于提供一种岩土边坡任意形状落石运动轨迹 三维分析方法，以解决现有技术中存在的落石运动轨迹预测结果与实际情况有较大差距 的问题。
为解决上述技术问题，本发明采用的技术方案为：
一种岩土边坡任意形状落石运动轨迹三维分析方法，包括以下步骤：
在系统坐标系下建立任意形状落石和坡面的三维模型；
通过接触搜索算法判断每一时刻落石是否与坡面发生接触；
通过接触碰撞模型计算落石与坡面发生接触时的接触力；
根据接触力和重力作用建立落石的运动方程；
采用时间步长法求解运动方程以获得落石的运动轨迹。
进一步的，所述落石的三维模型建立方法包括；
假设落石为刚体，对山区常见的落石形状进行量化表征；
将量化表征指标引入随机多面体生成算法中，生成形状任意且随机的多面体模型；
对所述多面体模型的表面进行三角剖分，生成由三角面元组成的落石模型。
进一步的，所述随机多面体生成算法包括：
根据基本延拓法和量化表征指标在正八面体上随机增加新的面，得到任意凹凸的多 面体；
对所述多面体进行拓扑变换得到多面体模型。
进一步的，所述三角剖分的过程如下：通过剖分方法不断地将多面体模型中最大的 三角面元剖分成两个较小的三角面元，直到所有的面元满足精度要求为止，得到由三角面元组成的落石模型。
进一步的，所述坡面的三维模型建立方法如下：
在整体坐标系中定义正交网格；
在所述正交网格的节点上输入坡面的高程信息得到坡面的点云模型；
将所述点云模型连接组成坡面的面元模型，即坡面的三维模型；
所述高程信息由自定义空间函数计算获取或地理信息系统导出的数字高程模型获 取。
进一步的，由数字高程模型获取的高程信息建立坡面的三维模型方法如下：
基于地理信息系统获取不规则的点云数据；
将不规则的点云数据转化为不规则的三角面元模型；
根据三角面元模型和正交网格将不规则的点云数据转化为规则分布的点云模型；
将规则分布的点云模型转化为规则化的三角面元模型，得到由多个三角面元组成的 坡面模型。
进一步的，通过接触搜索算法判断每一时刻落石是否与坡面发生接触包括；
根据落石面元的中心点位置获取其对应的坡面面元；
获取每个落石面元所有顶点的Z向高度和所有顶点在坡面面元上的投影点的Z向高 度；
根据顶点的Z向高度和投影点的Z向高度判断落石是否与坡面发生接触。
进一步的，所述接触碰撞模型为：
Figure BDA0002684080070000031
其中，
Figure BDA0002684080070000032
为落石质心受到坡面的接触力合力；
Figure BDA0002684080070000033
为落石质心受到坡面的总的法向阻力；
Figure BDA0002684080070000034
为落石质心受到坡面的总的切向阻力；
Figure BDA0002684080070000035
为落石质心受到坡面的总的阻 力矩；
Figure BDA0002684080070000036
为面元受到的沿坡面面元的法向接触力；
Figure BDA0002684080070000037
为面元受到的沿坡面面元的切向 接触力；
Figure BDA0002684080070000038
为面元对质心的力矩；
其中，所述面元受到的沿坡面面元的法向接触力的增量为：
Figure BDA0002684080070000039
其中，
Figure BDA00026840800700000310
是法向接触力
Figure BDA00026840800700000311
在时间步Δt内的增量，
Figure BDA00026840800700000312
是坡面法向方向的位移增量，k1、k2分别为每个面元的加载刚度和卸载刚度。
进一步的，所述落石的运动方程为：
Figure BDA00026840800700000313
Figure BDA00026840800700000314
其中，m表示落石质量，
Figure BDA00026840800700000315
表示质心的平动加速度，
Figure BDA00026840800700000316
为落石质心受到坡面的接触力合力，
Figure BDA00026840800700000317
表示重力矢量，
Figure BDA00026840800700000318
表示落石的惯性张量，
Figure BDA00026840800700000319
表示落石绕质心的转动加速 度，
Figure BDA00026840800700000320
为落石质心受到坡面的总的阻力矩。
进一步的，对所述运动方程进行求解包括：
获取落石的初始条件，所述初始条件包括落石质量、绕质心惯性张量、初始质心位置、初始质心平动速度、初始质心旋转角度、初始质心转动速度、每个面元中心点位置 和每个面元顶点位置；
根据初始条件计算t时刻落石所受合力和合力矩；
根据落石所受合力和合力矩计算时间增量Δt内落石的转角增量和质心平移增量；
在局部坐标系下，根据所述转角增量和质心平移增量计算每个面元中心点和顶点在 t+Δt时刻的坐标；
在局部坐标系下，根据落石所受合力和合力矩更新落石在t+Δt时刻的平动和转动速 度；
在局部坐标系下，更新落石在t+Δt时刻的惯性张量；
在整体坐标系下，更新所有面元中心点和顶点的坐标；
重复上述步骤直到落石运动达到终止条件。
与现有技术相比，本发明所达到的有益效果是：
(1)本发明通过建立任意形状的落石和坡面模型，根据落石与坡面间的接触搜索算法和接触碰撞模型，建立了落石运动方程，模拟了落石沿坡面的三维运动轨迹，缩小 了落石运动轨迹预测结果与实际情况的差距；
(2)提出了具有复杂形状的坡面建模方法，通过定义空间函数或从地理信息系统中导入坡面数字高程模型数据进行建模，坡面在整体坐标系下定义，而落石同时在整体 坐标系和局部坐标系下定义，更加方便地模拟了落石沿坡面的三维运动轨迹；
(3)通过接触搜索算法针对性解决了单个落石模型与大范围的坡面模型之间的接触问题，将落石与坡面模型之间的接触判定问题转换成两组三角面元之间的相交检测问题；通过将坡面模型按正交网格规则化后，可以直接检测落石面元与其所在坡面面元之 间位置关系，将传统的遍历搜索算法的计算次数由Np×Nt次降到Np次，极大地减少了计 算时间；
(4)为了反映落石任意形状的影响，根据接触刚度计算每个面元单位面积受力，相比于一些广泛应用的接触碰撞模型，该方法可以更精确地计算落石与坡面的实际接触面积、相互穿透深度和和永久变形，从而反映坡面的塑性变形和能量耗散；另外本发明 只需要根据试验确定坡面的加载和卸载刚度，不再依赖于经验性的回弹系数。
附图说明
图1为整体和局部坐标系示意图；
图2为典型落石形状示意图；
图3为任意形状落石球度和凹凸度示意图；
图4为基本延拓法示意图；
图5为典型的凸或凹多面体示意图；
图6为整体形状的拓扑变换示意图；
图7为落石表面的三角剖分示意图；
图8为坡面的建模方法示意图；
图9为典型坡面模型示意图；
图10为基于GIS的坡面模型建模方法示意图，其中，图10(a)为由GIS导出的 坡面点云模型及XOY面投影，图10(b)为不规则三角剖分模型和平面网格，图10(c) 正交规则化点云模型和平面网格，图10(d)规则三角剖分模型和平面网格。
图11为接触搜索算法示意图；
图12为落石面元与坡面面元的接触搜索和判定示意图；
图13为落石与坡面接触碰撞算法示意图；
图14为接触的两个阶段示意图；
图15为改进的双线性接触力模型示意图；
图16为落石飞落运动模式数值模拟结果与理论解对比示意图；
图17为落石反弹运动模式的数值模拟结果示意图；
图18为球形落石滚动运动模式的数值模拟结果图；
图19为球形落石沿不同坡度斜坡面滚动的数值模拟结果图；
图20为立方体落石滑动运动模式的数值模拟结果图；
图21为圆盘形落石沿不同坡度斜坡面滑动的数值模拟结果与理论解对比图；
图22为反弹与滑动运动模式之间转换图；
图23滚动与反弹运动模式之间转换图；
图24落石典型运动过程的数值模拟结果图。
具体实施方式
下面结合附图和具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本 发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。
一种岩土边坡任意形状落石运动轨迹三维分析方法，包括以下步骤：在系统坐标系 下生成任意形状的落石和坡面的三维模型；
通过接触搜索算法判断每一时刻落石是否与坡面发生接触；
通过接触碰撞模型计算落石与坡面之间的接触力(矩)；
在三维刚体动力学的框架内，根据落石与坡面的接触力(矩)和重力作用建立落石的运动方程，并采用时间步长法求解以获得落石的运动轨迹。具体步骤如下：
系统坐标系的建立
为了更加方便地描述落石与坡面的空间位置和相对运动，图1定义了两个正交坐标 系。整体坐标系G的原点设为O，而局部坐标系L的原点设为落石的质心O'。局部坐 标系通过将整体坐标系G的原点O平移至落石质心O'获得，其转换关系为
Figure BDA0002684080070000061
其中，
Figure BDA0002684080070000062
为局部坐标系L下落石表面任意一点坐标，
Figure BDA0002684080070000063
为整体坐标系G下落石表面任意一点坐标，
Figure BDA0002684080070000064
为整体坐标系G下落石质心O'的坐标。上角标L和G分别表示局 部和整体坐标系，下角标0表示落石的质心。
由于落石的飞落、反弹、滚动和滑动本质上均是落石在三维空间内平移运动和绕质 心旋转运动耦合组成的运动模式，因此，本发明采用6个自由度表征落石的运动，即质心的位移矢量
Figure BDA0002684080070000065
和转角矢量
Figure BDA0002684080070000066
另外，平动速度矢量
Figure BDA0002684080070000067
和转动速度矢量
Figure BDA0002684080070000068
分别定义为
Figure BDA0002684080070000069
和
Figure BDA00026840800700000610
的时间导数，以描述 落石运动速度的变化。
任意形状落石的建模
为了模拟山区落石形状的多样性，假设落石为刚体，不考虑落石的变形和可能发生 的劈裂。首先对山区常见的落石形状进行量化表征和分析，然后将量化表征指标引入改进的随机多面体生成算法中，生成形状任意且随机的多面体模型，最后将多面体模型表 面进行三角剖分，生成由三角面元组成的落石模型。
整体和局部形状的量化表征
在生成任意形状落石模型之前，选取一些典型的落石，观察和分析这些落石的形状 特征，如图2所示。落石的典型形状包括：近球形落石，通常由风化或雨水侵蚀而形成，具有光滑的曲面，因此滚动能力较强(图2a)；近长方体落石，通常由破碎岩体的裂缝 劈裂而产生，有平坦的表面和尖锐的棱角，所以滚动能力较差(图2b)；多面体形状落 石，这一类形状最常见，表现为落石面数在4～20左右，形状具有随机性，如图2c所示 的凸多面体形落石有很多平坦的面，图2d所示的凹多面体形落石有明显的凹面，而图 2e所示的落石有任意凹凸的表面。因此，凹凸性是落石重要的形状特征之一。下面从整 体和局部形状对落石形状进行量化表征。
为了描述落石的整体外观特征，本发明采用球度作为整体形状的量化表征指标。假 设某一任意凹凸的落石体积为Vp，三轴特征尺寸为Dx、Dy、Dz，其中Dz为最长尺寸， 如图3所示，球度Sp的计算方法为：
Figure BDA0002684080070000071
为了量化落石的局部形状特征，本发明采用凹凸度量化描述落石的凹凸性。如图3所示，假设某一椭球体与落石具有相同的特征尺寸，落石所有凸出椭球体范围的体积为Vout，落石所有凹进椭球体范围的体积为Vin。假定落石凸出和凹进的体积相等，所以有 Vout＝Vin。因此，凹凸度的量化指标Cp的计算方法为：
Figure BDA0002684080070000072
球度和凹凸度指标在一定程度上分别量化表征了落石的整体和局部形状特征，并且 可以作为模型参数引入任意形状落石的生成算法中。
任意形状落石的生成算法
根据任意形状落石生成算法，从正八面体开始，逐渐在原来的多面体基础上随机增 加新的面，直到根据整体和局部形状的量化指标(如Sp和Cp)生成任意凹凸的多面体， 具体步骤如下：
在多面体上随机增加新的面的方法称为基本延拓法，如图4所示。有两种从八面体向多面体的延拓方式，一种是向外延拓，可以使多面体向外凸起，如图4b所示，另一 种是向内延拓，使多面体向内凹陷，如图4c所示。
通过基本延拓法，可以根据落石的形状特征生成任意凹凸的多面体。多面体的形状 由延拓次数N和延拓系数Δ控制。多面体的面数取决于延拓次数N。每次延拓都会在原来的八面体上增加两个新的面。延拓系数控制落石形状的凹凸度，如图5a所示，展示 了凸十二面体，其中延拓次数N＝2，延拓系数Δ(n)＝0.20(n表示第n次延拓)。图5b 展示了凸多面体，其中N＝20，Δ(n)＝0.10，n＝1、2…20，而图5c展示了近球形凸多面体， 其中N＝20，Δ(n)＝0.20。通过比较以上例子可以发现，随着N和Δ的增加，多面体的形 状逐渐接近球形，特别是当N超过20，且Δ不超过0.20时。图5d展示了近球形的凹 多面体，其中N＝20，Δ(1)＝-0.20，Δ(n)＝0.20，n＝2、3…20，Δ(n)的变化范围是0.40。图 5e展示了任意凹凸的多面体，其中N＝20，Δ(n)在0～0.20范围内取随机值，Δ(n)的变化 范围是0.20。因此，N和Δ的取值方法会控制随机多面体的凹凸度Cp。
为了进一步控制落石的整体形状，需要在任意凹凸多面体的基础上进行拓扑变换。 如图6a所示，展示了一个近球形的落石(R＝1.41m，Sp＝1.00，Cp≈0)，图6(b-c)中落石则是从6a落石模型拓扑变换得到的，拓扑变换公式满足如下方程：
Figure BDA0002684080070000081
其中，ηx和ηy分别是x和y轴方向的比例系数，
Figure BDA0002684080070000082
和
Figure BDA0002684080070000083
是随机多面体每个顶点
Figure BDA0002684080070000084
的坐标值。图6b展示了近圆盘形的落石，其中，ηx＝0.30，ηy＝1.00，球度Sp＝0.67。而图 6c展示了近椭球形的落石，其中，ηx＝0.30，ηy＝0.30，球度Sp＝0.45。
此外，落石也可能会有一些特殊的形状，如长方体、圆柱体等，这些特殊形状很难从八面体中直接生成。因此，必须采用另一种拓扑变换方法将球形模型转换成长方体或 圆柱体模型。例如，将球形模型的所有面元均投影到其内接长方体或内接圆柱体中，如 图6d-e所示，这样就会获得新的长方体或圆柱体模型。
图5、图6中通过数值方法生成了多种形状的落石模型，与图2中展示的典型形状落石比较发现，本发明算法生成的任意形状落石模型可以从整体形状和局部凹凸性两方面，模拟大多数常见的山区落石形状。
落石表面的三角单元剖分方法
为了更加精确地模拟落石与坡面之间的面-面接触，需要将多面体的表面进一步剖 分成更精细的三角面元。
通过剖分方法不断地将最大的三角面元剖分成两个较小的三角面元，直到所有的面 元满足精度要求为止。如图7所示，首先搜索多面体的最长边PjPk，并生成PjPk的中点 E作为新的顶点；然后，三角面元PiPjPk和PjPlPk被分别剖分成两个更小的面元，如PiPjE、PiEPk、PjPlE和EPlPk；通过不断地重复这一过程，多面体的表面就可以被剖分成三角 面元，直到所有面元最长边小于自定义的长度限值dLmax为止。
假设落石的表面由Np个三角面元组成，每个顶点的坐标通过数组
Figure BDA0002684080070000085
记 录，其中下角标v表示顶点。每个三角面元的三个顶点按逆时针排列，由数组
Figure BDA0002684080070000086
)记录。除此之外，落石模型信息还包括所有三角面元的中心点，外法向量和面积，分别由数组
Figure BDA0002684080070000087
和dAi记录。
任意形状与特性的坡面模型的建立
坡面建模方法
首先在整体坐标系G中定义正交网格，然后在网格节点上输入坡面的高程坐标得到 坡面的点云模型，最后将点云模型连接组成坡面的面元模型，具体步骤介绍如下：
首先，在整体坐标系G下，X-Y平面上定义正交的格栅网格，所有的网格点用数组
Figure BDA0002684080070000088
记 录(下角标t表示坡面)，其中
Figure BDA0002684080070000091
Lx和Ly分别表示坡面 在X-Y平面上的矩形投影的长度和宽度。每个网格的尺寸可以设为一个固定值lx×ly(如图8a所示)， 网格尺寸的大小需要满足坡面模型的精度要求。
然后，确定每个网格节点的高程坐标
Figure BDA0002684080070000092
正交网格
Figure BDA0002684080070000093
将坡面的模型空间均匀划分，每 个网格节点都对应一个坡面上的三维空间点，可以通过数组
Figure BDA0002684080070000094
记录，由三维空间 点集描述的坡面模型可称为坡面的点云模型。点云的高程信息可以根据自定义的离散型二元函数
Figure BDA0002684080070000095
计算获得，或者由地理信息系统(GIS)导出的高程数据直接获得。
最后，将点云模型转换为坡面的面元模型。已知每个网格a'b'c'd'有四个节点，
Figure BDA0002684080070000096
Figure BDA0002684080070000097
如图8a所示。a'、b'、c'、d'分别是坡面上 的点a、b、c、d在X-Y平面上的投影。将点a、b、c、d依次连线可以得到空间四边形abcd，再连 接其对角线bd得到两个三角面元，如面abd和面bcd，三角面元的单位法向量可以分别用向量
Figure BDA0002684080070000098
和
Figure BDA0002684080070000099
表示，如图8c所示，所有的三角面元(可记为Ti)组成了坡面的面元模型。
由于坡面的不同区域的坡面材料可能不同，例如裸岩、覆土、植被、碎石层等，所以还需要定义坡面材料特性，可以通过数组
Figure BDA00026840800700000910
表示。坡面材料特性包括刚度kt(i,j)，摩 擦系数μt(i,j)。
综上所述，根据本发明坡面建模算法模拟坡面几何形状的前提是，获取并输入与格 栅网格点
Figure BDA00026840800700000911
对应的离散高程信息
Figure BDA00026840800700000912
而一般高程信息
Figure BDA00026840800700000913
的来源可分为两种途径：一种是通过各种自定义空间函数
Figure BDA00026840800700000914
计算获取，一般用于生成一些典型 形状的坡面模型；另一种是通过由地理信息系统(GIS)导出的数字高程模型获取，生 成的坡面模型形状更接近于实际的、形状复杂的边坡地形地貌，基于这两种途径生成坡 面模型的具体方法介绍如下。
典型坡面模型
一般地，在绕山公路或铁路建设过程中，会对部分道路沿线的山坡进行切削加工和 加固，依据山坡形状和走势选择构筑各种形式的人工边坡，如直线型、曲线型、阶梯型等。此类典型的坡面形状可以通过定义空间函数的方法模拟，以便于工程师比较分析落 石在各种形式坡面的运动特点。图9展示了几种常见的典型坡面坡型，图9a表示沿Y 轴方向倾斜的直线型坡面；图9b表示沿X和Y轴两个方向同时倾斜的直线型坡面；图 9c表示阶梯型人工坡面；图9d和图9e表示两种曲线型坡面；图9f展示的是，在沿Y 轴方向倾斜的直线型坡面上再加上坡面高度(Z轴)的随机变化
Figure BDA0002684080070000101
例如：
Figure BDA0002684080070000102
其中，
Figure BDA0002684080070000103
是线性方程，当
Figure BDA0002684080070000104
符合正态分布函数N(0,σ2)时，σ的大小可以反映坡面表面的局部起伏变化。
基于GIS的坡面模型的建立
基于地理信息系统(GIS)获取坡面的数字高程模型，其本质上是可以反映坡面形状的一组离散的点云数据，可用数组
Figure BDA0002684080070000105
表示，数据精度范围从几分米至几十米 不等。下面以某一地区道路边坡模型为例，介绍如何将由GIS导出的点云模型原始数据 转化为可用于落石运动轨迹分析的面元模型。
如图10a所示，通常情况下，基于GIS技术手段获取的点云数据在坡面上的位置是完全随机的，其在XOY面上的投影点呈现不规则分布的特点，这不利于落石与坡面模 型之间实现高效地接触判定，因此，需要预先将不规则分布的点云数据转化为坡面模型 中按正交格栅网格分布的点云数据。首先将原始不规则的点云数据转化不规则的三角面 元模型(图10b)，然后按照正交格栅网格点在对应面元上取对应的插值点，获得规则 分布的点云模型(图10c)，最后根据上述坡面建模算法将规则化的点云模型转化为规 则化的三角面元模型，即为最终的坡面模型(图10d)，具体算法介绍如下。
步骤1：将不规则的点云数据
Figure BDA0002684080070000106
转化为不规则的三角面元模型，其本质上是对三维空间点集的三角剖分。转化后的不规则三角剖分模型中每个三角面元TDj包含 三个顶点，可按逆时针排序记为
Figure BDA0002684080070000107
如图10b所示，三角面元的大小不 一，形状不规则，且空间分布不均匀。
步骤2：将不规则的点云数据转化为规则分布的点云模型，首先定义正交格栅网格点
Figure BDA0002684080070000108
假定其网格精度为lgrid(lgrid＝lx＝ly)，然后逐个计算网格点
Figure BDA0002684080070000109
对应的不规 则三角面元TDj上的投影点高度
Figure BDA00026840800700001010
获得规则分布的点云模型，用数组
Figure BDA00026840800700001011
表示，如图10c所示。
步骤3：将规则分布的点云模型转化为规则化的三角面元模型。具体方法同上，生成的坡面模型中每个三角面元Ti包含三个顶点，可按逆时针排序记为
Figure BDA00026840800700001012
对应面元的法向向量用
Figure BDA00026840800700001013
和
Figure BDA00026840800700001014
表示。观察图10d发现，转化后的坡面模型与原始 的不规则三角剖分模型的形状十分接近，但是三角面元的大小、形状和空间分布更加规 则和均匀。
落石与坡面模型的接触模型
在生成了落石和坡面模型之后，为了模拟落石沿坡面的运动过程，需要解决的关键 问题是如何判断落石与坡面之间是否发生接触，以及如何模拟落石与坡面之间的相互作 用。
接触搜索算法
接触搜索的目的是，判断落石沿坡面运动过程中的每一时刻是否与坡面接触。如果 发生接触，则根据接触碰撞模型在落石表面施加接触力，如果未发生接触，则不施加接触力。
根据落石与坡面模型建模方法，落石模型由Np个三角面元Ai组成(i＝1、2…Np)，坡面模型由Nt个三角面元Tj组成(j＝1、2…Nt)。落石与坡面均具有任意变化的复杂形状， 直接对两个任意形状几何体进行接触判定十分困难，因此，将落石与坡面间的接触判定 问题转化为两组三角面元之间的相交检测问题。对落石每个面元与坡面每个面元进行相 交检测，遍历一次需要Np×Nt次计算，十分耗时。为了提高计算效率，本发明将计算次 数降到Np次，其基本思路为，首先根据落石面元Ai的中心点
Figure BDA0002684080070000111
位置快速锁定其对应的 坡面面元Tj，然后再比较每个三角面元Ai所有顶点
Figure BDA0002684080070000112
与在Tj上的投影点
Figure BDA0002684080070000113
的高程坐标， 如果存在
Figure BDA0002684080070000114
的Z向高度低于
Figure BDA0002684080070000115
则表明落石与坡面发生了接触，如图11所示。具体算 法介绍如下：
首先，在整体坐标系内，每个三角面元Ai都有1个中心点
Figure BDA0002684080070000116
如图12a所示。中 心点
Figure BDA0002684080070000117
在X-Y平面上的投影点为
Figure BDA0002684080070000118
根据
Figure BDA0002684080070000119
的坐标
Figure BDA00026840800700001110
可以直接找到
Figure BDA00026840800700001111
所对应 的正交网格a'b'c'd'，方法为：
Figure BDA00026840800700001112
然后，由于每个网格a'b'c'd'包含两个三角面元，面元a'b'd'和面元b'c'd'，所以需要 进一步判定落石面元Ai具体位于哪个面元。如图12b所示，如果Δy≤1-Δx，那么投影点
Figure BDA00026840800700001113
位于上三角面元a'b'd'，而如果Δy>1-Δx，
Figure BDA00026840800700001114
位于下三角单元b'c'd'。
最后，落石的每个三角面元Ai有三个顶点
Figure BDA00026840800700001115
Figure BDA00026840800700001116
表示面元Ai每个顶点
Figure BDA00026840800700001117
的Z向高度，
Figure BDA0002684080070000121
表示
Figure BDA0002684080070000122
在对应的坡面面元Tj上的投影点，
Figure BDA0002684080070000123
则表示
Figure BDA0002684080070000124
的Z向高度， 如图12c所示。若Ai的任一顶点满足
Figure BDA0002684080070000125
则三角面元Ai位于坡面之下，说明该面 元Ai与坡面发生了接触；若Ai的每个顶点满足
Figure BDA0002684080070000126
则三角面元Ai位于坡面之上， 说明该面元Ai未与坡面发生接触。值得注意的是，每个面元Ai与坡面之间的接触面积 必须进行精确地计算，例如位于落石与坡面相交面上的面元Ai，部分侵入了坡面面元 Tj，则只计算并考虑面元接触部分的面积。
现有的一些广泛应用的接触搜索算法能够解决场景与场景中大量物体之间的相交 检测问题。但是这些算法花费较大代价用于计算具有复杂形状的物体之间的相交检测，而非考虑更加精细的结果，如相互穿透深度、材料的塑性和永久变形等。本发明提出的 接触搜索算法是为了针对性解决单个落石模型与大范围的坡面模型之间的接触问题而 设计的，可以将落石与坡面模型之间的接触判定问题转换成两组三角面元之间的相交检 测问题，不仅能得到上述更加精确的计算结果，而且改进后的接触搜索算法与传统的遍 历搜索方法相比极大地提高了计算效率。
接触碰撞模型
首先，面元Ai受到的接触力Fi L可以沿坡面面元Tj分解为法向接触力
Figure BDA0002684080070000127
和切向接触力
Figure BDA0002684080070000128
如图13所示。
Figure BDA0002684080070000129
的方向与坡面面元Tj的法向量
Figure BDA00026840800700001210
相同，而
Figure BDA00026840800700001211
的方向与落 石质心沿坡面面元Tj的切向速度
Figure BDA00026840800700001212
方向相反。因此面元Ai受到的接触力合力为
Figure BDA00026840800700001213
对质心O'的力矩
Figure BDA00026840800700001214
其中
Figure BDA00026840800700001215
为接触面元Ai中心到质心O'的 距离。
确定每个面元Ai受到接触力Fi L的大小。模拟落石与坡面间的面-面接触状态，假定接触力Fi L为面力，其大小取决于落石面元的面积dAi、侵入深度δi L和加载或卸载刚度。
首先，假定落石撞击垫层过程也可以分成两个阶段，在某一时刻t，当落石的位移δi L(t+Δt)＞δi L(t)时，落石对垫层处于压缩阶段或称加载阶段图14a；当 δi L(t+Δt)≤δi L(t)时，落石对垫层处于恢复阶段或卸载阶段(图14b，其中，Δt为时间 步长。
其次，引入加载刚度K1和卸载刚度K2，以分别计算加载阶段和卸载阶段球形物体间接触力，
Figure BDA00026840800700001216
其中，接触力
Figure BDA0002684080070000131
以加载刚度K1先加载至
Figure BDA0002684080070000132
然后以卸载刚度K2卸载至
Figure BDA0002684080070000133
因为 K1<K2，所以存在残余变形
Figure BDA0002684080070000134
为了在接触力模型中反映落石任意形状的影响，将根据 总的接触刚度K1和K2计算总接触力
Figure BDA00026840800700001325
改为根据每个面元Ai的接触刚度k1和k2计算 面元所受接触力
Figure BDA0002684080070000135
k1和k2可以看成是K1和K2折算到单位面积上的加卸载刚度，每 个面元Ai所受法向接触力的增量计算公式为：
Figure BDA0002684080070000136
其中，
Figure BDA0002684080070000137
是法向接触力
Figure BDA0002684080070000138
在时间步Δt内的增量，
Figure BDA00026840800700001324
是坡面法向方向的位移增量，dAi'是落石面元Ai沿坡面法向投影面积，k1和k2分别是坡面的加载和卸载刚度。一 般地，卸载段刚度大于加载段刚度，即k2＞k1。如图15所示，当
Figure BDA0002684080070000139
时，
Figure BDA00026840800700001310
说明该面元法向接触力增加，而当
Figure BDA00026840800700001311
时，
Figure BDA00026840800700001312
则该面元法向接触力减小。假 设当接触力
Figure BDA00026840800700001313
减小至零时，落石与坡面之间出现分离，此时的法向侵入深度
Figure BDA00026840800700001314
可视为 坡面的永久变形。图15中阴影面积反映了撞击过程中总能量的损失，等于坡面在加载 阶段的阻力所做的负功与卸载阶段恢复力所做的正功之和。双线性模型通过考虑塑性变 形反映能量的耗散，k1和k2是与坡面材料相关的力学参数，可以通过静力加载和卸载试 验获得。
切向接触力
Figure BDA00026840800700001315
是由Coulomb摩擦定律决定的，表达式为：
Figure BDA00026840800700001316
其中，μd为动摩擦系数，当vt L 0≠0时，落石发生滑动，且μd取决于落石与坡面的材料特性；当
Figure BDA00026840800700001317
降为零时，落石受到静摩擦力
Figure BDA00026840800700001318
静摩擦力
Figure BDA00026840800700001319
的大小是不确定的，与落 石所受外力保持平衡，方向与运动趋势相反，假设μs为最大静摩擦系数，则最大静摩擦 力等于
Figure BDA00026840800700001320
因为落石沿坡面的运动过程中很少会受到静摩擦力的作用，所以这方面 的影响可以忽略不计。
最后，将所有发生接触的面元Ai所受到的接触力和力矩累加在一起，即可得到落石 质心受到坡面的总的法向阻力
Figure BDA00026840800700001321
切向阻力
Figure BDA00026840800700001322
和阻力矩
Figure BDA00026840800700001323
计算公式为：
Figure BDA0002684080070000141
落石运动方程及求解
本发明将接触力(矩)施加在落石上建立落石的运动方程，采用时间步长法对方程进行显式求解，具体方法介绍如下。
运动方程
根据落石与坡面的接触力(矩)和重力作用，以及牛顿第二运动定律，建立落石的运动方程，
Figure BDA0002684080070000142
Figure BDA0002684080070000143
其中，m表示落石质量，
Figure BDA0002684080070000144
表示质心的平动加速度，
Figure BDA0002684080070000145
表示质心平动速度。
Figure BDA0002684080070000146
表示落石绕质心的转动加速度，而
Figure BDA0002684080070000147
表示绕质心的转动速度。
Figure BDA0002684080070000148
表示重力矢量，
Figure BDA0002684080070000149
表示重力的大小。接触力合力
Figure BDA00026840800700001410
和绕质心的合力矩
Figure BDA00026840800700001411
由公式(10)确定。
Figure BDA00026840800700001412
是落石的惯性张量，其 中
Figure BDA00026840800700001413
和
Figure BDA00026840800700001414
分别是绕局部坐标系L的x、y和z轴的转动惯量。值得注意的是，
Figure BDA00026840800700001415
Figure BDA00026840800700001416
和
Figure BDA00026840800700001417
在每个时间步内会发生改变，因为落石的形状是任意的，如图1所示。因此， 惯性张量
Figure BDA00026840800700001418
必须在每个时间步内重新计算和更新。
显式求解方法
通过对方程(11)和(12)求解，可以得到落石的运动状态，包括落石的质心位置
Figure BDA00026840800700001419
平动速度
Figure BDA00026840800700001420
和角速度
Figure BDA00026840800700001421
但是，由于落石与坡面几何形状的任意性和复杂 性导致运动方程没有理论解，因此，本发明根据时间步长法提出了运动方程的显式数值 计算方法，首先将落石连续的运动过程离散成众多时间步，根据落石在t时刻的运动状 态和合力，计算时刻t+Δt时刻的运动状态，然后在下一个时间步循环中重复上述计算， 直到落石运动达到终止条件为止。详细的求解步骤如下：
步骤1：给定落石的初始条件。这些条件包括落石质量m，绕质心惯性张量
Figure BDA00026840800700001422
初 始质心位置
Figure BDA00026840800700001423
初始质心平动速度
Figure BDA00026840800700001424
初始质心旋转角度
Figure BDA00026840800700001425
初始质心 转动速度
Figure BDA0002684080070000151
以及每个面元中心点位置
Figure BDA0002684080070000152
顶点位置
Figure BDA0002684080070000153
(为便于表达， 面元中心点
Figure BDA0002684080070000154
和顶点
Figure BDA0002684080070000155
在后文用
Figure BDA0002684080070000156
统一表示)。
步骤2：计算t时刻落石所受合力和合力矩。落石在整个运动过程中受到重力
Figure BDA0002684080070000157
作用。当与坡面接触时，落石受到接触力
Figure BDA0002684080070000158
和绕质心力矩
Figure BDA0002684080070000159
根据接触搜索算法 和接触碰撞模型，计算得到每个面元Ai上的面力Fi L(t)和力矩
Figure BDA00026840800700001510
然后，根据公式 (10)计算得到Np个面元的合力
Figure BDA00026840800700001522
和合力矩
Figure BDA00026840800700001511
步骤3：计算时间增量Δt内的转角增量
Figure BDA00026840800700001512
和质心平移增量
Figure BDA00026840800700001513
Figure BDA00026840800700001514
Figure BDA00026840800700001515
步骤4：在局部坐标系L下，通过先旋转
Figure BDA00026840800700001516
再平移
Figure BDA00026840800700001517
的方法，计算每个面元Ai中心点和顶点在t+Δt时刻的新坐标。
Figure BDA00026840800700001518
Figure BDA00026840800700001519
其中，
Figure BDA00026840800700001520
表示旋转矩阵，其表达式为，
Figure BDA00026840800700001521
步骤5：在局部坐标系L下，更新落石在t+Δt时刻的平动和转动速度。
Figure BDA0002684080070000161
Figure BDA0002684080070000162
步骤6：在局部坐标系L下，更新在t+Δt时刻的惯性张量
Figure BDA0002684080070000163
步骤7：在整体坐标系G下，更新所有面元中心点和顶点的坐标。
Figure BDA0002684080070000164
步骤8：重复步骤2至步骤7的求解过程，直到落石运动达到终止条件，例如，落 石是否停止运动，或撞到拦石墙等。如果满足其中之一的条件，则求解过程完成，如果 未满足条件，则继续求解。
步骤9：输出每个时间步的运动结果，包括空间位置
Figure BDA0002684080070000165
平动速度
Figure BDA0002684080070000166
和转动速度
Figure BDA0002684080070000167
实施例
落石典型运动模式的验证
落石的飞落运动模式一般出现于危岩从突出岩体中崩塌分离之后，或者滚动至断崖 等急剧下降的陡坡之后，或者与坡面碰撞反弹之后，是最重要的落石运动模式之一。飞落过程一般忽略空气阻力的影响，只考虑受到重力的作用，所以相较于其他运动模式更 易模拟，可以根据理论解对运动轨迹计算结果进行验证。下面分别通过二维平抛运动和 三维斜抛运动两个算例，比较理论解与数值模拟结果的运动轨迹。
假设某一球形落石的质量m＝10221kg，半径R＝1.0m，若将此落石从初始质心位置
Figure BDA0002684080070000168
Figure BDA0002684080070000169
以初始速度
Figure BDA00026840800700001610
作二维平抛运动，落石飞落过程仅受到重力作用(重力加速度记为g)，运动轨迹呈抛物线形式，其理论解为，
Figure BDA00026840800700001611
若将此落石从初始质心位置
Figure BDA00026840800700001612
以初始速度
Figure BDA00026840800700001613
作三维 斜抛运动，其理论解为，
Figure BDA00026840800700001614
根据发明提出的分析方法，并假设时间步长Δt＝0.10ms，代入上述初始条件，对运 动方程进行显式求解，得到的落石飞落运动轨迹的数值计算结果与理论解几乎一致(图16)，由此可见，当时间步长Δt足够小时，由显式求解方法所产生的误差基本可以忽 略不计。
反弹运动模式是指落石飞落至坡面时，以一定速度与其接触碰撞后，再以一定的速 度回弹的运动过程，在落石运动轨迹分析中，通常关注落石的反弹高度和反弹速度，这两者决定了落石运动轨迹的变化。下面通过球形落石与平坡面撞击反弹的典型算例对数值模拟结果进行验证。
假设某一球形落石的质量m＝10221kg，半径R＝1.0m，从初始质心位置
Figure BDA0002684080070000174
(t＝0)＝(6,6,21)作自由落体运动，落石与坡面碰撞后多次反弹。假定时间步长Δt＝0.10ms， 接触刚度k1＝6.0×109Pa/m和k2＝10.0×109Pa/m。其反弹高度和反弹速度随着每次碰撞而逐渐减小，如图17所示。
在不考虑落石的旋转和摩擦的情况下，标准的球形落石与平面正向撞击后的反弹高 度与反弹速度具有理论解。根据回弹系数法的定义，定义动能恢复性系数e2等于坡面对落石的恢复力与阻力的做功之比，所以根据考虑加载刚度k1和卸载刚度k2的双线性模 型，动能恢复性系数e2理论解的计算方法为，
Figure BDA0002684080070000171
根据算例中假定的k1和k2刚度值计算得到e2＝0.60，即为落石反弹运动的理论解。定义速度恢复性系数e等于反弹速度与入射速度之比，同理可求速度恢复性系数的理论 解
Figure BDA0002684080070000172
根据功能原理，坡面阻力做功大小等于落石入射时的动能E1，坡面恢复力做功等于 落石反射时的动能E2，由于算例中落石做自由落体运动，所以根据势能原理，动能恢复性系数e2等于反弹高度与下落高度之比，
Figure BDA0002684080070000173
因此，根据每次反弹前后的下落高度H1和反弹高度H2，可以获得动能恢复性系数e2的数值计算结果，根据每次反弹前后的入射速度v1和反射速度v2，可以获得速度恢复 性系数e的数值计算结果，数值解与理论解对比表1所示。结果表明，经过数值计算得 到的动能恢复性系数e2平均值等于0.60，速度恢复性系数e的平均值等于0.78，与理论 值之间的误差很小，验证了上述标准工况下数值计算方法的准确性。上述算例证明双线 性模型的刚度比值k1/k2可以控制落石撞击坡面过程中能量的耗散，在理想情况下落石速 度恢复性系数e等于理论解。但是如果考虑到实际情况的复杂性，如在落石随机形状、 撞击姿态、转动速度和摩擦作用等多种条件共同作用下模拟落石的反弹运动模式，其速 度恢复性系数e必然会产生一定的随机变化，这是当前基于回弹系数法的分析方法所难 以模拟的。
表1恢复性系数的数值计算结果与理论解对比
Figure BDA0002684080070000181
落石滚动是最难模拟的运动模式之一，大多数分析软件只能将落石运动假设为一系 列飞落和反弹的运动模式的集合，或者人为预设一些触发条件来假定落石进入滚动模式。 实际上，落石滚动是一种较为复杂的运动过程，从现象上观察为落石同时绕质心旋转和 沿坡面前进的复合运动，除了受到坡面法向接触力和重力的作用外，还受到滚动摩擦力的作用。因此，需要根据落石的转动速度
Figure BDA0002684080070000182
和落石与坡面之间的切向摩擦力
Figure BDA0002684080070000183
判断落 石是否发生了滚动。如图18a所示，当
Figure BDA0002684080070000184
且
Figure BDA0002684080070000185
时，落石发生了纯滚动， 根据本发明方法的模拟结果如图18b所示。
球形落石在不同坡度斜坡面顶端自由释放滚动的数值试验，通过定性比较的方法对 落石滚动模式进行验证。
假设落石质量m＝10221kg，半径R＝1.0m，分别置于坡角为θ＝30°和θ＝60°的坡面顶 端，斜坡的长度Ly均为40.0m，接触刚度k1＝6.0×109Pa/m和k2＝10.0×109Pa/m，动摩擦系数μd＝0.10，落石受到重力作用从静止开始滚动到坡面底端运动结束。如图19所示， 是球形落石沿不同坡度斜坡面滚动的数值模拟结果。结果表明，坡面越陡，落石的平动 速度和转动速度越快，与实际经验相符，证明本发明方法能够考虑落石绕质心转动和切 向摩擦力在滚动中起到的作用，验证了本发明方法在模拟落石滚动运动模式方面的可靠 性。
滑动模式也是常见的落石运动模式之一。滑动模式的判定条件如图20a所示，当
Figure BDA0002684080070000191
且
Figure BDA0002684080070000192
落石就会发生纯滑动，根据本发明方法的模拟结果如图20b 所示。
通常越扁平的落石越容易发生滑动，假设圆盘形落石质量m＝7845kg，圆盘半径 R＝1.0m，分别置于坡角为θ＝30°和θ＝60°的坡面上从静止开始滑动，接触刚度 k1＝6.0×109Pa/m和k2＝10.0×109Pa/m，动摩擦系数μd＝0.10，落石受到重力作用逐渐滑至 坡面底端。落石在坡面上的运动轨迹理论解为：
Figure BDA0002684080070000193
其中，s为落石沿坡面的滑动距离，g为重力加速度。图21中标有小正方形的实线和标有圆形的实线分别表示落石在坡度为30°和60°坡面上滑动距离s的理论解，正方形 和圆形的散点为按一定时间间隔输出的数值模拟结果。落石在坡面上滑动距离时程曲线 的数值解与理论解几乎完全一致，并且随着坡角的增加，滑动速度增大，证明本文方法 能够准确地模拟落石的滑动模式。
运动模式间的转换
落石的运动过程通常是一系列飞落、反弹、滚动和滑动四种运动模式的集合，因此必然涉及多种运动模式之间的转换，特别是反弹与滑动，滚动与反弹之间的转换。
当落石以一定速度斜向撞击坡面时，一般同时包含垂直于坡面的法向运动和平行于 坡面的切向运动，因此，落石每次反弹之前都有可能发生一小段距离的滑动，其触发条件如图22a所示，滑动距离取决于落石的冲击速度、坡面的法向刚度和摩擦系数，模拟 结果如图22b所示，图中c点和c点之间连接的弧线即为落石反弹过程中的滑动轨迹线。
对于正在滚动的落石也可能在特定的条件下发生反弹，如图23所示。此时正方体形落石由于滚动导致其与坡面发生了碰撞，而后反弹了出去，这导致部分的转动能量转 化为平动动能，所以才会发生如野外试验所观测到的反弹速度大于入射速度，回弹系数 超过1.0的情况。
上述结果表明本发明提出的分析方法能够模拟这些典型的运动模式及其之间的转 换。然后，我们从大量的数值算例中选取两组包含所有这些运动模式的典型例子，展示球形或者任意不规则形状的落石如何沿不规则的坡面进行运动。
第一个算例如图24a所示，一个球形的落石从坡面顶端自由释放，然后受重力作用开始沿斜坡面滚动。当它飞落到第二段坡面时，发生反弹再次飞到空中，最后落石继续 做飞落和反弹运动，直到到达坡面底端为止。运动过程包括滚动、飞落和反弹，没有发 生滑动，与球形落石的运动特征相一致。
第二个算例如图24b所示，一个任意形状落石由静止开始滑动而非滚动，然后飞落到第二段坡面。从坡面反弹之后，落石飞行了一小段距离，然后转换为滚动模式，直到 到达坡面底端为止。运动过程包括滑动、飞落、反弹和滚动，与球形落石的运动模式不 同。
综上所述，本文提出的落石运动轨迹分析方法，能够模拟落石的飞落、反弹、滚动和滑动等四种典型的运动模式，以及这些运动模式之间的转换，并得到了理论计算结果 的验证，具有更为可靠的物理依据，为进一步模拟落石沿坡面的运动轨迹奠定了基础。
以上仅为本发明的实施例而已，并不用于限制本发明，凡在本发明的精神和原则之 内，所做的任何修改、等同替换、改进等，均包含在申请待批的本发明的权利要求范围之内。

Claims (10)

1.一种岩土边坡任意形状落石运动轨迹三维分析方法，其特征在于，包括以下步骤：
在系统坐标系下建立任意形状落石和坡面的三维模型；
通过接触搜索算法判断每一时刻落石是否与坡面发生接触；
通过接触碰撞模型计算落石与坡面发生接触时的接触力；
根据接触力和重力作用建立落石的运动方程；
采用时间步长法求解运动方程以获得落石的运动轨迹。
2.根据权利要求1所述的一种岩土边坡任意形状落石运动轨迹三维分析方法，其特征在于，所述落石的三维模型建立方法包括；
假设落石为刚体，对山区常见的落石形状进行量化表征；
将量化表征指标引入随机多面体生成算法中，生成形状任意且随机的多面体模型；
对所述多面体模型的表面进行三角剖分，生成由三角面元组成的落石模型。
3.根据权利要求2所述的一种岩土边坡任意形状落石运动轨迹三维分析方法，其特征在于，所述随机多面体生成算法包括：
根据基本延拓法和量化表征指标在正八面体上随机增加新的面，得到任意凹凸的多面体；
对所述多面体进行拓扑变换得到多面体模型。
4.根据权利要求2所述的一种岩土边坡任意形状落石运动轨迹三维分析方法，其特征在于，所述三角剖分的过程如下：通过剖分方法不断地将多面体模型中最大的三角面元剖分成两个较小的三角面元，直到所有的面元满足精度要求为止，得到由三角面元组成的落石模型。
5.根据权利要求1所述的一种岩土边坡任意形状落石运动轨迹三维分析方法，其特征在于，所述坡面的三维模型建立方法如下：
在整体坐标系中定义正交网格；
在所述正交网格的节点上输入坡面的高程信息得到坡面的点云模型；
将所述点云模型连接组成坡面的面元模型，即坡面的三维模型；
所述高程信息由自定义空间函数计算获取或地理信息系统导出的数字高程模型获取。
6.根据权利要求5所述的一种岩土边坡任意形状落石运动轨迹三维分析方法，其特征在于，由数字高程模型获取的高程信息建立坡面的三维模型方法如下：
基于地理信息系统获取不规则的点云数据；
将不规则的点云数据转化为不规则的三角面元模型；
根据三角面元模型和正交网格将不规则的点云数据转化为规则分布的点云模型；
将规则分布的点云模型转化为规则化的三角面元模型，得到由多个三角面元组成的坡面模型。
7.根据权利要求1所述的一种岩土边坡任意形状落石运动轨迹三维分析方法，其特征在于，通过接触搜索算法判断每一时刻落石是否与坡面发生接触包括；
根据落石面元的中心点位置获取其对应的坡面面元；
获取每个落石面元所有顶点的Z向高度和所有顶点在坡面面元上的投影点的Z向高度；
根据顶点的Z向高度和投影点的Z向高度判断落石是否与坡面发生接触。
8.根据权利要求1所述的一种岩土边坡任意形状落石运动轨迹三维分析方法，其特征在于，所述接触碰撞模型为：
Figure FDA0002684080060000021
其中，
Figure FDA0002684080060000022
为落石质心受到坡面的接触力合力；
Figure FDA0002684080060000023
为落石质心受到坡面的总的法向阻力；
Figure FDA0002684080060000024
为落石质心受到坡面的总的切向阻力；
Figure FDA0002684080060000025
为落石质心受到坡面的总的阻力矩；
Figure FDA0002684080060000026
为面元受到的沿坡面面元的法向接触力；
Figure FDA0002684080060000027
为面元受到的沿坡面面元的切向接触力；
Figure FDA0002684080060000028
为面元对质心的力矩；
其中，所述面元受到的沿坡面面元的法向接触力的增量为：
Figure FDA0002684080060000029
其中，
Figure FDA00026840800600000210
是法向接触力
Figure FDA00026840800600000211
在时间步Δt内的增量，
Figure FDA00026840800600000212
是坡面法向方向的位移增量，k1、k2分别为每个面元的加载刚度和卸载刚度。
9.根据权利要求1所述的一种岩土边坡任意形状落石运动轨迹三维分析方法，其特征在于，所述落石的运动方程为：
Figure FDA00026840800600000213
Figure FDA00026840800600000214
其中，m表示落石质量，
Figure FDA00026840800600000215
表示质心的平动加速度，
Figure FDA00026840800600000216
为落石质心受到坡面的接触力合力，
Figure FDA0002684080060000031
表示重力矢量，
Figure FDA0002684080060000032
表示落石的惯性张量，
Figure FDA0002684080060000033
表示落石绕质心的转动加速度，
Figure FDA0002684080060000034
为落石质心受到坡面的总的阻力矩。
10.根据权利要求1所述的一种岩土边坡任意形状落石运动轨迹三维分析方法，其特征在于，对所述运动方程进行求解包括：
获取落石的初始条件，所述初始条件包括落石质量、绕质心惯性张量、初始质心位置、初始质心平动速度、初始质心旋转角度、初始质心转动速度、每个面元中心点位置和每个面元顶点位置；
根据初始条件计算t时刻落石所受合力和合力矩；
根据落石所受合力和合力矩计算时间增量Δt内落石的转角增量和质心平移增量；
在局部坐标系下，根据所述转角增量和质心平移增量计算每个面元中心点和顶点在t+Δt时刻的坐标；
在局部坐标系下，根据落石所受合力和合力矩更新落石在t+Δt时刻的平动和转动速度；
在局部坐标系下，更新落石在t+Δt时刻的惯性张量；
在整体坐标系下，更新所有面元中心点和顶点的坐标；
重复上述步骤直到落石运动达到终止条件。
CN202010976543.2A 2020-09-16 2020-09-16 岩土边坡任意形状落石运动轨迹三维分析方法 Pending CN112258643A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010976543.2A CN112258643A (zh) 2020-09-16 2020-09-16 岩土边坡任意形状落石运动轨迹三维分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010976543.2A CN112258643A (zh) 2020-09-16 2020-09-16 岩土边坡任意形状落石运动轨迹三维分析方法

Publications (1)

Publication Number Publication Date
CN112258643A true CN112258643A (zh) 2021-01-22

Family

ID=74231269

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010976543.2A Pending CN112258643A (zh) 2020-09-16 2020-09-16 岩土边坡任意形状落石运动轨迹三维分析方法

Country Status (1)

Country Link
CN (1) CN112258643A (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113188975A (zh) * 2021-05-07 2021-07-30 中南大学 基于图像处理技术的岩体裂隙及飞石运动分析系统及方法
CN113240735A (zh) * 2021-02-01 2021-08-10 北方工业大学 一种边坡位移活动性监测方法
CN113804166A (zh) * 2021-11-19 2021-12-17 西南交通大学 一种基于无人机视觉的落石运动参数数字化还原方法
CN116911000A (zh) * 2023-06-30 2023-10-20 中国科学院、水利部成都山地灾害与环境研究所 基于方位角的将岩块间角角接触转换为角边接触的方法
CN117743725A (zh) * 2023-11-30 2024-03-22 中国地质大学（北京） 地震作用下崩塌落石沿坡面运动轨迹计算方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103529455A (zh) * 2013-10-21 2014-01-22 中铁第四勘察设计院集团有限公司 一种基于机载激光雷达三维的危岩落石调查方法
CN107180450A (zh) * 2017-06-06 2017-09-19 广西师范学院 一种基于dem的河谷横断面形态的算法
US20180240262A1 (en) * 2015-08-17 2018-08-23 Side Effects Software Inc. Physically based simulation methods for modeling and animating two-and three-dimensional deformable objects
US20190392092A1 (en) * 2018-06-22 2019-12-26 Hao Huang Method and System for Generating Simulation Grids by Mapping a Grid from the Design Space

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103529455A (zh) * 2013-10-21 2014-01-22 中铁第四勘察设计院集团有限公司 一种基于机载激光雷达三维的危岩落石调查方法
US20180240262A1 (en) * 2015-08-17 2018-08-23 Side Effects Software Inc. Physically based simulation methods for modeling and animating two-and three-dimensional deformable objects
CN107180450A (zh) * 2017-06-06 2017-09-19 广西师范学院 一种基于dem的河谷横断面形态的算法
US20190392092A1 (en) * 2018-06-22 2019-12-26 Hao Huang Method and System for Generating Simulation Grids by Mapping a Grid from the Design Space

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
PENG YAN ET AL;: "Numerical simulation of rockfall trajectory with consideration of arbitrary shapes of falling rocks and terrain", 《COMPUTERS AND GEOTECHNICS》, pages 9 - 10 *
张蒙;王想红;徐胜华;肖冰;文化立;: "等高线约束的Delaunay三角网在土石方量计算中的应用", 地理与地理信息科学, no. 04, pages 1 - 10 *
龙满生;徐洪志;张道坤;陈同聚;: "侵蚀坡面数字高程模型重构算法研究", 计算机工程与应用, no. 36, pages 1 *

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113240735A (zh) * 2021-02-01 2021-08-10 北方工业大学 一种边坡位移活动性监测方法
CN113240735B (zh) * 2021-02-01 2023-06-23 北方工业大学 一种边坡位移活动性监测方法
CN113188975A (zh) * 2021-05-07 2021-07-30 中南大学 基于图像处理技术的岩体裂隙及飞石运动分析系统及方法
CN113188975B (zh) * 2021-05-07 2022-07-15 中南大学 基于图像处理技术的岩体裂隙及飞石运动分析系统及方法
CN113804166A (zh) * 2021-11-19 2021-12-17 西南交通大学 一种基于无人机视觉的落石运动参数数字化还原方法
CN113804166B (zh) * 2021-11-19 2022-02-08 西南交通大学 一种基于无人机视觉的落石运动参数数字化还原方法
CN116911000A (zh) * 2023-06-30 2023-10-20 中国科学院、水利部成都山地灾害与环境研究所 基于方位角的将岩块间角角接触转换为角边接触的方法
CN116911000B (zh) * 2023-06-30 2024-02-27 中国科学院、水利部成都山地灾害与环境研究所 基于方位角的将岩块间角角接触转换为角边接触的方法
CN117743725A (zh) * 2023-11-30 2024-03-22 中国地质大学（北京） 地震作用下崩塌落石沿坡面运动轨迹计算方法

Similar Documents

Publication Publication Date Title
CN112258643A (zh) 岩土边坡任意形状落石运动轨迹三维分析方法
Yan et al. Numerical simulation of rockfall trajectory with consideration of arbitrary shapes of falling rocks and terrain
Lu et al. Critical assessment of two approaches for evaluating contacts between super-quadric shaped particles in DEM simulations
Chen et al. Numerical simulation in rockfall analysis: a close comparison of 2-D and 3-D DDA
Lan et al. RockFall analyst: A GIS extension for three-dimensional and spatially distributed rockfall hazard modeling
Bourrier et al. Rockfall rebound: comparison of detailed field experiments and alternative modelling approaches
Markl et al. Powder layer deposition algorithm for additive manufacturing simulations
CN112052587B (zh) 砂土垫层的三维细观离散体模型密实方法
CN111737871A (zh) 一种结合岩土材料特性的落石轨迹三维预测分析方法
Coevoet et al. Adaptive merging for rigid body simulation
Garcia et al. Comparison of full-scale rockfall tests with 3D complex-shaped discrete element simulations
Lu et al. Modelling rockfall dynamics using (convex) non-smooth mechanics
Fu et al. Reproduction method of rockfall geologic hazards based on oblique photography and three-dimensional discontinuous deformation analysis
Lee Developments in large scale discrete element simulations with polyhedral particles
Chen et al. Discrete element simulation for polyhedral granular particles
Zheng et al. Investigating the role of earthquakes on the stability of dangerous rock masses and rockfall dynamics
Frery et al. Stochastic particle packing with specified granulometry and porosity
Zeng et al. A momentum‐based deformation system for granular material
Merrell et al. Constraint-based model synthesis
Gavrilova et al. Collision detection optimization in a multi-particle system
Feng et al. A GPU based Hybrid Material point and Discrete element method (MPDEM) algorithm and validation
CN115600449B (zh) 基于细长柔性物件组成的大规模颗粒系统的模拟预测方法
Numada et al. Material point method for run-out analysis of earthquake-induced long-traveling soil flows
Knowles et al. DEM modelling of 3D polyhedra with applications to gabion rockfall barriers
CN109492285B (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