岩土体细观结构精细化数值模拟方法分析

针对岩土材料不连续、各向异性和非线性等典型特征，本文以岩土材料的细观结构精细化建模方法为研究内容，采用离散元方法，基于计算几何理论，系统性提出一种包括二维不规则形状颗粒生成算法、三维不规则颗粒生成算法、可视化颗粒模型库管理系统开发、颗粒重叠检测以及颗粒投放算法在内的岩土体细观结构模型的方法，并将提出的方法应用于考虑块石形状的土石混合体开挖模拟中
1 绪 论

1.1 研究背景及研究意义

岩土材料是广泛分布在自然界的多相物质体系[1]，是人类赖以生存的重要物质基础。在人类社会的生产实践过程中，人们对岩土材料的工程地质性质认识不断深化，为改善居住条件，改造生活环境，或就地取材，或进行土质改良，修建了许多大型工程，其中有一些至今仍在发挥其功能。例如人们利用吐鲁番土质中砂砾和粘土胶结不易坍塌的工程地质特性，为解决干旱地区饮水灌溉问题而修建的坎儿井水利工程；人们“因地制宜，就地取材”，充分利用修建地的石头、沙砾、流沙、泥土等岩土材料修建的古代军事防御系统“万里长城”。这些工程的修建为解决古人的生产、生活以及军事等需求发挥了重要作用，对于人们深刻了解古代建造技术，传承历史文化精粹，达到古为今用的效果也有着深远意义。又例如古建筑的地基处理中，人们将黄土与白灰按一定比例掺和，分层夯实形成灰土地基[2]。现如今，这一地基处理技术经过改进后仍广泛应用于岩土工程实践中。近些年来，随着国民经济的迅速发展以及城镇化进程的加快，为改善居住环境，完善城市基础设施，解决水利、电力和能源短缺等问题，垃圾填埋场的修建、水利水电工程的修建、高边坡开挖、隧道开挖和基坑开挖等岩土工程项目数量日益增多，工程规模越来越大，所处的地质环境条件也越来越复杂。如白鹤滩水电站工程[3]、川藏铁路[4]和大柱山隧道等。这些工程中遇到的一些岩土体（如软岩、土石混合体、冻土等）由于地质成因复杂，地质环境复杂多变，岩土物理力学性质随着环境发生着时空变化，增加了工程建设中事故发生的不确定性，给工程建设带来了极大困扰。系统了解工程建设地区的岩土体类型及工程地质性质，对于减少地质环境破坏，避免和减少工程建设过程中地质灾害事故的发生，保障所建设施的后期正常运行具有重要意义。岩土材料是由颗粒、孔隙和裂隙等组合形成的多相结构系统，具有不连续、各向异性、不均匀和非线性等典型特征[5, 6]。而在传统的基于连续力学的岩土力学分析方法中，岩土材料常被假定为连续均匀各向同性的理想弹塑性材料。此外由于地质环境的复杂性，所建立的模型一定程度上都进行了简化，加之计算方法的简化，使得建立的岩土材料计算模型不能够准确反映岩土材料的应力应变关系，难以准确定位影响岩土结构稳定性的关键因素，甚至出现计算结果与实际情况严重偏离的情况，无法有效揭示岩土结构的破坏规律和机理。而一些复杂的岩土体本构模型，多是基于“现象”而建立，模型参数没有明确的物理意义，需要借助试验等手段进行数据拟合才能得到近似准确的结果。研究表明岩土材料的强度、变形和渗透性等宏观性质很大程度上受内部结构

如颗粒和孔隙的大小、分布和形态特征，颗粒的力学性质及接触方式等控制[7, 8]。以碎屑岩为例，图 1.1 所示为碎屑岩典型胶结类型示意图。根据碎屑颗粒和填隙物之间的关系，其胶结类型主要分为基底胶结、孔隙胶结、接触胶结和镶嵌胶结[9]。对于基底胶结的碎屑岩，其填隙物含量较高，颗粒之间基本不接触，其胶结强度较高，而渗透率和孔隙率低。而接触胶结类型的碎屑岩中碎屑颗粒多为点接触或线接触，呈颗粒支撑结构，其胶结强度低于基底胶结。此外，由于岩土材料赋存环境的复杂多变（如高地应力、高温、冻融循环等），其内部发生的力学、化学变化导致岩土材料的细观结构发生着时空变化，在宏观上表现为岩土材料物理力学性质的复杂性和多变性[10]。

