力学季刊

Parameter	Value
Particle radii r/mm	3~7
Normal stiffnesses k_n/(N/m)	3×10⁶
Tangential stiffnesses k_t/(N/m)	2.4×10⁶
Particle density ρ/(kg/m³)	2 650
Friction angle φ/rad	0.5

Parameter	Value
Particle radii r/mm	3~7
Normal stiffnesses k_n/(N/m)	3×10⁶
Tangential stiffnesses k_t/(N/m)	2.4×10⁶
Particle density ρ/(kg/m³)	2 650
Friction angle φ/rad	0.5

颗粒材料多尺度力学行为模拟及孔隙分形表征

PDF下载

仵仲阳 , 付龙龙 , 周顺华

力学季刊 | 2025,46(3): 649-661

收起

力学季刊 | 2025, 46(3): 649-661

颗粒材料多尺度力学行为模拟及孔隙分形表征

全屏

仵仲阳, 付龙龙, 周顺华

作者信息

同济大学　上海市轨道交通结构耐久与系统安全重点实验室，上海　201804

仵仲阳，博士生．研究方向：岩土颗粒材料的多尺度模拟．E-mail：wuzy@tongji.edu.cn

通讯作者:

付龙龙，博士，副教授．研究方向：交通岩土工程．E-mail：longlongfu@tongji.edu.cn

Modelling of Multi-Scale Mechanical Behaviours and Pore Fractal Characterization of Granular Materials

Zhongyang WU, Longlong FU, Shunhua ZHOU

Affiliations

Shanghai Key Laboratory of Rail Infrastructure Durability and System Safety, Tongii University, Shanghai 201804, China

出版时间: 2025-09-25 doi: 10.15959/j.cnki.0254-0053.2025.03.007

文章导航

摘要

收起

岩土颗粒材料的力学响应呈现多尺度特性，开展有限元法（Finite Element Method，FEM）和离散元法（Discrete Element Method，DEM）耦合的多尺度模拟可在反映多尺度响应的同时兼顾计算效率．为此，基于国产高性能离散元软件MatDEM开发了GPU并行FEM-DEM耦合计算程序，并结合颗粒孔隙分形特征对计算参数和结果进行分析．首先，采用多重分形理论研究了颗粒孔隙空间分布特征，确定了关键分形指标．其次，通过单元和双轴压缩试验验证了FEM-DEM耦合程序的可靠性．最后，基于双轴压缩试验研究了不同位置处表征元（Representative Volume Element，RVE）的细观响应．结果表明，RVE内孔隙空间分布具有多重分形特征，当颗粒数量超过400时，孔隙空间分布的自相似性确保了RVE内颗粒集合均匀化后输出应力-应变响应的稳定性．容量维数D₀和奇异指数α₀由于表征孔隙分布的平均信息而与RVE体应变具有线性相关性，可作为反映颗粒材料空间特征复杂性的简洁内变量．研究结果为分析工程尺度颗粒堆积体的宏-细观力学联系提供了一种方法探索．

关键词

岩土颗粒材料 / 多尺度 / 多重分形 / 有限元法 / MatDEM

Abstract

收起

The mechanical response of geotechnical granular materials exhibits multi-scale characteristics. The multi-scale simulations coupling the Finite Element Method (FEM) and Discrete Element Method (DEM) can effectively capture the multi-scale responses while maintaining high computational efficiency. A GPU-parallel FEM-DEM coupling code was developed based on the high-performance DEM software MatDEM, with computational parameters and results analyzed in conjunction with granular pore fractal characteristics. Firstly, multi-fractal theory was employed to investigate the spatial distribution characteristics of pores, identifying key fractal indices. Subsequently, the reliability of the FEM-DEM coupling code was verified through single element tests and biaxial compression tests. Finally, the meso-scopic responses of Representative Volume Elements (RVEs) at different locations were investigated based on biaxial compression tests. Results demonstrate that the pore spatial distribution within RVE exhibits multi-fractal characteristics. When the particle quantity exceeds 400, the self-similarity of pore spatial distribution ensures the stability of stress-strain responses output after homogenization of granular assemblies within RVE. The capacity dimension D₀ and singularity index α₀, which characterize the average information of pore distribution, show linear correlations with RVE volumetric strain. These indices can serve as internal variables reflecting the complexity of granular material spatial characteristics. This study provides an exploration for analyzing the macro-meso mechanical relationships in engineering-scale granular deposits.

Key words

geotechnical granular materials / multi-scale / multi-fractal / finite element method / MatDEM

引用本文

仵仲阳, 付龙龙, 周顺华. 颗粒材料多尺度力学行为模拟及孔隙分形表征. 力学季刊, 2025 , 46 (3) : 649 -661 . DOI: 10.15959/j.cnki.0254-0053.2025.03.007

