Loading [MathJax]/jax/element/mml/optable/GreekAndCoptic.js
  • 全国中文核心期刊
  • 中国科技核心期刊
  • 美国工程索引(EI)收录期刊
  • Scopus数据库收录期刊

考虑黏土片不规则形状的黏土离散元模拟方法

申志福, 张栩银, 高峰, 王志华, 高洪梅

申志福, 张栩银, 高峰, 王志华, 高洪梅. 考虑黏土片不规则形状的黏土离散元模拟方法[J]. 岩土工程学报, 2022, 44(9): 1654-1662. DOI: 10.11779/CJGE202209010
引用本文: 申志福, 张栩银, 高峰, 王志华, 高洪梅. 考虑黏土片不规则形状的黏土离散元模拟方法[J]. 岩土工程学报, 2022, 44(9): 1654-1662. DOI: 10.11779/CJGE202209010
SHEN Zhi-fu, ZHANG Xu-yin, GAO Feng, WANG Zhi-hua, GAO Hong-mei. Discrete element method for clay considering irregular planar shape of clay platelets[J]. Chinese Journal of Geotechnical Engineering, 2022, 44(9): 1654-1662. DOI: 10.11779/CJGE202209010
Citation: SHEN Zhi-fu, ZHANG Xu-yin, GAO Feng, WANG Zhi-hua, GAO Hong-mei. Discrete element method for clay considering irregular planar shape of clay platelets[J]. Chinese Journal of Geotechnical Engineering, 2022, 44(9): 1654-1662. DOI: 10.11779/CJGE202209010

考虑黏土片不规则形状的黏土离散元模拟方法  English Version

基金项目: 

江苏省科技厅自然科学基金青年基金项目 BK20190667

国家自然科学基金青年基金项目 51908284

江苏省高等学校自然科学研究面上项目 19KJB560015

江苏省高等学校自然科学研究重大项目 18KJA560002

江苏省高校“青蓝工程”中青年学术带头人资助项目 

详细信息
    作者简介:

    申志福(1988—),男,四川绵阳人,工学博士,副教授,主要从事土的宏微观土力学特性与本构关系研究。E-mail: zhifu.shen@njtech.edu.cn

    通讯作者:

    高洪梅, E-mail: hongmei54@163.com

  • 中图分类号: TU431; TU470