1.2 国内外研究现状

1.2.1 岩土材料细观结构试验研究方法

岩土材料细观结构的研究方法主要包括岩土细观结构试验和数值模拟。在岩土试验方面，常规的宏观试验方法侧重于对岩土材料宏观物理力学性质的测定。自 20 世纪 70 年代以来，随着计算机技术与先进物理探测技术的发展，扫描电子显微镜 （Scanning Electron Microscopy, SEM）、CT 技术、核磁共振技术等逐渐用于对岩土细观结构的识别和获取，且其精度逐渐提高。岩土细观结构试验主要用于获取岩土体的细观结构、颗粒形态特征、颗粒间的接触等，并将其结果应用于岩土体损伤破坏细观机理研究以及细观结构模型的模拟重建。对于岩土微观结构实验中常用的技术方法总结如下。① 扫描电子显微镜（SEM）扫描电子显微镜是岩土细观结构研究中较为常用的一种研究手段。其主要通过获取颗粒表面形态特征，分析比较在不同条件下岩土体的细观结构变化，以探究岩土体破坏的细观机理。在土体细观结构研究中，许多学者采用 SEM 技术研究了不同土体（如砂土、粘土、黄土、膨胀土、污染土等）的细观结构特征[20, 21]，以及在施加物理作用（如降雨[22, 23]等）、力学作用[24]、化学作用[25-29]条件下土体细观结构特征变化。在宏观现象的细观机理揭示方面，Barden 等[24]采用 SEM 技术对风成砂土、粘土和黄土的亚稳定结构进行了研究，并探讨了不同土的施加荷载条件下的湿陷机理。陈亚明和印艳华[30]采用SEM技术分析了黄土的细观结构并以此分析了黄土湿陷性的细观机理。Baracos[31]利用 SEM 和 X 射线衍射技术从细观层面分析了不同土体剪切强度和渗透率各向异性的原因。Puppala 等[32]采用 X 射线衍射和 SEM 技术研究了钙矾石在土体中形成过程及引土体隆起的机理。Gratche 等[33]采用 SEM 技术研究了粘土的细观结构对液化能力的影响，并在此基础上分析了粘土液化机理。唐朝生等[34]使用 SEM 分析纤维表面和土体之间的界面上的相互作用，总结了 PP 纤维增强粘土强度和力学行为的细观机制。张先伟和孔令伟[35]采用扫描电镜、压汞法、

氮气吸附法对近海地区黏土中的孔隙特征进行了定量分析，建立了其与宏观物理力学性质之间的关系，并对其高灵敏性的原因进行了细观分析。

2 基于转向角理论的二维不规则颗粒生成方法

2.1 引言

目前已有的二维不规则形状生成方法主要有基于图像处理的岩土材料整体建模方法、基于图像处理的模板库构建方法、基于标准形状扩展的方法（如三角形）、泰森多边形方法、R(θ)方法以及反向傅里叶变换的方法等。在采用数字图像技术进行岩土材料细观结构模拟时，能够生成与实际颗粒形状以及岩土体真实细观结构一致的模型。然而在实际模拟中，很难获得整个研究区域内岩土体的数字图像，也很难获得用于特定研究目的的理想化岩土体截面的数字图像。基于颗粒图像信息重建生成颗粒模板库，进而生成具有各种颗粒形状和尺寸的岩土试样，弥补了采用数字图像信息直接重建岩土体细观结构的这一不足。但是建立含有各种颗粒形状的模板库需要付出巨大的努力。另外采用基于数字图像生成的颗粒模型库中的颗粒进行岩土体细观结构建模时，不能准确控制试样中颗粒的形状。基于 R(θ)以及反向傅里叶变换方法广泛应用于不规则颗粒的生成，然而这些方法不适用于非星形颗粒的生成。针对以上方法在模拟二维不规则颗粒方面的不足，本节基于 Z-R ShapeFunction 理论（Zahn 和 Roskies[144]于 1972 年提出的一种分析描述封闭多边形轮廓特征的方法），提出了新的不规则颗粒轮廓生成方法—随机转向角算法（Randomangular bend algorithm, RAB）。通过将提出的二维颗粒生成算法采用 Python 语言编程以生成单个的不规则颗粒轮廓，并从整体形状特征、圆度和表面纹理三个方面对生成的颗粒形状进行定量描述，并建立 RAB 方法中特征输入参数与颗粒形状描述符之间的定量关系，以精确控制颗粒的生成。最后以土石混合体的生成为例验证所提出模拟方法的可行性。本章研究路线如图 2.1 所示。