Zhongyang WU, Longlong FU, Shunhua ZHOU. Modelling of Multi-Scale Mechanical Behaviours and Pore Fractal Characterization of Granular Materials[J]. Chinese Quarterly of Mechanics, 2025 , 46 (3) : 649 -661 . DOI: 10.15959/j.cnki.0254-0053.2025.03.007

正文

收起

岩土颗粒材料（如粗粒土）因具有压实性能好、抗剪强度高等特点，被广泛用作铁路路基填料、道砟等关键承力层^[1]．在荷载作用下岩土颗粒材料表现出复杂的宏观力学特性，如宏观指标量测结果的离散性^[2]、高频振动下的流化^[3]等．现有研究表明，这些宏观现象本质上取决于颗粒材料中接触组构^[4]、力链^[5]等多样细观结构的演化．因此，建立岩土颗粒材料细观结构与宏观力学行为的关联，始终是岩土力学领域研究的热点与难点^[6]．

近年来，离散元法^[7]（Discrete Element Method，DEM）被广泛应用于揭示颗粒材料的宏-细观力学参数之间的联系，成为研究其多尺度力学行为的有效工具^[1]．然而，要实现大尺度岩土工程的数值模拟，即便在算力极大提升的当下，DEM也因计算效率低而难以胜任^[8]．因此，一种行之有效的方法是对颗粒材料进行均匀化处理后在连续介质力学框架下进行计算，其中FEM（Finite Element Method）-DEM耦合方法结合FEM在宏观尺度上的成熟高效性与DEM在细观尺度上的完善准确性建立颗粒材料宏-细观尺度间的联系，受到了国内外学者的广泛关注^[6,9]．基于FEM平台Escript和DEM平台Yade，Guo等^[10]开发了颗粒材料的分级多尺度模拟方法，并将其应用于三轴试验^[11]、边坡^[12]、挡土墙和基础^[13]等问题．基于FEM平台Lagamine，Desrues等^[14]开发了二阶梯度增强FEM×DEM框架，并探讨了数值解的变化性及非唯一性^[15]．刘其鹏^[16]在宏观上采用Cosserat连续体构建了FEM-DEM耦合多尺度模型并将其应用于颗粒破碎问题^[17]．以上研究表明，FEM-DEM耦合多尺度方法能够准确表征颗粒材料的多尺度力学行为．同时可以发现，该方法在每一个FEM积分点位置嵌入一个颗粒集合体表征元（RVE）进行DEM计算，导致计算成本随FEM积分点数量的增加而快速增长^[11]．因此以上FEM-DEM耦合程序主要依赖基于分布式内存架构MPI（Message Passing Interface）的CPU并行化进行加速，对于工程尺度问题，通常需要在多核服务器或超级计算机平台上运行^[12,14]．而GPU相较CPU具有并行计算能力强、性价比高等特点，逐渐被用于DEM并行计算^[18]，因此亟需开发一种基于GPU并行的颗粒材料FEM-DEM耦合多尺度计算程序^[19]．

岩土颗粒材料由颗粒和颗粒间孔隙构成，其中细观孔隙结构深刻影响颗粒材料的宏观物理力学特性，如渗透性^[12]、应力波传播^[20]等，因此国内外学者采用Voronoi分割、Delaunay分割和图像分割等多种手段定量表征细观孔隙结构和孔隙拓扑形态^[21]．然而颗粒材料中孔隙的形状和大小各异、数量众多、尺度上跨越多个数量级，以上基于欧氏几何的方法难以完整刻画颗粒材料孔隙的复杂性^[22]，而分形/多重分形理论^[23]以其自相似性与尺度不变性特征，在量化复杂不规则结构方面展现出独特优势，自创立以来被广泛应用于岩石孔隙^[24]、掌纹纹理^[22]等复杂几何形态的精细化表征．现有颗粒材料FEM-DEM耦合多尺度模拟中通常将孔隙率作为初始输入条件，主要表征接触组构^[10]、力链^[19]、孔隙比^[11]等细观特征在加载过程中的演化规律，鲜有量化孔隙结构的指标，进而较少探究细观孔隙结构特征与颗粒材料宏观响应的关系．

本文首先基于国产高性能离散元软件MatDEM开发了RVE颗粒运动的GPU并行DEM计算程序，并根据RVE内颗粒孔隙空间分布多重分形指标α₀/D₀演化确定了RVE的合理尺寸，从而实现了岩土颗粒材料FEM-DEM耦合的高效多尺度模拟．然后通过单元压缩和双轴压缩试验验证所开发程序的有效性，为后期进行大尺度工程问题的FEM-DEM耦合多尺度计算奠定了基础．最后基于双轴压缩试验研究了颗粒材料在加载过程中不同位置处RVE的应力-应变响应和孔隙结构分形特征的演化规律，为建立岩土颗粒材料细观结构与宏观力学行为的定量化关系提供新的思路．