Discrete element method for clay considering irregular planar shape of clay platelets

  • 摘要: 离散元模拟是研究黏土复杂宏观力学特性微观机理的有力工具。当前黏土离散元模拟方法中,黏土片一般假设为矩形、椭圆形等规则的平片状,与真实的黏土片形态不符,也尚无文献报道黏土片形态对黏土力学特性的影响规律。提出了不规则凸多边形黏土片间双电层排斥作用和范德华吸引作用的计算方法,发展了相应的微观接触模型;探讨了黏土片间相对排列对相互作用的影响,据此总结了大幅提高模拟效率的途径,最终形成了基于不规则平片状颗粒的黏土离散元模拟方法。模拟表明,此方法能很好再现黏土的一维压缩特性,黏土片的平面形状、厚度对黏土一维压缩曲线和微观组构演化都有明显影响,有必要在未来的离散元模拟中考虑黏土片不规则形状的影响。
    Abstract: The discrete element method (DEM) is a powerful tool to study the microscopic mechanics behind the complicated macroscopic mechanical behaviors of clay. In the current DEM simulations of clay, the planar shape of clay platelets is commonly assumed to be regular, such as rectangle and ellipse, which are quite different from the true planar shapes of clay. There is no research report about the effects of shape of clay platelets on the mechanical behaviors of clay. In this study, the methods to calculate the Van der Waals force and the double-layer repulsive force between two convex-shaped clay platelets are proposed, and an inter-clay platelet contact model is developed. Then, the effects of relative alignment of two clay platelets on the interactions are studied, and the approaches to significantly improve the simulation efficiency are summarized accordingly. The DEM simulation for clay based on irregular convex planar shape of clay platelets is thus developed. The simulated results indicate that the developed DEM simulation scheme can well reproduce the mechanical behaviors of clay in oedometer tests. The planar shape and thickness of clay platelets have significant effects on the one-dimensional compressive curves and the evolution of platelet arrangement fabric. It is necessary to incorporate the irregular planar shape of clay platelets in future DEM simulation of clay.
  • 黏土复杂的宏观力学特性源自其复杂的微细观结构及其在外力作用下的变化过程。不同放大倍数的电镜扫描图像反映出黏土微细观结构的多层级特征[1-3],软黏土受力变形过程中黏土片间边–边、边–面不稳定接触趋向于转变为面–面接触,进而导致大孔隙结构受损坍塌[4-6]。然而,电镜扫描图像仅能反映某一时刻、某一极小区域的黏土微细观结构,且微观测试数量有限,无法全面反映其完整演变过程,这些不足在短时间内尚难以通过现有技术突破。

    离散单元法将土体视为大量离散微观单元的集合体,可方便地连续、全面揭示土体微细观结构及其演变规律,有助于揭示土体复杂宏观力学特性的微细观机理。早期学者采用二维离散元模拟黏土,将黏土片简化为线状单元,能定性反映黏土的主要宏观力学特性[7-9]。随着计算能力的增强,将黏土片简化为矩形、椭圆形平片状单元的三维离散元法成为发展的主流[10-15]。也有部分研究将黏土视为球体集合体(可认为对应黏土微观结构中的团簇集合)用于模拟黏土颗粒沉积[16]、失水干缩[17]等物理过程。然而,实际的黏土片形状远比上述理想假定复杂得多,目前尚无文献系统探讨黏土片形态对黏土宏微观力学特性的影响。

    黏土片之间存在机械相互作用和物理–化学相互作用(包括范德华相互作用(引力)和双电层相互作用(斥力))。其中,双电层斥力和机械作用计算较为成熟。而对于范德华引力,Ebrahimi等[12]采用Gay-Berne势函数描述黏土片间范德华作用,该方法数学形式简单,计算效率高,但其将黏土片的形状限定为薄椭球体、且势函数对各种排列方式的黏土片的适用性也未进行系统验证;Yao等[10]用有限尺寸薄板与无限大薄板间范德华引力(存在解析解)近似作为两有限尺寸正方形黏土片间的范德华力,存在较大误差;Guo等[14]采用规则排列的球体团簇模拟黏土片,将球体之间的范德华引力之和(存在解析解)作为两黏土片间范德华引力,该方法忽略了球体间空隙部分(特别是两黏土片表面最近点附近空隙部分)对范德华引力的重要贡献,将产生较大误差。

    综上,当前离散元模拟中对黏土片形状假定过于理想化,范德华引力计算方法又与黏土片形状假定密切关联,使得离散元模拟黏土方法本身不具有普适性,其准确性也受限,也限制了离散元方法本身应有优势的充分发挥。为此,本文在平衡计算效率与准确性的双重需求下,提出一种基于不规则平片状颗粒的黏土离散元模拟方法,通过算例论证了其必要性、可靠性和高效性。

    将黏土视为由黏土片堆叠而成的散粒体材料,黏土片可简化为由不规则凸多边形沿厚度方向垂直拉伸而成的薄片状多面体。因而本文研究适用于空间翘曲不明显的黏土矿物(如高岭石)构成的黏土片。本文研究对象为饱和黏土,虽在离散元数值模拟中不直接模拟水,但以下双电层排斥作用是基于黏土片周围完全填充水溶液的情形得到的,相当于间接反映了孔隙水的存在。

    与粗粒土不同,黏土片间存在两种特殊的物理–化学相互作用:双电层排斥作用和范德华吸引作用。这两种作用与黏土片的形状、相对排列都密切相关,本节介绍其计算方法。

    图 1(a)所示,对任意排列的两黏土片1和2,其单位法向量分别为n1n2。定义中平面为两黏土片夹角(锐角)的角平分面。以中平面法线方向为z轴,以矢量n2×n1方向为y轴,按照右手定则建立局部坐标系oxyz,该局部坐标系原点可在中平面上任意位置。黏土片1的下表面与黏土片2的上表面构成一对倾斜排列的带电平板,两带电平板在中平面上的投影相交部分如图 1(b)中实线包围的区域,称为交错区域。该区域为双电层排斥作用的有效范围,交错区域外的黏土片因相互“错开”,对双电层作用贡献不计[10, 15]

    图  1  黏土片间双电层排斥作用计算方法示意图
    Figure  1.  Schematics of calculation method for double-layer repulsive interaction between two clay platelets

    图 1(b)(c)中,将双电层排斥作用有效范围沿局部坐标系的x轴方向等间距划分为若干宽度为w的条带,并将有效作用范围近似等效为若干矩形条带的集合。采用Yao等[10]的简化方法,将两倾斜排列的带电平板间双电层排斥作用近似为若干对平行板之间的双电层排斥作用之和。以矩形条带i为例,其x方向中心位置处对应的y向重叠范围为条带长度li,如图 1(b)所示。在图 1(c)所示的oxz平面投影上,两带电平板在条带i处的z向间距为hi。根据Yao等[10]的工作,矩形条带i处的双电层斥力为

    Fe,i=2nkTliw[coshϕ(hi)1], (1)

    式中,n为带电平板间溶液(即黏土的孔隙水)中的离子浓度,k=1.38×10-23 J/K为玻尔兹曼常数,T为开氏温度,ϕ(hi)为间距为hi的两平行板的中平面归一化电势。式(1)中Fe,i方向垂直于中平面。

    图 2为通过数值方法求解泊松-玻尔兹曼方程得到的不同黏土片表面电势ψ0下,中平面归一化电势ϕ与归一化距离ζ的关系[10]。此处,归一化电势定义为

    ϕ=veψkT, (2)
    图  2  归一化中平面电势与归一化距离的关系
    Figure  2.  Relationship between normalized mid-plane electric potential and normalized distance

    式中,v为离子带电荷数,e=1.6×10-19C为单电子电荷数,ψ为中平面电势。

    归一化距离定义为

    \zeta = \frac{{{h_i}}}{2}\sqrt {\frac{{8{{\rm{ \mathsf{ π} }}}n{e^2}{v^2}}}{{\varepsilon {\varepsilon _0}kT}}} , (3)

    式中, \varepsilon 为水的相对介电常数, {\varepsilon _0} =8.8542×10-12 C/ (V·m)为真空的介电常数。

    根据图 2,通过插值可获得某表面电势下,间距为hi的两平行条带间中平面归一化电势 \phi ({h_i}) ,并根据式(1)得到条带间双电层斥力。黏土片1受到的双电层斥力和相应的力矩在局部坐标系oxyz中可表示为

    {\boldsymbol{F}}_{\text{e}, 1}\text{=}\left[0, 0, {\displaystyle \sum\limits_{i=1}^{{N}_{x}}{F}_{\text{e}, i}}\right] , (4)
    {\boldsymbol{M}}_{\text{e}, 1}=\left[{\displaystyle \sum\limits_{i=1}^{{N}_{x}}{F}_{\text{e}, i}{y}_{i}}, -{\displaystyle \sum\limits_{i=1}^{{N}_{x}}{F}_{\text{e}, i}{x}_{i}}, 0\right] , (5)

    式中,Nx图 1(b)中矩形条带划分数,xiyi为矩形条带i中心在局部坐标系oxyz中的坐标。该计算方法的本质是将倾斜的黏土片等效为台阶形排列的一系列与中平面平行的条带,条带划分越多越准确。因两个正对平行的条带间双电层斥力远大于两交错平行条带间斥力[10, 15],故式(4),(5)中求和只需进行Nx次,而非Nx2次。

    本文发展的黏土片间相互作用计算方法适用于任意凸多边形。为便于后续算例使用,定义如图 3所示的不同交错区域面积的3个黏土片。不失一般性,黏土片A、B、C的平面形状分别为长方形、七边形、五边形,平面边长0.5 μm到0.6 μm,厚度均为0.06 μm(高岭石矿物构成的黏土片的典型尺寸范围)。oxy平面位于下部黏土片厚度方向的中分面(见后文各图中的黏土片排列)。后文不加说明时,各算例计算参数如下:温度293K,表面电势50 mV,离子浓度0.0015 mol/L,相对介电常数80.1,离子带电荷数1.0。

    图  3  算例中黏土片平面几何形状
    Figure  3.  Planar geometries of clay platelets in analyses

    图 4给出了不同条带划分数Nx下,黏土片B与C之间的双电层排斥作用(用|Fe, z|、|Me, y|表征)随黏土片间夹角的变化,夹角变化过程中保持两黏土片间最小距离dn恒为1 nm。图 4表明,双电层排斥作用随黏土片间夹角增大而快速减小,当夹角超过7°后,

    图  4  不同条带划分数下双电层排斥作用随片间夹角的变化
    Figure  4.  Variation of double-layer repulsive interaction with angle between two clay platelets calculated with different slicing numbers

    双电层排斥作用低于两黏土片平行时排斥作用的1‰,可以忽略,与Jaradat等[15]对矩形片的分析结果定性一致。当Nx达到8后,双电层排斥作用不再随条带划分数量增加而发生明显变化。

    图 5给出了黏土片A与B在不同黏土表面电势下,双电层斥力随黏土片间夹角和距离的变化,计算中取Nx=8。图 5(a)表明,在各表面电势下,当黏土片间夹角超过10°后,双电层斥力低于两黏土片平行时斥力的1%,可以忽略。图 5(b)表明,双电层斥力随距离增大快速减小,当间距超过10 nm后,双电层斥力低于间距1 nm时斥力的1%,可以忽略。

    图  5  不同黏土片表面电势下双电层斥力随黏土片间夹角和距离的变化
    Figure  5.  Variation of double-layer repulsive interaction with angle and distance between two clay platelets with different surface potentials

    黏土片间范德华作用是黏土片1中分子与黏土片2中分子间相互作用的总和。两分子间范德华吸引能表达式最初由London[18]给出

    u(\vec r) = - \frac{B}{{{{\left| {\vec r} \right|}^6}}}, (6)

    式中,B是伦敦常数,\vec r是两分子中心连线矢量,\left| {\vec r} \right|为两分子中心距离。

    实际上,完整的范德华作用包括吸引和排斥两部分,分别与距离的6次方和12次方成反比,都是分子内部电荷分布不均匀在不同尺度的表现。对于离散元模拟的黏土,仅式(6)给出的吸引作用具有实际意义,后文将进一步说明。两分子间范德华作用可通过能量对距离求导得到

    \vec f(\vec r) = - \frac{{\partial u(\vec r)}}{{\partial \vec r}}。 (7)

    图 6(a)中,将两黏土片与中平面投影重叠区域对应的部分分别称为黏土片1'、黏土片2'(分别是黏土片1、黏土片2的一部分,图中用实体示出)。因范德华作用随距离衰减非常快,问题可简化为求解图 6(a)所示的黏土片1'、黏土片2'之间的范德华吸引作用。借助图 1(b)中划分的条带,将黏土片1'划分为若干条块,分别求得每个条块与黏土片2'之间的范德华作用再求和即可。以黏土片1'上的条块i为例,其与黏土片2'的作用可近似为与黏土片2'向各方向无限延伸形成的无穷大黏土片间的范德华作用(无限延伸部分对范德华作用贡献很小,但却极大地方便了解析求解)。对式(7)积分可得一个点与厚度为t的无穷大薄片间范德华力为

    {F_{\text{v}}}(h) = \frac{{{{\rm{ \mathsf{ π} }}}{\rho _1}B}}{2}\left[ {\frac{1}{{h_{}^4}} - \frac{1}{{(h + t)_{}^4}}} \right], (8)
    图  6  黏土片间范德华吸引作用计算方法示意图
    Figure  6.  Schematic diagram of calculation method for Van der Waals attractive interaction between two clay platelets

    式中,h为点与薄片的最小距离,{\rho _1}为薄片中的分子密度。

    条块i与无穷大黏土片间范德华力可通过在条块i上对式(8)进行体积积分求得,一般情形解为

    \begin{array}{l} {F_{{\rm{v}}, i}} = \frac{{A{l_i}\cos 2\theta }}{{24 {\rm{ \mathsf{ π}}} }}\left[ {\frac{1}{{d_1^2}} - \frac{1}{{d_2^2}} - \frac{1}{{d_3^2}} + \frac{1}{{d_4^2}} - } \right.\\ \qquad\quad \left. {\frac{1}{{{{({d_1} + t)}^2}}} + \frac{1}{{{{({d_2} + t)}^2}}} + \frac{1}{{{{({d_3} + t)}^2}}} - \frac{1}{{{{({d_4} + t)}^2}}}} \right], \end{array} (9)

    式中,A={{{\rm{ \mathsf{ π} }}}^2}B{\rho _1}{\rho _2}为哈梅克常数,θ为两黏土片间夹角,d1-d4为条块ioxz平面投影的4个角点到黏土片2'上表面的距离,t为黏土片2'的厚度,详见图 6(b)

    由范德华作用的可叠加性,黏土片1受到的范德华力和相应的力矩为

    {F}_{\text{v}, 1}={\displaystyle \sum\limits_{i=1}^{{N}_{x}}{F}_{\text{v}, i}}\left(\text{sin}\left(\frac{\theta }{2}\right), 0, \text{ cos}\left(\frac{\theta }{2}\right)\right) , (10)
    \begin{array}{l} {M_{{\rm{v}}, 1}} = \left( {{\rm{cos}}\left( {\frac{\theta }{2}} \right) \displaystyle \sum\limits_{i = 1}^{{N_x}} {{F_{{\rm{v}}, i}}} {y_i}} \right., \\ \qquad\quad\;\; - {\rm{cos}}\left( {\frac{\theta }{2}} \right) \displaystyle \sum\limits_{i = 1}^{{N_x}} {{F_{{\rm{v}}, i}}} {x_i} + {\rm{sin}}\left( {\frac{\theta }{2}} \right) \displaystyle \sum\limits_{i = 1}^{{N_x}} {{F_{{\rm{v}}, i}}} {z_i}, \\ \qquad\quad\;\; \left. { - {\rm{sin}}\left( {\frac{\theta }{2}} \right) \displaystyle \sum\limits_{i = 1}^{{N_x}} {{F_{{\rm{v}}, i}}{y_i}} } \right) \end{array} 。 (11)

    图 7给出了黏土片B与C之间在不同条带划分数Nx下,范德华吸引作用随黏土片间夹角的变化,哈梅克常数A取2×10-20J(黏土典型值[10])。当Nx达到8后,范德华吸引作用不再随条带划分数量增大而发生明显变化。图 7还表明,范德华吸引作用随黏土片间夹角增大而快速减小;当黏土片间夹角超过4°后,范德华作用低于黏土片平行时的1‰,可以忽略。

    图  7  不同条带划分数下范德华吸引作用随片间夹角的变化
    Figure  7.  Variation of Van der Waals attractive interaction with angle between two clay platelets calculated with different slicing numbers

    当两黏土片间夹角接近垂直时,上述基于图 6的计算方法可用以下更简便的算法替代。如图 8所示情形,定义以黏土片2所在平面为局部坐标系的oxy平面,z正向从黏土片2指向黏土片1。设黏土片1上距黏土片2最近的点为E;以EA为边的两个侧面EABFEADH中,设面EABFoxy平面夹角更小。两黏土片间的范德华引力可近似为黏土片1'(EABF-IJKL)与无穷大黏土片(近似代替黏土片2)之间的范德华力;此处,黏土片1'是以面EABF为基础沿该面法向远离黏土片2方向扩展而形成的长方体(图 8中以实体示出),扩展距离可取为黏土片1平面尺寸最大值的一半。实际上,由于黏土片1'与黏土片2接近垂直,而范德华作用随距离增加而急剧减小,试算表明扩展距离再大对范德华力影响极小。此外,试算表明,该简便算法与基于图 6的计算方法结果相差极小,但效率更高。

    图  8  两黏土片接近垂直时范德华吸引作用计算方法示意图
    Figure  8.  Schematics diagram of method for Van der Waals attractive interaction between nearly perpendicular clay platelets

    在黏土片1'上对式(8)进行体积积分可得范德华力的一般情形解为

    \left. \begin{array}{l}{F}_{\text{v}}=\frac{A({F}_{\text{a}}-{F}_{\text{b}})}{12{\rm{ \mathsf{ π} }}{n}_{EA, 3}{n}_{EF, 3}{n}_{EH, 3}}, \\ {F}_{\text{a}}=\left(\frac{1}{{h}_{C}}-\frac{1}{{h}_{B}}-\frac{1}{{h}_{D}}+\frac{1}{{h}_{A}}-\frac{1}{{h}_{G}}+\frac{1}{{h}_{F}}+\frac{1}{{h}_{H}}-\frac{1}{{h}_{E}}\right), \\ {F}_{\text{b}}=\left(\frac{1}{{h}_{C}+t}-\frac{1}{{h}_{B}+t}-\frac{1}{{h}_{D}+t}+\frac{1}{{h}_{A}+t}- \right.\\ \qquad\quad \left.\frac{1}{{h}_{G}+t}+\frac{1}{{h}_{F}+t}+\frac{1}{{h}_{H}+t}-\frac{1}{{h}_{E}+t} \right), \end{array} \right\} (12)

    式中,hA~hH为黏土片1'的点A~G到黏土片2上表面的距离,t为黏土片2的厚度, {n}_{EA, 3} {n}_{EF, 3} {n}_{EH, 3} 为局部坐标系中 \overrightarrow {EA} \overrightarrow {EF} \overrightarrow {EH} 方向单位矢量的z方向分量。同理可积分求得范德华力矩。

    图 9给出了不同哈梅克常数下、黏土片A与B之间的范德华引力随黏土片间夹角和距离的变化。当两黏土片间夹角小于73°时用式(10)计算范德华引力,而大于73°时用式(12)计算(此处仅以73°为例,如采用该值附近的其他分界角度,下列结论不变)。由图 9(a)可知,范德华引力随两黏土片夹角增大先快速降低,而后变化较小,随后在夹角较大时再次增大,与Jaradat等[15]对矩形片的分析结果定性一致;式(10),(12)计算的结果过渡平稳,说明两种方法在各自范围内合理适用。图 9(a)还表明,当夹角θ在10°~80°之间时,范德华作用低于黏土片平行时的1%,可以忽略。图 9(b)表明,范德华引力随距离增大也快速减小,当距离超过10 nm后,范德华引力小于间距1 nm时的1‰,可以忽略。

    图  9  不同哈梅克常数下范德华引力随黏土片间夹角和距离的变化
    Figure  9.  Variation of Van der Waals attractive interaction with angle and distance between two clay platelets with different Hamaker constants

    在离散元模拟中,接触模型用于描述颗粒之间的相互作用。黏土片间相互作用包括双电层排斥作用、范德华吸引作用和机械接触作用,三者叠加即可得到黏土片间接触模型。前两者可由本文第2节方法计算。机械接触作用可沿用经典离散元中的抗转动接触模型[19],其法向、切向、转动向作用分别表达为

    {\boldsymbol{F}_{\rm{n}}} = \left\{ {\begin{array}{*{20}{l}} {{K_{\rm{n}}}{u_{\rm{n}}}{\boldsymbol{n}_{\rm{c}}} \quad \left( {{u_{\rm{n}}}{\rm{ > 0}}} \right)}\\ {0 \qquad\qquad \left( {{u_{\rm{n}}} \le {\rm{0}}} \right)} \end{array}} \right. (13)
    {\boldsymbol{F}}_{\text{s}}={\boldsymbol{F}}_{\text{s}}+{K}_{\text{s}}\Delta{\boldsymbol{u}}_{\text{s}}, \left|{\boldsymbol{F}}_{\text{s}}\right|\le \mu \left|{\boldsymbol{F}}_{\text{n}}\right| , (14)
    {\boldsymbol{M}}_{\text{r}}={\boldsymbol{M}}_{\text{r}}+{K}_{\text{r}}\Delta{\boldsymbol{\theta} }_{\text{r}}, \left|{\boldsymbol{M}}_{\text{r}}\right|\le {\mu }_{\text{r}}R\left|{\boldsymbol{F}}_{\text{n}}\right| , (15)

    式中, {\boldsymbol{F}_{\text{n}}} {\boldsymbol{F}_{\text{s}}} {\boldsymbol{M}_{\text{r}}} 为法向、切向、转动向接触力(矩),KnKsKr为相应接触刚度,un为接触重叠量,nc为垂直接触面的单位法向量, {{\Delta }}{\boldsymbol{u}_{\text{s}}} {{\Delta }}{\boldsymbol{\theta} _{\text{r}}} 为相对切向位移和转动角增量,μμr为摩擦系数和抗转动系数,R为离散元接触算法算得的接触半径。

    式(13)~(15)中的接触刚度可表示为

    {K_{\text{n}}}{\text{ = }}2ER , {K_{\text{s}}}{\text{ = }}{K_{\text{n}}}/{\xi _{\text{k}}} , {K_{\text{r}}}{\text{ = }}{K_{\text{s}}}{R^2} , (16)

    式中,E为接触模量, {\xi _{\text{k}}} 为法、切向刚度比。

    值得注意的是,根据式(9),(12),当黏土片表面非常接近(如净距dn小于1 nm)时,范德华引力将非常大,此时范德华斥力(此处未考虑)也将开始发挥作用。此外,黏土片表面存在一层强结合水膜,当两黏土片表面强结合水膜接触后,其相互作用也尚无成熟理论描述。综合上述两个因素,当前常见的做法是设置力学截断距离dm,当黏土片间净距dn小于dm后范德华作用和双电层作用不再随距离而变化,而黏土片固相部分(包括强结合水)开始发生机械接触[10]

    图 10给出了平行黏土片A与C之间的相互作用随净距的变化,图中取力学截断距离dm=1 nm,机械接触法向刚度Kn=60 N/m。图 10中的合力由双电层斥力、范德华引力和机械接触斥力合成,正值表示排斥。两黏土片间相互作用合力在距离较大时表现为斥力,随距离靠近,斥力逐渐增大并在dn=2.0 nm达到局部峰值,黏土片间相互作用由双电层斥力主导;当净距进一步减小,合斥力逐渐降低并在dn=1.25 nm处转变为引力,在dn=dm时引力最大,黏土片间相互作用由范德华引力主导;当净距进一步减小时,合力再次表现为斥力,黏土片间相互作用由机械接触斥力主导。

    图  10  黏土片间法向接触响应曲线
    Figure  10.  Normal contact response curves between clay platelets

    图 11进一步给出了黏土片A与C之间在平行和和成7.5°夹角时,哈梅克常数(控制范德华引力大小)和表面电势(控制双电层斥力大小)对黏土片间相互作用合力的影响。其曲线形态与图 10中的合力曲线一致,但斥力峰值、引力峰值及峰值发生的位置,引力、斥力转变点的位置随具体参数取值变化而有所区别。

    图  11  哈梅克常数和表面电势对法向接触响应的影响
    Figure  11.  Influences of Hamaker constant and surface electric potential on normal contact responses

    将上述接触模型植入PFC3D自定义接触模型中以开展离散元模拟[19]。本节离散元模拟中,为比较黏土片形状对黏土数值试样一维压缩特性的影响,构建了3个黏土数值试样。在试样1中,黏土片平面形状为正方形,边长为0.5 μm、厚度为0.025 μm。在试样2中,参考图 12(a)中黏土电镜扫描图像[4],不妨设定黏土片形状为凸八边形,共有9种不同边长组合,如图 12(b)所示,构造方法如下:①在圆形轮廓上随机选取8个点,依次连接构成长短轴比近似为1.0的初始凸八边形,其面积与试样1中边长0.5 μm的正方形面积相等,两者厚度也相等;②以此为模板,分别沿两个垂直方向拉伸形成长短轴比分别为1.1,1.2,1.3和1.4的凸八边形,对这8个凸八边形进行缩放使其面积与初始凸八边形面积相等。在试样3中,黏土片平面形状与试样2中相同,但厚度为0.05 μm。上述尺寸范围在高岭石矿物构成的黏土片的典型尺寸范围内。

    图  12  黏土片形态
    Figure  12.  Shapes of clay platelets

    为构建黏土的离散元数值试样,首先在由刚性墙边界围合的模型空间中随机生成500个黏土片构成初始试样,试样2和3中的9种黏土片所占数量比例相同。数值模型中的黏土片逐一生成,位置和方向随机,初始试样孔隙比为21(黏土片间无重叠),相当于模拟静水环境中的极稀泥浆。随后,黏土片在自重作用下达到平衡状态,孔隙比降低至某一稳定值(为提高计算效率,此过程中不引入范德华作用和双电层作用)。随后引入范德华作用和双电层作用,再次平衡至稳定。以此为初始状态进行一维压缩模拟。参考Yao等[10]的工作,接触模型参数如表 1所示。

    表  1  离散元接触模型参数
    Table  1.  Contact model parameters in DEM simulation
    参数 取值 参数 取值
    哈梅克常数A 10.2×10-20J 离子浓度n 0.015 mol/L
    表面电势ψ0 193.5 mV 开氏温度T 293 K
    离子带电荷数v 1.0 接触模量E 10 MPa
    溶液相对介电常数ε 80.1 摩擦系数μ 0.5
    法、切向刚度比ξk 2.0 抗转动系数μr 0.15
    力学截断距离dm 0.9 nm
    下载: 导出CSV 
    | 显示表格

    图 13给出了试样1在不同条带划分数Nx下的一维压缩e–logp曲线,其中e为孔隙比、p为竖向有效应力。可以看到,当Nx超过6后不再影响模拟结果,结合图 47结论,综合考虑计算效率和精度的平衡,本文后续离散元模拟中取Nx=8。结合图 5(a)9(a),通过离散元模拟验证表明,若仅在黏土片间夹角θ < 10°或θ > 80°时计算范德华引力、仅在夹角θ < 10°时计算双电层斥力,与所有角度下都计算范德华引力和双电层斥力相比,模拟结果几乎不受影响,但综合节省计算时间约67%。结合图 5(b)9(b)1011,通过离散元模拟验证表明,相互作用计算的起始距离取为10 nm对计算结果几乎无影响。通过上述方法可极大提高离散元模拟黏土力学特性的计算效率,但几乎不影响精度。

    图  13  条带划分数对模拟的一维压缩曲线的影响
    Figure  13.  Influences of slicing number on simulated oedometer test curves

    图 14为试样2在不同孔隙比下的黏土片排列情况。图 15(a)给出了3个黏土数值试样的一维压缩响应。试样1、2和3的压缩指数分别为1.27,1.04和0.92,在黏土压缩指数的常见范围。比较试样1和试样2可知,由于试样1中黏土片平面形状为正方形,其棱角相比于试样2中的凸八边形更加分明,使得试样1中黏土片更倾向于在空间中搭接交错,故而试样1的初始孔隙比更大,压缩性更高。比较试样2和试样3可知,由于试样2中黏土片更薄,比表面积更大,试样中的双电层斥力作用更加明显,故试样2的初始孔隙比更大,压缩性更高。

    图  14  不同孔隙比下的黏土离散元数值试样(试样2)
    Figure  14.  Numerical samples of clay in DEM simulation with different void ratios (sample 2)
    图  15  离散元模拟的黏土一维压缩响应
    Figure  15.  Responses of clay in oedometer tests simulated by DEM

    为表征黏土试样的微观结构,定义试样中黏土片排列的二阶组构张量Fij

    {F_{ij}} = \frac{1}{N}\sum\limits_{k = 1}^N {{n_{k, i}}{n_{k, j}}} , (17)

    式中,N为黏土片的数量,nk, ink, j为黏土片k的单位法向量分量,ij的取值范围为1,2,3。定义组构各向异性系数{I_{\text{F}}} = \sqrt {3{J_{2, F}}} {J_{2, F}}为组构张量Fij的第二不变量)。图 15(b)表明,IF随着孔隙比e减小而增大,这是由于随着一维压缩的进行,黏土片倾向于水平排列,各向异性程度增大。3个数值试样的eIF关系相差较大。

    值得注意的是,此处算例仅为定性反映黏土片形态的影响规律,不针对特定黏土进行定量分析。可以看到,离散元试样中的黏土片形状、厚度对一维压缩模拟结果有明显影响。随着离散元模拟在黏土力学特性研究中的深入应用,有必要考虑黏土片形状这一细致且重要的影响因素。

    本文围绕不规则形状的黏土片间相互作用计算模型和黏土离散元数值模拟,开展了理论分析与模拟工作,得到以下结论:

    (1)提出了不规则凸多边形黏土片间双电层排斥作用和范德华吸引作用的计算方法,其本质是将黏土片等效为台阶形排列的一系列与中平面平行的条带,将条带间相互作用解析解叠加而成。双电层排斥作用和范德华吸引作用都与黏土片的相对排列(以夹角、净距、交错区域面积等表征)密切相关。

    (2)将双电层排斥作用、范德华吸引作用和经典离散元中的抗转动接触模型三者叠加即得黏土片间接触模型。离散元模拟表明,采用以下方法可极大提高模拟效率而几乎不影响准确性:①划分条带数Nx取8;②相互作用计算的起始距离取10 nm;③仅在黏土片间夹角θ < 10°或θ > 80°时计算范德华引力、仅在夹角θ < 10°时计算双电层斥力。

    (3)黏土片形状、厚度对黏土一维压缩模拟的e–lgp曲线、试样微观组构演化都有明显影响,后续研究中有必要考虑黏土片形状这一重要的影响因素。

    本文对黏土片的平片状假设使得研究结果只适用于高岭石等矿物构成的空间翘曲不明显的黏土片。但通过将平片进行空间组合可构建蒙脱石等矿物构成的、具有空间翘曲特征的黏土片,则本文模型可方便地扩展并适用于不同类型黏土的离散元模拟。

    本文及已有相关研究都未考虑黏土片边缘可能带正电荷的影响,笔者将在后续研究中对此深入探讨。后续研究还需考虑黏土片尺寸级配以定量模拟特定黏土。在此基础上,将宏微观土力学的研究思想[20]应用于黏土力学特性的研究将更加便利。

  • 图  1   黏土片间双电层排斥作用计算方法示意图

    Figure  1.   Schematics of calculation method for double-layer repulsive interaction between two clay platelets

    图  2   归一化中平面电势与归一化距离的关系

    Figure  2.   Relationship between normalized mid-plane electric potential and normalized distance

    图  3   算例中黏土片平面几何形状

    Figure  3.   Planar geometries of clay platelets in analyses

    图  4   不同条带划分数下双电层排斥作用随片间夹角的变化

    Figure  4.   Variation of double-layer repulsive interaction with angle between two clay platelets calculated with different slicing numbers

    图  5   不同黏土片表面电势下双电层斥力随黏土片间夹角和距离的变化

    Figure  5.   Variation of double-layer repulsive interaction with angle and distance between two clay platelets with different surface potentials

    图  6   黏土片间范德华吸引作用计算方法示意图

    Figure  6.   Schematic diagram of calculation method for Van der Waals attractive interaction between two clay platelets

    图  7   不同条带划分数下范德华吸引作用随片间夹角的变化

    Figure  7.   Variation of Van der Waals attractive interaction with angle between two clay platelets calculated with different slicing numbers

    图  8   两黏土片接近垂直时范德华吸引作用计算方法示意图

    Figure  8.   Schematics diagram of method for Van der Waals attractive interaction between nearly perpendicular clay platelets

    图  9   不同哈梅克常数下范德华引力随黏土片间夹角和距离的变化

    Figure  9.   Variation of Van der Waals attractive interaction with angle and distance between two clay platelets with different Hamaker constants

    图  10   黏土片间法向接触响应曲线

    Figure  10.   Normal contact response curves between clay platelets

    图  11   哈梅克常数和表面电势对法向接触响应的影响

    Figure  11.   Influences of Hamaker constant and surface electric potential on normal contact responses

    图  12   黏土片形态

    Figure  12.   Shapes of clay platelets

    图  13   条带划分数对模拟的一维压缩曲线的影响

    Figure  13.   Influences of slicing number on simulated oedometer test curves

    图  14   不同孔隙比下的黏土离散元数值试样(试样2)

    Figure  14.   Numerical samples of clay in DEM simulation with different void ratios (sample 2)

    图  15   离散元模拟的黏土一维压缩响应

    Figure  15.   Responses of clay in oedometer tests simulated by DEM

    表  1   离散元接触模型参数

    Table  1   Contact model parameters in DEM simulation

    参数 取值 参数 取值
    哈梅克常数A 10.2×10-20J 离子浓度n 0.015 mol/L
    表面电势ψ0 193.5 mV 开氏温度T 293 K
    离子带电荷数v 1.0 接触模量E 10 MPa
    溶液相对介电常数ε 80.1 摩擦系数μ 0.5
    法、切向刚度比ξk 2.0 抗转动系数μr 0.15
    力学截断距离dm 0.9 nm
    下载: 导出CSV
  • [1]

    YONG R N. Overview of modeling of clay microstructure and interactions for prediction of waste isolation barrier performance[J]. Engineering Geology, 1999, 54(1/2): 83–91.

    [2] 唐朝生, 施斌, 王宝军. 基于SEM土体微观结构研究中的影响因素分析[J]. 岩土工程学报, 2008, 30(4): 560–565. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC200804018.htm

    TANG Chao-sheng, SHI Bin, WANG Bao-jun. Factors affecting analysis of soil microstructure using SEM[J]. Chinese Journal of Geotechnical Engineering, 2008, 30(4): 560–565. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC200804018.htm

    [3] 尹小涛. 基于微结构量化分析的软土各向异性特征研究[J]. 地下空间与工程学报, 2015, 11(增刊2): 486–490. https://www.cnki.com.cn/Article/CJFDTOTAL-BASE2015S2021.htm

    YIN Xiao-tao. Study on anisotropy of soft soil based on quantitative analysis of microstructure[J]. Chinese Journal of Underground Space and Engineering, 2015, 11(S2): 486–490. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-BASE2015S2021.htm

    [4]

    HICHER P Y, WAHYUDI H, TESSIER D. Microstructural analysis of inherent and induced anisotropy in clay[J]. Mechanics of Cohesive-Frictional Materials, 2000, 5(5): 341–371. doi: 10.1002/1099-1484(200007)5:5<341::AID-CFM99>3.0.CO;2-C

    [5]

    WANG Y H, SIU W K. Structure characteristics and mechanical properties of kaolinite soils: Ⅱ effects of structure on mechanical properties[J]. Canadian Geotechnical Journal, 2006, 43(6): 601–617. doi: 10.1139/t06-027

    [6] 刘治清, 宋晶, 杨玉双, 等. 饱和细粒土固结过程的三维孔隙演化特征[J]. 工程地质学报, 2016, 24(5): 931–940. https://www.cnki.com.cn/Article/CJFDTOTAL-GCDZ201605024.htm

    LIU Zhi-qing, SONG Jing, YANG Yu-shuang, et al. Three-dimensional pores evolution characteristics during consolidation process of saturated fine-grained soil[J]. Journal of Engineering Geology, 2016, 24(5): 931–940. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-GCDZ201605024.htm

    [7]

    ANANDARAJAH A. Discrete element modeling of leaching-induced apparent overconsolidation in kaolinite[J]. Soils and Foundations, 2003, 43(6): 1–12. doi: 10.3208/sandf.43.6_1

    [8]

    BAYESTEH H, MIRGHASEMI A A. Numerical simulation of pore fluid characteristic effect on the volume change behavior of montmorillonite clays[J]. Computers and Geotechnics, 2013, 48: 146–155. doi: 10.1016/j.compgeo.2012.10.007

    [9] 商翔宇, 鲁巨明, 杨晨, 等. 考虑黏土特性的离散元程序开发[J]. 防灾减灾工程学报, 2016, 36(4): 657–663. https://www.cnki.com.cn/Article/CJFDTOTAL-DZXK201604023.htm

    SHANG Xiang-yu, LU Ju-ming, YANG Chen, et al. Development of discrete element code considering the characteristics of clay[J]. Journal of Disaster Prevention and Mitigation Engineering, 2016, 36(4): 657–663. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-DZXK201604023.htm

    [10]

    YAO M, ANANDARAJAH A. Three-dimensional discrete element method of analysis of clays[J]. Journal of Engineering Mechanics, 2003, 129(6): 585–596. doi: 10.1061/(ASCE)0733-9399(2003)129:6(585)

    [11]

    KATTI D R, MATAR M I, KATTI K S, et al. Multiscale modeling of swelling clays: a computational and experimental approach[J]. KSCE Journal of Civil Engineering, 2009, 13(4): 243–255. doi: 10.1007/s12205-009-0243-0

    [12]

    EBRAHIMI D, WHITTLE A J, PELLENQ R J M. Mesoscale properties of clay aggregates from potential of mean force representation of interactions between nanoplatelets[J]. The Journal of Chemical Physics, 2014, 140(15): 154309. doi: 10.1063/1.4870932

    [13]

    SJOBLOM K J. Coarse-grained molecular dynamics approach to simulating clay behavior[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2016, 142(2): 06015013. doi: 10.1061/(ASCE)GT.1943-5606.0001394

    [14]

    GUO Y, YU X. A holistic computational model for prediction of clay suspension structure[J]. International Journal of Sediment Research, 2019, 34(4): 345–354. doi: 10.1016/j.ijsrc.2018.12.002

    [15]

    JARADAT K A, ABDELAZIZ S L. On the use of discrete element method for multi-scale assessment of clay behavior[J]. Computers and Geotechnics, 2019, 112: 329–341. doi: 10.1016/j.compgeo.2019.05.001

    [16]

    LU N, ANDERSON M T, LIKOS W J, et al. A discrete element model for kaolinite aggregate formation during sedimentation[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2008, 32(8): 965–980. doi: 10.1002/nag.656

    [17]

    JUN S M, JIANG M J, ZHOU C B. Numerical simulation of desiccation cracking in a thin clay layer using 3D discrete element modeling[J]. Computers and Geotechnics, 2014, 56: 168–180. doi: 10.1016/j.compgeo.2013.12.003

    [18]

    LONDON F. Zur theorie und systematik der molekularkräfte[J]. Zeitschrift Für Physik, 1930, 63(3/4): 245–279. http://www.mendeley.com/research/zur-theorie-und-systematik-der-molekularkrfte/

    [19]

    LONDON F. On the theory and systematics of the molecular forces[J]. Journal of Physics, 1930, 63(3/4): 245–279. (in German) doi: 10.1142/9789812795762_0023

    [20]

    ITASCA CONSULTING GROUP, INC. Documentation of Particle Flow Code 3D V6.0[M]. Minneapolis, 2019.

    [21] 蒋明镜. 现代土力学研究的新视野: 宏微观土力学[J]. 岩土工程学报, 2019, 41(2): 195–254. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201902002.htm

    JIANG Ming-jing. New paradigm for modern soil mechanics: geomechanics from micro to macro[J]. Chinese Journal of Geotechnical Engineering, 2019, 41(2): 195–254. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201902002.htm

图(15)  /  表(1)
计量
  • 文章访问数:  258
  • HTML全文浏览量:  66
  • PDF下载量:  62
  • 被引次数: 0
出版历程
  • 收稿日期:  2021-07-21
  • 网络出版日期:  2022-09-22
  • 刊出日期:  2022-08-31

目录

/

返回文章
返回