2.2 Z-R Shape Function 理论

形状函数 φ(l) (Z-R shape function)最早是由 Zahn 和 Roskies 于 1972 年提出的一种根据累计边长 l 和累计方向角 φ 描述多边形轮廓特征的方法。对于一个由 n条边构成的任意形状多边形（如图 2.2 所示），假定其由 n 个顶点分别为 P0, P1,…, Pn(P0 与 Pn 为同一个顶点)。其中每个顶点对应的方向角分别为 Δφ0, Δφ1,…, Δφn，对应的边长为 Δl0, Δl1, …, Δln。为构建形状函数 φ(l)，首先沿着该多边形的轮廓以某一点为起始点顺时针搜集顶点 Pi的笛卡尔坐标(xi, yi)，并基于这些顶点坐标计算多边形的边长 Δli以及角弯曲值 Δφi。其中从 Pi-1(xi-1, yi-1)到 Pi(xi, yi)的边长 Δli计算公式为

3 基于二维颗粒轮廓的三维不规则颗粒生成方法............................................33

3.1 前言 ........................................................................................................................................33

3.2 三维不规则颗粒生成算法.....................................................................................................34

3.3 三维不规则颗粒形状描述 .................................................................................................... 40

3.4 三维颗粒生成方法验证及讨论 ............................................................................................ 50

3.5 本章小结 ................................................................................................................................ 53

4 基于图形界面的岩土颗粒模型生成及管理系统构建 ................................. 55

4.1 前言........................................................................................................................................ 55

4.2 岩土颗粒模型库管理系统结构框架 .................................................................................... 55

4.3 岩土颗粒模型管理系统各模块功能 .................................................................................... 58

4.4 岩土颗粒模型库系统的应用实例 ........................................................................................ 63

5 基于离散元的岩土体精细化数值建模方法....................................................... 67

5.1 前言........................................................................................................................................ 67

5.2 基于离散元方法的岩土体细观结构建模过程 .................................................................... 67

5.3 基于颗粒簇单元的颗粒重叠检测算法 ................................................................................ 68

5.4 颗粒投放算法 ........................................................................................................................ 80

5.5 本章小结 ................................................................................................................................ 89

6 岩土体精细化模拟方法应用算例—以土石混合体边坡开挖为例....... 91

6.1 前言 ........................................................................................................................................91

6.2 边坡开挖应力分析 ................................................................................................................92

6.3 土石混合体边坡模型及模拟方案.........................................................................................93

6 岩土体精细化模拟方法应用算例—以土石混合体边坡开挖为例

6.1 前言