1　岩土颗粒材料FEM-DEM耦合多尺度模拟方法

收起

本文构建的岩土颗粒材料FEM-DEM耦合多尺度模拟方法在宏观上采用FEM算法，将边值问题的连续几何域离散化为FEM网格以求解控制方程，在细观上构建颗粒集合体表征元（RVE）镶嵌于FEM网格中的高斯积分点位置，以高斯点应变确定RVE边界条件，采用GPU并行技术进行RVE内颗粒运动的DEM模拟，并将基于均匀化方法确定的RVE应力响应传递至FEM中进行全局求解．如图1所示，本文多尺度方法不使用唯象本构模型模拟颗粒材料应力-应变响应，能够合理表征颗粒材料的非线性、状态相关特征和多尺度效应．

1.1　有限元法

有限元中准静态力学问题的平衡方程为

(1)

式中，σ为应力张量，b为体积力向量．在FEM求解过程中，将几何域（Ω）离散化为合适的网格（如图1所示），并采用虚功原理将式（1）转化为以下线性方程组

(2)

式中，u为未知节点位移向量，f^ext为由边界条件引起的等效节点力向量，K为全局刚度矩阵

(3)

式中，B为应变矩阵，D为弹性矩阵．对于非线性颗粒材料响应，采用Newton-Raphson迭代法进行求解^[10]，该方法在每一次迭代中，使用切线模量矩阵D_t代替弹性矩阵D并更新应力，以最小化残余力R为收敛条件

(4)

本文采用MATLAB语言编写了以上FEM过程的计算代码，其中切线模量矩阵D_t和应力σ通过高斯积分点位置处颗粒集合体RVE的DEM计算获得，从而代替了传统唯象本构模型．

1.2　表征元GPU并行化离散元计算

图1所示宏观几何域离散化后生成的FEM网格中存在大量高斯积分点，每一个高斯点位置对应一个颗粒集合体RVE需要进行DEM计算，本文基于国产高性能离散元平台MatDEM^[18]实现了RVE的GPU并行化计算．为了满足宏-细观尺度转换过程中的Hill-Mandel能量准则^[25]，RVE采用周期性边界条件^[26]，即位于RVE边界处的颗粒通过虚拟颗粒与对面边界附近颗粒进行接触，以颗粒直接接触消除边界效应．RVE内颗粒接触相互作用以线弹性模型表征，同时遵从摩尔-库伦破坏准则^[10]．在本文开发的RVE并行DEM计算程序中，首先采用静态分配策略为每一个RVE分配一个GPU的CUDA（Compute Unified Device Architecture）核心，同时考虑周期性边界条件，基于窗口法实现了所有RVE在GPU上同时进行真实与虚拟颗粒接触检测^[18]，从而获得了所有RVE颗粒的邻居矩阵．根据颗粒邻居矩阵和颗粒属性（坐标、半径等）计算颗粒相对位移矩阵，并与颗粒接触刚度矩阵相乘获得颗粒受力矩阵．然后对颗粒受力矩阵求和计算颗粒所受合力，并结合颗粒质量矩阵计算颗粒加速度矩阵，最后结合颗粒速度矩阵更新颗粒坐标^[18]．以上矩阵运算过程可在MatDEM中调用GPU进行并行计算，详细过程可参考文献[18]．

1.3　表征元均匀化

表征元（RVE）内颗粒集合体的应力张量σ通过Love公式计算^[27]

(5)

式中，V为RVE总体积，N_c为RVE内颗粒接触点的数量，f^c为两接触颗粒的接触力向量，d^c为连接两接触颗粒中心的枝向量（以压为正）．对于本文二维多尺度计算过程，平面外方向的应力和应变均设定为0．基于均匀应变场假设估计RVE内颗粒集合体的切线模量矩阵^[28]

(6)

式中，k_n和k_t分别为颗粒间法向及切向等效刚度，n^c和t^c分别为接触点处法向和切向方向上的单位向量．

1.4　FEM-DEM耦合多尺度计算过程

本文构建的岩土颗粒材料FEM-DEM耦合多尺度模拟方法的详细计算过程如下

（1）输入数据：将FEM和DEM计算所需信息输入至耦合程序中，其中FEM信息包括网格信息（节点坐标等）和边界条件（位移/应力边界）等，DEM信息包括颗粒单元性质（粒径，单元法向/切向刚度，摩擦系数和单元密度等），RVE内颗粒数量和RVE内颗粒集合体的初始孔隙比等．

（2）进行FEM预处理：根据采用的单元类型（本文采用4节点4高斯点四边形等参元）和输入的网格信息计算单元形函数、雅可比矩阵和应变矩阵B，同时确定高斯点数量．

