Calculation of loosening earth pressure in deep trapdoor tests on unsaturated sandy soils
-
摘要: 合理确定深埋活动门上方松动土压力对分析土与结构物之间的相互作用具有重要意义。传统的土拱模型和松动土压力计算方法主要针对干燥和完全饱和土体，较少涉及到非饱和土的研究。基于ABAQUS有限元软件建立不同饱和度下深埋活动门数值模型，根据数值试验结果，通过考虑非饱和土中基质吸力的影响，建立了非饱和砂土中深埋活动门上方松动土压力分析模型，随后结合大主应力圆弧拱理论和非饱和土抗剪强度理论，给出了松动区高度和松动土压力计算方法。通过与非饱和活动门试验结果和数值模拟结果对比，验证了模型的合理性。参数结果分析表明：松动土压力随着饱和度的增加呈现先减小后增加的趋势，在临界饱和度处，松动土压力值最小；随着活动门埋深的增加，归一化松动土压力呈现缓慢减小的趋势。Abstract: The loosening earth pressure acting on the deep trapdoor test is of great significance for the interaction analysis between soils and structures. The traditional model for soil arching and the solutions for loosening earth pressure mainly focus on dry and fully saturated soils. The unsaturated conditions commonly encountered in engineering practice are seldom considered. Based on the ABAQUS finite element software, a series of numerical deep trapdoor tests are conducted to study the model for soil failure by considering the effects of soil saturation. Then, an analytical model for the calculation of the loosening earth pressure is established according to the influences of soil saturation. By combining the major principal arc arch theory with the shear strength theory of unsaturated soils, the analytical solutions for the height of the loose zone and the loosening earth pressure are proposed. The rationality of the proposed model is verified by comparing the results with those of unsaturated trapdoor tests, existing analytical solutions and numerical simulations. The results of parameter analysis show that with the increase of the soil saturation, the loosening earth pressure exhibits a decreasing trend first and then an increasing trend at the critical saturation point, and the value of the loosening earth pressure reaches its minimum value. With the increase of the relative depth ratio, the loosening earth pressure shows a slowly decreasing trend.
-
0. 引言
土拱效应是土与结构物相互作用的重要现象，在工程实践中普遍存在。由于土体内部不均匀变形产生的剪切作用，使得土体内部应力发生转移和重分布，导致作用在结构物上的土压力大幅减小。活动门试验是通过以一定速率下降的活动门模拟实际土层中出现的沉降差，通过观测下降过程中土层的变形和活动门上方松动土压力的变化研究土拱效应的演化过程。合理考虑土拱效应并准确计算作用在结构物上的土压力对研究土与结构物相互作用具有重要的意义。
Terzaghi[1]通过活动门试验证实了土拱效应的存在，并在此基础上，给出了平面应变条件下的松动土压力解析解，被广泛应用于岩土工程实践中，如盾构开挖面稳定性分析、桩承式路堤桩间以及管道上方土压力计算等[2-6]。在Terzaghi模型基础上，许多学者通过理论分析、数值模拟和模型试验等方式对土拱效应和松动土压力展开研究。试验方面，Costa等[7]研究了主动和被动模式下深埋矩形活动门上方土体压力的变化情况。试验结果表明，随着活动门位移，土体中出现了一个拱形卸载区域和一个拱形加载区域。拱形卸载区域主要由活动门上方坍塌土体和周围的部分土体组成，而拱形加载区域主要由稳定回填区组成。Liang等[8]采用经典活动门试验，结合PIV技术研究了土拱效应随相对位移的演变。研究表明松动土压力的最小值和残余值分别对应三角形和垂直滑移面。Liu等[9]通过模型试验，研究了上软下硬复合浅层地铁隧道的渐进坍塌机理。理论研究方面，Wan等[10]考虑土体相对位移与拱形变形规律的相关性，提出了一种新的深埋隧道下松动土压力计算模型。Lin等[11-12]基于极限平衡法，分别给出了矩形隧道和圆形隧道上覆松动土压力计算方法。宫全美等[13]结合颗粒椭球体理论推导了椭圆形滑移面松动土压力计算公式。数值模拟方面，Lai等[14]和Chen等[15]基于有限元数值模拟分析，分别对黏性土在浅埋和深埋情况下的破坏模式进行了研究，提出了考虑破坏机制和荷载传递机制的解析解。Liang等[16]采用PFC2D颗粒流方法研究了受约束颗粒物质作用下的平均松动土压力和活动门松动区的应力分布。
以上关于活动门试验松动土压力计算研究大多针对干燥和饱和土体，较少涉及在自然界中普遍存在的非饱和状况。崔蓬勃等[17]通过开展低含水率砂土活动门试验发现相对于干砂，非饱和工况土压力极值明显减小，低含水率下的松动区破坏模式呈现椭圆状，且松动区边缘出现明显裂缝，上部形成自然平衡拱。Song等[18]通过活动门试验研究了水位变化对土拱效应的影响，地下水位对土拱效应的极限影响范围是活板门宽度的两倍。蔺港等[19]在Terzaghi模型的基础上，引入非饱和土抗剪强度准则，推导了考虑基质吸力影响的非饱和土松动土压力计算公式。崔蓬勃等[20]基于大主应力迹线假设，建立了考虑土拱效应发挥过程的非饱和砂土松动土压力计算方法。
可以看出，目前对于非饱和土中松动土压力已有一定的研究，然而由于非饱和问题的复杂性，该问题仍需进一步深入研究。非饱和土中存在的基质吸力对土拱效应的演化具有重要影响。降雨、地下水位升降和管道渗漏等因素引起基质吸力发生变化，会导致松动区范围和土体强度的改变，从而给工程带来潜在的安全隐患。本文基于ABAQUS有限元软件建立不同饱和度下深埋活动门数值模型，根据数值试验结果，考虑非饱和土中基质吸力的影响，建立了非饱和砂土中深埋活动门松动土压力分析模型，随后结合大主应力圆弧拱理论和非饱和土抗剪强度理论，给出了松动区高度和松动土压力计算方法。最后与活动门试验结果和数值模拟结果对比验证了本文模型的合理性。
1. 非饱和土抗剪强度理论
非饱和土力学特性和活动门试验较为复杂，目前相关试验和理论研究尚处于起步阶段。随着活动门的移动，活动门上方土体发生应力重分布，该问题本质上属于土力学中的应力-应变关系和强度问题。目前，非饱和土抗剪强度表达式主要有两种。
（1）Bishop[21]提出的强度公式：
τ=c′+(σ−ua+χ(ua−uw))tanφ′。 (1) 式中：τ为破坏时剪切面上的剪应力；c′为有效黏聚力；σ−ua为破坏时剪切面上的净法向应力；ua为气压力；uw为孔隙水压力；ua−uw为基质吸力；χ为与饱和度相关的参数；φ′为有效内摩擦角。
（2）Fredlund等[22]提出的双变量公式：
τ=c′+(ua−uw)tanφb+(σ−ua)tanφ′。 (2) 式中：φb为抗剪强度随基质吸力变化的参数；(ua−uw)tanφb又被称为表观黏聚力。
从式（1），（2）可以看出，非饱和土中的抗剪强度主要由土体本身的有效黏聚力、外力相关的摩擦强度和吸力相关的黏聚强度。若令：
χtanφ′=tanφb。 (3) 则两者可以统一写成扩展的莫尔-库仑形式：
τ=c+(σ−ua)tanφ′。 (4) 式中：c为总黏聚力，c′+χ(ua−uw)tanφ′。
含水率或者饱和度对非饱和土的力学参数有较大影响，反映到强度参数中主要体现在内摩擦角和表观黏聚力或总黏聚力上，这两个土体参数和活动门本身的几何参数控制着活动门上方土体的破坏模式，因此，相对于其他复杂的模型，采用内摩擦角和表观黏聚力或总黏聚力来反映非饱和土的强度特性抓住了问题的主要矛盾，避免了问题复杂化，具有较好的合理性和适用性。
2. 非饱和砂土深埋活动门数值试验
对于活动门或地下结构，“深埋”和“浅埋”并没有严格的定义，当前研究中通常认为埋深比小于等于2的时候为浅埋，大于2为深埋。崔蓬勃等[17]在非饱和砂土中进行的活动门试验表明，对于埋深比小于等于2的“浅埋”活动门，随着饱和度的增加，滑移面表现为梯形-三角形-梯形的变化规律，滑移面并非都发展到地表。笔者认为“深埋”和“浅埋”应该以滑移面是否发展到地表为准，对于深埋活动门，随着活动门下移，形成稳定的土拱，松动区高度和滑移面形状随着饱和度的变化发生改变。本节首先建立深埋活动门数值模型，研究不同饱和度下活动门破坏模式。由于非饱和问题的复杂性，为定量评价不同饱和度对土拱效应演化和松动土压力的影响，试验土体假定为含水率均匀分布的均质土，同时，由于活动门试验时间较短，则在试验过程中可以不考虑重力作用引起的水分迁移和吸力再平衡过程，此外，忽略渗流对土拱效应的影响。
2.1 数值模型
如图 1所示，采用ABAQUS有限元软件建立活动门数值分析模型。该模型尺寸长为30 m，活动门宽D=6 m，上覆地层高度C取为2D~5D。模型两侧和底部分别进行法向约束和固定约束，上部为自由边界。
土体采用莫尔-库仑塑性本构模型，陈若曦等[23]研究结果表明弹性模量对土体破坏模式和最终的松动土压力影响较小。不改变其他参数，本文对弹性模量E=20~100 MPa的土体进行数值模拟发现：弹性模量的变化对地层位移曲线中的弹性阶段和最大拱阶段产生的影响较大，对于极限拱（稳定阶段）产生的影响较小，弹性模量越小，土压力达到各阶段的峰值所需要的位移也就越大，但其破坏模式和达到稳定状态时的所需的松动土压力几乎是没有变化的。因此本文算例中不考虑饱和度对土体弹性模量的影响，弹性模量E设为30 MPa，泊松比v为0.3，其他材料参数如表 1所示。
饱和度/% 干密度ρd/(g·cm-3) 有效内摩擦角φ′/(°) 总黏聚力c/kPa 0 1.70 36.0 0 9 1.70 35.3 1.0 14 170 35.0 2.0 22 1.70 34.3 6.2 36 1.70 33.8 9.4 45 1.70 32.9 8.3 63 1.70 31.4 5.5 81 1.70 30.4 3.2 90 1.70 30.1 2.4 100 1.70 29.6 0 活动门数值试验计算流程如下：首先，建立活动门试验数值模型，设置好边界条件，对模型施加初始应力状态，消除因重力载荷产生的竖向位移；然后，采用位移法控制活动门下降，其下降速率为-1×10-5 m/步，同时观察上覆土体的位移和应力状态；最后，当活动门上覆荷载不再随活动门位移发生变化后，即完成计算。
2.2 破坏模式的确定
图 2给出了C/D=2时，上覆土体的塑性应变云图和位移云图随饱和度Sr的变化。由该图可以确定上覆土体的破坏模式[14-15]，在非饱和工况中，基质吸力提供一定的表观黏聚力，该表观黏聚力随着饱和度增加呈现先增加后减小的趋势[19, 22]，表观黏聚力达到最大值时对应的饱和度本文称为临界饱和度。不同饱和度下活动门上覆土体的破坏模式随着基质吸力的变化呈现出不同的变化形式。低饱和度下，图 2所示位移云图和松动土压力结果与崔蓬勃等[17]结果较为一致，可以说明本文数值模拟结果的可靠性。当饱和度处于临界饱和度（Sr=36%）时，非饱和土中基质吸力引起的表观黏聚力最大，活动门上方产生的破坏剪切带从两侧开始发展，逐渐向内弯曲，土体滑移面最终形成一个松动区高度较低的半椭圆形破坏模式，随着基质吸力的减弱（Sr=22%， 45%），即表观黏聚力降低，新的破坏剪切带逐渐向上延伸，呈现出近似椭圆弧状，松动区高度相较于临界饱和度是有明显的增加，破坏模式仍为半椭圆形。随着表观黏聚力进一步降低（Sr=81%），半椭圆滑移面不断上移，松动区高度不断增加，底部等应变线和等位移线与垂直方向夹角逐渐减小。破坏剪切带沿活动门两侧垂直发展，随后逐渐向内弯曲，最终形成一个矩形加半椭圆的破坏模式。需要指出的是，此时上方半椭圆区域塑性应变相对较小，但已处于塑性破坏状态，同时该区域仍具有不小的位移，考虑到结构的安全性和可靠性，该部分仍需视为松动区，应予以考虑。在干砂和饱和工况中（Sr=0%， 100%），表观黏聚力为0，破坏剪切带呈垂直发展，并延伸至地表，最终呈垂直状。由上述可知，在非饱和工况中，随着饱和度的变化，活动门上覆土体的破坏模式主要分为半椭圆形和矩形加半椭圆组合模型。
![]()
图 2 不同饱和度时塑性应变云图和位移云图Figure 2. Plastic strain and displacement clouds at different saturations3. 分析模型
对于深埋活动门而言，由于土拱效应，土体塌陷往往难以到达地表。在干砂中，土拱并不是真实存在的自然平衡拱，需要松动区提供一定的支撑力，而在非饱和土中，基质吸力提供一定的表观黏聚力，随着活动门向下位移，将会形成自然平衡拱承载上部荷载，松动区高度和滑移面形状随着饱和度的变化发生改变。结合上述数值模型试验结果，假设活动门上方非饱和砂土的破坏模式为以下两种，即半椭圆模型和矩形+半椭圆模型，如图 3所示。图 3中，活动门埋深为C，宽度为D，V1是由于活动门下移造成地层损失的土体体积。图 3（a）中该部分由上方土体形成的半椭圆滑移面体积V2填充，图 3（b）由半椭圆滑移面体积V′2和矩形滑移面体积V3填充。对于半椭圆破坏模式和矩形+半椭圆破坏模式，其松动区高度分别为H1和H。
基于上述模型，本文假设如下：①活动门上方的非饱和土体符合扩展的莫尔-库仑准则；②活动门向下移动时，滑移面内包含的土体随着活动门整体下移；③活动门下移引起的地层损失完全由滑面内的土体回弹填充；④活动门移动一定距离后（松动土压力趋于稳定状态），滑面内的土体处于极限平衡状态；⑤考虑到弹性模量对最终的松动土压力几乎没有影响，本文假设弹性模量不随土体饱和度变化；⑥由于活动门试验时间较短，忽略渗流对含水率分布的影响。
3.1 松动区高度确定
随着活动门下降，活动门上方土体随之下降，并发生回弹变形，填充由活动门移动造成的土层损失体积V1。假设活动门下移位移为s，则损失的体积V1为
V1=sD。 (5) 该损失体积由上方松动区土体填充。对于半椭圆破坏模式，该部分体积为V2：
{V_2} = {{\mathtt{π} D{H_1}} \mathord{\left/ {\vphantom {{\mathtt{π} D{H_1}} 4}} \right. } 4} 。 (6a) 对于半椭圆+矩形破坏模式，该部分体积为 {V'_2} 和 {V_3} ，其中， {V'_2} 则为
{V'_2} = {\mathtt{π}}{{{D^2}} \mathord{\left/ {\vphantom {{{D^2}} {(8\sqrt {1 - {\varepsilon ^2}} }}} \right. } {(8\sqrt {1 - {\varepsilon ^2}} }}) 。 (6b) 式中： \varepsilon 为椭球体的偏心率，Janelid等[25]提出的椭球体理论表明， \varepsilon 取值一般为0.90~0.98，最常用取值为0.92~0.96。
{V_3} 可由下式求得
{V_3} = D(H - {D \mathord{\left/ {\vphantom {D {(2\sqrt {1 - {\varepsilon ^2}} )}}} \right. } {(2\sqrt {1 - {\varepsilon ^2}} )}}) 。 (6c) 根据假设（3）可得
{V_1} = \delta {V_2} （半椭圆破坏模式）\text{，} (7a) {V_1} = \delta ({V'_2} + {V_3}) （半椭圆+矩形破坏模式）。 (7b) 式中： \delta 为上方土体的体积膨胀系数，可通过坍塌试验来确定，取值范围为0.02~0.24[12, 26]。
将式（5），（6a）代入式（7a），将式（5），（6a），（6c）代入式（7b），可得归一化后的松动区高度：
{H_1} = {{4s} \mathord{\left/ {\vphantom {{4s} {(\delta \mathtt{π} )}}} \right. } {(\delta \mathtt{π} )}} \text{，} (8a) \frac{H}{D} = \frac{s}{{\delta D}} + \left( {\frac{1}{2} - \frac{\mathtt{π} }{8}} \right)\frac{1}{{\sqrt {1 - {\varepsilon ^2}} }} 。 (8b) 由式（8）可知，松动土压力的计算高度与活动门归一化的位移 s{\text{/}}D 和土体的体积膨胀系数 \delta 有关。
3.2 侧压力系数
图 4为活动门下移过程中，松动区土单元体在扩展的莫尔应力圆 \tau {\text{ - }}(\sigma - {u_\text{a}}) 内的应力状态变化。当活动门未发生位移时，土体处于静止土压力 {K_0} 状态，随着活动门下移，地层损失增大，土拱效应逐渐发挥并达到极限平衡状态，图中A点和B点分别为滑移面处和中心线处点的应力状态，N为水平微分土条上任意一点。基于大主应力圆弧拱理论， {\sigma _\text{v}} - {u_\text{a}} 和 {\sigma _\text{h}} - {u_\text{a}} 分别为竖向净应力和水平向净应力， \theta 为大主应力方向与水平方向的夹角，其范围为[0， {\theta _0} ]。当滑移面为矩形时， {\theta _0} = \mathtt{π} /4 + \varphi '/2 ；当滑移面为“椭圆状”时， {\theta _0} = \mathtt{π} /4 + \varphi '/2 - \alpha ， \alpha 为滑移面与垂直面的夹角。
考虑水平土条微分单元内任意一点，竖向净应力、水平净应力和剪应力可分别表示为
{\sigma _\text{v}} - {u_\text{a}} = ({\sigma _3} - {u_\text{a}}){\cos ^2}\theta + ({\sigma _1} - {u_\text{a}})\text{sin}{^2}\theta \text{，} (9a) {\sigma _\text{h}} - {u_\text{a}} = ({\sigma _1} - {u_\text{a}}){\cos ^2}\theta + ({\sigma _3} - {u_\text{a}})\text{sin}{^2}\theta \text{，} (9b) \tau = (({\sigma _1} - {u_\text{a}}) - ({\sigma _3} - {u_\text{a}}))\cos \theta \sin \theta 。 (9c) 式中： {\sigma _1} - {u_\text{a}} 和 {\sigma _3} - {u_\text{a}} 分别为最大和最小主应力。
主动土压力系数 {K_\text{a}} 为
{K_\text{a}} = {\tan ^2}({45^ \circ } - \varphi '/2) = \frac{{{\sigma _3} - {u_\text{a}} + c\cot \varphi '}}{{{\sigma _1} - {u_\text{a}} + c\cot \varphi '}} 。 (10) 则水平向和竖向净应力由如下比值关系 \bar K ：
\bar K = \frac{{{\sigma _\text{h}} - {u_\text{a}} + c\cot \varphi '}}{{{\sigma _\text{v}} - {u_\text{a}} + c\cot \varphi '}} = \frac{{{{\cos }^2}\theta + {K_\text{a}}\text{sin}{^2}\theta }}{{{K_\text{a}}{{\cos }^2}\theta + \text{sin}{^2}\theta }} 。 (11) 对于圆弧拱，其轨迹线控制方程为
{x^2} + {z^2} = {{{D^2}} \mathord{\left/ {\vphantom {{{D^2}} {(4{{\sin }^2}{\theta _0})}}} \right. } {(4{{\sin }^2}{\theta _0})}} 。 (12) 则坐标 x 可写为
x = {{D\sin \theta } \mathord{\left/ {\vphantom {{D\sin \theta } {(2\sin {\theta _0})}}} \right. } {(2\sin {\theta _0})}} 。 (13) 作用在水平微分土条上的平均竖向净应力，可由下式求得
\begin{gathered} \left(\sigma_{\mathrm{av}}-u_{\mathrm{a}}\right)+c \cot \varphi^{\prime}=\frac{2}{D} \int_0^{D / 2}\left(\left(\sigma_3-u_{\mathrm{a}}+c \cot \varphi^{\prime}\right) \cos ^2 \theta+\right. \\ \left.\left(\sigma_1-u_{\mathrm{a}}+c \cot \varphi^{\prime}\right) \sin ^2 \theta\right) \mathrm{d} x \\ =\frac{2}{D} \int_0^{\theta_0}\left(\left(\sigma_1-u_{\mathrm{a}}+c \cot \varphi^{\prime}\right)\left(K_{\mathrm{a}} \cos ^2 \theta+\sin ^2 \theta\right)\right) \\ \frac{D}{2 \sin \theta_0} \cos \theta \mathrm{~d} \theta 。 \end{gathered} (14) 积分后可得
{\sigma _\text{av}} - {u_\text{a}} + c\cot \varphi '\\ =\left(\frac{1}{3}(1-{K}_\text{a}){\mathrm{sin}}^{2}{\theta }_{0}+{K}_\text{a}\right)({\sigma }_{1}-{u}_\text{a}+c\mathrm{cot}{\varphi }^{\prime })\text{ }。 (15) （1）椭圆滑移面侧压力系数
将椭圆松动体划分 n 个无限小的水平微分土条，每个土条形状近似为梯形，其中滑移面与垂直面的夹角为 \alpha 。图 5所示为水平微分土条在倾斜滑移面上的应力状态以及主应力轴偏转状态示意图，根据单元体在水平方向和竖直方向的受力平衡条件可得
{\sigma }_\text{n}-{u}_\text{a}=({\sigma }_{1}-{u}_\text{a}){\mathrm{cos}}^{2}({\theta }_{0}+\alpha )+({\sigma }_{3}-{u}_\text{a}){\mathrm{sin}}^{2}({\theta }_{0}+\alpha )\text{，} (16) {\tau }_{A}=(({\sigma }_{1}-{u}_\text{a})-({\sigma }_{3}-{u}_\text{a}))\mathrm{sin}({\theta }_{0}+\alpha )\mathrm{cos}({\theta }_{0}+\alpha )\text{ }。 (17) 式中： {\tau _A} 为倾斜滑移面上的剪应力； {\sigma _\text{n}} 为倾斜滑移面上的法向应力，由于滑移面倾斜角度的原因，此时 {\theta _0} 取值为 \mathtt{π} {\text{/4}} + \varphi '{\text{/2}} - \alpha 。经整理，滑移面上的剪应力和法向应力可分别表示为
{\sigma _\text{n}} - {u_\text{a}} + c\cot \varphi ' = ({\sigma _1} - {u_\text{a}} + c\cot \varphi ')({\cos ^2}({\theta _0} + \alpha ) + \\ {K_\text{a}}{\sin ^2}({\theta _0} + \alpha )) \text{，} (18) {\tau _A} = ({\sigma _\text{n}} - {u_\text{a}} + c\cot \varphi ')\tan \varphi ' 。 (19) 由式（15），（18）可得椭圆滑移面A点处侧压力系数 {K_{A1}} 为
{K}_{A1}=\frac{{\sigma }_\text{n}-{u}_\text{a}+c\mathrm{cot}{\varphi }^{\prime }}{{\sigma }_\text{av}-{u}_\text{a}+c\mathrm{cot}{\varphi }^{\prime }}=\frac{{\mathrm{cos}}^{2}({\theta }_{0}+\alpha )+{K}_\text{a}{\mathrm{sin}}^{2}({\theta }_{0}+\alpha )}{\frac{1}{3}(1-{K}_\text{a}){\mathrm{sin}}^{2}{\theta }_{0}+{K}_\text{a}}。 (20) （2）矩形滑移面侧压力系数
对于竖直滑移面，任意点处的竖向净应力和竖向平均净应力间关系 m 可表达为
m = \frac{{{\sigma _\text{v}} - {u_\text{a}} + c\cot \varphi '}}{{{\sigma _\text{av}} - {u_\text{a}} + c\cot \varphi '}} = \frac{{{K_\text{a}}{{\cos }^2}\theta + \text{sin}{^2}\theta }}{{\frac{1}{3}(1 - {K_\text{a}}){{\sin }^2}{\theta _0} + {K_\text{a}}}} 。 (21) 则有：
{\sigma _\text{av}} - {u_\text{a}} = {{({\sigma _{\text{v}A}} - {u_\text{a}} + c\cot \varphi ')} \mathord{\left/ {\vphantom {{({\sigma _{\text{v}A}} - {u_\text{a}} + c\cot \varphi ')} {{m_A}}}} \right. } {{m_A}}} - c\cot \varphi ' 。 (22) 式中： {\sigma _{\text{v}A}} - {u_\text{a}} 和 {m_A} 分别为竖直滑移面A点处竖向净应力和系数 m 值。
由式（11），（22）可得竖直滑移面A处等效侧压力系数 {K_{A2}} 为
{K}_{A}{}_{2}=\frac{{\sigma }_{\text{h}A}-{u}_\text{a}+c\mathrm{cot}{\varphi }^{\prime }}{{\sigma }_{\text{v}A}-{u}_\text{a}+c\mathrm{cot}{\varphi }^{\prime }}=\overline{K}{m}_{A}=\frac{{\mathrm{cos}}^{2}{\theta }_{0}+{K}_\text{a}{\mathrm{sin}}^{2}{\theta }_{0}}{\frac{1}{3}(1-{K}_\text{a}){\mathrm{sin}}^{2}{\theta }_{0}+{K}_\text{a}}。 (23) 3.3 椭圆体松动土压力计算
将式（18），（20）相结合可得
{\sigma }_{1}-{u}_\text{a}+c\mathrm{cot}{\varphi }^{\prime }=\frac{{K}_{A}{}_{1}({\sigma }_\text{av}-{u}_\text{a}+c\mathrm{cot}{\varphi }^{\prime })}{{\mathrm{cos}}^{2}({\theta }_{0}+\alpha )+{K}_\text{a}{\mathrm{sin}}^{2}({\theta }_{0}+\alpha )}。 (24) 将式（24）代入式（18），（19）可得
{\sigma _\text{n}} - {u_\text{a}} + c\cot \varphi ' = {K_{A1}}({\sigma _\text{av}} - {u_\text{a}} + c\cot \varphi ') \text{，} (25) {\tau _A} = {K_A}_1({\sigma _\text{av}} - {u_\text{a}} + c\cot \varphi ')\tan \varphi ' 。 (26) 对图 6所示的梯形水平微分土条进行受力平衡分析可得
{F}_\text{v}+ \text{d}W={F}_\text{v}+d{F}_\text{v}+2({F}_{\tau }\mathrm{cos}\alpha -{F}_{\text{n}}\mathrm{sin}\alpha )\text{ }。 (27) 式中： \text{d}W 为土层单元的重力； {F_\text{v}} 和 {F_\text{v}} + \text{d}{F_\text{v}} 分别为作用于土层单元顶部和底部的竖向力； {F_\tau } 为滑移面的摩擦力， {F_{\text{n}}} 为滑移面的法向力。其中：
\left.\begin{array}{l} F_{\mathrm{v}}=\left(\sigma_{\mathrm{av}}-u_{\mathrm{a}}\right) D_z, \\ \mathrm{~d} F_{\mathrm{v}}=\mathrm{d}\left(\sigma_{\mathrm{av}}-u_{\mathrm{a}}\right) D_z+2\left(\sigma_{\mathrm{av}}-u_{\mathrm{a}}\right) \tan \alpha \mathrm{d} z, \\ \mathrm{~d} W=\gamma D_z \mathrm{~d} z, \\ F_{\mathrm{n}}=\sigma_{\mathrm{n}} \mathrm{~d} z / \cos \alpha, \\ F_\tau=\tau_A \mathrm{~d} z / \cos \alpha, \\ D_z=D+2\left(z-H_{\mathrm{r}}\right) \tan \alpha 。 \end{array}\right\} (28) 将式（25）、（26）代入式（28）可得：
\frac{{\text{d}({\sigma _\text{av}} - {u_\text{a}})}}{{dz}} + \frac{P}{{{D_z}}}({\sigma _\text{av}} - {u_\text{a}} + c\cot \varphi ') = \gamma 。 (29) 式中： P = 2[(1 - {K_{A1}})\tan \alpha + {K_{A1}}\tan \varphi '] 。
在对式（29）进行求解时，其椭圆面上任意位置夹角 \alpha 不相同且难以确定。因此，本文采用文献[13]提出的等效角度 {\alpha _1}^\prime 来代替椭圆面上任意位置与水平方向的夹角 {\alpha _1} 。通过对椭圆面整个弧长上的角度进行积分求得等效角度，如下所示：
\left.\begin{array}{l} \alpha_1^{\prime}=\operatorname{arcoc}\left(\frac{\int_{A^{\prime}}^{B^{\prime}} \cos \alpha_1 \mathrm{~d} l}{\int_{A^{\prime}}^{B^{\prime}} \mathrm{d} l}\right), \\ \alpha^{\prime}=90^{\circ}-\alpha_1^{\prime} \text { 。 } \end{array}\right\} (30) 式中： A' 点为椭圆的顶点； B' 点为椭圆体与矩形体的交界点。
将式（30）代入式（29）求解可得
{\sigma _\text{av}} - {u_\text{a}} = {C_1}{{\text{(}}{D_z}{\text{)}}^{ - \frac{P}{{2\tan \alpha }}}} + \frac{{\gamma {D_z}}}{{2\tan \alpha + P}} - c\cot \varphi ' 。 (31) 根据现有模型试验[17-18]研究发现，对于深埋活动门而言，松动区往往难以达到地表。由于非饱和土中基质吸力的作用，未扰动区域与松动区之间存在明显的“开裂”现象，因此，上部土体自重并未对下部松动区产生影响。根据边界条件，当 z = 0 时， {\sigma _\text{av}} - {u_\text{a}} = 0 ，将边界条件代入式（31）可得 {C_1} 的表达式为
{C}_{1}=-\frac{\gamma {(D-2{H}_{1}\mathrm{tan}\alpha )}^{1+\frac{P}{2\mathrm{tan}\alpha }}}{2\mathrm{tan}\alpha +P}+c\mathrm{cot}{\varphi }^{\prime }{(D-2{H}_{1}\mathrm{tan}\alpha )}^{\frac{P}{2\mathrm{tan}\alpha }}。 (32) 当式（31）中 z = {H_1} 时，可得作用在矩形松动区顶部的压力 {q_1} 为
{q_1} = {\sigma _\text{av}} = {C_1}{D^{ - \frac{P}{{2\tan \alpha }}}} + \frac{{\gamma D}}{{2\tan \alpha + P}} - c\cot \varphi ' + {u_\text{a}} 。 (33) 3.4 矩形体松动土压力计算
图 7给出水平微分土条受力示意图，由受力平衡条件可得
D \text{d}({\sigma _\text{av}} - {u_\text{a}}) + 2\tau \text{d}z = D\gamma \text{d}z 。 (34) 滑移面处的剪应力可表示为
\tau = \bar K({\sigma _{\text{v}A}} - {u_\text{a}} + c\cot \varphi ')\tan \varphi ' 。 (35) 将式（22），（23），（35）代入（34）得
\frac{{\text{d}({\sigma _\text{av}} - {u_\text{a}})}}{{\text{d}z}} + \frac{{2{K_{A2}}({\sigma _\text{av}} - {u_\text{a}})\tan \varphi '}}{D} + \frac{{2c{K_A}_2}}{D} = \gamma 。 (36) 求解式（36）可得
{\sigma _\text{av}} - {u_\text{a}} = C{\text{e}^{ - \frac{{2z{K_{A2}}\tan \varphi '}}{D}}} + \frac{{D\gamma }}{{2{K_{A2}}\tan \varphi '}} - \frac{c}{{\tan \varphi '}} 。 (37) 将边界条件 {\left. {({\sigma _\text{av}} - {u_\text{a}})} \right|_{z = 0}} = {q_1} 代入式（37）可得
{\sigma _\text{av}} = \left( {\frac{{D\gamma }}{{2{K_{A2}}\tan \varphi '}} - \frac{c}{{\tan \varphi '}}} \right)(1 - {\text{e}^{ - \frac{{2z{K_A}_2\tan \varphi '}}{D}}}) + \\ {q_1}{\text{e}^{ - \frac{{2z{K_{A2}}\tan \varphi '}}{D}}} + {u_\text{a}} 。 (38) 4. 模型验证
为验证理论模型的合理性，将本文数值模拟和理论模型计算结果与崔蓬勃等[17]活动门模型试验和现有理论解答进行对比验证。崔蓬勃等[17]针对中粗砂，配置了饱和度为0%~36%的4组试样，开展了埋深比 C{\text{/}}D 为1~3的低含水率活动门试验，分别通过粒子图像测速（PIV）技术获得活动门下移过程中的位移场云图和埋置在活动门上方的土压力盒获得作用在活动门上的土压力。黏聚力c和内摩擦角 \varphi ' 的取值参考表 1，参数 \delta 的经验系数难以确定，本文通过式（8）推算不同含水率下的膨胀系数。首先根据崔蓬勃等[17]文中给出的位移云图和地层位移曲线确定不同含水率情况下的松动区高度H和活动门移动距离 \varepsilon ，活动门的宽度D已知，椭球体的偏心率 \varepsilon 取为0.92，即可由式（8）反算出膨胀系数 \delta 。当饱和度 {S_r} 分别为0%，14%，22%，36%时，膨胀系数 \delta 为0.03，0.05，0.11，0.18。
图 8给出了不同饱和度 {S_{\text{r}}} 下归一化松动土压力 {\sigma _av}/({\gamma _{\text{d}}}C) 的比较，从该图可以看出，除 C{\text{/}}D = 2 时饱和度0.36的点外，本文理论解答和数值解与活动门试验数据相对误差在分别在1%~30%和3%~27%，其原因可能有：①理论解和数值解计算的是活动门上松动土压力的平均值，而试验数据给出的是具体点位的土压力，未考虑松动土压力在活动门上呈非均匀分布状态。②本文理论解答和数值模型均以莫尔-库仑理论为本构模型，作为一种理想弹塑性模型，该模型简单易用，抓住了土体呈现剪切破坏的本质，然而，该模型尚不能全面反映土体的力学特性，尤其是对于复杂的非饱和土体。整体而言，本文理论和数值模型的计算结果与文献[17]的试验值无论在量级还是变化趋势上都表现出了较好的一致性，相对误差在工程尺度可接受的范围，证明本文的理论和数值模型具有较好的可靠性。已有理论模型的计算结果要大于试验结果，尤其在临界饱和度时，计算结果偏差较大。本文的数值模型结果要略小于理论计算结果，更接近模型试验的最小值，而理论模型的计算结果要更接近模型试验的最终值。活动门模型试验和数值模型可以很好验证本文理论模型的合理性。
![]()
图 8 不同饱和度下松动土压力结果对比Figure 8. Comparison of loosening pressures under different saturations5. 结果分析
5.1 松动土压力
图 9给出不同 C{\text{/}}D 下，归一化松动土压力 {\sigma _{{\text{av}}}}/({\gamma _{\text{d}}}C) 随饱和度 {S_{\text{r}}} 的变化。由该图可以看出，随着饱和度的增加，松动土压力先缓慢减小，在某一饱和度时松动土压力达到最小值，随后开始逐渐增加。由非饱和抗剪强度公式（4）可知，非饱和土中的抗剪强度主要由土体的抗剪强度指标（黏聚力c和内摩擦角 \varphi ' ）所决定，其中，黏聚力包括有效黏聚力和表观黏聚力，而表观黏聚力是随基质吸力发生变化的。对于本文研究的土体，如表 1所示，其表观黏聚力随着基质吸力的增加呈现先增加后减小的趋势，即存在某一临界饱和度下对应的表观黏聚力最大。在临界饱和度之前，表观黏聚力随着饱和度的增加逐渐增大，松动土压力逐渐减小。当土体饱和度超过临界点时，表观黏聚力随饱和度的增加开始减小，松动土压力随之增大。当饱和度处于临界点时，表观黏聚力达到最大值，活动门上方受到的松动土压力最小。由上述可知，饱和度对活动门上方受到的松动土压力有较大影响。
![]()
图 9 不同C/D下松动土压力随饱和度的变化Figure 9. Loosening pressures with saturation under different values of C/D与干砂和饱和土相比，非饱和土中存在的基质吸力对活动门的破坏模式具有较大影响。由表观黏聚力的表达式 \chi (u{}_\text{a} - {u_\text{w}})\tan \varphi ' 可知，土体饱和时，基质吸力为0，此时表观黏聚力为0，土体干燥时，饱和度为0，非饱和参数 \chi 值较小，当假定非饱和参数 \chi 取值等于饱和度时[26-29]，此时表观黏聚力同样为0，如图 2所示，两者在滑移面的形态和破坏模式上具有较多类似之处。临界饱和度处土体的松动土压力分别为干砂和饱和情况的1.83倍和3.83倍，对于砂性土，由于其渗透系数较大，极易受到外界环境如降雨、水位升降和渗漏等因素影响，从而造成土体饱和度的改变，如果不能及时处理，极易诱发潜在的工程风险。
5.2 松动区高度
图 10给出不同 C{\text{/}}D 下，归一化的松动区高度 H{\text{/}}D 随饱和度 {S_{\text{r}}} 的变化。由该图可知，随着饱和度的增加，松动区高度先降低后增加，松动区高度的变化与活动门归一化的位移 s{\text{/}}D 和土体的体积膨胀系数 \delta 有关。在临界饱和处，松动区失稳高度最低，在临界饱和度之前，松动区高度随着饱和度的增加逐渐降低，当土体饱和度超过临界点时，松动区高度随饱和度的增加逐渐增高。
![]()
图 10 不同C/D下松动区高度随饱和度的变化Figure 10. Loosening heights with saturation under different values of C/D6. 结论
根据数值模拟结果，提出了非饱和砂土中深埋活动门松动土压力分析模型，通过与现有模型试验、理论方法和数值模拟结果比较，验证了本文模型的合理性。最后通过系统分析土体参数对松动土压力的影响，得到3点结论。
（1）根据数值试验结果，考虑到非饱和中基质吸力的影响，将非饱和砂土深埋活动门松动土压力分析模型分为半椭圆形和矩形加半椭圆组合模型。
（2）推导了松动区归一化高度H/D的预测公式和松动土压力计算方法。其中，松动土压力的计算高度与归一化活动门位移 s{\text{/}}D 和土体的体积膨胀系数 \delta 有关。
（3）松动土压力随着饱和度的增加呈现先减小后增加的趋势，在临界饱和度处，松动土压力值最小。
-
![]()
图 2 不同饱和度时塑性应变云图和位移云图
Figure 2. Plastic strain and displacement clouds at different saturations
![]()
图 8 不同饱和度下松动土压力结果对比
Figure 8. Comparison of loosening pressures under different saturations
![]()
图 9 不同C/D下松动土压力随饱和度的变化
Figure 9. Loosening pressures with saturation under different values of C/D
![]()
图 10 不同C/D下松动区高度随饱和度的变化
Figure 10. Loosening heights with saturation under different values of C/D
饱和度/% 干密度 {\rho _d} /(g·cm-3) 有效内摩擦角 \varphi ' /(°) 总黏聚力c/kPa 0 1.70 36.0 0 9 1.70 35.3 1.0 14 170 35.0 2.0 22 1.70 34.3 6.2 36 1.70 33.8 9.4 45 1.70 32.9 8.3 63 1.70 31.4 5.5 81 1.70 30.4 3.2 90 1.70 30.1 2.4 100 1.70 29.6 0 
下载: 导出CSV
-
[1] TERZAGHI K. Theoretical Soil Mechanics[M]. New York: Wiley, 1943.
[2] 张宇, 陶连金, 刘军, 等. 考虑主应力偏转和土拱效应的干砂盾构隧道掌子面极限支护力计算方法研究[J]. 岩土工程学报, 2023, 45(3): 530-540. doi: 10.11779/CJGE20211349 ZHANG Yu, TAO Lianjin, LIU Jun, et al. Method for calculating limit support pressure of face of shield tunnels considering principal stress axis rotation and soil arching effects in dry sand[J]. Chinese Journal of Geotechnical Engineering, 2023, 45(3): 530-540. (in Chinese) doi: 10.11779/CJGE20211349
[3] 张常光, 吴凯, 隋建浩. 基于小主应力轨迹的上埋式涵管竖向土压力非线性描述[J]. 岩土工程学报, 2021, 43(12): 2200-2208. doi: 10.11779/CJGE202112006 ZHANG Changguang, WU Kai, SUI Jianhao. Nonlinear descriptions of vertical earth pressure against positive buried pipelines based on minor principal stress trajectory[J]. Chinese Journal of Geotechnical Engineering, 2021, 43(12): 2200-2208. (in Chinese) doi: 10.11779/CJGE202112006
[4] 王涛, 姬建. 砂土边坡桩间水平土拱机理与演变规律离散元分析[J]. 岩土工程学报, 2024, 46(8): 1742-1752. doi: 10.11779/CJGE20230309 WANG Tao, JI Jian. DEM analysis of the mechanism and evolution of horizontal soil arching between piles in sandy slope[J]. Chinese Journal of Geotechnical Engineering, 2024, 46(8): 1742-1752. (in Chinese) doi: 10.11779/CJGE20230309
[5] WANG H L, CHEN R P, CHENG W, et al. Full-scale model study on variations of soil stress in geosynthetic-reinforced pile-supported track bed with water level change and cyclic loading[J]. Canadian Geotechnical Journal, 2019, 56(1): 60-68. doi: 10.1139/cgj-2017-0689
[6] LAI F W, YANG D Y, LIU S Y, et al. Towards an improved analytical framework to estimate active earth pressure in narrow c–ϕ soils behind rotating walls about the base[J]. Computers and Geotechnics, 2022, 141: 104544. doi: 10.1016/j.compgeo.2021.104544
[7] COSTA Y D J, ZORNBERG J G. Active and passive arching stresses outside a deep trapdoor[J]. Acta Geotechnica, 2020, 15(11): 3211-3227. doi: 10.1007/s11440-020-00969-x
[8] LIANG L J, XU C J, CHEN Q Z, et al. Experimental and theoretical investigations on evolution of soil-arching effect in 2D trapdoor problem[J]. International Journal of Geomechanics, 2020, 20(6): 06020007. doi: 10.1061/(ASCE)GM.1943-5622.0001643
[9] LIU C, ZHANG S L, ZHANG D L, et al. Model tests on progressive collapse mechanism of a shallow subway tunnel in soft upper and hard lower composite strata[J]. Tunnelling and Underground Space Technology, 2023, 131: 104824. doi: 10.1016/j.tust.2022.104824
[10] WAN T, LI P F, ZHENG H, et al. An analytical model of loosening earth pressure in front of tunnel face for deep-buried shield tunnels in sand[J]. Computers and Geotechnics, 2019, 115: 103170. doi: 10.1016/j.compgeo.2019.103170
[11] LIN X T, SU D, SHEN X, et al. Effect of tunnelling-induced ground loss on the distribution of earth pressure on a deep underground structure[J]. Computers and Geotechnics, 2022, 147: 104781.
[12] LIN X T, CHEN R P, WU H N, et al. Calculation of earth pressure distribution on the deep circular tunnel considering stress-transfer mechanisms in different zones[J]. Tunnelling and Underground Space Technology, 2022, 119: 104211. doi: 10.1016/j.tust.2021.104211
[13] 宫全美, 张润来, 周顺华, 等基于颗粒椭球体理论的隧道松动土压力计算方法[J]. 岩土工程学报, 2017, 39(1): 99-105. doi: 10.11779/CJGE201701008 GONG Quanmei, ZHANG Runlai, ZHOU Shunhua, et al. Method for calculating loosening earth pressure around tunnels based on ellipsoid theory of particle flows[J]. Chinese Journal of Geotechnical Engineering, 2017, 39(1): 99-105. (in Chinese) doi: 10.11779/CJGE201701008
[14] LAI F W, CHEN S X, XUE J F, et al. New analytical solutions for shallow cohesive soils overlying trench voids under various slip surfaces[J]. Transportation Geotechnics, 2020, 25: 100411.
[15] CHEN F Q, LUO S C, LAI F W. New analytical solutions for cohesive–frictional soils above deep active trapdoors[J]. International Journal of Geomechanics, 2022, 22(12): 04022235.
[16] LIANG L J, XU C J. Numerical and theoretical research on stress distribution in the loosening zone of the trapdoor problem[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2019, 43(7): 1426-1447. doi: 10.1002/nag.2906
[17] 崔蓬勃, 朱永全, 刘勇, 等. 非饱和砂土隧道土拱效应模型试验及颗粒流数值模拟研究[J]. 岩土力学, 2021, 42(12): 3451-3466. CUI Pengbo, ZHU Yongquan, LIU Yong, et al. Model test and particle flow numerical simulation of soil arch effect for unsaturated sandy soil tunnel[J]. Rock and Soil Mechanics, 2021, 42(12): 3451-3466. (in Chinese)
[18] SONG J H, CHEN K F, LI P, et al. Soil arching in unsaturated soil with different water table[J]. Granular Matter, 2018, 20(4): 78. doi: 10.1007/s10035-018-0849-3
[19] 蔺港, 孔令刚, 詹良通, 等. 基于太沙基土拱效应考虑基质吸力影响的松动土压力计算模型[J]. 岩土力学, 2015, 36(7): 2095-2104. LIN Gang, KONG Linggang, ZHAN Liangtong, et al. An analytical model for loosening earth pressure considering matric suction based on Terzaghi soil arch effect[J]. Rock and Soil Mechanics, 2015, 36(7): 2095-2104. (in Chinese)
[20] 崔蓬勃, 朱永全, 刘勇, 等. 考虑土拱发挥过程的非饱和砂土盾构隧道极限支护力计算方法研究[J]. 岩土工程学报, 2020, 42(5): 873-881. doi: 10.11779/CJGE202005009 CUI Pengbo, ZHU Yongquan, LIU Yong, et al. Calculation of ultimate supporting forces of shield tunnels in unsaturated sandy soils considering soil arching effects[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(5): 873-881. (in Chinese) doi: 10.11779/CJGE202005009
[21] BISHOP A W. Factors controlling the shear strength of partly saturated soils[J]. ASCE Res Conf Shear Strength of Cohesive Soils, 1960.
[22] FREDLUND D G, MORGENSTERN N R, WIDGER R A. The shear strength of unsaturated soils[J]. Canadian Geotechnical Journal, 1978, 15(3): 313-321. doi: 10.1139/t78-029
[23] 陈若曦, 朱斌, 陈云敏, 等. 基于主应力轴旋转理论的修正Terzaghi松动土压力[J]. 岩土力学, 2010, 31(5): 1402-1406. CHEN Ruoxi, ZHU Bin, CHEN Yunmin, et al. Modified Terzaghi loozening earth pressure based on theory of main stress axes rotation[J]. Rock and Soil Mechanics, 2010, 31(5): 1402-1406. (in Chinese)
[24] 崔蓬勃. 低含水率砂层盾构隧道施工土拱效应研究[D]. 石家庄: 石家庄铁道大学, 2022. CUI Pengbo. Research on Soil Arching Effect of Shield Tunneling Construction in Low Moisture Content Sand Layer[D]. Shijiazhuang: Shijiazhuang Tiedao University, 2022. (in Chinese)
[25] JANELID I, KVAPIL R. Sublevel caving[J]. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1966, 3(2): 129-132.
[26] ZHANG H F, ZHANG P, ZHOU W, et al. A new model to predict soil pressure acting on deep burial jacked pipes[J]. Tunnelling and Underground Space Technology, 2016, 60: 183-196.
[27] ALONSO E E, PINYOL N M, GENS A. Compacted soil behaviour: initial state, structure and constitutive modelling[J]. Géotechnique, 2013, 63(6): 463-478.
[28] BOLZON G, SCHREFLER B A, ZIENKIEWICZ O C. Elastoplastic soil constitutive laws generalized to partially saturated states[J]. Géotechnique, 1996, 46(2): 279-289.
[29] BORJA R I. On the mechanical energy and effective stress in saturated and unsaturated porous continua[J]. International Journal of Solids and Structures, 2006, 43(6): 1764-1786.
[30] TARANTINO A. A possible critical state framework for unsaturated compacted soils[J]. Géotechnique, 2007, 57(4): 385-389.
-
其他相关附件
计量
- 文章访问数: 256
- HTML全文浏览量: 23
- PDF下载量: 76
