MatDEM离散元：岩土体振动信号分布式光纤监测和离散元建模研究

─ ☕ ─

导读：MatDEM作为一款国产的离散元模拟软件，凭借其高效计算、功能完善和自主可控等特点，在岩土工程领域展现了强大的应用潜力。通过不断的研发和推广，MatDEM有望为我国在岩土工程领域的数值分析提供更加先进和可靠的工具。

本文是MatDEM官方团队毕业生闵寅通同学关于《岩土体振动信号分布式光纤监测和离散元建模研究》对建模过程与模拟过程中思考和总结，希望能引发读者朋友的共鸣。诚邀读者朋友积极报名6月15日19时MatDEM官方团队成员将在仿真秀官网和APP带来《岩土体振动信号分布式光纤监测和离散元建模研究》和《星球车爬坡过程轮壤作用多体动力学—离散元数值模拟》公开课。详情见后文：

一、建模过程

二维模型的尺寸为60× 5 m，颗粒的平均半径为0.02m，模型的颗粒总数为19.80万。颗粒生成后施加重力进行堆积，平衡颗粒之间的作用力。随后在土体上表面增加压力板，后续重锤夯击和平板夯的振动作用都直接作用于压力板再传递至土体。在模型中压力板的作用可视作为“薄膜”，在起到防止颗粒飞溅的作用的同时，对振源作用于土体的影响小。在离散元模拟中过强的瞬时外力会导致颗粒的飞溅，使得模拟计算量不收敛。

Fig. 4-1Two-dimensional discrete element model

2、材料赋值

在完成土体模型堆积后需要构建DAS模型的两个重要组件：振源和光纤。对于振源，试验中振源分为铁球下落和平板夯振动两类，因而模拟中也要构造两种不同的振源模型。铁球下落振源：试验中所用的铁球的直接分别为80mm，120mm和150mm，分别对应2kg，5kg和10kg的铁球。平板夯振源：平板夯振动是其底部的平板受到发动机的偏心距影响从而产生振动，试验中平板夯底部尺寸为400mm × 400 mm × 100 mm。在二维的模型中振源的尺寸由三维转化为二维，铁球转为铁圆盘，长方体平板夯转为长方形的平板夯，并且对质量进行折算，得到二维模拟下的质量：

式中，V2d是二维圆盘的体积，V3d是三维球体的体积，m3d是三维振源的质量，位于二维模型中振源距光纤水平距离0m（正上方）。

试验所用的单模光纤由石英裸纤和外部的紧包套等组成，当外力作用光纤时会产生应变，随后由光电接收器获得相关数据。光纤的直径为2mm，而大部分材质都为涂料、保护套等，中间裸芯的半径极小。若要在离散元模拟中构建相同结构的光纤模型，则会造成颗粒数量多，计算效率下降的后果。因此为提高计算效率，本文使用半径为0.01m的弹性Clump模型来组成光纤，设置过程和原理如下：

• ① 设定光纤的坐标，可以根据工况任意改变光纤的形态和位置，本文中的光纤的铺设为水平铺设（图4-1）。通过MatDEM自带的make3DCurve函数，带入XYZ坐标参数和颗粒半径得到光纤结构体。

• ② 将颗粒用弹性Clump的方式聚合到一起，其颗粒的重叠度为0.5即有一半的重叠度，此时颗粒之间的球心距为0.01m。光纤位于土体埋深0.5m位置，位于模型中间，删去模型中光纤位置处的颗粒，将光纤“埋入”土体内。

• ③将光纤的三个坐标方向的自由度锁定，将颗粒间摩擦系数设为零，并对模型进行平衡计算，让颗粒能够与光纤相接触，减少模型中土体与光纤的空隙。随后分别撤销XYZ三个方向的自由度锁定，连接土体颗粒，并消散由于自由度锁定取消后光纤与土体颗粒之间的预压力。

• ④光纤的轴向应变和相对应变计算。公式（3-5）给出了DAS的相位差信号与应变率之间的理论关系，因而可以简化复杂的光学解调过程，通过记录应变率随时间的变化来反映模拟的振动信号，通过模型中设置的光纤结构体计算其发生的轴向应变。

首先根据光纤相邻颗粒直接的坐标，可以计算出相邻颗粒之间的初始距离，如图4-2(a)所示：

式中x，y，z为颗粒的坐标。当光纤受到外力作用后发生位移，其球心距发生改变（图4-2(b)），此时的距离用坐标差用公式（4-2）计算可得到。其发生的距离变化