（3）进行DEM预处理：采用周期性边界条件，制备出满足给定孔隙比条件的DEM颗粒集合体，并根据高斯点数量使用制备的试样建立GPU并行化DEM模型．

（4）进行全局加载，对每一加载步进行如下流程

（a）根据FEM边界条件计算当前加载步的等效节点力f^ext．

（b）使用式（6）计算每一个高斯点位置处RVE的切线模量矩阵D_t．

（c）使用式（3）计算单元刚度矩阵并集成至全局刚度矩阵K中．

（d）求解式（2）获得单元节点位移，并根据节点位移计算出每一个高斯点位置处的应变．

（e）根据高斯点应变更新RVE边界位置，生成边界处颗粒对应的虚拟颗粒，将所有RVE的真实与虚拟颗粒信息输入至GPU的不同CUDA核心中同时进行颗粒接触检测，然后计算颗粒合力并更新颗粒位置，从而实现所有RVE的并行化DEM计算．

（f）使用式（5）计算出每一个RVE的应力σ．

（g）根据高斯点处RVE的应力使用式（4）计算FEM内的残余力R，并进行误差判断

(7)

（h）当计算获得的误差大于容许误差（TOL_f，本文取0.005^[10]）时，以残余力R代替步骤（a）中的等效节点力f^ext，并继续迭代计算；反之则跳出本次Newton-Raphson迭代过程．

（5）当所有加载步完成后，保存数据并结束程序．

2　基于孔隙分形特性的表征元尺寸确定

收起

本文构建的颗粒材料多尺度模拟方法中以RVE内颗粒集合体经均匀化后获得的应力-应变响应代替唯象本构模型，因此确定合理的RVE尺寸成为后续多尺度力学行为分析的关键^[10]．前人研究表明^[10]，RVE内应包含足够的颗粒数量从而保证能够充分表征颗粒细观结构并精确模拟材料的整体力学响应；同时，RVE内颗粒数量不宜过多以避免加载过程中的局部化现象．为了确定合理的RVE尺寸，本文构建了11个RVE模型，分别包含10²、12²、14²、16²、18²、20²、22²、24²、26²、28²、30²个颗粒，每个RVE模型分别采用5个随机数生成颗粒．颗粒细观参数列于表1．对于每个RVE试样，先将其压缩至初始孔隙比e₀=0.18（密实状态），后采用p₀=100 kPa的压力进行等压固结，从而获得了55个具有相同初始状态但尺寸和随机数不同的RVE试样．

2.1　颗粒孔隙空间分布多重分形分析

本文引入多重分形理论量化RVE内颗粒集合体孔隙在空间中的分布特征．首先将RVE试样通过图像处理获得其孔隙空间分布的二值化图像（分辨率取0.3 mm）^[21]，然后采用计盒法^[23]计算RVE孔隙图像的多重分形分维数．如图2所示，观测盒子大小ε_f通过以下公式确定

(8)

式中，k为正整数（本文取[1，8]^[24]），L为孔隙图像大小．

对于第i个大小为ε_f的盒子，其孔隙局部密度为

(9)

式中，N_i（ε_f）为第i个盒子内所包含的像素数量，N_T为整个孔隙图像中的所有像素数量，p_i（ε_f）称为孔隙质量密度．不同观测尺度ε_f下孔隙质量密度满足如下关系^[22]

(10)

式中，α_i称为奇异指数，体现区域孔隙生成的概率大小．具有相同α的盒子数量N（α）与ε_f满足如下关系^[23]

(11)

式中，f（α）为具有相同奇异指数α的盒子集的分形维数，α~f（α）称为奇异谱，可通过以下公式计算^[29]

(12)

(13)

(14)

式中，q为矩序（本文取[-10，10]^[24]）．通过信息论可以获得广义分维数D_q^[30]

(15)

图3给出了包含不同颗粒数量（10²、16²、20²、26²、30²）的RVE孔隙空间分布广义分维谱q~D_q和奇异谱α~f（α），从图中可以看出，q~D_q谱为非对称Sigma形单调递减曲线，α~f（α）谱呈非对称凸抛物线形状，表明RVE内颗粒孔隙空间分布满足多重分形特征^[31]．

矩序q=0时，可获得容量维数D₀和奇异谱峰值对应的奇异指数α₀，它们是量化孔隙空间分布非均匀性的重要参数．前人研究表明^[31]，α₀/D₀表征孔隙空间分布的缩放指数范围，反映孔隙结构的自相似特征．图4给出了55个RVE试样的α₀/D₀值，从图中可以看出：当RVE内颗粒少于20²个时，α₀/D₀值随颗粒数量增加剧烈波动，表明RVE内颗粒孔隙空间分布结构差异明显；当RVE内颗粒超过20²个时，α₀/D₀值随颗粒数量增加而趋于稳定，表明RVE内颗粒孔隙空间分布满足自相似性，即可以认为包含20²、22²、24²、26²、28²和30²个颗粒的RVE内具有基本一致的孔隙空间分布特征，预示着20²是RVE宏观行为稳定的最小颗粒数．

