土体干缩开裂过程的边界效应试验与离散元模拟

林朱元; 唐朝生; 曾浩; 王怡舒; 程青; 施斌

doi:10.11779/CJGE202002019

土体干缩开裂过程的边界效应试验与离散元模拟 English Version

林朱元^1,,
唐朝生^1, ,,
曾浩¹,
王怡舒²,
程青¹,
施斌¹

1.
南京大学地球科学与工程学院,江苏　南京210023
2.
河海大学水利水电学院,江苏　南京210024

基金项目:

国家杰出青年科学基金项目 41925012

国家自然科学基金项目 41572246

国家自然科学基金项目 41772280

国家自然科学基金项目 41902271

江苏省自然科学基金项目 BK20171228

江苏省自然科学基金项目 BK20170394

中央高校基本科研业务费专项资金项目

详细信息

作者简介:
林朱元（1997— ）,男,硕士研究生,主要从事地质工程及环境岩土工程方面的研究工作。E-mail：crecyslin@smail.nju.edu.cn

通讯作者:
唐朝生, E-mail：tangchaosheng@nju.edu.cn

中图分类号: TU43
计量
- 文章访问数: 411
- HTML全文浏览量: 43
- PDF下载量: 268
出版历程
- 收稿日期: 2019-03-27
- 网络出版日期: 2022-12-07
- 刊出日期: 2020-01-31

Laboratory characterization and discrete element modeling of desiccation cracking behavior of soils under different boundary conditions

1.
School of Earth Sciences and Engineering, Nanjing University, Nanjing 210023, China
2.
College of Water Conservancy and Hydropower Engineering, Hohai University, Nanjing 210024, China

摘要

摘要: 为了探究土体干缩开裂过程的边界效应问题,采用不同底面粗糙度的容器开展了多组干燥试验,发现干缩裂隙存在从顶面向下和从底面向上两种典型的发育形式。并且,裂隙发育程度与土样/容器界面接触条件密切相关,从而验证了裂隙发育过程的边界效应。通过理论分析,阐明了上边界的蒸发条件及下边界的接触条件对裂隙发育形式的控制作用。为了能更深入地理解土体干缩开裂边界效应的内在机制,在试验的基础上建立离散元模型,创新性地引入了沿深度的失水速率梯度参数,模拟土样上边界的蒸发条件变化。通过设置底面摩擦系数,模拟土样下边界的接触条件变化。将模拟结果与试验结果进行了对比分析,发现二者具有较好的吻合度。总体上,土体干缩裂隙的发育过程是顶面蒸发失水与底面摩擦两种边界条件共同作用的结果。当底面摩擦系数相对较小时,裂隙发育由蒸发失水主导,大部分裂隙由顶面向下发育。随着底面摩擦系数的增加,底面接触条件对裂隙发育过程的主导作用逐渐增强,由底面向上发育的裂隙数量所占比重也相应增加。
- 干缩开裂 /
- 边界条件 /
- 离散元模拟 /
- 底面摩擦 /
- 蒸发失水
Abstract: In order to explore the boundary effect of the desiccation cracking process, multiple sets of drying tests are carried out using the containers with different bottom roughnesses. Two different forming patterns can be observed in the laboratory tests, initiating from the top/bottom, and there is propagation closely related to the sample/container interface contact conditions. This verifies the boundary effect of the crack propagation. In order to understand the internal mechanism of the desiccation boundary effect of the soils more deeply, a discrete element model is established based on the drying tests. A water loss rate gradient parameter along the depth is introduced innovatively to simulate the change of the evaporation condition of the upper boundary of the soil samples. By setting the friction coefficient of the bottom surface, the contact condition of the lower boundary of the sample is simulated. The simulated results are compared with the experimental ones and found to have good agreement. In general, the initiation and propagation of desiccation cracks are the result of the combination of water loss due to surface evaporation and bottom friction. When the coefficient of friction of the bottom surface is relatively small, the development of the fracture is dominated by water loss, and most of the fractures develop from the top surface. With the increase of the friction coefficient of the bottom surface, the effect of the contact condition of the bottom surface on the development of the crack gradually increase, and the proportion of the number of cracks developed from the bottom surface increases accordingly.
- desiccation cracking /
- boundary conditions /
- DEM /
- bottom friction /
- evaporation

HTML全文

0. 引言

自然界中的岩体常年受复杂水文地质条件等环境因素影响，因此岩体内部结构通常存在着裂缝等不连续结构面，表现出较强的非连续性。岩体裂缝的力学效应作为岩体力学的中心内容，且岩体裂缝扩展和贯通的原理是预防和控制重大工程地质灾害的基础，所以研究岩体裂缝的力学效应具有十分重要的现实意义和理论价值。然而，基于连续介质的方法只能够模拟裂缝形成之前的断裂过程，无法模拟裂缝形成后岩石块体间的碰撞、运动过程^[1]。相反，基于不连续介质的数值方法，例如非连续变形分析方法(discontinuous deformation analysis，简称DDA)^[2]、数值流行方法(numerical manifold method，简称NMM)^[3]、离散单元法(discrete element method，简称DEM)^[4-5]和光滑粒子流体动力学方法(smoothed particle hydro- dynamics，简称SPH)^[6]可以作为有效的替代方法弥补此类不足，已广泛应用于分析岩体之间的相互作用^[7]。其中，离散单元法(DEM)因其理论的成熟性，已广泛应用于研究离散颗粒的力学行为^[8-9]。

球体单元因其简单有效的接触模式常被用于离散单元法的数值模拟。但是，它不能准确模拟不规则形状的块体。多面体单元因其形状和岩体具有相似性，在离散元中常被用于模拟不规则块体。但是，对于不同的接触模式，多面体单元缺少一个简单统一的数学模型。在计算块体单元之间的接触力之前，需要确定块体之间的接触模式，例如点-面接触、点-边接触、点-点接触、面-面接触、边-边接触和边-面接触^[10]。此外，由于顶点处法向接触力的不确定性，它不能解决点-点接触模式的问题。这也导致了法向接触力的不一致，进而带来能量失衡和数值不正确的问题。

Munjiza等^[11]提出了有限-离散单元耦合分析方法(combined finite-discrete element method，简称FDEM)，与传统离散单元法相比，能够在不考虑多变的接触模式的情况下实现接触力的计算。在该方法中，接触力是由分布力代替中心接触力进行计算，通过对嵌入区域的势函数进行积分，实现接触力的计算，既方便又统一，并保证了能量平衡和动量守恒。大量的应用实践已经证实了FDEM的优越性^[12-13]。但FDEM的计算过程是基于复杂的布尔运算，这会大幅度提高计算成本。此外，它对单元类型有严格的限制，复杂模型只能离散成四面体单元而不能离散成其他多面体单元，这会导致单元数量巨增^[14]，大幅度降低计算效率。