如下所示：

式中，和是相邻颗粒的半径，Ci是相邻颗粒之间的重叠量。通过公式（4-3）计算得到每个颗粒之间发生的距离变化，除以其之间初始距离，就能得到相邻颗粒之间的应变：

对于每一个颗粒而言，其记录的轴向应变数据为左右两侧颗粒间相邻应变的平均值:

公式（4-5）即为每个颗粒所在每一时刻所记录的应变量大小。将得到的相对应变作差并除以时间间隔，得到光纤相对应变变化率，可通过转化公式（3-5）得到试验结果物理量纲（相变率）。

Fig. 4-2 (a)Axial distance between adjacent particles of optical fiber (b) Axial relativestrain calculation of optical fiber

试验时光纤埋入土体，光纤表面会与土体耦合，振动作用力传递过程中需要对介质表面进行求解。本文所用的离散元接触模型为线弹性体接触模型，为简化模拟复杂程度，不进行光纤土体耦合的接触模型开发。

3、振动信号模拟

模拟中需要实现两种形式的振动信号，平板夯振动与铁球下落。平板夯的振动原理和第二章中振动动力作用模型相同，即通过施加周期变化的激振力来实现，其振动参数与试验相同：振动频率100Hz，振幅2mm，压力板质量100kg。具体原理在第二章已经详细阐述，不在此赘述。

为节省计算效率，在离散元模拟中不实现铁球在空中下落的阶段，因而铁球下落的作用可抽象为小球与地面接触时的状态计算。通过将铁球砸至地面的状态计算出来并进行赋值的方式实现该振动信号的模拟，其中铁球改变的物理量为下落时的速度v：

式中，g为重力加速度，h为小球下落的高度。

公式（4-2）～（4-5）给出了离散元模拟中记录振动信号数据的原理，在模拟中设置保存数据的时间为0.5ms，对应的采样频率为2000Hz和试验相同。在模拟开始后会施加动力作用，并记录开始后十秒的光纤结构体的相对应变。综上，模拟最终保存的数据是整条光纤上所有颗粒在10s内，2000Hz采样频率下的轴向相对应变数据。

二、模拟过程相关

在试验中，DAS的标距为1m，每一个通道的数据信号代表的是1m的光纤振动信号，但是在上述建模过程中记录光纤的分辨率为0.01m，因此需要对1m长度上同时刻的100个数据进行积分，最终得到和试验一致的通道应变率数据。本文进行不同质量、下落高度的铁球下落模拟和平板夯模拟，振动参数与试验保持一致，具体模型参数如下表4-1所示。

取下落高度为1m，铁球质量为5kg的模拟信号与试验对比，如图4-5(a) 和(b) 所示：通道信号的模拟结果与试验的结果相似，通道信号振幅从中心向两边衰减，表明了离散元模型的合理性和有效性。

取该工况下中心通道的数据，对试验的结果进行相变率到应变率的转换，使得模拟和试验的量纲相同（图4-5(c) (d)）：其波形特征都是单个峰值信号并快速衰减，表明模拟的有效性。但其信号衰减时间约为1s，高于试验中振动信号的衰减时长0.25～0.3s，该现象是因为边界阻尼层未能完全实现对边界反射的吸收，机械波多次反射作用于光纤上，使的光纤信号衰减值背景噪音的时间增加。

图4-5 下落高度为1 m，铁球质量为5 kg，振源距离0 m: (a) 二维模拟通道信号 (b) 试验通道信号 (c) 二维模拟中心通道 (d) 试验中心通道

Fig. 4-5 Falling height is 1 m, iron ball mass is 5kg, vibration source distance is 0 m: (a) two-dimensional simulation channelsignal (b) test channel signal (c) two-dimensional simulation center channel(d) test center channel

三、2024MatDEM国产离散元线上报告会

为了帮助读者朋友更好的理解《岩土体振动信号分布式光纤监测和离散元建模研究》，促进MatDEM用户技术交流。6月15日19时，2024MatDEM国产离散元用户案例线上分享会将邀请MatDEM官方团队毕业生闵寅通与曹政同学分别带来《岩土体振动信号分布式光纤监测和离散元建模研究》和《星球车爬坡过程轮壤作用多体动力学—离散元数值模拟》公开课。以下是直播安排：

MatDEM离散元：岩土体振动信号分布式光纤监测与星球车爬坡过程轮壤作用模拟-仿真秀直播

扫码观看直播和回放