2.2　表征元应力-应变曲线

选择5种颗粒数量（10²、16²、20²、26²、30²）的RVE试样进行双轴加载（z方向加载，x方向保持σ_x=p₀=100 kPa），获得的应力-应变曲线如图5所示．从图中可以看出：相较于包含10²和16²个颗粒的RVE试样，包含20²个颗粒的RVE试样能够与26²和30²个颗粒RVE试样产生基本一致的宏观应力-应变响应．前人研究中通常将不同颗粒数下应力-应变曲线达到稳定作为判断RVE最小颗粒数的依据^[10]，图5结果表明本文从孔隙空间分布分形自相似性的角度确定选择20²作为RVE的最小颗粒数是合理有效的，也说明了孔隙空间分布分形指标α₀/D₀表征的颗粒材料细观结构能够有效关联宏观应力-应变响应．值得一提的是，基于孔隙分形特征的结果也与前人基于颗粒接触法向分布的研究结果一致^[10]．结合DEM计算效率，本文确定20²个颗粒数量作为RVE的合理大小，并用于后续颗粒材料的多尺度模拟过程．

3　FEM-DEM耦合多尺度程序验证

收起

3.1　单元压缩试验

采用本文FEM-DEM耦合多尺度计算程序进行单元压缩试验，边界条件如图6（a）所示，单元包含4个节点和4个高斯点，每个高斯点对应一个包含400个颗粒的RVE，颗粒细观参数与表1一致．图6（b）和（c）给出了单元压缩试验的应力-应变曲线，从图中可以看出：多尺度程序计算结果不仅与纯DEM模拟结果基本一致，且与文献[10]中结果较为吻合．此外，图6中应力-应变曲线展现了颗粒材料在加载过程中的非线性、应变软化和剪胀现象，验证了本文多尺度程序在模拟岩土颗粒材料力学行为方面的合理性．

3.2　双轴压缩试验

采用本文FEM-DEM耦合多尺度计算程序进行双轴压缩试验，如图7（a）所示，长和高分别为50 mm和100 mm的试样被离散化为150个FEM单元，单元类型采用4节点4高斯点四边形等参元，因此共有600个节点和600个高斯点．每个高斯点对应的RVE包含400个颗粒，因此共有24万个DEM颗粒．颗粒细观参数与表1一致．图7（b）给出了双轴压缩试验获得的应力-应变曲线，从图中可以看出：本文多尺度程序获得的应力-应变响应与纯DEM模拟结果基本一致，同时也与文献[10]中结果较为吻合，验证了本文多尺度程序在模拟岩土颗粒材料力学行为方面的可靠性．比较图6（b）和图7（b）发现，达到峰值应力状态后，相较于单元压缩试验获得的应力-应变曲线（图6（b）），图7（b）中的应力-应变曲线更为光滑．进一步将图7（b）中的应力-应变曲线划分为Ⅰ弹性状态、Ⅱ峰值状态和Ⅲ临界状态三个阶段用于后续进一步分析．

值得一提的是，本文单元和双轴压缩试验获得的应力-应变结果与文献[10]结果的差异，主要来自本文所采用MatDEM平台与文献[10]采用Yade平台在初始阶段生成RVE试样时颗粒位置分布和排列方式的差异，后续将进行岩土颗粒材料的室内试验，并将试验结果与本文开发的耦合计算程序进行对比验证．

4　双轴压缩条件下颗粒材料的宏细观响应

收起

4.1　宏观变形特征

图8给出了峰值应力状态下（2%轴向应变）的试样变形云图，从图中可以看出：x和z方向的位移场呈现非对称特征并向右下角方向倾斜，试样剪应变云图呈现出明显的贯穿剪切带，展现了颗粒材料在荷载作用下的内在非均匀性．值得一提的是，传统FEM在模拟颗粒材料时需人为添加缺陷点才能模拟应变局部化现象，本文构建的FEM-DEM耦合多尺度模拟方法不需添加任何缺陷即可合理表征颗粒材料的内在各向异性，这是因为RVE内的颗粒集合体天然具有颗粒材料的非共轴特性^[10]．

4.2　细观特征演化规律