土石混合体边坡是一种由石块和土颗粒等组成的典型非均质边坡，在全国各地有着广泛分布，例如在西南地区残坡积边坡以及古滑坡等。由于其形成原因以及工程地质性质的复杂性与特殊性，这类边坡的变形以及失稳破坏机理与土质边坡和岩质边坡显著不同。研究土石混合体边坡失稳破坏机理，提出相应的防治措施一直是工程人员与研究人员需要面对的工程问题。在边坡稳定性分析方法方面，极限平衡法（LEM）被广泛用于边坡稳定性的定量评估[168, 169]。然而极限平衡法的结果取决于对滑移面、条间力方向及力位置等的假定，且忽略了土体的变形。有限元方法（FEM）克服了 LEM 方法在这些方面的不足，近些年来被众多研究人员应用于各种类型边坡的稳定性分析[111, 170, 171]。但是由于网格畸变引起的计算结果不稳定，使得该方法在大变形计算中受到限制，另外这种方法将土石混合体作为连续介质，不能有效反映土石混合体的结构特征以及变形不连续特性。作为替代数值方法，Lucy 和 Gingold 等于 1977 年提出的光滑粒子流方法[172, 173]（Smoothed Particle Hydrodynamics , SPH）和 Cundall 于 1971年提出的离散元方法[93]（Discrete Element Method, DEM）以及石根华和 Goodman[84]在 1984 年提出非连续变形分析方法（Discontinuous Deformation Analysis, DDA）等在模拟岩土体的大变形方面，与有限元/有限差分方法相比，显示出一些明显的优势。尽管如此，SPH 方法存在着边界处理不够精确等问题[174]。DDA 方法主要应用于岩体方面的分析计算。DEM 方法允许岩土颗粒之间的接触断开和岩土体中裂隙的产生，且其原理相对简单，便于进行计算机编程实现，能够对从多尺度对岩土体变形问题进行分析，逐渐成为一种研究复杂条件下边坡稳定性问题的强有力工具。

7 结论及展望

针对岩土材料不连续、各向异性和非线性等典型特征，本文以岩土材料的细观结构精细化建模方法为研究内容，采用离散元方法，基于计算几何理论，系统性提出一种包括二维不规则形状颗粒生成算法、三维不规则颗粒生成算法、可视化颗粒模型库管理系统开发、颗粒重叠检测以及颗粒投放算法在内的岩土体细观结构模型的方法，并将提出的方法应用于考虑块石形状的土石混合体开挖模拟中。本文的主要结论如下：

① 提出了一种基于随机转向角理论的二维不规则颗粒形状生成方法。将提出的二维颗粒生成算法采用 Python 语言编程以生成指定形状的不规则颗粒几何模型，并从整体形状特征、圆度和表面纹理三个方面对生成的颗粒形状进行了定量描述，基于统计分析结果，建立了形状描述符与算法中所需输入参数之间的定量关系，以精确控制颗粒的生成。以指定块石形状的二维土石混合体试样生成为例验证所提出方法的可行性。

② 以二维不规则颗粒生成方法（RAB 方法）为基础，提出了一种基于二维颗粒轮廓叠加的三维颗粒轮廓生成方法。通过控制层与层之间二维颗粒轮廓上点的转向角大小范围，生成由一层层二维轮廓点组成的三维颗粒点云数据。基于计算几何理论，计算点云法向量，并采用泊松重建的方法重构出三维不规则颗粒几何模型。随后针对三维颗粒，选取合适的形状描述符对生成的三维颗粒几何模型进行定量描述，建立了形状描述符与颗粒生成算法之间的经验关系。最后以土石混合体的三维模拟为例，验证了所提出方法的可行性。结果表明，所提出算法中的四个控制性输入参数（Rab、Rbc、Rl/b 和 Δφmax）能够定量控制生成的三维颗粒形状，且具有明确的几何意义，便于理解。

③ 结合本文提出的二维和三维岩土颗粒形状生成方法，开发构建了由人机交互界面、颗粒生成导入模块、颗粒分析管理模块和颗粒查询调用模块四个基本结构模块组成的基于图形界面的岩土颗粒模型库系统。采用开源编程语言 Python 中的图形界面开发库 Tkinter 开发了岩土颗粒模型库管理系统图形用户界面，采用嵌入式 SQLite 数据库作为数据库平台,进行数据的存取、增删以及查询等操作。该系统可以实现颗粒生成、第三方颗粒数据导入、颗粒形状分析、颗粒模型信息数据库构建、颗粒查询调用等功能的可视化操作。最后通过具体的实例应用展示了该系统的使用过程。