工程地质学报  2017, Vol. 25 Issue (6): 1430-1437   (3437 KB)    
Article Options
  • PDF (3437 KB)
  • Full Text HTML
  • Abstract
  • Figures
  • References
  • History
  • 收稿日期:2016-09-14
  • 收到修改稿日期:2017-06-04
  • 扩展功能
    把本文推荐给朋友
    加入引用管理器
    Email Alert
    文章反馈
    浏览反馈信息
    本文作者相关文章
    张晓宇
    许强
    刘春
    施斌
    黏性土失水开裂多场耦合离散元数值模拟
    张晓宇, 许强, 刘春①②③, 施斌    
    ① 南京大学地球科学与工程学院 南京 210023;
    ② 地质灾害防治与地质环境保护国家重点实验室(成都理工大学) 成都 610059;
    ③ 南京大学(苏州)高新技术研究院 苏州 215123
    摘要:黏性土在失水过程中逐渐变形开裂,裂隙相互交错形成网络,这一过程涉及水热力等多场耦合的作用。本文基于南京大学自主研发的三维离散元软件MatDEM,采用紧密堆积离散元模型对土体开裂进行模拟。在此模型中,每个离散元单元代表一定体积的土颗粒、孔隙和孔隙水的集合体。单元具有含水量属性,并采用有限差分思想计算水分运移量,实现了水分场模拟。同时,考虑水分场对土体抗拉强度等力学性质的影响,建立水分场和应力场的耦合。在数值模拟中,假定土体表面水分以一定速率蒸发,根据试验数据由含水量计算单元直径和力学参数,从而实现黏土蒸发失水、收缩和开裂变形过程模拟。数值模拟与前人室内试验结果基本一致,能较好地再现开裂过程中的各个阶段。本文为研究多场作用下土体变形破坏模拟提供了一个新的思路。
    关键词黏土    土体开裂    离散元    数值模拟    
    NUMERICAL SIMULATION OF DRYING CRACKING USING MULTI-FIELD COUPLING DISCRETE ELEMENT METHOD
    ZHANG Xiaoyu, XU Qiang, LIU Chun①②③, SHI Bin    
    ① School of Earth Sciences and Engineering, Nanjing University, Nanjing 210023;
    ② State Key Laboratory of Geo-hazard Prevention and Geo-environment Protection, Chengdu University of Technology, Chengdu 610059;
    ③ Nanjing University(Suzhou) High-Tech Institute, Suzhou 215123
    Abstract: Desiccation cracks can be formed in clay soil during drying, which connects with each other to generate a network. This process involves multi-field coupling effect, including water, thermal and stress. A close-packed model is used to simulate the formation of desiccation cracks. The simulation is based on the 3D discrete element modeling software MatDEM that is independently developed by Nanjing University. In the model, a discrete unit represents an assembly of soil particles, pore and pore water. Elements in the model are assigned moisture properties, where the transport of moisture can be simulated based on the Finite Difference Method. Furthermore, the coupling of moisture field and stress field is built, which considers the influence of moisture on the tensile strength of soil. In numerical simulation, assuming the water of soil surface evaporates at a certain rate, element diameters and mechanical parameters are calculated from moisture according the experimental data, so as to simulate evaporation, shrinkage and cracking processes of clay soil. Numerical simulation results coincide with the laboratory results, and each step of cracking process can be well simulated. This paper provides an alternative way for the numerical simulation of soil deformation and failure under multi-field effect.
    Key words: Clay soil    Desiccation crack    DEM    Numerical simulation    

    0 引言

    黏性土分布极其广泛,常作为建筑物地基或用作堤坝、路堤填土等材料。由于干燥失水等原因,黏性土表面与内部常常会出现收缩裂缝,提高了土体的渗透性降低了整体强度以及稳定性。在边坡工程中,坡面的裂隙加速了土体的风化,加剧了坡体的水土流失,雨水能沿着裂隙渗入坡体内部,诱导滑坡的发生。在环境岩土工程中,黏性土常作为一种缓冲材料,用于垃圾卫生填埋场的衬垫层和核废料地质处置库的工程屏障中,在温度梯度和地下水的相互作用下,土裂隙的产生将极大缩短这一类隔离材料的使用年限,增加污染物泄露的风险(施斌等,2009)。

    土体干缩裂隙的机制理论和数值模拟是研究和预测土体裂隙形成和发展的重要手段。国际上的相关研究较早,主要基于弹性理论和断裂力学理论来研究和解释土体干缩裂隙的一些现象和规律:Morris et al.(1992)指出土体的干缩裂隙受土吸力,土的弹性模量,泊松比,抗拉强度和表面能等的控制。在线弹性理论基础上,Konrad et al.(1997)提出了黏土干缩裂隙间距理想化预测模型。基于断裂力学的方法可以分析裂隙在发展过程中的应力分布,但是这种方法是建立于连续均匀介质的基础上,无法模拟土体裂隙的形成。Chertkov(2012)从理论上研究了不同因素对裂隙网络形成的影响,这些工作对裂隙的模拟和预测有指导意义。在国内,沈珠江等(2004)提出结合非饱和土简化固结理论和有限元法来模拟裂隙形成过程。高志朋(2007)提出了膨胀土路堑边坡开裂的有限元分析方法,获得了膨胀土路基边坡在自然干燥过程中产生收缩开裂及其裂缝发展的规律。这些模型还相对理想化,无法动态地模拟整个裂隙网络的形成过程。

    最近,离散元数值法开始应用于土体干缩裂隙的模拟中(Peron et al., 2009Amarasiri et al., 2011),这种方法将土体离散成一系列具有弹性的颗粒,颗粒间通过可断裂的连接联系并相互作用,这种数值方法可以很方便地模拟土体的动态开裂过程。通过定义颗粒单元不同的力学参数,可模拟土体的不均匀性等各种性质。离散元法是模拟和研究土体裂隙形成的一个可选方法。近几年,具有弹塑性、黏聚力等特性的土体离散元模型有了较快的发展(El Shamy et al., 2012)。在国内,离散元法也在土体变形和破坏机理研究方面得到应用(赵学亮等,2012)。

    可见,离散元法在研究黏性土变形开裂和模拟黏性土力学特性上具有突出的优势。然而,三维离散元多场耦合(如:温度场,水分场,应力-应变场等)理论和应用研究还处于起步阶段,相应的数值模拟工具还有所欠缺。目前,离散元数值模拟研究多基于一些商业软件(如PFC等),很难用于解决多场作用问题。刘春(2012)开发了三维离散元软件MatDEM,对三维离散元多场耦合做了初步的理论研究,并通过实验和模拟研究了黏性土抗剪强度受温度和含水量等因素的影响规律及工程效应(Liu et al., 2013),这些工作为进一步开展多场作用下黏性土变形开裂三维离散元建模和模拟提供了基础。

    本文基于自主研发的三维离散元软件MatDEM,建立紧密堆积离散元土体模型,模拟研究在水分场、应力场和土体耦合的情形下,黏土体的失水收缩开裂过程。随着土体水分的蒸发,根据试验数据建立含水量和单元直径的关系,并由含水量推导计算土体单元力学参数变化,实现黏土蒸发失水、收缩和开裂变形过程模拟。并对模拟过程中的一些结果和现象进行初步探讨,验证离散元法模拟土体多场作用下变形开裂的可行性。

    1 离散元基本模型

    黏土的微观结构非常复杂,颗粒之间还存在着各种物理化学作用,前人研究表明(Tessier,1984; Ferber et al., 2008),微观层面上大量黏土颗粒会聚集在一起构成一个团粒状结构,大量团粒彼此连结又构成黏土体,水和空隙赋存其中。因此,为了简化计算,本文采用的离散元基本单元代表一系列黏土颗粒和水组成的团粒。

    1.1 离散元单元力学性质

    本文采用紧密规则堆积的离散元模型(图 1a)(Liu et al., 2013),模型由一系列服从牛顿的运动方程的弹性单元(即圆球)组成。单元由可断裂的弹簧相互连接,作用力只能出现在相邻单元之间的接触点(图 1b图 1c)。

    图 1 (a) 三维紧密堆积网格、(b)单元通过弹簧力与邻近颗粒连接和(c)单元间的剪切力通过剪切方向的弹簧表示 Fig. 1 (a)The three-dimensional close-packed lattice. (b)Two particles are bonded by breakable elastic spring along the normal direction. (c)Two particles are bonded by breakable elastic spring along the tangential direction to simulate the shear force

    在模型中(图 1a),相同直径的球单元紧密堆积。单元间的法向力用弹簧力Fn来表示,并由以下公式给出:

    $ {F_n} = {K_n}{\cdot}{X_n} $ (1)

    式中,Kn为弹簧法向刚度;Xn为法向相对位移。单元相互连接时,可产生拉力的作用,当Xn超过断裂位移Xb时,连接断开。因此,单元间最大法向力(Fnmax)为:

    $ {F_{n{\rm{max}}}} = {K_n}{\cdot}{X_b}\left({连接完整时} \right) $ (2)

    当连接在力的作用下断开,单元间拉力不再存在。而当两单元回到压缩接触状态时,单元间的压力仍然存在:

    $ {F_n} = {K_n}{\cdot}{X_n}({X_n} < 0)\left({连接断裂时} \right) $ (3)

    同时,模型中还考虑了剪切力(切向)FS的作用。如图 1c所示,当两个单元表面相互接触并发生相对滑动时,产生滑移方向相反滑动摩擦力。两个单元在切线方向也由可断裂的弹簧连接,切向弹簧力(FS)由下式定义:

    $ {F_s} = {K_s}{\cdot}{X_s} $ (4)

    式中,Ks为剪切刚度;Xs为两单元之间的切向相对位移。对于完整的单元连接,其可承受的最大剪切力(Fsmax)由库仑准则确定:

    $ {F_{s{\rm{max}}}} = {F_{{s_0}}} - {\mu _p}{\cdot}{F_n}\left({连接完整时} \right) $ (5)

    式中,Fs0为单元间抗剪力;μp为单元间摩擦系数;Fn是法向力(压力为负)。当单元所受外力超过式(5)确定的最大剪切力时,单元间的完整连接断开,单元间的剪切阻力(Fs0)消失。此时,切向力(Fs)小于或等于断裂连接时的最大剪切力:

    $ {F_{s{\rm{max}}}} = - {\mu _p}{\cdot}{F_n}\left({连接断裂时} \right) $ (6)

    在连接断裂的情况下,当外力超过由方程(6)确定的最大剪切力时,两个单元开始相对滑动;当两单元彼此分开时(Xn>0),单元间法向力和切向力均为零。

    1.2 离散元法数值模拟方法

    在物理世界中,由于岩土体的颗粒并非完全弹性,存在流体或者微裂隙等原因,动能会逐渐转化成热能,颗粒振动很快停止。在数值模拟中,通过引入了黏滞阻力(Fv)减弱来模拟机械能的衰减。黏滞阻力与单元速度成正比,并由下式给定:

    $ {F_v} = - \eta {\cdot}V $ (7)

    式中,η为阻尼系数;V为单元的速度。

    作用于粒子上的合力是所有连接上的力、重力和黏滞力等的总和。颗粒堆积体的动态演化可以用时间步算法来模拟和观察(Cundall,1979)。在一个非常小的时间步内,假设球单元的运动是线性的,可以计算出在这个时间步内球单元的受力,加速度,速度和位移。通过迭代运算实现颗粒的运动变化,并反映模型土体的宏观变化。

    2 土体离散元多场耦合理论

    为有效地模拟土体的失水收缩开裂过程,本模型对离散单元附加含水量属性。计算水分运移时,首先进行蒸发计算,然后计算水分传导。根据含水量的变化,计算土体体积收缩和强度计算由此实现水分场、应力场的耦合,最后根据牛顿运动定律计算颗粒位移。具体描述如下:

    2.1 蒸发过程和水分传导

    土体水分蒸发过程是孔隙中的水分通过蒸发面(土-气界面)逸入大气,从而使土体含水量逐渐降低而变干的过程。在此模型中,每个离散元单元代表一定体积的土颗粒、孔隙和孔隙水的集合体。把水分赋给每个单元,即单元的含水量表示土颗粒及其周围区域的平均含水量。单元具有一定的初始含水率,且表面单元含水率以一定速度减少,以模拟蒸发作用。经过一个时间步,土体表面单元的含水量(ω)降为:

    $ \omega = {\omega _0}{\cdot}\left({1 - \beta } \right) $ (8)

    式中,ω0为上一时间步结束时单元的含水量;β为水分减少率,其大于并接近0。

    单元具有特定的含水量和扩散率,采用有限差分法思想,实现水分由高含水量单元向低含水量单元运移,以此来模拟土中水的运移。假定单位时间内,相邻单元的水分运移量与两者含水量的差值成正比,单位时间dt内,单元1向单元2的水分运移量为:

    $ {\rm{d}}{M_{12}} = {\alpha _\omega }{\cdot}({\omega _1} - {\omega _2}){\cdot}A{\cdot}{\rm{d}}t $ (9)

    其中,αω为单元间水分扩散率,ω为含水量,A为接触面积。通过此方法实现了土中水的传导模拟。

    2.2 单元收缩和强度计算

    随着土中水分的蒸发,土体的含水量由高到低,黏性土颗粒表面的双电层中结合水膜厚度逐渐变薄并且导致吸力的产生。在吸力作用下,土颗粒会发生重排而逐渐靠拢,孔隙不断减小,宏观上表现为土体收缩。由于收缩变形的不均匀性或边界的约束,在土体中会形成张拉应力场。当拉应力的大小超过土体的抗拉强度时,土体破坏,裂隙便会产生。在失水过程中,拉应力和抗拉强度均是不断增大的,含水量和抗拉强度密切相关,通过建立含水量和抗拉强度等力学性质之间的对应关系,在数值模拟中需根据单元的含水量设定其抗拉力、初始抗剪力等力学参数。

    为了实现土体收缩模拟,假定土颗粒半径与含水量成正相关:

    $ R = {R_0}{\cdot}(1 - \gamma {\cdot}\frac{{{\omega _i} - \omega }}{{{\omega _i}}}) $ (10)

    式中,R0为颗粒初始半径;γ为收缩相关的系数;ωi为初始含水率;ω为当前含水率。根据冉龙洲等(2011)的研究,Romainville黏土的单轴抗拉强度Tu与含水率ω的关系为ω与ln(Tu)线性相关:

    $ {T_u} = 3729{e^{-\frac{{ \omega }}{{12}}}} - 25.21 $ (11)

    在模拟中,利用上式含水率和抗拉强度的关系,由单元含水率来计算单元力学性质,实现水分场对土体强度作用的模拟。

    3 土体多场耦合模型建立
    3.1 参数设置及测试

    本研究基于以上原理改进自主研发的三维离散元软件MatDEM(刘春等,2014),实现黏土开裂多场耦合模拟。本次模拟构建了一个16cm×16cm×2.4cm的土体模型(图 2),单元直径为d=4mm,个数为13078个,其中边界的单元数为3067个,为隔水光滑边界。土的性质以Romainville土为参考,取:ρ=1.8g ·cm-3,杨氏模量E=12MPa,泊松比v=0.18,抗拉强度Tu=60kPa,抗压强度Cu=200kPa,内摩擦系数μ=0.7。参数设置完毕后,运用Liu et al.(2015)给出的转换公式,可以计算出模型单元的力学参数。

    图 2 预平衡前后模型垂直方向位移变化 Fig. 2 Variation of vertical displacement of model after the balance

    根据前人研究(Liu et al., 2013),由转换公式得到紧密堆积模型,其整体力学性质通常会比设定值小。因此,需要通过数值模拟测试得到模型实际力学性质。测试时,模型选用如图 1a所示紧密堆积立方体,模型材料性质如上文所述。

    测试中,对模型底部和上部施加垂向压力,模型产生水平和垂向变形量,根据应力和应变得到模型的杨氏模量和泊松比,抗拉强度和抗压强度。模拟中构建了8个不同大小模型,其水平方向数分别为1,2,6,8,12,16,20和24,相应的颗粒数则为14,24,334,717,2467,6052,11884和19350。图 3图 4为不同颗粒数下模型的弹性模量,泊松比,抗拉强度和抗压强度的对应值。从图 3图 4可以看出,弹性模量的值在颗粒较小时会有一定程度的波动,当颗粒数较多、模型较大基本和设定值一致,误差比较小,约为0.3%。泊松比与设定值较好符合,误差在2%左右。对于抗压强度,模型颗粒数较少时,模型运算中的实际强度和设定强度在颗粒数量较多时,由200kPa减小到110kPa,大约折减为原来的55%。而抗拉强度60kPa减小到了53kPa,约为原来的88%(图 4)。而在模拟中所用模型颗粒数为13078,此时土的强度基本稳定下来,所以模拟中实际的土体抗压强度和抗拉强度分别为110kPa和53kPa。

    图 3 弹性模量和泊松比随模型变大时的变化情况 Fig. 3 The change of tested Young's modulus and Poisson's ratio with increasing element number

    图 4 抗拉强度和抗压强度随模型变大时的变化情况 Fig. 4 The change of compressive strength and tensile strength with increasing element number

    3.2 预平衡

    在颗粒堆积好后,需要通过预平衡使单元在重力作用下运动至稳定状态。如图 2所示是模拟中建立的16cm×16cm×2.4cm土体模型。图 2a是未平衡时颗粒位移情况,图 2b是平衡后颗粒的位移情况,由黄色到蓝色向下的位移逐渐增大。可见在重力作用下颗粒明显下沉收缩,中间部分下沉最为剧烈,而靠近边界部分由于和稳定不动的边界拖拽,所以下沉量小于中心部分。多次平衡直至土体的位移量不再发生变化时,说明模型已经稳定下来。

    3.3 数值模拟

    采用MatDEM进行一系列离散元多场耦合数值模拟。在一个时间步内,模型上表面水分以一定速率蒸发(式8);同时,用有限差分思想计算每个颗粒新的含水量(式9),并根据含水率与抗拉强度(式11)得到新的抗拉强度,由转换公式(Liu et al., 2015)计算颗粒的抗拉力,和其他力学参数。最后通过常规的离散元法计算颗粒受力和位移,一直循环重复,直到开裂基本结束,土体保持稳定。

    4 模拟结果和分析
    4.1 裂隙形成和发展分析

    模拟结果初始开裂位置总是发生在边缘处(图 5a),这与前人的试验结果一致(Tang et al., 2008)。图 6是缩小一倍的模型的主裂隙形成前的含水量分布情况和连接关系情况。由于模拟中边界隔水,当表面蒸发失水后,与边界单元相连的土体单元得到的水分补给少,含水率较低,由于收缩量较大,易于开裂,故周边总是先产生裂隙。同样的,由于表面蒸发,表面单元的含水量比底部低,开裂时土体会先从表面开裂,随后裂隙逐渐往下加深,这也符合土体开裂的实际情况。

    图 5 裂隙发育演变过程 Fig. 5 Development and evolution processes of cracks

    图 6 (a) 主裂隙形成前土体含水量;(b)单元间的连接网 Fig. 6 (a)Distribution of soil water content before the formation of cracks; (b)The intact bonds between elements

    土体的裂隙随含水率的演变过程(图 5),与室内试验结果和其他数值模拟结果(唐朝生等,2012; 司马军等,2013Hirobe et al., 2016)一致,数值模拟的裂隙发育过程同样可以分为3个阶段:(1)从土体中随机产生的薄弱点处逐渐产生裂缝,随着时间的发展,将土块分割成几个大块,这些裂缝即为主裂缝(图 5a图 5b)。(2)在初级块区形成的同时,主裂缝上开始出现分叉,这种自主衍生的分支裂缝被称为子裂缝,若干子裂缝将初级块区细分成更多的次级块区,裂缝数量长度等进一步增加(图 5c图 5d)。(3)裂缝网格基本稳定,土样表面不再出现新生裂缝,裂缝逐渐变宽直到最终稳定(图 5e图 5f)。

    4.2 含水率变化分析

    模型表面由于蒸发失水,表面含水率低于底部含水率,四周含水率低于中部含水率。随着土体表面继续蒸发,在边界和土体中强度较低位置的单元附近开裂。裂隙向下发展,并形成新的表面。在模拟中,这种新形成的表面也进行蒸发失水。由于上表面和裂隙侧表面的共同蒸发作用,裂隙处的蒸发速率较大,含水量较低,单元收缩较快,并进一步促进了裂隙向前发展和加宽。图 7a给出了裂隙网络形成时的模型含水量分布图。可以看到,在裂隙附近单元的含水率明显较低,说明水分的运移和含水量的变化对裂隙形成和分布有重要作用。图 7b给出了图 7a 3个位置处的纵剖面(0mm,16mm,48mm),单元含水量由下向上减少。

    图 7 (a) 土体含水量分布;(b)含水量分布切片 Fig. 7 (a)The soil moisture distribution; (b)The slices of soil moisture distribution map

    4.3 模拟中能量变化

    MatDEM软件实现了能量转化和能量守恒计算(刘春等,2014)。通过离散元模拟里的能量监测模块,我们可以得到模拟过程中的各种能量变化。图 8是模拟的16cm×16cm×2.4cm的土体在开裂过程的能量变化图。在0.5时刻前为预平衡阶段,几乎没有能量变化。在0.5~1.5时刻间为土体收缩(但无裂隙产生)阶段,可以看到模型的弹性应变能不断增加(图 8b),而几乎没有断裂热等热量产生。土体裂隙产生后(1.5时刻),弹性势能部分在阻尼过程中转化为热量(图 8a)。同时,土体表面水分蒸发促使单元不断收缩,整个系统能量不断增加。在约2.6时刻处机械能突然减少,热量相应增加,是由于生成了大规模的新的内部裂隙,使部分弹性势能在阻尼过程中转化为了热量。

    图 8 (a) 系统热量曲线;(b)能量转化曲线 Fig. 8 (a)Heat generation curve; (b)Energy conversion curve

    5 结论

    土体开裂受到多场耦合作用,影响因素十分复杂,加上土颗粒的不均匀性和离散性,使得对土体开裂的模拟异常困难。本文提出一种三维离散元多场耦合模型,实现了水分场和应力场的耦合情况下土体失水开裂的模拟。通过模拟可以得出以下几个结论:

    (1) 模型的模拟结果基本符合实验室内土体失水开裂试验的规律。模型在蒸发失水时,首先从边界和缺陷位置处开裂,并能较好地模拟黏土失水开裂的3个阶段。

    (2) 含水量分析表明,裂隙的形成和发展与水分的运移和分布有密切关系。裂隙的产生又加剧了水的分布的不均匀性。

    (3) 模型实现了黏土表面蒸发失水、水分运移、单元收缩和强度变化、模型变形开裂的多场耦合。能够有效地进行土体开裂问题模拟,为岩土体的多场耦合模拟研究提供一个新的思路。

    本文通过实例模拟,主要验证了黏土失水开裂多场耦合模拟方法的有效性。而土体裂隙的形成和发展受不同土体性质,不同的蒸发速度,不同厚度等复杂因素的影响,这些因素的作用机制和规律还有待于进一步的研究。本文在考虑水分场对应力场的影响时,仅探究了水分场影响抗拉强度的影响,还缺乏抗剪强度,抗压强度,内摩擦系数等其他参数与含水量之间的关系,需进一步开展工作。

    参考文献
    Amarasiri A L, Kodikara J K, Susanga C. 2011. Numerical modelling of desiccation cracking[J]. International Journal for Numerical & Analytical Methods in Geomechanics, 35(1): 82~96.
    Chertkov V Y. 2012. An integrated approach to soil structure, shrinkage, and cracking in samples and layers[J]. Geoderma, (s173~174): 258~273.
    Cundall P A, Strack O D L. 1979. A discrete numerical mode for granular assemblies[J]. Géotechnique, 29(1): 47~65. DOI:10.1680/geot.1979.29.1.47
    El Youssoufi M S, Delenne J Y, Radjai F. 2004. Self-stresses and crack formation by particle swelling in cohesive granular media[J]. Physical Review E, 71(5Pt1).
    El Shamy U, Zamani N. 2012. Discrete element method simulations of the seismic response of shallow foundations including soil-foundation-structure interaction[J]. International Journal for Numerical & Analytical Methods in Geomechanics, 36(10): 1303~1329.
    Ferber V, Auriol J C C, Cui Y J J, et al. 2008. Wetting-induced volume changes in compacted silty clays and high-plast[J]. Canadian Geotechnical Journal, 45(45): 252~265.
    Gao Z P. 2007. The research on the cracking property of expansive soils in wet-dry cycling of natural weather and cracking evolution of cutting slope[D]. Changsha:Changsha University of Science & Technology.
    Hirobe S, Oguni K. 2016. Coupling analysis of pattern formation in desiccation cracks[J]. Computer Methods in Applied Mechanics & Engineering, 307: 470~488.
    Konrad J M, Ayad R. 1997. A idealized framework for the analysis of cohesive soils undergoing desiccation[J]. Canadian Geotechnical Journal, 34(4): 477~488. DOI:10.1139/t97-015
    Liu C, Pollard D D, Gu K, et al. 2015. Mechanism of formation of wiggly compaction bands in porous sandstone:2.Numerical simulation using discrete element method[J]. Journal of Geophysical Research Solid Earth, 120(12): 8153~8168. DOI:10.1002/2015JB012374
    Liu C, Pollard D D, Shi B. 2013. Analytical solutions and numerical tests of elastic and failure behaviors of close-packed lattice for brittle rocks and crystals[J]. Journal of Geophysical Research Solid Earth, 118(1): 71~82. DOI:10.1029/2012JB009615
    Liu C, Shi B, Gu K, et al. 2014. Development and application of large-scale discrete element simulation system for rock and soil[J]. Journal of Engineering Geology, 22(S): 551~557.
    Liu C, Shi B, Shao Y, et al. 2013. Experimental and numerical investigation of the effect of the urban heat island on slope stability[J]. Bulletin of Engineering Geology & the Environment, 72(3-4): 303~310.
    Liu C, Wang B J, Shi B, et al. 2008. Analytic method of morphological parameters of cracks for rock and soil based on image processing and recognition[J]. Chinese Journal of Geotechnical Engineering, 30(9): 1383~1388.
    Liu C. 2012. Image analysis and numerical simulation of soil temperature field, strength and deformation in urban heat island environment[D]. Nanjing:Nanjing University.
    Morris P H, Graham J, Williams D J. 1992. Cracking in drying soils[J]. Canadian Geotechnical Journal, 29(2): 263~277. DOI:10.1139/t92-030
    Peron H, Delenne J Y, Laloui L, et al. 2009. Discrete element modelling of drying shrinkage and cracking of soils[J]. Computers and Geotechnics, 36(1-2): 61~69. DOI:10.1016/j.compgeo.2008.04.002
    Ran L Z, Song X D, Tang C S. 2011. Laboratorial investigation on tensile strength of expansive soil during drying[J]. Journal of Engineering Geology, 19(4): 620~625.
    Shen Z J, Deng G. 2004. Numerical simulation of crack evolution in clay during drying and wetting cycle[J]. Rock and Soil Mechanics, 25(S2): 1~6.
    Shi B, Tang C S, Wang B J, et al. 2009. Development and mechanism of desiccation cracking of clayey soil under different Temperatures[J]. Geological Journal of China Universities, 15(2): 192~198.
    Shuichiro Yoshida, Kazuhide Adachi. 2001. Effects of cropping and puddling practices on the cracking patterns in paddy fields[J]. Soil Science & Plant Nutrition, 47(3): 519~532.
    Si M J, Jiang M J, Zhou C B. 2013. Numerical simulation of desiccation cracking of clay soils by DEM[J]. Chinese Journal of Geotechnical Engineering, 35(S2): 286~291.
    Tang C S, Cui Y J, Anh-minh Tang, et al. 2012. Shrinkage and desiccation cracking process of expansive soil and its temperature-dependent behaviour[J]. Chinese Journal of Geotechnical Engineering, 34(12): 2181~2187.
    Tang C S, Shi B, Liu C. 2012. Study on desiccation cracking behaviour of expansive soil[J]. Journal of Engineering Geology, 20(5): 663~673.
    Tang C, Shi B, Liu C, et al. 2008. Influencing factors of geometrical structure of surface shrinkage cracks in clayey soils[J]. Engineering Geology, 101(3-4): 204~217. DOI:10.1016/j.enggeo.2008.05.005
    Tessier D. 1984. Étude expérimentale de l'organisation des matériaux argileux. Hydratation, gonflement et structuration au cours de la dessiccation et de la réhumectation(France)[D]. Paris, France:Université de Paris Ⅶ.
    Zhao X L, He J M, Dong G F, et al. 2012. New developments of microscale study on granular soil using discrete element method[J]. Journal of Engineering Geology, 20(4): 607~613.
    高志朋. 2007. 膨胀土干湿循环开裂特性及路堑边坡开裂演化过程研究[D]. 长沙: 长沙理工大学. http://d.wanfangdata.com.cn/Thesis/Y1188879
    刘春, 施斌, 顾凯, 等. 2014. 岩土体大型三维离散元模拟系统的研发与应用[J]. 工程地质学报, 22(S): 551~557.
    刘春, 王宝军, 施斌, 等. 2008. 基于数字图像识别的岩土体裂隙形态参数分析方法[J]. 岩土工程学报, 30(9): 1383~1388.
    刘春. 2012. 城市热岛环境中土体温度场、强度和变形的数值模拟与图像分析[D]. 南京: 南京大学. http://d.g.wanfangdata.com.cn/Thesis_Y2941922.aspx
    冉龙洲, 宋翔东, 唐朝生. 2011. 干燥过程中膨胀土抗拉强度特性研究[J]. 工程地质学报, 19(4): 620~625.
    沈珠江, 邓刚. 2004. 黏土干湿循环中裂缝演变过程的数值模拟[J]. 岩土力学, 25(S2): 1~6.
    施斌, 唐朝生, 王宝军, 等. 2009. 黏性土在不同温度下龟裂的发展及其机理讨论[J]. 高校地质学报, 15(2): 192~198.
    司马军, 蒋明镜, 周创兵. 2013. 黏性土干缩开裂过程离散元数值模拟[J]. 岩土工程学报, 35(S2): 286~291.
    唐朝生, 崔玉军, Anh-minhTang, 等. 2012. 膨胀土收缩开裂过程及其温度效应[J]. 岩石工程学报, 34(12): 2181~2187.
    唐朝生, 施斌, 刘春. 2012. 膨胀土收缩开裂特性研究[J]. 工程地质学报, 20(5): 663~673.
    赵学亮, 赫建明, 董高峰, 等. 2012. 离散单元法对粒状土的微观特性研究探讨[J]. 工程地质学报, 20(4): 607~613.