选择如图7（a）所示的三个不同位置处表征元（RVE51，RVE339和RVE429），其中RVE51位于剪切带以外区域，RVE339位于剪切带内，RVE429位于剪切带边缘位置．图9给出了这三个RVE在双轴压缩过程中应力、体应变、容量维数D₀和奇异指数α₀演化规律，从图中可以看出：位于剪切带外的RVE51在加载过程中仅发生剪缩，偏应力逐渐增加；位于剪切带区域的RVE339和RVE429在加载过程中先剪缩后剪胀，并在进入峰值状态后（阶段Ⅱ和Ⅲ）产生应力软化现象；进入阶段Ⅲ临界状态后，位于剪切带边缘的RVE429体应变超过剪切带内的RVE339，导致其偏应力趋近于0，表明试样在阶段Ⅲ临界状态时的宏观响应主要受剪切带边缘RVE的应力状态影响；位于剪切带外的RVE51在加载过程中容量维数D₀和奇异指数α₀逐渐减小；位于剪切带区域的RVE339和RVE429在加载过程中容量维数D₀和奇异指数α₀先减小后增大，同时位于剪切带边缘的RVE429进入阶段Ⅲ临界状态后，容量维数D₀和奇异指数α₀逐渐超过剪切带内的RVE339．对比图9（b）、（c）和（d）可以发现，容量维数D₀和奇异指数α₀与体应变曲线具有一致的演化规律，即RVE发生剪缩时D₀和α₀减小，发生剪胀时D₀和α₀增大．

图10分别给出了三个RVE在双轴压缩过程中容量维数D₀和奇异指数α₀与体应变的关系，从图中可以看出：容量维数D₀和奇异指数α₀均与体应变满足线性关系，表明量化孔隙分布平均程度的D₀和量化孔隙分布集中程度的α₀都能有效关联RVE的宏观体应变响应．此外，剪切带区域的RVE339和RVE429的线性拟合直线斜率基本一致，并与非剪切带区域的RVE51不同，表明容量维数D₀和奇异指数α₀与体应变的关联性主要受RVE应力状态的影响．

值得一提的是，本文采用的RVE内颗粒细观参数（表1）是真实岩土颗粒级配和形状的适当简化，导致RVE内形成的颗粒孔隙大小和分布较为均匀（如图2所示），同时容量维数D₀和奇异指数α₀均为表征孔隙分布平均信息的分形指标，因此在本文相对简单的孔隙形态中D₀和α₀主要通过揭示孔隙占比（类似于孔隙比）而能够与体应变建立良好的线性关系．这意味着分形参数D₀和α₀可作为反映颗粒材料空间特征复杂性的简洁内变量：颗粒形状、堆积形式对颗粒材料力学行为有直接影响，但因具有空间复杂性而难以统一有效量化，分形参数则提供了一种量化的思路．此外，对于真实岩土颗粒材料中的宽级配、复杂颗粒形状和复杂应力路径问题，由于多重分形理论在表征复杂孔隙形态方面的优势，本文提出的基于孔隙空间分布自相似性来确定RVE合理尺寸的思路同样适用，而此时孔隙分形指标与RVE宏观响应的关联性需要进一步研究．

4.3　表征元GPU并行化离散元计算程序效率

为了验证本文1.2节开发的RVE并行DEM计算程序的效率，采用包含3种颗粒数量的RVE（10²、20²和30²），并对不同RVE数量的DEM试样进行双轴加载，分别计算串行模拟与并行模拟10⁴次DEM迭代步所需要的时间，并定义串行模拟与并行模拟时间之比为加速系数．采用AMD 5900X和NVIDIA GeForce GTX 3080 12 GB计算平台，获得了加速系数与RVE数量的关系．从图11中可以看出：对于包含10²和20²个颗粒的RVE，加速系数随RVE数量增加近似线性增长；对于包含30²个颗粒的RVE，加速系数随RVE数量增加先线性增长后趋于稳定．这是由于对于12 GB显存的NVIDIA GPU，MatDEM软件进行高效模拟的颗粒数量为40~80万^[18]．稳定状态下本文开发的GPU并行DEM计算程序较串行计算可实现约150倍的计算速度提升，相较于前人基于CPU并行的加速方案^[10,14]，本文开发的岩土颗粒材料FEM-DEM耦合多尺度计算程序在实际大尺度岩土工程中具有更广阔的应用前景．

5　结论

收起

本文针对岩土颗粒材料的多尺度力学特性，构建了GPU并行FEM-DEM耦合多尺度模拟方法，并探究了表征元（RVE）内颗粒孔隙空间分布的多重分形特征．本文主要结论如下

（1）基于MatDEM平台灵活的GPU调用和强大的矩阵运算能力实现了RVE的GPU并行化计算，所开发的FEM-DEM耦合多尺度计算程序能够实现岩土颗粒材料多尺度力学行为的高效模拟，在所采用的计算平台中，并行计算较串行计算的计算速度提升约150倍，展现出多尺度计算程序在大尺度工程适用性方面的显著优势．

（2）RVE内颗粒孔隙空间分布满足多重分形特征，通过多重分形指标α₀/D₀量化孔隙空间分布自相似性，得出少于400个颗粒的RVE无法满足孔隙空间分布自相似性而难以产生稳定的应力-应变响应，并通过以往直接分析不同颗粒数RVE的应力-应变响应来确定合理颗粒数的方法证实了上述思路的有效性．