相比于FDEM，圆化多面体离散单元法具有计算效率高、接触力计算形式简单等优点。圆化多面体是由多面体和球体的闵可夫斯基和生成^[15]，是一个和多面体具有同等尺寸但角点已经圆化的多面体^[16]，因此两个多面体的接触类似于两球体之间的接触，简化接触模式的判断，提高了计算效率。虽然圆化多面体离散单元法是可行的，但由于该方法认为离散单元服从刚体假设，因此它并不能模拟软颗粒或弹性颗粒的变形。此外，该方法还缺少一个合适的模型来模拟断裂的全过程，包括岩体的变形、断裂和接触相互作用。颗粒的变形会影响岩体的应力分布，进而影响断裂破坏过程。准确模拟岩体的断裂全过程需要既能反映颗粒的相互作用，又能反映颗粒自身的变形。

为解决上述问题，提出一种可用于模拟断裂全过程的三维可变形圆化多面体离散单元法，并建立了一个三维断裂模型用于该方法预测断裂的演化过程。该方法将岩体离散成有限单元，并沿着有限单元的边界嵌入无厚度节理单元。采用虚拟裂缝模型^[17]和Mohr-Coulomb准则来判断节理单元的破坏状态。圆化多面体离散单元法解决了离散岩体之间的接触力问题，单元的应力应变问题则由有限单元法求解。与现有的圆化多面体离散单元法相比，该方法能够反映颗粒的应力应变，并能够模拟岩体的断裂全过程，包括变形、断裂和接触相互作用。与FDEM相比，该方法在计算接触力时避免了复杂的布尔运算，大幅度提高了计算效率。此外，它还突破四面体单元的限制，能够应用于更多类别的单元类型。

1. 三维可变形圆化多面体离散单元法

在本文提出的方法中，每个岩体用一个离散块体表示，并离散成多个有限单元来捕捉弹性响应。如图 1所示，在有限单元的边界嵌入无厚度节理单元用来模拟岩石材料的破坏过程。一旦节理单元处的应力达到破坏条件，就会在有限单元的表面产生新的裂缝。在此情况下，采用圆化多面体计算颗粒和断裂面之间的接触相互作用。

图 1 有限单元边界嵌入节理单元

Figure 1. Boundary of finite elements inserted into joint elements

下载: 全尺寸图片幻灯片

1.1 连续体到非连续体之间的过渡

断裂破碎通常被认为是从连续体到非连续体之间的过渡过程。节理单元和有限单元的结点相互重合，节理单元表面分离会产生应力 $\boldsymbol\sigma$ ，它是节理单元在两边界之间相对位移 $\boldsymbol\delta$ 的函数。节理单元的位移 $\boldsymbol\delta$ 和应力 $\boldsymbol\sigma$ 定义为

$\boldsymbol\delta ={({\boldsymbol\delta }_{\text{n}},{\boldsymbol\delta }_{\text{s}})}^{\text{T}}\text{，}$

(1)

$\boldsymbol\sigma =({\boldsymbol\sigma }_{\text{n}},{\boldsymbol\sigma }_{\text{s}})。$

(2)

式中： ${\boldsymbol\sigma _n}$ 为法向应力分量，与法向位移 ${\boldsymbol\delta _{\text{n}}}$ 相对应； ${\boldsymbol\sigma _{\text{s}}}$ 为切向应力分量，与切向位移 ${\boldsymbol\delta _{\text{s}}}$ 相对应。

根据断裂力学可知，节理单元的破坏模式可分为3种：模式Ⅰ(拉伸破坏)，模式Ⅱ(剪切破坏)，复合模式Ⅰ/Ⅱ(拉剪破坏)。本文研究的断裂模型是Hillerborg等所提出的虚拟裂缝模型^[17]，如图 2所示，节理单元的应力计算是基于应力-位移曲线完成的。如图 2(a)，(b)所示，在应力峰值后的曲线形状有相同的变化过程，而且峰后曲线可以用任意脆性材料的具体数据集来定义。

图 2 节理单元的本构特性

Figure 2. Constitutive characteristics of joint elements

下载: 全尺寸图片幻灯片

拉伸荷载作用下，每个节理单元的变化状态是由抗拉强度 ${f_t}$ 和模式Ⅰ的断裂能 $G{f_{\text{I}}}$ 决定，断裂能的值等于中曲线阴影部分的面积。法向应力 ${\sigma _{\text{n}}}$ 与法向位移 ${\delta _{\text{n}}}$ 之间的关系表示为

$\sigma_{\mathrm{n}}= \begin{cases}k_{\mathrm{n}} \delta_{\mathrm{n}} & \left(\delta_{\mathrm{n}}<w_{\mathrm{ft}}\right) \\ f_{\mathrm{t}}\left(1-\frac{\delta_{\mathrm{n}}-w_{\mathrm{ft}}}{w_{\mathrm{GI}}}\right) & \left(w_{\mathrm{ft}}<\delta_{\mathrm{n}}<w_{\mathrm{ft}}+w_{\mathrm{GI}}\right) 。 \\ 0 & \left(\delta_{\mathrm{n}}>w_{\mathrm{ft}}+w_{\mathrm{GI}}\right) 。\end{cases}$

(3)

式中： ${k_{\text{n}}}$ 为法向接触因子的值； ${w_{{\text{ft}}}}$ 为对应于抗拉强度 ${f_t}$ 法向位移的值； ${w_{{\text{ft}}}} + {w_{{\text{GI}}}}$ 为法向位移的最大值。

如所示，当法向应力 ${\sigma _n}$ 未达到材料的抗拉强度 ${f_{\text{t}}}$ ，它与法向位移之间的关系是线性的。当 ${\sigma _{\text{n}}}$ 的值达到抗拉强度 ${f_{\text{t}}}$ 时，它随着法向位移 ${\delta _{\text{n}}}$ 的增大而减小。如所示，一旦节理单元的法向位移超过最大法向位移 ${w_{{\text{ft}}}} + {w_{{\text{GI}}}}$ ，相互连接的两个有限单元就会完全分开，由此产生裂缝，块体间的相互作用由离散单元法计算。

图 3 岩体破裂过程的应力和位移

Figure 3. Stress and displacements of rock mass during failure process

下载: 全尺寸图片幻灯片

剪切载荷作用下，各个节理单元的状态是由剪切强度f_s和模式Ⅱ的断裂能 $G{f_{\text{\Pi }}}$ 决定，后者等于中曲线阴影部分的面积。模式Ⅱ断裂的法向应力 ${\sigma _{\text{s}}}$ 和切向位移 ${\delta _{\text{s}}}$ 之间的关系为

