DEM simulation of desiccation cracking in clay considering capillarity and adsorption
-
摘要: 黏土干缩开裂过程中土颗粒间的作用非常复杂，导致难以从颗粒层面对干缩裂缝的形成机理开展定量研究。基于非饱和土中的吸应力曲线，明确区分黏土中的毛细与吸附作用，在此基础上建立考虑黏土颗粒间吸引力随含水率变化的离散元接触模型，开展黏土干缩开裂的数值模拟，并与室内试验结果进行对比。结果表明，数值模拟得到的裂缝形态、裂缝发育历程、土样应变均与室内试验较为吻合，验证了离散元接触模型的可靠性。对模拟结果的进一步分析发现：①在黏土干缩开裂过程中，毛细与吸附作用共同发挥效果。随着含水率的降低，吸附作用的影响逐渐超过毛细作用。在吸附作用占优阶段，模拟土样的土颗粒沿裂缝分布方向的平均位移为最终位移值的73%。②土颗粒间总接触数量随含水率的降低先减少后增加。③在毛细/吸附作用占优阶段，土颗粒间的物理接触力会迅速增加并引起应力集中，土颗粒间接触因此集中断开进而形成裂缝。Abstract: The interaction between clay particles during the process of desiccation cracking is highly complex, making it challenging to conduct quantitative researches on the formation mechanism of cracks at the particle level. The capillarity and adsorption in clay is distinguished based on the suction stress characteristic curve (SSCC) in unsaturated soils. A contact model of discrete element method (DEM) that accounts for the change of attraction between clay particles with water content is then established, and the numerical simulation of desiccation cracking in clay is carried out. The simulated results are compared with those of laboratory tests, and the findings indicate that the crack morphology, crack development history and strain of the soil sample obtained by the DEM simulation are in good agreement with the laboratory results, verifying the reliability of the DEM model. Further analysis of the simulated results reveals that: (1) The capillarity and adsorption both play a role in the desiccation cracking process of clay. With the decrease of water content, the effects of adsorption gradually exceed those of capillarity. At the dominant stage of adsorption, the average displacement of simulated soil particles along the crack distribution direction accounts for 73% of its final value. (2) The total contact number between soil particles initially decreases and then increases as water content decreases. (3) The physical contact force between soil particles will increase rapidly at the dominant stage of capillarity/adsorption, resulting in stress concentration, and the contact between soil particles will be centralized and cracks will be formed. The proposed DEM contact model is of significant physical implications and can offer valuable insights into the underlying mechanisms of desiccation cracking in clay at the particle level.
-
Keywords:
- clay /
- unsaturated soil /
- desiccation cracking /
- discrete element method /
- capillarity /
- adsorption
-
0. 引言
基坑的开挖会导致周围土层发生位移，极易引起临近地铁隧道的变形破坏。实际工程中常常通过在基坑和隧道之间设置隔离桩来阻隔部分基坑开挖产生的影响[1-2]，从而对既有隧道形成保护作用。因此，对隔离桩的保护效果进行研究与评估，从而对隔离桩的设计进行优化，具有较大的现实意义和研究价值。
在针对隔离桩保护效果展开的国内外研究中，主要的研究方法包括数值模拟[2]、实测数据分析[2]和室内模型试验[3]等，其中以数值模拟为主，实测数据作为验证的研究较多，理论解[1]较少。其中，郑刚等[2]通过有限元模拟和实测数据验证的方法展开，提出隔离桩对隧道水平位移的影响具体表现为“阻隔作用”和“牵引作用”。在直接相关的理论解中，Wei等[1]建立了基坑-隔离桩-隧道三维力学模型，通过基坑开挖后隔离桩外隧道水平位移的三阶段计算方法得到隧道位移，并与实测对比验证了该模型的可靠性。相类似的单桩理论解研究中，张治国等[5]则是在二维坐标系上提出了基于地层损失的二阶段分析法，并给出了单桩水平位移解析解。综上所述，目前直接相关的理论解研究较少，大多数相似研究多是建立在二维坐标系上，但在实际工程中，隔离桩布置时大多是多根隔离桩间隔布置，作用在单桩上的基坑卸荷应力在平面坐标体系中难以满足研究要求。
实际工程中，基坑变形会呈现出明显的空间效应，二维坐标系上的研究也很难考虑到隔离桩的纵向排布情况。Wei等[1]虽给出三维坐标系上的理论解，但未进行影响因素上的后续研究。本文建立基坑-隔离桩-隧道计算模型，在三维坐标系上考虑基坑侧壁、隔离桩变形区域上的实际地层损失，并跟进研究了隔离桩长度、间距、位置分布三类影响因素。
目前已有的计算模型大多采用应力释放法计算基坑卸荷应力，该方法计算过程中需要较多的土体参数，如土体弹性模量、泊松比、重度、黏聚力、内摩擦角、土体主动极限位移等。此外，在预测基坑侧壁变形曲线计算上，目前的应力释放法仅适用于一定范围内的围护结构变形量，无法考虑大变形或破坏情况。本文根据Zhang等[7]提出的基坑侧壁变形曲线也可结合现场监测数据，由围护结构变形曲线或侧壁破坏土体损失量计算得到基坑旁侧隔离桩桩身附加应力，此过程仅需土体弹性模量以及泊松比，所需土体参数少，且可考虑大变形或破坏情况。
在基坑-隔离桩-隧道三维力学模型中，考虑地层损失计算的结果与地层损失区域形状和区域与计算点的相对位置有关，故在不同研究隔离桩挠曲变形区需要更精确的拟合区域以表示地层的损失，隔离桩地基模型的不同也可能导致周边土体或桩基的变形差异，需要更符合实际工程的地基模型。本文根据隔离桩长度对Kerr地基模型中的地基系数进行调整，计算出的隔离桩挠曲变形区域更加精确，故得到的地层损失区域的形状、位置更加精确，隧道水平位移更加符合工程实际。
1. 建立力学模型
图 1，2为本文计算模型图。由图可知，基坑位于既有隧道左侧，隔离桩位于两者之间，以基坑中心为原点建立三维坐标系。基坑为长方形，平行于隧道一侧的侧壁长为L，另一侧壁长为Lb，基坑开挖深度为d，坑底以下围护结构埋深为d0，围护结构总长为H（H=d+d0）。多根隔离桩组成的围护结构坐标位于xg处，总宽度为Lg，每一根隔离桩的桩径为Dg，长度为Hg，相邻间隔距离为lg。既有隧道直径为Ds，中心轴线埋深为Hs，位于xs处。
本部分分为3个步骤：①首先忽略隧道和隔离桩的存在，根据相应基坑变形曲线，考虑地层损失计算基坑开挖在隔离桩处的附加应力；②基于Kerr地基模型，根据桩身参数对地基系数进行调整，计算隔离桩挠曲变形；③基于地层补偿法，根据地层损失求得隔离桩位移引起的既有隧道位移值。
隧道水平位移小于20 mm时可不用设置隔离桩。隔离桩施工位置与隧道结构外边线水平投影间净距[21]最小为3 m，本模型中隔离桩施工均应位于隧道安全范围内。
2. 计算公式
2.1 基坑侧壁位移计算公式
如图 3所示，Zhang等[7]以分段余弦函数拟合围护结构变形增量，提出“内凸型”围护结构变形曲线
δi(η,Hei)=δmaxi2[1−cos(πηHei)](0≤η≤Hei)，δi(η,Hei)=δmaxi2{1−cos{π[η+(H−2Hei)]H−Hei}}(Hei≤η≤H)。} (1) 式中：δmaxi为第i层开挖引起的围护结构最大变形值；δi(η,Hei)为第i层开挖引起η深度处基坑侧壁围护结构的变形增量；Hei为第i层开挖后的开挖面深度。
第i层开挖引起的最大变形增量为
δmaxi=vmaxHe⋅Hei(i=1)，δmaxi=vmaxHe⋅Hei−i−1∑j=1δj(Hei)(i≥2)。} (2) 式中：vmax为围护结构积累最大变形值，实际工程中通常以围护结构的累计最大变形作为围护结构变形的控制要求，每一层开挖累积变形都应满足控制值。
考虑基坑变形的空间效应，由于围护结构系统刚度和土体应力状态不同，靠近基坑边角区域基坑变形较小，较长围护结构中心区域变形较大。刘念武等[8]采用围护结构的PSR变化趋势[9]有边角附近处PSR为0.72；当PSR接近于1.00时，对应的λ/Hei值为4.00。
PSR(λ,Hei)=1.671−e−0.1λ/Hei(0≤λ<4Hei)，PSR(λ,Hei)=1(λ≥4Hei) 。} (3) 式中：λ为计算位移点与基坑宽度边缘的距离。
整合式（1）~（3），即可获得考虑空间效应后基坑开挖时侧壁任一位置处的围护结构变形量计算公式：
v(λ,η)=PSR(λ,Hei)⋅δ(η,Hei)。 (4) 2.2 桩身附加应力计算公式
Sagasteta[15]在1987年提出地层位移变形计算模型，得到弹性半空间中任意点地层损失引起的周围土体位移场。该方法对土体位移作了一定假设：①土体不可压缩，地层损失是导致土体变形的原因；②忽略土体固结、流变以及孔隙水压力等影响。
图 4为地层位移变形计算模型，在图示的直角坐标系中，点F(Lb/2,η)处半径为a的间隙在隧道轴线处上的点P(x1,z1)处产生的位移沿x轴分量为
Sx1=−a22x1−Lb/2r′12。 (5) 式中：r′1=√(Lb/2−x1)2+(η−z1)2，表示点F与点P的距离。
式（5）是基于全空间条件下得到的，而实际中为半空间体。因此，需要转化为半空间体的求解问题。于是将F(Lb/2,η)镜像为F′(Lb/2,−η)，在该点发生大小相等的体积膨胀，F′(Lb/2,−η)在点P(x1,z1)处产生的沿x轴位移分量为
Sx2=a22x1−Lb/2r′22。 (6) 式中：r′2=√(Lb/2−x1)2+(η + z1)2，表示点F′与点P的距离。
在式（5）建立过程中，将在地面产生剪应变γ，
γ = −4a2η(x1−Lb/2)[(x1−Lb/2)2+η2]2。 (7) 则在地表产生的附加剪应力为
τ=Gγ=−4Ga2η(x1−Lb/2)[(x1−Lb/2)2+η2]2。 (8) 式中：G为土体剪切模量，G=Est6(1+νs)。
考虑地表附加剪应力影响，则剪应力在P点处产生的位移沿x方向分量为
Sx3=−a2(x1−Lb/2)r′22[1−2z1(z1+η)r′22]。 (9) 综上半无限空间内点F(Lb/2,η)处半径为a的空隙在P(x1,z1)处产生的土体位移沿x轴分量为
Sx(x1,z1)=Sx1+Sx2+Sx3。 (10) 为了计算基坑变形引起的既有隧道隔离桩附加应力值，本文引入地层位移变形计算模型，王涛等[6]利用该模型推导了隧道开挖产生的地层损失引起的土层中任意一点沿x轴、y轴、z轴3个方向上的附加应力计算公式。
在点（x0, y0, z0）处半径为a的球形空隙在水平x方向产生的应力分量σx - s1，σx - s2分别为
σx - s1=a3E3(1−2μ)(1r32−1r31)+a3E(1−μ)(1+μ)(1−2μ)(x−x0)2(1r51−1r52)+a3Eμ(1+μ)(1−2μ)[(y−y0)2(1r51−1r52)+(z−z0)21r51−(z+z0)21r52]， (11) σx-s2=a3πE1+μlimm→∞limn→∞∫y0+my0−m∫x0+nx0−nz0(x−u)(u−x0)[(u−x0)2+(t−y0)2+z20]5/2{1−2μ(r0+z)2[1r0−(y−t)2r03−2(y−t)2r02(r0+z)]−3(x−u)2r05}dudt+a3πE1+μlimm→∞limn→∞∫y0+my0−m∫x0+nx0−nz0(y−t)(t−y0)[(u−x0)2+(t−y0)2+z20]5/2⋅{1−2μ(r0+z)2[1r0−(x−u)2r03−2(x−u)2r02(r0+z)]−3(y−t)2r05}dudt 。 (12) 单位体积空隙引起的x方向的附加应力σ′x - s为
σ′x - s=σx - s1+σx - s243π a3。 (13) 基坑围护结构变形引起的水平附加应力σ1为
σ1=∭ (14) 式中：{r_0} = \sqrt {{{(x - u)}^2} + {{(y - t)}^2} + {z^2}} ；
{r_1} = \sqrt {{{(x - {x_0})}^2} + {{(y - {y_0})}^2} + {{(z - {z_0})}^2}} \text{；} {r_2} = \sqrt {{{(x - {x_0})}^2} + {{(y - {y_0})}^2} + {{(z + {z_0})}^2}} \text{；} E为土体弹性模量；V1为基坑围护结构变形形成的空隙区域，可由现场监测手段获得，也可通过式（4）计算区域各点坐标获得。
2.3 隔离桩挠曲变形计算公式
（1）隔离桩位移模型
参考冯国辉等[10]的方法求解隔离桩的挠曲变形。基于假定的开挖土体自由场位移模式，采用Kerr地基梁模型模拟隔离桩。图 5为Kerr地基桩土相互作用计算模型。
如图 6所示，假设基坑开挖引起的隔离桩处的自由土体位移为hsz（曲线AD），但由于隔离桩的存在，土体自由位移会受到约束，设桩身的最终水平位移为w2（曲线AC）。根据桩-土变形协调，桩身水平位移即为该部分土体的最终侧向位移。土体的约束位移，即桩-土的相对位移为\Delta =hsz-w2（阴影区域ACD），在附加荷载q的作用下，桩基水平位移w2的平衡微分方程为
\frac{{EIG}}{{{D_{\text{g}}}c}}\frac{{{{\text{d}}^6}{w_2}}}{{{\text{d}}{z^6}}} - \frac{{EI(c + k)}}{{{D_{\text{g}}}c}}\frac{{{{\text{d}}^4}{w_2}}}{{{\text{d}}{z^4}}} + G\frac{{{{\text{d}}^2}{w_2}}}{{{\text{d}}{z^2}}} - k{w_2} = - p。 (15) 式中：p为作用在桩基上的附加应力；EI为抗弯刚度；k，c分别为上层弹簧和下层弹簧的地基反力，k满足地基系数取值，参照简化弹性空间法[11]取c=3k；t为剪切层厚度。
\begin{array}{l} A{({w_2})_i} + B{({w_2})_i} - 3 + C{({w_2})_{i - 1}} + D{({w_2})_i} +\\ \qquad C{({w_2})_{i + 1}} + B{({w_2})_{i + 3}} + A{({w_2})_{i + 3}} = - {p_i}。 \end{array} (16) 式中：m=0，1，2，···，n，n+1。
根据差分法知
\frac{{{\text{d}}{w_2}}}{{{\text{d}}z}} = \frac{{{{({w_2})}_{i + 1}} + {{({w_2})}_{i - 1}}}}{{2l}}\text{，} (17) \frac{{{{\text{d}}^2}{w_2}}}{{{\text{d}}{z^2}}} = \frac{{{{({w_2})}_{i + 1}} - 2{{({w_2})}_i} + {{({w_2})}_{i - 1}}}}{{2l}}\text{，} (18) \frac{{{{\text{d}}^3}{w_2}}}{{{\text{d}}{z^3}}} = \frac{{{{({w_2})}_{i + 2}} - 2{{({w_2})}_{i + 1}} + 2{{({w_2})}_{i - 1}} - {{({w_2})}_{i - 2}}}}{{2{l^3}}}\text{，} (19) \frac{{{{\text{d}}^4}{w_2}}}{{{\text{d}}{z^4}}} = \frac{{{{({w_2})}_{i + 2}} - 4{{({w_2})}_{i + 1}} + 6{{({w_2})}_i} - 4{{({w_2})}_{i - 1}} + {{({w_2})}_{i - 2}}}}{{{l^4}}}\text{，} (20) \begin{array}{c} \frac{{{{\text{d}}^5}{w_2}}}{{{\text{d}}{z^5}}} = \left[ {{{({w_2})}_{i + 3}} - 5{{({w_2})}_{i + 2}} + 7{{({w_2})}_{i + 1}} - 7{{({w_2})}_{i - 1}}} \right. +\\ {{\left. {5{{({w_2})}_{i - 2}} - {{({w_2})}_{i - 3}}} \right]} \mathord{\left/ {\vphantom {{\left. {5{{({w_2})}_{i - 2}} - {{({w_2})}_{i - 3}}} \right]} {2{l^5}}}} \right. } {2{l^5}}}\text{，} \end{array} (21) \begin{array}{l}\frac{{{{\text{d}}^6}{w_2}}}{{{\text{d}}{z^6}}} = \left[ {{{({w_2})}_{i + 3}} - 6{{({w_2})}_{i + 2}} + 15{{({w_2})}_{i + 1}} - 20{{({w_2})}_i}} \right. +\\ \qquad \qquad \qquad \qquad {{\left. {15{{({w_2})}_{i - 1}} - 6{{({w_2})}_{i - 2}} + {{({w_2})}_{i - 3}}} \right]} \mathord{\left/ {\vphantom {{\left. {15{{({w_2})}_{i - 1}} - 6{{({w_2})}_{i - 2}} + {{({w_2})}_{i - 3}}} \right]} {{l^6}}}} \right. } {{l^6}}}\text{，} \end{array} (22) \left( \begin{gathered} {A_i} \hfill \\ {B_i} \hfill \\ {C_i} \hfill \\ {D_i} \hfill \\ \end{gathered} \right) = \left( {\begin{array}{*{20}{c}} 1&0&0&0 \\ { - 6}&{ - 1}&0&0 \\ {15}&4&1&0 \\ { - 20}&{ - 6}&{ - 2}&1 \end{array}} \right) \cdot \left( {\begin{array}{*{20}{c}} {EIG/{D_{\text{g}}}c{l^6}} \\ {EIW/{D_{\text{g}}}c{l^4}} \\ {G/{l^2}} \\ k \end{array}} \right)。 (23) 式中：W=c+k；l=Hg/n，Hg为隔离桩桩长。
（2）桩身位移矩阵表示
结合桩基边界条件，得到单桩位移方程为
\left\{ {{w_2}} \right\} = {[K]^{ - 1}} \cdot \left\{ { - {p_i}} \right\}。 (24) 式中：
\begin{array}{l}\left\{{w}_{2}\right\}=\left\{{\left\{{w}_{2}\right\}}_{0},{\left\{{w}_{2}\right\}}_{1},\cdots \cdots ,{\left\{{w}_{2}\right\}}_{n-1},{\left\{{w}_{2}\right\}}_{n}\right\}\\ \left\{-{p}_{i}\right\}=\left\{-{p}_{0},-{p}_{1},\cdots \cdots ,-{p}_{n-1},-{p}_{n}\right\}\text{ }。\end{array} \text{，} 假设桩端两处均为摩擦桩时，桩端剪力和弯矩为零，剪切层弯矩也为零，有
\left.\begin{array}{c}{M}_{S0}={M}_{Sn}=-EI{\frac{{\text{d}}^{2}{w}_{2}}{\text{d}{z}^{2}}|}_{i=0/n}=0\text{，}\\ {M}_{0}={M}_{n}=-EI{\frac{{\text{d}}^{2}w}{\text{d}{z}^{2}}|}_{i=0/n}=0\text{ }\text{，}\\ {Q}_{0}={Q}_{n}=-EI{\frac{{\text{d}}^{3}w}{\text{d}{z}^{3}}|}_{i=0/n}=0\text{ }。\end{array}\right\} (25) 式中：MS0和MSn分别为隔离桩上下两端剪切层的弯矩；M0和Mn为隔离桩两端弯矩；Q0和Qn为隔离桩上下两端剪力，那么[K]为
[K]=\left[\begin{array}{cccccccccc}{D}_{i}+2{C}_{i}+4{D}_{i}+8{A}_{i}& -4{B}_{i}-10{A}_{i}& 2{B}_{i}+2{A}_{i}& 2{A}_{i}& & & & & & \\ {C}_{i}+2{B}_{i}+2{A}_{i}& {D}_{i}-{B}_{i}& {C}_{i}-{A}_{i}& {B}_{i}& {A}_{i}& & & & & \\ {B}_{i}+2{A}_{i}& {C}_{i}-{A}_{i}& {D}_{i}& {C}_{i}& {B}_{i}& {A}_{i}& & & & \\ {A}_{i}& {B}_{i}& {C}_{i}& {D}_{i}& {C}_{i}& {B}_{i}& {A}_{i}& & & \\ & {\bf{0}}& {\bf{0}}& {\bf{0}}& {\bf{0}}& {\bf{0}}& {\bf{0}}& {\bf{0}}& & \\ & & & {A}_{i}& {B}_{i}& {C}_{i}& {D}_{i}& {C}_{i}& {B}_{i}& {A}_{i}\\ & & & & {A}_{i}& {B}_{i}& {C}_{i}& {D}_{i}& {C}_{i}-{A}_{i}& {B}_{i}+2{A}_{i}\\ & & & & & {A}_{i}& {B}_{i}& {C}_{i}-{A}_{i}& {D}_{i}-{B}_{i}& {C}_{i}+2{B}_{i}+2{A}_{i}\\ & & & & & & 2{A}_{i}& 2{B}_{i}+2{A}_{i}& -4{B}_{i}-10{A}_{i}& {D}_{i}+2{C}_{i}+4{B}_{i}+8{A}_{i}\end{array}\right]。 (26) （3）地基系数k调整
鲁子爱[13]考虑了地基系数k值沿深度分布发生变化，引入参数 {S_{\text{a}}} 考虑土质对地基系数的影响，
{S_{\text{a}}} = {a_z}({c_0} + \gamma \overline x {\text{tan}}\varphi ) 。 (27) 式中：c0为地基土黏聚力， \gamma 为地基土重度， \varphi 为地基土内摩擦角，\overline x 为桩身中心点所在土层平均深度， {a_z} 为计算参数，黏土为1，砂土可取3~5。
给出不同深度下地基系数计算公式：
k=\left\{\begin{array}{c}{c}_{1}{\left(\frac{{H}_{p}}{{D}_{\text{g}}^{2}}\right)}^{-{d}_{1}}+{c}_{2}{\left(\frac{{H}_{\text{p}}}{{D}_{\text{g}}^{2}}\right)}^{-{d}_{2}}\cdot {x}^{{c}_{3}\frac{{H}_{\text{p}}}{{D}_{\text{g}}{}^{2}}+{d}_{3}}\text{，}\text{ }表层土\text{ }\text{；}\\ {c}_{4}{\left(\frac{{H}_{\text{p}}}{{D}_{\text{g}}^{2}}\right)}^{-{d}_{4}}+\left({d}_{5}-{c}_{5}\frac{{H}_{\text{p}}}{{D}_{\text{g}}^{2}}\right)\cdot {x}^{{c}_{6}\frac{{H}_{\text{p}}}{{D}_{\text{g}}{}^{2}}+{d}_{8}}\text{，}深层土。\end{array}\right. (28) 式中： {H_{\text{p}}} 为桩身横向荷载； {c}_{1}=1.3{S}_{\text{a}}{}^{4.88}{D}_{\text{g}}{}^{-0.85}； {c}_{2}=248914{S}_{\text{a}}{}^{-0.58}{D}_{\text{g}}{}^{-8}； {c_3} = {10^{ - 8}}(193 + 3.312{S_{\text{a}}}){D_{\text{g}}}^{ - 0.2}； {c}_{4}=(331216+43463{S}_{\text{a}}){D}_{\text{g}}{}^{-0.15}； {c}_{5}=313.13{S}_{\text{a}}{}^{-2}{D}_{\text{g}}{}^{-2.25}； {c}_{6}=0.0227{S}_{\text{a}}{}^{-1.158}{D}_{\text{g}}{}^{-0.25}；{d}_{1}=(0.8485+0.026{S}_{\text{a}}){D}_{\text{g}}{}^{0.24}； {d}_{2}=(0.98-0.00787{S}_{\text{a}}){D}_{\text{g}}{}^{-0.15}；{d}_{3}=(1.782-0.0043{S}_{\text{a}}){D}_{\text{g}}{}^{0.28}； {d}_{4}=(0.9827-0.00243{S}_{\text{a}}){D}_{\text{g}}{}^{-0.05}；{d}_{5}=12820{S}_{\text{a}}{}^{-1.268}{D}_{\text{g}}{}^{-8}； {d}_{6}=(2.318+0.00886{S}_{\text{a}}){D}_{\text{g}}{}^{0.85}。
2.4 隔离桩侧移引起既有隧道位移计算公式
隔离桩的侧向挠曲变形将会影响既有隧道一侧的土体，引起土体位移，从而造成既有隧道的变形。成怡冲等[4]在研究中根据地层补偿原理[14]计算了隔离桩侧移引起的桩外土体水平位移，证明了采用该方法的合理性。如图 7所示，由于隔离桩的向左偏移，将会在ABC区域产生一个空隙，该空隙在实际工程中会被隔离桩右侧的周围土体所填充，因此，隔离桩右侧的土体会发生位移，从而引起隧道的位移。
![]()
图 7 基坑-隔离桩-隧道三者位移变形示意图Figure 7. Schematic diagram of displacement deformation of foundation pit-isolation pile-tunnel（1）隔离桩侧移后产生空隙引起的附加应力计算
结合式（11）~（13），计算得出隔离桩侧向偏移引起的水平附加应力{\sigma _2}为
{\sigma _2} = \iiint\limits_{{V_2}} {{{\sigma '}_{x{\text{ - }}s}}(x,y,z)}{\text{d}}x{\text{d}}y{\text{d}}z。 (29) 式中：V2为隔离桩挠曲变形形成的空隙区域。
经过上述分析，为了计算既有隧道轴线处的附加应力{\sigma _{x{\text{ - }}s}}，首先需要明确隔离桩挠曲变形后形成的空隙区域，并对空隙区域的边界进行函数表达。
如图 8（a）所示，隔离桩能够起到阻隔作用，减小所在位置的土体位移。当各个隔离桩相互紧靠，即间距为零时，隔离桩就能起到最好的阻隔作用。但在实际工程中，相邻隔离桩之间往往存在空隙，周围土体在空隙处的移动并未受到隔离桩的阻隔，故桩间土体会偏移更多，形成“小鼓包”，如图 8（b）阴影处所示。这一现象对空隙区域的函数表达极为不利，为了有效解决上述问题，本文首先以图 8（a）所示的紧密隔离桩模型（间距为零）进行计算。
![]()
图 8 紧密隔离桩与非紧密隔离桩的对比Figure 8. Comparison between tight piles and non-tight isolation piles如图 9（a）所示，各面围成的部分即为隔离桩挠曲变形后形成的空隙区域，将位移后的隔离桩桩心依次连接，即形成图 9（b）中的ABC区域，曲线ABC即为同一水平面（同一埋深）下隔离桩对应的位移函数xy。则式（29）可以改为
{\sigma _2}{\text{ = }}\int_{ - 0.5{L_g}}^{0.5{L_{\text{g}}}} {\int_{{x_g} - {x_y}}^{{x_g}} {\int_{{z_{xy}}}^0 {{{\sigma '}_{x{\text{ - s}}}}(x,y,z)} } } {\text{d}}z{\text{d}}x{\text{d}}y 。 (30) ![]()
图 9 隔离桩挠曲变形形成的空隙区域示意图Figure 9. Schematic diagram of void area formed by deflection deformation of isolation piles式中：zxy为空隙部分底部的弧面所对应的函数表达式。
（2）隧道位移计算
本文引用剪切错台和刚体转动协同变形模型[16]对既有隧道水平位移进行计算。该模型已在多篇论文[1, 20]中进行了应用并得到了验证。本文不再对公式推导过程进行赘述，具体方法可参照原文。
既有隧道水平位移W的最终计算公式为
W{\text{ = }}\left\{ {{T_n}(y)} \right\}A_a^{\text{T}}。 (31) 式中： {T_n}(y) = \left\{ {1{\text{,}}\cos \frac{{{\text{π }}y}}{{N{D_{\text{t}}}}}{\text{,}}\cos \frac{{2{\text{π }}y}}{{N{D_{\text{t}}}}}{\text{,}} \cdots ,\cos \frac{{{n_0}{\text{π }}y}}{{N{D_{\text{t}}}}}} \right\} ； {A_a} = {({a_1}{\text{, }}{a_2}{\text{, }}{a_3}{\text{,}} \cdots ,{a_{{n_0}}})^{\text{T}}} ；ai为矩阵元素；N为既有隧道单侧受影响管片环数； {D_{\text{t}}} 为管片环环宽；n0为计算阶数，隧道位移计算时n0取10即可满足精度要求。
根据式（31），无隔离桩时既有隧道的水平位移量W0，设置无间隙隔离桩时既有隧道的水平位移量W1均可计算得到。
紧密隔离桩阻隔效果w1为设置隔离桩前后既有隧道水平变形量的减小值，即满足
{w_1} = {W_0} - {W_1}。 (32) 而在有间距的隔离桩布置方案下，隔离桩的整体阻隔效果会有所下降。假设隔离桩桩径为Dg，间隔为lg，每一根隔离桩等距间隔布置，隔离桩墙的总体阻隔效果w0-1为紧密隔离桩阻隔效果w1的\frac{{{D_{\text{g}}}}}{{{D_{\text{g}}} + {l_{\text{g}}}}}倍[1]。既满足
{w_{0{\text{ - }}1}} = \frac{{{D_{\text{g}}}}}{{{D_{\text{g}}} + {l_{\text{g}}}}}{w_1}。 (33) 故设置有间隙隔离桩时的既有隧道水平位移量W0-1为
{W_{0{\text{ - }}1}}{\text{ = }}{W_0} - {w_{0{\text{ - }}1}} = {W_0} - \frac{{{D_{\text{g}}}}}{{{D_{\text{g}}} + {l_{\text{g}}}}}({W_0} - {W_1})。 (34) 3. 工程案例及分析验证
3.1 案例一
宁波市轨道交通1号线世纪大道站—海晏北路站区间隧道周围需进行大型基坑开挖，工程中通过设置隔离桩来对既有隧道进行保护，相关工程计算参数[17]如下：L=70 m，d=10 m；隔离桩布置平面位于xg=5 m，Lg=70 m，Dg=1.2 m，Hg=50 m，lg=1.6 m，隔离桩等效抗弯刚度EpIp=2.85×106 kN·m2；既有隧道直径Ds=6.2 m，Hs=17.1 m，间距xs=40.1 m，隧道采用C50混凝土管片，厚度ts=0.35 m，环宽Dt=1.2 m，管片环之间由16根M30纵向螺栓连接，根据计算得到环间剪切强度ks= 7.45×105 kN/m，环间抗拉刚度kt=1.94×106 kN/m，隧道等效抗弯刚度EtIt=1.1×108 kN·m2，计算精度N=100，比例系数j=0.3。μ=0.35， \gamma =18.09 kN/m3，静止土压力系数{K_0}=0.753，土的黏聚力{c_0}=12.9 kPa，φ= 11.37°，地基土的压缩模量Es=3.99 MPa，E=9.975 MPa。
如图 10所示为本文计算方法所得隧道水平位移与文献[17]所得实测值及数值模拟结果的对比。在图中，工况1为无隔离桩工况，工况2为设置紧密隔离桩工况，工况3为工程案例1中实际的隔离桩布置工况。按照本文方法首先计算工况1和工况2下的隧道水平位移曲线，再通过式（33）即可求得实际工况下既有隧道的水平位移曲线。由图可知，工程案例1对应的隧道水平位移曲线与实测值及仿真值的变化规律基本一致，大致呈现中心大，两端小的变化趋势。本文计算方法所得隧道最大水平位移值发生在y=0 m处，其值大小为1.38 mm，对应实测值为1.5 mm，对应仿真值为1.42 mm，差值均较小，满足准确性要求。
![]()
图 10 工程案例1隧道水平位移对比Figure 10. Comparison of horizontal displacements of tunnels with different spacings in project case 13.2 工程案例二
武汉市汉阳区地铁3号线旁侧深基坑，基坑与隧道最小水平距离为20 m。L=62 m，Lb=95 m，d=13 m；xg=60 m，Lg=77 m；Dg=1 m，Hg=16.5 m，lg依次取为0.5，1，1.5，2 m；EpIp=1.47×106 kN·m2；Ds=6 m，Hs=13 m，xs=67.5 m，ks=7.45×105 kN/m，ts=0.35 m，Dt=1.5 m，kt=1.94×106 kN/m，EtIt=1.1×108 kN·m2，N=60，j=0.5，μ=0.31， \gamma =19.42 kN/m3，K0=0.795，c0=20.26 kPa，φ=8.89°，Es=4.28 MPa，E=10.7 MPa。
图 11为本文计算方法所得隧道水平位移曲线与文献[18]数值模拟结果的对比。如图 11所示，隔离桩间距lg=1，1.5，2 m时，本文计算方法所得曲线与数值模拟所得曲线吻合度较好，近似呈正态分布。本文方法所得隧道最大水平位移发生在中心处，最大值约为4.00 mm，而数值模拟所得最大值为4.25 mm，差值仅为0.25 mm，满足准确性要求。如图 11（b）所示，隔离桩间距lg=1 m，两条曲线总体变化趋势相同，数值模拟值略小于本文计算所得，中心最大值分别为4.47 mm和4.79 mm，差值较小，满足准确性要求。而图 11（c）为隔离桩间距lg=1.5 m，图 11（d）为隔离桩间距lg=2 m，其变化规律相同，差值较小，均满足计算精度要求。
![]()
图 11 工程案例2隧道水平位移对比Figure 11. Comparison of horizontal displacements of tunnels in project case 23.3 工程案例三
浙江杭州地铁1号线一侧需要开挖一个大型基坑。由于最小距离仅为24 m。在对数据进行整理和简化后，有关工程计算参数如下[19]：L=186 m，d=5.7 m，xg=20 m，Lg=120 m，Dg=1.2 m，Hg=18 m，lg=2.4 m，EpIp=3.52×105 kN·m2；ks=7.45×105 kN/m，Ds=6.2 m，Hs=14 m，kt=1.94×106 kN/m，xs=30 m，ts=0.35 m，Dt=1.2 m，EtIt=1.1×108 kN·m2，计算精度N=100，比例系数j=0.5。土层及厚度从上至下依次为等效粉砂层15.9 m，等效黏土层27.4 m，μ=0.31， \gamma =18.9 kN/m3。
如图 12所示，工况1、工况2分别为隔离桩间距2.4，1.2 m时隧道水平位移计算值。工况3为隔离桩间距1.6 m时隧道水平位移计算值，其y=0处计算值为1.78 mm与原文的数值模拟值为1.66 mm对比仅相差0.12 mm，满足计算精度要求。
![]()
图 12 工程案例3隧道水平位移对比Figure 12. Comparison of horizontal displacements of tunnels in project case 34. 隔离桩保护效果影响因素分析
4.1 隔离桩桩长Hg
以案例一为基本工况，在其他参数不变情况下，仅改变隔离桩桩长，研究不同长度下隔离桩对隧道变形的影响。考虑实际工程中隔离桩长度、基坑开挖深度、隧道埋深分别为50，10，17.1 m，故按照50，40，30，20 m数值选取隔离桩桩长Hg。
图 13为隔离桩长度分别在50，40，30，20 m时，邻近隧道水平位移计算值分布曲线。如图 13所示，桩长为20至40 m范围内隧道水平位移随着隔离桩桩长的增加而减少，但40至50 m间隧道位移不降反增。
![]()
图 13 不同桩长隔离桩影响下的水平位移Figure 13. Horizontal displacements of isolation piles under influences of different pile lengths原因分析：①在相同的附加应力下，桩挠曲变形区域起始会随着桩长度的增加而减小，故桩与隧道间的土体损失量会减小，进而减小了作用在隧道上的附加荷载，减小隧道变形；②但是当桩达到一定长度时，出现临界点，桩的挠曲变形区域开始增大，具有一定“牵引作用”[2]，导致桩后土体损失量增加，从而增加了隧道水平位移。根据理论计算规律，原案例工程建议可选择设置长度40 m的隔离桩。
4.2 隔离桩桩间距lg
以案例二为基本工况，在其他参数不变情况下，仅改变隔离桩桩间距，研究不同间距下隔离桩对隧道变形的影响。图 14为隔离桩桩间距分别为紧密隔离桩间距0.5，1，1.5，2 m以及无隔离桩时隧道水平位移。由图 14可知，根据本文计算方法，隧道水平位移分布曲线均呈现中心大，两端小的趋势。隧道整体水平位移随桩间距的减小而减小；隧道最大处水平位移与桩间距呈非线性增大趋势。
![]()
图 14 桩间距影响下隧道的水平位移Figure 14. Horizontal displacements of tunnels under influences of pile spacing4.3 隔离桩桩位置xg
以案例二为基本工况，记基坑与隧道水平距离为Ls=20 m，隔离桩与基坑水平距离计算值Lz=L/2-xg。
图 15为Lz/Ls分别为0.375，0.5，0.625，0.75时隧道的水平位移计算值分布曲线。由图 15可知，在不同位置隔离桩保护下，隧道水平位移分布曲线均呈现正态分布趋势；当隔离桩靠近基坑时，隧道水平位移分布曲线逐渐变缓。y=0对应隧道最大水平位移处，随隔离桩位置逐渐靠近基坑一侧逐渐减小；0.375处与0.5处计算曲线近似重叠，0.375处曲线两侧计算值大于0.5处，但考虑到隔离桩施工会对基坑维护结构产生扰动，故建议原工程将隔离桩设置在基坑与隧道水平间隔中心处，即xg=57.5 m。
![]()
图 15 设置不同位置桩时隧道水平位移分布曲线Figure 15. Distribution curves of horizontal displacement of tunnels under different positions of piles5. 结论
（1）本文计算结果与实测值及数值模拟结果较吻合，设置隔离桩前后既有隧道的水平位移曲线均为正态分布，且隔离桩的设置能够有效减小既有隧道的水平位移值，起到很好的阻隔保护作用。
（2）本文在三维力学模型上考虑地层损失推导出基坑侧壁变形在坑外任意一点处产生的附加应力计算公式，比应力释放法计算参数少且更符合工程实际。
（3）在计算隔离桩挠曲变形区域过程中，本文方法综合考虑地基系数随桩参数和地层深度的变化规律，对地基系数作出调整，发现不同长度隔离桩下变形区域的差异规律，从理论上解释了隔离桩对桩后土体的“牵引作用”。
（4）隔离桩桩间距越小，其对基坑旁侧既有隧道的保护作用就越好。隔离桩长度低于临界值时，保护作用随着桩长度的增加而减小，高于临界值时规律相反。隔离桩设置在靠近基坑一侧时，保护效果较好，但考虑到隔离桩施工也会对基坑围护结构产生影响，案例二中隔离桩设置在xg=57.5 m时效果最好。
本文在分析过程中做了较多简化，如没有考虑隔离桩施工对基坑围护结构产生的风险以及基坑施工前隧道就已经产生位移的情况，理论结果会有一定误差，可在本文基础上作进一步研究。
-
![]()
![]()
![]()
图 8 离散元模拟得到的不同含水率情况下土颗粒的轴向位移
Figure 8. Axial displacements of soil particles under different water contents obtained by DEM simulation
![]()
图 9 室内试验[29]与离散元模拟的土样应变-含水率关系
Figure 9. Strain-water content relationship of soil samples in laboratory tests and DEM simulation
![]()
图 10 土颗粒轴向位移绝对值的最大值、平均值的归一化百分比
Figure 10. Normalized percentages of maximum value and average value of absolute axial displacement of soil particles
![]()
图 11 土颗粒间总接触数量与含水率的关系
Figure 11. Relationship between total contact number of DEM soil samples and water content
![]()
图 12 干缩过程吸应力与含水率的关系
Figure 12. Relationship between suction stress and water content during desiccation
![]()
图 13 土颗粒间接触力与含水率的关系
Figure 13. Relationship between contact force of soil particles and water content
![]()
图 14 土颗粒间物理接触力力链与含水率的关系
Figure 14. Relationship between physical contact force chain of soil particles and water content
表 1 离散元模拟参数取值
Table 1 Values of parameters for DEM simulation
参数 变量 取值 初始颗粒半径/m {R_0} 1.0 \times {10^{-3}} 阻尼系数 \beta 0.7 颗粒法向刚度/(N·m-1) k_{\text{n}}^{{\text{ball-ball}}} 1.0 \times {10^6} 颗粒切向刚度/(N·m-1) k_{\text{s}}^{{\text{ball-ball}}} 8.0 \times {10^5} 颗粒-墙法向刚度/
(N·m-1)k_{\text{n}}^{{\text{ball-wall}}} 1.0 \times {10^7} 颗粒-墙切向刚度/
(N·m-1)k_{\text{s}}^{{\text{ball-wall}}} 8.0 \times {10^6} 颗粒摩擦系数 {\mu ^{{\text{ball-ball}}}} 0.20 颗粒-墙摩擦系数 {\mu ^{{\text{ball-wall}}}} 0.0 吸引力/N {F_{\text{a}}} d_{50}^2{\sigma _{\text{s}}} 吸引力作用范围/m {d_{\text{a}}} 1.0{R_{\text{w}}} 
下载: 导出CSV
-
[1] 唐朝生, 施斌, 崔玉军. 土体干缩裂隙的形成发育过程及机理[J]. 岩土工程学报, 2018, 40(8): 1415-1423. doi: 10.11779/CJGE201808006 TANG Chaosheng, SHI Bin, CUI Yujun. Behaviors and mechanisms of desiccation cracking of soils[J]. Chinese Journal of Geotechnical Engineering, 2018, 40(8): 1415-1423. (in Chinese) doi: 10.11779/CJGE201808006
[2] TANG C S, ZHU C, CHENG Q, et al. Desiccation cracking of soils: A review of investigation approaches, underlying mechanisms, and influencing factors[J]. Earth-Science Reviews, 2021, 216: 103586. doi: 10.1016/j.earscirev.2021.103586
[3] 殷宗泽, 袁俊平, 韦杰, 等. 论裂隙对膨胀土边坡稳定的影响[J]. 岩土工程学报, 2012, 34(12): 2155-2161. http://cge.nhri.cn/cn/article/id/14958 YIN Zongze, YUAN Junping, WEI Jie, et al. Influences of fissures on slope stability of expansive soil[J]. Chinese Journal of Geotechnical Engineering, 2012, 34(12): 2155-2161. (in Chinese) http://cge.nhri.cn/cn/article/id/14958
[4] ALBRECHT B A, BENSON C H. Effect of desiccation on compacted natural clays[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2001, 127(1): 67-75. doi: 10.1061/(ASCE)1090-0241(2001)127:1(67)
[5] RAYHANI M H T, YANFUL E K, FAKHER A. Physical modeling of desiccation cracking in plastic soils [J]. Engineering Geology, 2008, 97(1-2): 25-31. doi: 10.1016/j.enggeo.2007.11.003
[6] TANG C S, SHI B, LIU C, et al. Influencing factors of geometrical structure of surface shrinkage cracks in clayey soils [J]. Engineering Geology, 2008, 101(3/4): 204-217.
[7] TANG C S, CUI Y J, TANG A M, et al. Experiment evidence on the temperature dependence of desiccation cracking behavior of clayey soils[J]. Engineering Geology, 2010, 114(3/4): 261-266.
[8] TRABELSI H, JAMEI M, ZENZRI H, et al. Crack patterns in clayey soils: Experiments and modeling[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2012, 36(11): 1410-1433. doi: 10.1002/nag.1060
[9] 林朱元, 唐朝生, 曾浩, 等. 土体干缩开裂过程的边界效应试验与离散元模拟[J]. 岩土工程学报, 2020, 42(2): 372-380. doi: 10.11779/CJGE202002019 LIN Zhuyuan, TANG Chaosheng, ZENG Hao, et al. Laboratory characterization and discrete element modeling of desiccation cracking behavior of soils under different boundary conditions[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(2): 372-380. (in Chinese) doi: 10.11779/CJGE202002019
[10] GUO Y, HAN C J, YU X. Laboratory characterization and discrete element modeling of shrinkage and cracking in clay layer[J]. Canadian Geotechnical Journal, 2018, 55(5): 680-688. doi: 10.1139/cgj-2016-0674
[11] 沈珠江, 邓刚. 黏土干湿循环中裂缝演变过程的数值模拟[J]. 岩土力学, 2004, 25(增刊2): 1-6, 12. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX2004S2000.htm SHEN Zhujiang, DENG Gang. Numerical simulation of crack evolution in clay during drying and wetting cycle[J]. Rock and Soil Mechanics, 2004, 25(S2): 1-6, 12. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX2004S2000.htm
[12] WANG X N, YU P, YU J L, et al. Simulated crack and slip plane propagation in soil slopes with embedded discontinuities using XFEM[J]. International Journal of Geomechanics, 2018, 18(12): 04018170. doi: 10.1061/(ASCE)GM.1943-5622.0001290
[13] SANCHEZ M, MANZOLI O L, GUIMARAES L J N. Modeling 3-D desiccation soil crack networks using a mesh fragmentation technique[J]. Computers and Geotechnics, 2014, 62: 27-39. doi: 10.1016/j.compgeo.2014.06.009
[14] YAN C Z, WANG T, KE W H, et al. A 2D FDEM-based moisture diffusion-fracture coupling model for simulating soil desiccation cracking [J]. Acta Geotechnica, 2021, 16(8): 2609-2628. doi: 10.1007/s11440-021-01297-4
[15] PERON H, DELENNE J Y, LALOUI L, et al. Discrete element modelling of drying shrinkage and cracking of soils [J]. Computers and Geotechnics, 2009, 36(1-2): 61-69. doi: 10.1016/j.compgeo.2008.04.002
[16] EL YOUSSOUFI M S, DELENNE J Y, RADJAI F. Self-stresses and crack formation by particle swelling in cohesive granular media [J]. Physical Review E, 2005, 71(5): 051307. doi: 10.1103/PhysRevE.71.051307
[17] 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. doi: 10.1016/j.compgeo.2013.12.003
[18] 司马军, 蒋明镜, 周创兵. 黏性土干缩开裂过程离散元数值模拟[J]. 岩土工程学报, 2013, 35(增刊2): 286-291. http://cge.nhri.cn/cn/article/id/15396 SIMA Jun, JIANG Mingjing, ZHOU Chuangbing. Numerical simulation of desiccation cracking of clay soils by DEM[J]. Chinese Journal of Geotechnical Engineering, 2013, 35(S2): 286-291. (in Chinese) http://cge.nhri.cn/cn/article/id/15396
[19] LIN Z Y, WANG Y S, TANG C S, et al. Discrete element modelling of desiccation cracking in thin clay layer under different basal boundary conditions[J]. Computers and Geotechnics, 2021, 130: 103931. doi: 10.1016/j.compgeo.2020.103931
[20] LE T C, LIU C, TANG C S, et al. Numerical simulation of desiccation cracking in clayey soil using a multifield coupling discrete-element model[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2022, 148(2): 04021183. doi: 10.1061/(ASCE)GT.1943-5606.0002747
[21] MITCHELL J K, SOGA K. Fundamentals of Soil Behaviour[M]. Hoboken, New Jersey: John Wiley & Sons, Inc, 2005.
[22] 施斌, 唐朝生, 王宝军, 等. 黏性土在不同温度下龟裂的发展及其机理讨论[J]. 高校地质学报, 2009, 15(2): 192-198. doi: 10.3969/j.issn.1006-7493.2009.02.007 SHI Bin, TANG Chaosheng, WANG Baojun, et al. Development and mechanism of desiccation cracking of clayey soil under different temperatures[J]. Geological Journal of China Universities, 2009, 15(2): 192-198. (in Chinese) doi: 10.3969/j.issn.1006-7493.2009.02.007
[23] LU N, LIKOS W J. Suction stress characteristic curve for unsaturated soil[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2006, 132(2): 131-142. doi: 10.1061/(ASCE)1090-0241(2006)132:2(131)
[24] ZHANG C, LU N. Unified effective stress equation for soil [J]. Journal of Engineering Mechanics, 2020, 146(2): 04019135. doi: 10.1061/(ASCE)EM.1943-7889.0001718
[25] ZHANG C, LI J Z, ZHANG Y X, et al. Experimental and discrete element modeling study on suction stress characteristic curve and soil-water characteristic curve of unsaturated reticulated red clay [J]. Bulletin of Engineering Geology and the Environment, 2022, 81(9): 363. doi: 10.1007/s10064-022-02834-5
[26] GUO L, CHEN G, DING L, et al. Numerical simulation of full desiccation process of clayey soils using an extended DDA model with soil suction consideration [J]. Computers and Geotechnics, 2023, 153: 105107. doi: 10.1016/j.compgeo.2022.105107
[27] VAN GENUCHTEN M T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils [J]. Soil Science Society of America Journal, 1980, 44(5): 892-898. doi: 10.2136/sssaj1980.03615995004400050002x
[28] LU N, GODT J W, WU D T. A closed-form equation for effective stress in unsaturated soil [J]. Water Resources Research, 2010, 46(5): 567-573.
[29] PERON H, LALOUI L, HUECKEL T, et al. Experimental study of desiccation of soil [M]. Unsaturated Soils 2006. ASCE geotechnical special publication 147. 2006: 1073-1084.
[30] YAO M. Three-Dimensional Discrete Element Method Analysis of Cohesive Soil[D]. Baltimore: The Johns Hopkins University, 2002.
[31] 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)
[32] SOULIE F, CHERBLANC F, EL YOUSSOUFI M S, et al. Influence of liquid bridges on the mechanical behaviour of polydisperse granular materials[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2006, 30(3): 213-228. doi: 10.1002/nag.476
-
期刊类型引用(1)
1. 赵飞涛. 基于锚固界面力学特性的拉压型锚杆承载特性研究. 长沙理工大学学报(自然科学版). 2025(02): 99-109 . 
百度学术
其他类型引用(0)
-
其他相关附件
计量
- 文章访问数: 331
- HTML全文浏览量: 84
- PDF下载量: 82
- 被引次数: 1