（3）容量维数D₀和奇异指数α₀分别量化孔隙分布的平均程度和集中程度，由于其反映RVE内孔隙分布的平均信息而与体应变具有高度线性相关性，RVE剪缩时D₀和α₀减小，剪胀时D₀和α₀增大，线性拟合系数主要受剪切带内、外RVE应力状态的影响．

基于上述研究，后续可考虑岩土颗粒材料的真实级配特性、复杂形状特征和复杂应力路径，系统开展孔隙结构分形特征在FEM-DEM耦合多尺度计算方法中的应用探究，并结合室内试验与多尺度仿真结果交叉验证，以研究实际工程中广泛存在的多尺度多物理场耦合问题．

基金

收起

国家自然科学基金(52378457)
国家留学基金委资助(CSC202206260139)

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

[1]

黄志刚, 王轩, 傅力, 等. 加载速率和摩擦系数对颗粒材料系统剪切强度的影响研究[J]. 力学季刊, 2024, 45(4): 1032-1042.

Huang

Z G

, Wang

, Fu

, et al. Influence of loading rate and friction coefficient on the shear strength of granular material system[J]. Chinese Quarterly of Mechanics, 2024, 45(4): 1032-1042 (in Chinese).

[2]

孙思齐, 肖军华, 刘志勇, 等. 高频振动下岩土颗粒材料动力行为及细观机理[J]. 力学季刊, 2023, 44(4): 835-844.

Sun

S Q

, Xiao

J H

, Liu

Z Y

, et al. Dynamic behavior of granular materials under high frequency loading and the meso-mechanism[J]. Chinese Quarterly of Mechanics, 2023, 44(4): 835-844 (in Chinese).

[3]

肖军华, 白英琦, 张骁, 等. 考虑应力波透反射作用的分层颗粒材料细观动力响应分析[J]. 力学季刊, 2023, 44(3): 620-632.

Xiao

J H

, Bai

Y Q

, Zhang

, et al. Analysis of mesoscopic dynamic response for layered granular materials considering transmission and reflection of stress wave[J]. Chinese Quarterly of Mechanics, 2023, 44(3): 620-632 (in Chinese).

[4]

Zhou

S H

, Jiang

H B

, Fu

L L

, et al. Experimental study on deformation and strength characteristics of granular soil-structure interface under coupled monotonic shear and vibration using a modified direct shear apparatus[J]. Acta Geotechnica, 2023, 18: 2899-2913.

[5]

L L

, Zhou

S H

, Guo

P J

, et al. Induced force chain anisotropy of cohesionless granular materials during biaxial compression[J]. Granular Matter, 2019, 21: 1-16.

[6]

吴焕然, 刘汉龙, 赵吉东, 等. 高孔隙率砂岩中破坏模式演化的多尺度分析[J]. 岩土工程学报, 2020, 42(12): 2222-2229.

H R

, Liu

H L

, Zhao

J D