$\sigma_{\mathrm{s}}=\left\{\begin{array}{cc} k_{\mathrm{s}} \delta_{\mathrm{s}} & \delta_{\mathrm{s}}<w_{\mathrm{fs}}, \\ f_{\mathrm{s}}\left(1-\frac{\delta_{\mathrm{s}}-w_{\mathrm{fs}}}{w_{G \mathrm{GI}}}\right) & w_{\mathrm{fs}}<\delta_{\mathrm{s}}<w_{\mathrm{fs}}+w_{G \Pi} ，\\ -\sigma_{\mathrm{n}} \tan \varphi & \delta_{\mathrm{s}}>w_{\mathrm{fs}}+w_{G \Pi}。 \end{array}\right.$

(4)

式中：k_s为剪切接触因子的值； ${w_{{\text{fs}}}}$ 为对应于剪切强度 ${f_\text{s}}$ 的切向位移；剪切强度 ${f_\text{s}}$ 为由发生拉伸截断破坏的Mohr-Coulomb准则来定义的：

$f_{\mathrm{s}}= \begin{cases}-\sigma_{\mathrm{n}} \tan \varphi+c & \left(\sigma_{\mathrm{n}}<f_{\mathrm{t}}\right)， \\ -f_{\mathrm{t}} \tan \varphi+c & \left(\sigma_{\mathrm{n}} \geqslant f_{\mathrm{t}}\right)。\end{cases}$

(5)

式中： $c$ 为内部凝聚力； $\varphi$ 为材料的摩擦角。

当剪切应力 ${\sigma _{\text{s}}}$ 达到剪切强度f_s时，裂缝就会产生。随着切向位移 ${\delta _{\text{s}}}$ 的增加，切向应力 ${\sigma _{\text{s}}}$ 就会逐渐减小，直到切向位移超过最大剪切位移 $\begin{equation} w_{\mathrm{fs}}+w_{\mathrm{G} \Pi} \end{equation}$ 。如所示，切向应力 ${\sigma _{\text{s}}}$ 逐渐减小到与纯摩擦阻力相对应的残余值f_μ= $-\sigma_{\mathrm{n}} \tan \varphi$ 。

岩体的破坏是复合型破坏模式而不是单一的破坏模式。当发生复合型破坏时，法向位移或者切向位移可能还没达到单个模式断裂的破坏值，但断裂已经发生。图 2(c)为复合模式(Ⅰ/Ⅱ)的断裂模型，两个有限单元之间的最大分裂距离可以分解为水平轴和垂直轴上的法向位移和切向位移。法向位移和切向位移之间的关系若符合以下条件，复合模式的节理单元就会发生断裂：

${\left(\frac{{\delta }_{\text{n}}-{w}_{\text{ft}}}{{w}_{\text{GΙ}}}\right)}^{2}+{\left(\frac{{\delta }_{\text{S}}-{w}_{\text{ft}}}{{w}_{\text{GΠ}}}\right)}^{2}\ge 1。$

(6)

节理单元的应力沿着它的边界是呈非线性变化的，节理单元变形所产生的结点力 ${f_\text{joint}}$ 需进行数值积分：

${f}_{\text{joint}}={\displaystyle \int {N}^{\text{T}}}\sigma d\Omega 。$

(7)

式中：N为对应节理单元的形函数。

1.2 接触力计算

圆化多面体是由多面体和球体的闵可夫斯基和构成。圆化多面体生成过程如所示，多面体 $P$ 先收缩成多面体骨架 $H$ ，然后经球体 $S$ 膨胀得到圆化多面体 $SP$ ，由于收缩和膨胀过程都是基于相同的球半径开展，故圆化多面体与原始多面体具有相同尺寸，只是角点已被圆化^[16]。

图 4 多面体

$H$ 和球体

$S$ 的闵可夫斯基和构成

Figure 4. Polyhedral

$H$ and sphere

$S$ under Minkowski sum

下载: 全尺寸图片幻灯片

圆化多面体 $S{P_1}$ 有一组顶点 $\left\{ {V_1^i} \right\}$ ，边 $\left\{ {E_1^i} \right\}$ ，面 $\left\{ {F_1^i} \right\}$ 。圆化多面体 $S{P_2}$ 有一组 $\left\{ {V_2^i} \right\}$ ，边 $\left\{ {E_2^i} \right\}$ ，面 $\left\{ {F_2^i} \right\}$ 。假定 $S{P_1}$ 和 $S{P_2}$ 相互接触，如所示，它们之间的接触方式有点-点接触，点-边接触，点-面接触，边-边接触。 $S{P_1}$ 和 $S{P_2}$ 之间的法向接触力 $\begin{equation} \boldsymbol{P}_{\mathrm{n}}\left(S P_1, S P_2\right) \end{equation}$ ，它是考虑了所有接触方式的法向接触力的矢量和，可以由以下公式表示：

$\begin{aligned} \boldsymbol{P}_{\mathrm{n}}\left(S P_1, S P_2\right)= & \sum\limits_{\substack{i=1, N V_1 \\ j=1, N V_2}} \boldsymbol{P}_{\mathrm{n}}\left(V_1^i, V_2^j\right)+\sum\limits_{\substack{i=1, N E_1 \\ j=1, N E_2}} \boldsymbol{P}_{\mathrm{n}}\left(E_1^i, E_2^j\right)+ \\ & \left(\sum\limits_{\substack{i=1, N V_1 \\ j=1, N V_2}} \boldsymbol{P}_{\mathrm{n}}\left(V_1^i, E_2^j\right)+\sum\limits_{\substack{i=1, N E_1 \\ j=1, N V_2}} \boldsymbol{P}_{\mathrm{n}}\left(E_1^i, V_2^j\right)\right)+ \\ & \left(\sum\limits_{\substack{i=1, N V_1 \\ j=1, N F_2}} \boldsymbol{P}_{\mathrm{n}}\left(V_1^i, F_2^j\right)+\sum\limits_{\substack{i=1, N F_1 \\ j=1, N V_2}} \boldsymbol{P}_{\mathrm{n}}\left(F_1^i, E_2^j\right)\right) 。 \end{aligned}$

(8)

图 5 接触对的几何模型

Figure 5. Geometrically-shaped contact pair model

下载: 全尺寸图片幻灯片

式中： $N{V_1},N{E_1},N{F_1}$ 分别为圆化多面体 $S{P_1}$ 顶点、边和面的数量； $N{V_2},N{E_2},N{F_2}$ 分别为 $S{P_2}$ 顶点、边和面的数量。

此外， $S{P_1}$ 由法向接触力引起的力矩 ${\boldsymbol{M}_n}$ 表示为

$\boldsymbol{M}_{\text{n}}^{}=\boldsymbol{P}_{\text{n}}^{}\times \boldsymbol{r}_{\text{n}}^{}。$

(9)

式中： ${r_n}$ 为法向接触力 ${P_n}$ 的作用点位置到 $S{P_1}$ 质心的矢量，可以由下式计算得到

$\boldsymbol{r}_{n}=\frac{{\displaystyle \sum _{i=1}^{m}\left({r}_{n}^{i}\left|{P}_{n}^{i}\right|\right)}}{{\displaystyle \sum _{i=1}^{m}\left|{P}_{n}^{i}\right|}}。$

(10)

式中： $m$ 为 $S{P_1}$ 接触对的数量， $P_n^i$ 为 $S{P_1}$ 一个接触对得到的法向接触力。

法向接触力的作用点 ${p_0}({x_0},{y_0},{z_0})$ 位置是根据力矩平衡原理来确定的

${p}_{0}({x}_{0},{y}_{0},{z}_{0})={r}_{n}+{p}_{\text{c}}({x}_{\text{c}},{y}_{\text{c}},{z}_{\text{c}})。$

(11)

式中： ${p_0}({x_0},{y_0},{z_0})$ 和 ${p_c}({x_c},{y_c},{z_c})$ 分别为法向接触力作用点的位置和 $S{P_1}$ 质心的位置。

切向接触力作用点的位置与法向接触力的相同。通常情况下，法向接触力合力的矢量方向在每个时间步长与切向接触力是相互垂直的。在 $t + \Delta t$ 时间内 $S{P_1}$ 的切向接触力 $\boldsymbol{P}_{\text{s}}^{t + \Delta t}$ 可以表示为

$\boldsymbol{P}_{\text{s}}^{t+\Delta t}=\boldsymbol{R}_{\text{ro}}\times \boldsymbol{P}_{\text{s}}^{t}+{K}_{\text{s}}\cdot \Delta {\delta }_{\text{s}}^{t+\Delta t}。$

(12)

式中： ${K_{\text{s}}}$ 为切向接触刚度； $\Delta \boldsymbol{\delta} _{\text{s}}^{t + \Delta t}$ 为切向位移增量； ${\boldsymbol{R}_{{\text{ro}}}}$ 为旋转矩阵； $\boldsymbol{P}_{\text{s}}^t$ 为 $t$ 时刻的切向力。

切向力 $\boldsymbol{P}_s^{t + \Delta t}$ 用Mohr-Coulomb准则校核，表示为

$\boldsymbol{P}_{\text{s}}^{t+\Delta t}=\mathrm{min}\left(\left|\boldsymbol{P}_{\text{s}}^{t+\Delta t}\right|,\mu \left|\boldsymbol{P}_{\text{n}}^{t+\Delta t}\right|\right)\cdot \boldsymbol{n}_{\text{s}}。$

(13)

式中： $\mu$ 为摩擦系数；n_s为 $t + \Delta t$ 时刻切向方向的单位向量。

由切向接触力产生的力矩 ${\boldsymbol{M}_{\text{s}}}$ 表示为

$\boldsymbol{M}_{\text{s}}=\boldsymbol{P}_{\text{s}}\times \boldsymbol{r}_{\text{s}}。$

(14)

式中： ${\boldsymbol{r}_{\text{s}}}$ 为切向作用力P_s作用点到 $S{P_1}$ 质心的矢量。

$S{P_1}$ 接触力的结果表示为

$\boldsymbol{P}_{\text{contact}}=\boldsymbol{P}_{\text{n}}+\boldsymbol{P}_{\text{s}}。$

(15)

1.3 等效结点力

本文提出的方法中，每个块体都离散成有限单元来模拟块体的变形情况，为更好捕捉各离散块体之间的力学响应，需要将接触单元间的接触力转换成有限单元的等效结点力。

如所示，计算 $S{P_1}$ 和 $S{P_2}$ 两个圆化多面体接触单元之间的接触力时，会忽略角点处的接触力。故在求解等效结点力时，如所示，接触单元应当恢复原始形状，等效结点力采用虚功原理进行换算， $S{P_1}$ 单元的等效结点力，

$\boldsymbol{f}_{\text{contact}}^{i}={N}_{S{P}_{1}}^{i}\boldsymbol{P}_{\text{contact}},\text{ }i=1,2,\cdots ,{n}_{S{P}_{1}}。$

(16)

式中： $N_{S{P_1}}^i$ 为 $S{P_1}$ 单元的形函数， ${n_{S{P_1}}}$ 为接触单元 $S{P_1}$ 的结点数。

图 6 虚功法的换算方法

Figure 6. Conversion method of virtual work principle

下载: 全尺寸图片幻灯片

1.4 应力应变计算

几何非线性，例如大位移和大应变，常发生于单个块体的运动中。大应变等几何非线性问题通过非线性有限单元法求解，本节用改进的拉格朗日公式(U.L.)来解释说明可变形块体的增量平衡方程^[18-19]。所有单元 $e$ 的内部结点力 $\boldsymbol{f}_{{\text{in}}}^{t + \Delta t}$ 可由下式计算：

$\boldsymbol{f}_{\text{in}}^{t+\Delta t}={\displaystyle \int {}_{{V}^{e}}}\boldsymbol{B}^{\text{T}}\boldsymbol{\tau }^{t+\Delta t}\text{d}v。$

(17)

式中：B为应变-位移变换矩阵。 $t + \Delta t$ 时刻的柯西应力 ${\tau ^{t + \Delta t}}$ 可以表示为

${\tau }^{t+\Delta t}=({\tau }^{t}+\Delta \boldsymbol{S})。$

(18)

式中： ${\tau ^t}$ 表示在 $t$ 时刻的应力。 $\Delta S$ 表示在 $t$ 到 $t + \Delta t$ 时刻应力增量，可按下式计算：

$\begin{equation} \Delta S=\boldsymbol{D} \Delta \varepsilon \end{equation}$

(19)

式中：D为弹性矩阵， $\Delta \varepsilon$ 表示在 $t$ 到 $t + \Delta t$ 时刻应变增量。

$t + \Delta t$ 时刻应变 ${\varepsilon ^{t + \Delta t}}$ 可表示为

$\begin{equation} \varepsilon^{t+\Delta t}=\varepsilon^t+\Delta \varepsilon \end{equation}$

(20)

式中： ${\varepsilon ^t}$ 表示在 $t$ 时刻的应变。

1.5 控制方程

由1.4节提到的内部节点力 ${\boldsymbol{f}_{{\text{in}}}}$ 和外力 ${\boldsymbol{f}_{{\text{ext}}}}$ 的计算方法可知，单个颗粒的运动是由作用在它上面的不平衡力决定的^[11]。单个颗粒的运动方程为

$\boldsymbol{Ma}=\boldsymbol{f}_{\text{ext}}-\boldsymbol{f}_{\text{in}}-\boldsymbol{Cv}。$

(21)

式中：M和C分别为质量矩阵和瑞利阻尼矩阵；v和a分别为速度和加速度的矢量； ${\boldsymbol{f}_{{\text{in}}}}$ 和 ${\boldsymbol{f}_{{\text{ext}}}}$ 分别为内力和外荷载的矢量。

外力 ${\boldsymbol{f}_{{\text{ext}}}}$ 的矢量包括接触力矢量 ${\boldsymbol{f}_{contact}}$ ，外荷载矢量 ${\boldsymbol{f}_{{\text{load}}}}$ 和节理单元的结点力 ${\boldsymbol{f}_{{\text{joint}}}}$ ：

$\boldsymbol{f}_{\text{ext}}=\boldsymbol{f}_{\text{joint}}+\boldsymbol{f}_{\text{contact}}+\boldsymbol{f}_{\text{load}}。$

(22)

运动方程采用显式求解格式，结合集中质量法进行求解，可以避免整体刚度矩阵的计算。采用Velocity Verlet算法对运动方程进行积分，能够稳定求解可变形块体的位移和速度^[20]。假定在时间步 $t$ 完成计算，算法将会执行 $t + \Delta t$ 时间步的计算， $t + \Delta t$ 时刻的运动方程可表示为

$\boldsymbol{M{a}}^{t+\Delta t}=\boldsymbol{f}_{\text{ext}}^{t+\Delta t}-\boldsymbol{f}_{\text{in}}^{t+\Delta t}-\boldsymbol{C{v}}^{t+\Delta t}。$

(23)

本文所提方法的技术路线如图 7所示，下面列出了进行模拟计算的步骤：

图 7 本文方法的技术方案

Figure 7. Technical proposal of proposed method

下载: 全尺寸图片幻灯片

(1) 建立几何模型。将块体分解为有限单元并在其边界嵌入节理单元。输入外荷载和约束条件。

(2) 计算节理单元的应力。采用虚拟裂缝模型和Mohr-Coulomb准则来判断节理单元的破坏状态。当应力达到破坏条件时，块体发生断裂。

(3) 对分离的块体和裂缝表面进行接触检测。圆化多面体用来计算作用在它们之间的接触力。

(4) 将接触力转换成有限单元的等效结点力。在此步骤中，单元不再是圆化多面体，而是恢复到单元的初始形状。

(5) 采用Velocity Verlet算法求解运动方程。

(6) 更新所有有限单元的结点坐标。

2. 数值试验

为验证本文提出的模拟断裂全过程的三维可变形圆化多面体离散单元法的准确性和可行性，本节给出5个算例加以佐证。首先，通过对颗粒堆积试验的模拟，验证了该方法计算的高效性；然后通过巴西圆盘劈裂试验验证了所提出的方法能够模拟脆性材料的断裂过程；随后，通过单缺口梁裂缝扩展试验的数值模拟验证了该方法能够捕捉裂缝的萌生和扩展；最后，通过四点弯曲梁试验和球体冲击试验验证了该方法可以模拟岩石承受冲击载荷后断裂的过程。

2.1 颗粒堆积试验模拟

本节通过模拟四面体颗粒在16.0 m×13.8 m×15.2 m的刚性盒子里的堆积过程，来验证本文提出方法的高效性。为便于比较，采用FDEM对同一个算例进行了模拟。四面体颗粒材料属性为：密度 $\rho$ 为2000 kg/m³，阻尼比 $\xi$ 为0.0005。两种方法均采用相同的计算时间，时间步长和迭代次数，所有结果均在同一台电脑上获得。

采用上述两种方法模拟5346个四面体颗粒的最终堆积形式如图 8所示，FDEM模拟结果如图 8(a)所示，本文提出方法的模拟结果图 8(b)所示，两种方法模拟出来的结果相似。本文提出的方法模拟消耗的时间为62060 s，FDEM模拟消耗的时间为117813 s，节省将近一半的时间。本文提出的方法在颗粒间的接触采用非均匀颗粒快速线性接触检测算法，该算法的复杂度为O(n)，当颗粒增加时，会节省大量的接触检测时间，计算效率线性增加。因此，该方法在计算效率方面大幅度提高。

图 8 5346个四面体颗粒最终堆积形式

Figure 8. Final packing morphology exhibited by 5346 tetrahedral particles

下载: 全尺寸图片幻灯片

2.2 巴西圆盘劈裂试验模拟

巴西圆盘劈裂试验是一个间接测量脆性材料抗拉强度简单有效的方式。试验的三维模型是由一个圆盘与两块板组成。如所示，在顶板和底板施加恒定的垂直速度 ${v_y} = 0.02$ m/s，该圆盘被离散成12335个有限单元。该试验的材料属性为：杨氏模量 $E =$ 25.8 GPa，泊松比 $\nu =$ 0.15，抗拉强度 ${f_\text{t}} =$ 2.3 MPa，Ⅰ型断裂能 $G{f_{\text{I}}} =$ 10 N/m, 凝聚力 $c =$ 20 MPa，内摩擦角 $\varphi =$ 45°，Ⅱ型断裂能 $\begin{equation} G f_{\Pi}=1500 \mathrm{~N} / \mathrm{m} \end{equation}$ 。

图 9 模型网格剖分图

Figure 9. Mesh configuration of model

下载: 全尺寸图片幻灯片

本文提出的方法模拟圆盘的断裂过程如所示。在断裂发生之前，拉应力沿着圆盘的垂直直径分布，加载板附近有小范围的压应力分布。在 $t$ =1 ms，如所示，未观察到任何裂缝。在 $t$ =57 ms，如所示，裂缝出现在圆盘的纵轴上。在 $t$ =58 ms，如所示，圆盘处于完全破坏的状态。为更清晰地了解圆盘沿直径方向上的应力分布，进一步分析试样表面和中心的两条压缩径线(AB和CD)上的拉应力 ${\sigma _{xx}}$ 和压应力 ${\sigma _{yy}}$ 。直径线AB上的拉应力 ${\sigma _{xx}}$ 和压应力 ${\sigma _{yy}}$ 的理论值可按下式计算^[21]：

图 10 本文提出方法预测圆盘的断裂演化过程

Figure 10. Prediction of fracture evolution process of disk by proposed method

下载: 全尺寸图片幻灯片

${\sigma }_{xx}=\frac{2P}{\text{π}DL}\text{，}$

(24)

${\sigma }_{yy}=\frac{2P}{\text{π}DL}\left(1-\frac{4{D}^{2}}{{D}^{2}-4{y}^{2}}\right)。$

(25)

圆盘试验在33.0 kN荷载下沿轴线上的应力分布如所示， ${\sigma _{xx}}$ 和 ${\sigma _{yy}}$ 的数值计算和解析解的结果是一致的。外荷载施加在上板处和试件中心处拉应力 ${\sigma _{xx}}$ 的结果如图 12所示，上板处的最大荷载为40.6 kN，用公式求得的抗拉强度的数值结果是2.54 MPa，与输入值相比减小了1.57%，结果差别相对较小，说明该方法模拟脆性材料的断裂过程是可行的。

图 11 本文提出方法进行巴西圆盘数值试验

Figure 11. Numerical tests of Brazilian disk by proposed method

下载: 全尺寸图片幻灯片

图 12 上板处的外荷载和试件中心处的拉应力

Figure 12. External loads on upper plate and tensile stresses at center of specimen

下载: 全尺寸图片幻灯片

2.3 单缺口梁裂缝扩展模拟

本节对单缺口梁的裂缝扩展试验进行模拟，验证本文提出的方法能够捕捉裂缝的萌生与扩展。Gálvez等^[22]进行的试验研究被认为是基准性研究，梁的几何形状、边界条件和网格结构如所示。梁的下部由两个支座点支撑，当上部施加恒定速度 ${v_y}$ =0.005 m/s，梁会发生变形。单缺口位于梁的底边中心，裂缝将在这块区域内产生。梁的材料属性为：杨氏模量 $E$ =41.8 GPa，泊松比 $\nu$ =0.2，抗拉强度 ${f_\text{t}}$ =30.0 MPa，凝聚力 $c$ =15.59 MPa，内摩擦角 $\varphi$ =30^°，Ⅰ型断裂能 $G{f_{\text{I}}} =$ 69 N/m，Ⅱ型断裂能 $G{f_{{\Pi }}} =$ 1380 N/m。

图 13 梁的网格结构

Figure 13. Mesh configuration of beam

下载: 全尺寸图片幻灯片

单缺口梁断裂之后的几何形状如图 14所示。如图 15所示，数值模拟得到的裂缝轨迹与试验结果是一致的。如图 16所示，荷载-裂缝口张开位移曲线表示该方法的结果与试验记录的结果相似。因此，该方法可以有效模拟拉剪混合应力作用下的裂缝扩展。

图 14 单缺口梁的变形结构(放大100倍)

Figure 14. Deformation struture of single notched beam (Amplified 100 times)

下载: 全尺寸图片幻灯片

图 15 裂缝轨迹数值模拟结果与试验结果的对比

Figure 15. Comparison between numerical and experimental results of crack trajectory

下载: 全尺寸图片幻灯片

图 16 荷载与CMOD数值模拟与试验结果的对比

Figure 16. Comparison of numerical simulation and experimental results of load and CMOD

下载: 全尺寸图片幻灯片

2.4 四点弯曲梁的冲击试验模拟

以四点弯曲梁为例进行冲击试验，验证本文提出来的方法能够模拟承受冲击载荷后混凝土材料的断裂状态。实验室试验是由Murray等^[23]开展进行的，在该试验中，两块金属块以 ${v_y}$ =5 m/s的速度撞击混凝土简支梁。如所示，梁试件的尺寸为1524 mm× 114.3 mm×100 mm，梁的顶部和底部的块体简化为半径为16 mm的圆柱体，混凝土梁被离散成4080个六面体单元。梁的材料属性为：杨氏模量 $E$ =25.8 GPa，泊松比 $\nu$ =0.15，抗拉强度 ${f_\text{t}}$ =2.7 MPa，Ⅰ型断裂能 $G{f_{\text{I}}} =$ 98 N/m。为确保发生Ⅰ型断裂，将剪切强度设置得足够大。

图 17 模型的几何形状

Figure 17. Geometry of model

下载: 全尺寸图片幻灯片

模拟混凝土梁的破坏过程如图 18所示，裂缝产生于冲击块正下方梁的下表面，然后梁端1/4处出现裂缝。最后，梁断裂成5个部分，这与试验结果完全一致，如图 19所示。如图 20所示，梁在断裂过程中最大挠度曲线数值模拟的计算结果与试验结果是一致的，表明数值模拟结果是可靠的。

图 18 本文提出方法预测混凝土梁的断裂变化过程

Figure 18. Change process of fracture of concrete beam by proposed method

下载: 全尺寸图片幻灯片

图 19 素混凝土最终断裂成5个部分^[23]

Figure 19. Eventually fractured five parts of plain concrete beams

下载: 全尺寸图片幻灯片

图 20 梁的最大挠度数值模拟与试验结果的对比

Figure 20. Comparison between numerical and experimental results of maximun deflection

下载: 全尺寸图片幻灯片

2.5 球体冲击模拟

该算例涉及到岩石系统的瞬态动力学，将直径为200 mm的岩石球体放置在刚性板附近，如所示，在负Y方向以不同的速度与刚性板发生碰撞，将球体离散为4539个四面体单元。球体的材料属性：杨氏模量 $E$ =50.0 GPa，密度 $\rho$ =2500 kg/m³，泊松比 $\nu$ =0.2，抗拉强度 ${f_\text{t}}$ =3.0 MPa，凝聚力 $c$ =55 MPa，内摩擦角 $\varphi$ =30°，Ⅰ型断裂能 $G{f_{\text{I}}} =$ 80 N/m，Ⅱ型断裂能 $G{f_{\Pi}} =$ 200 N/m。

图 21 模型的几何形状

Figure 21. Geometry of model

下载: 全尺寸图片幻灯片

球体在不同冲击速度下的破裂形态如图 22所示，当冲击速度小于5 m/s时，观察不到明显裂缝。当速度为7.5 m/s时，可以识别出一些裂缝，此时球体被分割成几个比较大的区域。一旦冲击速度达到20 m/s，球体底部破碎为多个块体。随着冲击速度的进一步增大，球体的破坏则更加明显，在球体底部可以观察到更多的破碎块体。

图 22

$t$ =0.75 ms时刻不同冲击速度下的断裂形态

Figure 22. Fracture patterns under different impact velocities when

$t$ =0.75 ms

下载: 全尺寸图片幻灯片

如图 23所示，采用FDEM获得的混凝土球体破碎模式和Khanal等^[24]观察到的实验室结果与本文提出方法模拟的结果相似。结果表明，该方法能够成功获得承受冲击载荷后岩体裂缝的主要特征。

图 23 破碎模式

Figure 23. Broken petterns

下载: 全尺寸图片幻灯片

3. 结论

(1) 本文提出了一种用于模拟断裂全过程的三维可变形圆化多面体离散单元法。该方法将圆化多面体离散单元法与有限单元法相结合，有限单元法计算颗粒的应力应变，圆化多面体离散单元法计算颗粒之间的接触相互作用，在有限单元的边界嵌入无厚度节理单元，通过判断节理单元的破坏状态来捕捉裂缝萌生与扩展，可以有效模拟岩体断裂的全过程。

(2) 通过5个算例来验证该方法的有效性和可行性，其中颗粒堆积试验模拟结果表明，该方法比FDEM更高效；巴西圆盘劈裂试验验证了该方法在处理拉伸问题的准确性；单缺口梁试验证明了该方法能够捕捉到拉剪复合应力下的断裂裂缝扩展；最后，四点弯曲梁试验和球体冲击试验表明，该方法能够准确捕捉岩体受冲击载荷后裂缝的萌生和扩展。因此，该方法为模拟准脆性材料的断裂全过程提供了一个简单高效的工具。

图 1 裂隙从顶面向下发育过程（E1组侧面局部放大）

Figure 1. Downward propagation of crack from top (side view of E1, enlarged)

下载: 全尺寸图片幻灯片

图 2 裂隙从底面向上发育过程（E1组侧面局部放大）

Figure 2. Upward propagation of crack from bottom (side view of E1, enlarged)

下载: 全尺寸图片幻灯片

图 3 室内试验最终裂隙图像（侧面）

Figure 3. Final crack patterns in laboratory tests (side view)

下载: 全尺寸图片幻灯片

图 4 contact bond模型的力学性质

Figure 4. Mechanical performance of contact bond model

下载: 全尺寸图片幻灯片

图 5 离散元模型初始状态

Figure 5. Original state of DEM model

下载: 全尺寸图片幻灯片

图 6 模拟结束时模型中的土样开裂情况

Figure 6. Final crack patterns of DEM samples

下载: 全尺寸图片幻灯片

图 7 3组摩擦系数对应不同失水速率梯度条件下的裂隙发育情况

Figure 7. Variation of crack development with water loss gradient under three different friction coefficients

下载: 全尺寸图片幻灯片

图 8 （a）模拟中裂隙从顶面向下发育情况；（b）模拟中裂隙从底面向上发育情况（线段越粗表示颗粒间连接受力越.大,应力集中越明显）

Figure 8. Simulated results of crack initiating from (a) top and crack initiating from bottom (b) simulation results (greater line weight indicates stronger force in contact)

下载: 全尺寸图片幻灯片

表 1 室内试验参数

Table 1 Parameters of laboratory tests

试样编号砂纸规格/目砂颗粒粒径/mm 摩擦系数初始厚度/cm

E1 80 0.180 大 1
E2 120 0.120 中 1
E3 240 0.063 小 1

下载: 导出CSV

表 2 室内试验裂隙条数统计

Table 2 Number of cracks in laboratory tests

试样编号摩擦系数总裂隙条数顶面发育裂隙条数底面发育裂隙条数

E1 大 9 5 4
E2 中 5 3 2
E3 小 3 2 1

下载: 导出CSV

表 3 数值模拟裂隙条数统计（λ=0.05）

Table 3 Number of cracks in numerical tests (λ=0.05)

编号摩擦系数总裂隙条数顶面发育裂隙条数底面发育裂隙条数

① 0.3 8 5 3
② 0.2 5 4 1
③ 0.1 3 3 0

下载: 导出CSV

参考文献(30)

[1]	LOZADA C, THOREL L, CAICEDO B. Effects of cracks and desiccation on the bearing capacity of soil deposits[J]. Géotechnique Letters, 2015, 5(3): 112-117. doi: 10.1680/jgele.15.00021
[2]	袁俊平, 殷宗泽. 膨胀土裂隙的量化指标与强度性质研究[J]. 水利学报, 2004, 35(6): 108-113. https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB200406019.htm YUAN Jun-ping, YIN Zong-ze. Quantitative index of fissure and strength characteristics of fissured expansive soils[J]. Journal of Hydraulic Engineering, 2004, 35(6): 108-113. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB200406019.htm
[3]	姚海林, 郑少河, 葛修润, 等. 裂缝膨胀土边坡稳定性评价[J]. 岩石力学与工程学报, 2002, 21(增刊2): 2331-2335. YAO Hai-lin, ZHENG Shao-he, GE Xiu-run, et al. Assessment on slope stability in cracking expansive soils[J]. Chinese Journal of Rock Mechanics and Engineering, 2002, 21(S2): 2331-2335. (in Chinese)
[4]	MITCHELL J K. Fundamentals of Soil Behavior[M]. New York: Wiley, 1993.
[5]	LAKSHMIKANTHAM. R, PRATPERE C, ALBERTO LEDESMA. Experimental evidence of size effect in soil cracking[J]. Canadian Geotechnical Journal, 2008, 49(3): 264-284.
[6]	唐朝生, 施斌, 顾凯. 土中水分的蒸发过程试验研究[J]. 工程地质学报, 2011, 19(6): 875-881. https://www.cnki.com.cn/Article/CJFDTOTAL-GCDZ201106011.htm TANG Chao-sheng, SHI Bin, GU Kai. Experimental investigation on Evaporation process of water in soil during drying[J]. Journal of Engineering Geology, 2011, 19(6): 875-881. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-GCDZ201106011.htm
[7]	唐朝生, 施斌, 刘春. 膨胀土收缩开裂特性研究[J]. 工程地质学报, 2012, 20(5): 663-673. https://www.cnki.com.cn/Article/CJFDTOTAL-GCDZ201205004.htm TANG Chao-sheng, SHI Bin, LIU Chun. Study on desiccation cracking behaviour of expansive soil[J]. Journal of Engineering Geology, 2012, 20(5): 663-673. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-GCDZ201205004.htm
[8]	唐朝生, 王德银, 施斌, 等. 土体干缩裂隙网络定量分析[J]. 岩土工程学报, 2013, 35(12): 2298-2305. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201312026.htm TANG Chao-sheng, WANG De-yin, SHI Bin, et al. Quantitative analysis of soil desiccation crack network[J]. Chinese Journal of Geotechnical Engineering, 2013, 35(12): 2298-2305. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201312026.htm
[9]	唐朝生, 施斌, 崔玉军. 土体干缩裂隙的形成发育过程及机理[J]. 岩土工程学报, 2018, 40(8): 1415-1423. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201808008.htm TANG Chao-sheng, SHI Bin, CUI Yu-jun. Behaviors and mechanisms of desiccation cracking of soils[J]. Chinese Jouranl Geotechnical Engineering, 2018, 40(8): 1415-1423. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201808008.htm
[10]	PERON H, HUECKEL T, LALOUI L, et al. Fundamentals of desiccation cracking of fine-grained soils: experimental characterisation and mechanisms identification[J]. Canadian Geotechnical Journal, 2009, 46(10): 1177-1201. doi: 10.1139/T09-054
[11]	WEINBERGER R. Initiation and growth of cracks during desiccation of stratified muddy sediments[J]. Journal of Structural Geology, 1999, 21(4): 379-386. doi: 10.1016/S0191-8141(99)00029-2
[12]	曾浩, 唐朝生, 林銮. 土体干缩裂隙发育方向及演化特征的层间摩擦效应研究[J]. 岩土工程学报, 2019, 41(6): 1172-1180. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201906026.htm ZENG Hao, TANG Chao-sheng, LIN Luan. Interfacial friction dependence of propagation direction and evolution characteristics of soil desiccation cracks[J]. Chinese Journal of Geotechnical Engineering, 2019, 41(6): 1172-1180. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201906026.htm
[13]	CUNDALL P A, STRACK O D L. A discrete numerical mode for granular assemblies[J]. Géotechnique, 1979, 29(1): 47-65.
[14]	蒋明镜, 胡海军, 彭建兵. 结构性黄土一维湿陷特性的离散元数值模拟[J]. 岩土力学, 2013, 34(4): 1121-1130. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201304036.htm JIANG Ming-jing, HU Hai-jun, PENG Jian-bing. Simulation of collapsible characteristics of structural loess under one-dimensional compression condition by discrete element method[J]. Rock and Soil Mechanics, 2013, 34(4): 1121-1130. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201304036.htm
[15]	张程林, 周小文. 砂土颗粒三维形状模拟离散元算法研究[J]. 岩土工程学报, 2015, 37(增刊1): 115-119. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC2015S1024.htm ZHANG Cheng-lin, ZHOU Xiao-wen. Algorithm for modelling three-dimensional shape of sand based on discrete element method[J]. Journal of Engineering Geology, 2015, 37(S1): 115-119. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC2015S1024.htm
[16]	PERON H, DELENNE J Y, LALOUI L, et al. Discrete element modelling of drying shrinkage and cracking of soils[J]. Computers & Geotechnics, 2009, 36(1/2): 61-69.
[17]	司马军, 蒋明镜, 周创兵. 黏性土干缩开裂过程离散元数值模拟[J]. 岩土工程学报, 2013, 35(增刊2): 286-291. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC2013S2049.htm SIMA Jun, JIANG Ming-jing, ZHOU Chuang-bing, Numerical simulation of desiccation cracking of clay soils by DEM[J]. Chinese Journal of Geotechnical Engineering, 2003, 35(S2): 286-291. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC2013S2049.htm
[18]	张晓宇, 许强, 刘春, 等. 黏性土失水开裂多场耦合离散元数值模拟[J]. 工程地质学报, 2017, 25(6): 1430-1437. ZHANG Xiao-yu, XU Qiang, LIU Chun, et al. Numerical simulation of drying cracking using multifield coupling discrete element method[J]. Journal of Engineering Geology, 2017, 25(6): 1430-1437. (in Chinese)
[19]	COSTA S, KODIKARA J, SHANNON B. Salient factors controlling desiccation cracking of clay in laboratory experiments[J]. Géotechnique, 2013, 63(1): 18-29.
[20]	KODIKARA J, CHOI X. A simplified analytical model for desiccation cracking of clay layers in laboratory tests[C]//Proceedings of the 4th International Conference on Unsaturated Soil. New York, 2006: 2558-2569.
[21]	SIMA J, 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.
[22]	GUO M Y, HAN M C, YU X. Laboratory characterization and discrete element modeling of shrinkage[J]. Canadian Geotechnical Journal, 2016, 14: 5-13.
[23]	冉龙洲, 宋翔东, 唐朝生. 干燥过程中膨胀土抗拉强度特性研究[J]. 工程地质学报, 2011, 19(4): 620-625. https://www.cnki.com.cn/Article/CJFDTOTAL-GCDZ201104029.htm RAN Long-zhou, SONG Xiang-dong, TANG Chao-sheng. Laboratorial investigation on tensile strength of expansive soil during drying[J]. Journal of Engineering Geology, 2011, 19(4): 620-625. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-GCDZ201104029.htm
[24]	El YOUSSOUFI M S, DELENNE J Y, RADIAI F. Self-stresses and crack formation by particle swelling in cohesive granular media[J]. Physical Review E Statistical Nonlinear & Soft Matter Physics, 2005, 71(51): 5-7.
[25]	刘昌黎, 唐朝生, 李昊达, 等. 界面粗糙度对土体龟裂影响的试验研究[J]. 工程地质学报, 2017, 25(5): 1314-1321. https://www.cnki.com.cn/Article/CJFDTOTAL-GCDZ201705018.htm LIU Chang-li, TANG Chao-sheng, LI Hao-da, et al. Experimental study on eddect of interfacial roughness on desiccation cracking behavior of soil[J]. Journal of Engineering Geology, 2017, 25(5): 1314-1321. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-GCDZ201705018.htm
[26]	曾浩, 唐朝生, 刘昌黎, 等. 控制厚度条件下土体干缩开裂的界面摩擦效应[J]. 岩土工程学报, 2019, 41(3): 544-553. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201903021.htm ZENG Hao, TANG Chao-sheng, LIN Chang-li, et al. Effect of boundary friction and layer thickness on soil desiccation cracking behavior[J]. Chinese Journal of Geotechnical Engineering, 2018, 41(3): 544-553. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201903021.htm
[27]	ITASCA. PFC3D (Particle Flow Code in 3 Dimensions), Version 3.00[R]. Minneapolis: Itasca Consulting Group, Inc.; 2003.
[28]	AMARASIRI A L, KODIKARA J K, SUSANGA C. Numerical modelling of desiccation cracking[J]. International Journal for Numerical & Analytical Methods in Geomechanics, 2011, 35(1): 82-96.
[29]	JIANG M J, KONRAD J M, LEROUEIL S. An efficient technique to generate homogeneous specimens for DEM studies[J]. Computers and Geotechnics, 2003, 30(7): 579-597.
[30]	HILLIL D. Introduction to Environmental Soil Physics[M]. San Diego, CA: Elsevier Science, 2004.