, et al. Multiscale analyses of failure pattern transition in high-porosity sandstones[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(12): 2222-2229 (in Chinese).

[7]

Cundall

P A

, Strack

O D L

. A discrete numerical model for granular assemblies[J]. Géotechnique, 1979, 29(1): 47-65.

[8]

付龙龙, 叶伟涛, 周顺华, 等. 高频低幅振动下颗粒材料抗剪强度衰减细观机理研究[J]. 中国科学: 技术科学, 2023, 53: 1687-1698.

L L

, Ye

W T

, Zhou

S H

, et al. Shear resistance reduction of granular materials subjected to high-frequency and low-intensity vibration: A micro-mechanical study[J]. Scientia Sinica Technologica, 2023, 53: 1687-1698 (in Chinese).

[9]

Nitka

, Tejchman

. A two-scale numerical approach to granular systems[J]. Archives of Civil Engineering, 2011, 57(3): 313-330.

[10]

Guo

, Zhao

. A coupled FEM/DEM approach for hierarchical multiscale modelling of granular media[J]. International Journal for Numerical Methods in Engineering, 2014, 99(11): 789-818.

[11]

Guo

, Zhao

. 3D multiscale modeling of strain localization in granular media[J]. Computers and Geotechnics, 2016, 80: 360-372.

[12]

, Yang

Z X

, Guo

, et al. Multiscale modeling of seepage-induced suffusion and slope failure using a coupled FEM-DEM approach[J]. Computer Methods in Applied Mechanics and Engineering, 2022, 398: 115177.

[13]

Guo

, Zhao

. Multiscale insights into classical geomechanics problems[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2016, 40(3): 367-390.

[14]

Desrues

, Argilaga

, Caillerie

, et al. From discrete to continuum modelling of boundary value problems in geomechanics: An integrated FEM-DEM approach[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2019, 43(5): 919-955.

[15]

Nguyen

T K

, Vo

T T

, Nguyen

N H T

, et al. Variability and loss of uniqueness of numerical solutions in FEM×DEM modeling with second gradient enhancement[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2024, 48(9): 2381-2407.

[16]

刘其鹏. 基于平均场理论的颗粒材料离散颗粒集合-Cosserat连续体模型多尺度模拟[D]. 大连: 大连理工大学, 2010.

Liu

Q P

. Multi-scale modeling of discrete particle assembly-Cosserat continuum model for granular materials based on the average-field theory[D]. Dalian: Dalian University of Technology, 2010 (in Chinese).

[17]

王增会. 基于多尺度方法的颗粒材料破碎行为与损伤—愈合—塑性表征研究[D]. 大连: 大连理工大学, 2018.

Wang

Z H

. Study on the particles breakage and damage-healing-plasticity characterization for granular materials based on multiscale approach[D]. Dalian: Dalian University of Technology, 2018 (in Chinese).

[18]

刘春, 乐天呈, 施斌, 等. 颗粒离散元法工程应用的三大问题探讨[J]. 岩石力学与工程学报, 2020, 39(6): 1142-1152.

Liu

, Le

T C

, Shi

, et al. Discussion on three major problems of engineering application of the particle discrete element method[J]. Chinese Journal of Rock Mechanics and Engineering, 2020, 39(6): 1142-1152 (in Chinese).

[19]

Zhou

, Xu

W J

, Lubbe

. Multi-scale mechanics of sand based on FEM-DEM coupling method[J]. Powder Technology, 2021, 380: 394-407.

[20]

Seo

, Luscher

D J

, Scovel

, et al. Mesoscale simulation of granular materials under weak shock compaction-pore size distribution effects[J]. Journal of Applied Physics, 2024, 136: 155902.

[21]

Zhang

, Zhao

, et al. Three-dimensional Voronoi analysis of realistic grain packing: An XCT assisted set Voronoi tessellation framework[J]. Powder Technology, 2021, 379: 251-264.

[22]

李彤. 多重分形原理及其若干应用[D]. 北京: 北京交通大学, 2008.

. Multifractal theory and some applications[D]. Beijing: Beijing Jiaotong University, 2008 (in Chinese).

[23]

Halsey

T C

, Jensen

M H

, Kadanof

L P

, et al. Fractal measures and their singularities: the characterization of strange sets[J]. Physical Review A, 1986, 33(2): 1141-1151.

[24]

Xia

, Cai

, Perfect

, et al. Fractal dimension, lacunarity and succolarity analyses on CT images of reservoir rocks for permeability prediction[J]. Journal of Hydrology, 2019, 579: 124198.

[25]

Miehe

, Dettmar

, Zäh

. Homogenization and two-scale simulations of granular materials for different microstructural constraints[J]. International Journal for Numerical Methods in Engineering, 2010, 83: 1206-1236.

[26]

Miehe

, Dettmar

. A framework for micro-macro transitions in periodic particle aggregates of granular materials[J]. Computer Methods in Applied Mechanics and Engineering, 2004, 193: 225-256.

[27]

Christoffersen

, Menrabadi

M M

, Nemat

. A micromechanical description of granular material behavior[J]. Journal of Applied Mechanics, 1981, 48: 339-344.

[28]

Kruyt

, Rothenburg

. Statistical theories for the elastic moduli of two-dimensional assemblies of granular materials[J]. International Journal of Engineering Science, 1998, 36: 1127-1142.

[29]

Rényi

. On a new axiomatic theory of probability[J]. Acta Mathematica Hungarica, 1995, 6(3): 285-335.

[30]

Chhabra

, Jensen

R V

. Direct determination of the f(α) singularity spectrum[J]. Physical Review Letters, 1989, 62(12): 1327-1330.

[31]

Gould

D J

, Vadakkan

T J

, Poche

R A

, et al. Multifractal and lacunarity analysis of microvascular morphology and remodeling[J]. Microcirculation, 2011, 18(2): 136-151.

2025年第46卷第3期

PDF下载

187

引用本文

BibTeX

文章信息

doi: 10.15959/j.cnki.0254-0053.2025.03.007

接收时间：2025-05-11
首发时间：2026-03-24
出版时间：2025-09-25

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2025-05-11

基金

国家自然科学基金(52378457)

国家留学基金委资助(CSC202206260139)

作者信息

同济大学　上海市轨道交通结构耐久与系统安全重点实验室，上海　201804

通讯作者:

付龙龙，博士，副教授．研究方向：交通岩土工程．E-mail：longlongfu@tongji.edu.cn

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/lxjk/CN/10.15959/j.cnki.0254-0053.2025.03.007

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT