考虑温度与围压影响的冻结砂土非正交弹塑性本构模型

梁靖宇; 齐吉琳; 张跃东; 路德春; 李昊雯

doi:10.11779/CJGE20230455

考虑温度与围压影响的冻结砂土非正交弹塑性本构模型 English Version

1.
北京建筑大学土木与交通工程学院, 北京 100044
2.
北京工业大学城市与工程安全减灾教育部重点试验室, 北京 100124

基金项目:

国家自然科学基金项目 52108294

国家自然科学基金项目 41972279

国家自然科学基金项目 52025084

北京建筑大学金字塔人才培养工程项目 JDYC20220812

详细信息

作者简介:
梁靖宇(1988—)，男，博士，副研究员，主要从事岩土材料工程特性及其本构理论方面的教学与科研工作。E-mail: liangjy@bucea.edu.cn

通讯作者:
齐吉琳, E-mail: jilinqi@bucea.edu.cn

中图分类号: TU43
计量
- 文章访问数: 469
- HTML全文浏览量: 74
- PDF下载量: 115
出版历程
- 收稿日期: 2023-05-22
- 网络出版日期: 2024-03-24
- 刊出日期: 2024-08-31

Non-orthogonal elastoplastic model for frozen sand incorporating effects of temperature and confining pressure

1.
School of Civil and Transportation Engineering, Beijing University of Civil Engineering and Architecture, Beijing 100044, China
2.
Key Lab of Urban Security and Disaster Engineering, Ministry of Education, Beijing University of Technology, Beijing 100124, China

摘要

摘要: 冻土的受力变形特性受温度与围压的影响显著。为描述温度的影响，通过建立冻土三向拉伸强度与温度的非线性关系，并采用平移变换的方法将其考虑到屈服函数中。为描述围压的影响，通过建立潜在强度衰减因子，并将其引入统一硬化参量中，发展得出了考虑围压影响的硬化参量。最终，基于非正交弹塑性模型框架，在平移变换应力空间内建立了能够考虑温度和围压影响的冻土非正交弹塑性本构模型。模型预测结果与冻结粉砂三轴压缩试验结果的对比表明，所建立模型能够模拟不同温度和围压下冻结砂土应力–应变关系。其不仅能够反映冻土峰值抗剪强度随温度降低而增大的温度效应，还能描述围压增大条件下应力应变曲线由剪胀与软化型逐渐转变为剪缩与硬化型的规律。
- 冻土 /
- 本构模型 /
- 非正交塑性流动法则 /
- 温度效应 /
- 围压效应
Abstract: The mechanical properties of frozen soils are significantly affected by temperature and confining pressure. To characterize the effects of temperature, a nonlinear relationship between the three-dimensional tensile strength and the temperature is established, which is incorporated into the yield function based on the coordinate transformation method. To characterize the effects of confining pressure, a potential strength degradation factor is established, which is used to develop the hardening parameters that effectively account for the effects of the confining pressure. Finally, based on the framework of the non-orthogonal elastoplastic model, a non-orthogonal elastoplastic constitutive model for frozen soils that can consider the effects of temperature and confining pressure is developed in the coordinate transformation space. Comparisons between the model predictions and the triaxial compression test results of the frozen silty sand demonstrate that the developed constitutive model can simulate the stress-strain relationship of frozen silty sand under different temperatures and confining pressures. The developed constitutive model characterizes the temperature effects, i.e., the increase in the peak shear strength with decreasing temperature, and the confining pressure effects, i.e., the transition from shear dilation and softening to shear contraction and hardening as reflected by the stress-strain curve under the increasing confining pressure.
- frozen soil /
- constitutive model /
- non-orthogonal plastic flow rule /
- temperature effect /
- confining pressure effect

HTML全文

0. 引言

相较于传统的单台盾构掘进方式，盾构地中对接法的最大优点是无需设置中间到达井，两台盾构机相向掘进有效解决了盾构机掘进能力不足的问题，并且可大幅降低施工工期。1997年建成的东京湾海底隧道是早期实现盾构对接解体的经典工程案例^[1]，中国则是在2011年建成的广东狮子洋隧道上首次应用该方法^[2-3]。盾构对接法通常施工流程如下：先行盾构掘进至对接点、对接点注浆加固、后行盾构掘进至对接点、洞内拆解刀盘等设备、浇筑二次衬砌。刀盘拆卸后，对接段临空面拱顶发生下沉，围岩应力重新分布。为了保证地层稳定性，隧道围岩有限的变形并不能使压力拱充分发展，由于抗剪应力未达到极限条件，此时的成拱效应称为“不完全拱效应”，探究不完全拱的上覆围岩压力计算方法有助于更好指导实际工程。

明确压力拱拱体的范围是准确计算隧道围岩压力的必要前提。傅鹤林等^[4]根据隧道围岩弹塑性应力分布特点，提出以切向应力高于原岩应力区域作为拱体范围。王家全等^[5]提出以最大主应力最大值为内边界，以大小主应力偏转点为外边界。李奎^[6]认为以最大主应力最大值作为内边界是不充分的，提出根据水平应力减小区域判定内边界，以最大主应力矢量流线确定拱体外边界。Kong等^[7]提出以水平应力相较于初始状态进行判定隧道上方压力拱内外边界，但缺少对两侧边界的判定。Chen等^[8]通过计算水平应力和垂直应力在扰动后与初始状态之比区分对拱顶、拱脚范围。从上述研究可知对于压力拱范围的确定主要有围岩弹塑性分区法、主应力偏转法和主应力矢量流线法^[9-10]，但大多都是单独讨论拱圈范围，如何划分起转移上部垂直荷载作用的主拱圈范围并不明确。

针对地层成拱作用下的围岩上覆围岩压力计算，叶飞等^[11]针对软弱围岩中隧道掌子面前方形成的不完全拱，假设水平微分土条上的垂直应力为定值，推导了砂土和低黏聚力岩体中拱顶上方垂直应力的计算方法。Chen等^[12]指出由于隧道拱顶沉降由中心向两端递减，实际拱体范围内的垂直应力呈现为向隧道内凹的形式，极限状态下的砂土大主应力拱迹线为圆弧形，其计算所得的侧压力系数和同一土条上的垂直应力分布系数更符合实际。徐长杰等^[13]进一步探讨了圆弧、抛物线和悬链线3种不同大主应力拱迹线的区别。上述关于主应力拱迹线的研究均是基于水平微分土条上各点应力莫尔圆相同的假设^[14-15]，这主要适用于浅埋隧道。但是，随着隧道埋深的增大，水平微分土条中点与两端的主应力大小差距越加明显。

本文依托青岛盾构地中对接工程，通过COMSOL Multiphysics建立三维数值模型，基于地层主应力矢量流线偏转和各向应力变化确定盾构地中对接段压力拱的主拱圈范围，并通过考虑大主应力不完全偏转和内摩擦角不完全调用推导上覆围岩压力计算公式。

1. 计算模型的建立

青岛地铁8号线工程全长约3.26 km，隧道外直径为6.7 m，内直径为6.0 m，盾构机直径6.98 m，对接段隧道埋深为17.5 m。该区域上覆地层为沿海地区常见的土岩复合地层，地层参数见表 1。三维计算模型如图 1所示。模型外部尺寸：长80.0 m，宽73.7 m，高44.3 m。内部尺寸：盾构主机长9.0 m，相当于6环管片宽，刀盘宽度1.0 m，管片外直径6.7 m，内直径6.0 m。结合青岛地铁8号线的现场施工结果，先、后行盾构机实际对接偏差小于10 mm，刀盘最小间隔10 cm，假设模拟开挖过程中先、后行盾构机轴线水平和高程位置始终保持一致，盾构对接间距设置为0。刀盘仅模拟其侧面，通过施加梯度水平压力模拟掌子面支撑作用^[16]。

表 1 材料物理力学参数

Table 1. Physical and mechanical parameters of materials

地层	H/m	g/(kg·m^-3)	E/MPa	μ	c/kPa	$\varphi$ /(°)
回填土	2.2	1750	20	0.35	15	10
含有机质粉质黏土	4.5	1800	15	0.36	9	7
粉质黏土	4.1	2000	28	0.34	17	13
强风化安山岩	2.1	2200	1600	0.33	45	21
微风化安山岩	32.3	2500	10500	0.30	300	36
管片	—	2500	27600	0.20	—	—
盾壳	—	7800	206000	0.25	—	—
注浆层	—	2000	100	0.25	—	—
注：H为层厚，g为密度，E为弹性模量，μ为泊松比，c为黏聚力， $\varphi$ 为内摩擦角。

下载: 导出CSV

| 显示表格

图 1 三维计算模型

Figure 1. Three-dimensional model

下载: 全尺寸图片幻灯片

盾构模型设计3个研究步骤：①步骤一为初始地应力平衡，模型侧面采用辊支承约束水平位移，底面采用固定约束，施加重力；②步骤二为先行隧道和后行隧道相向掘进开挖，通过循环移除两环开挖土体并施加盾尾管片及水平推力实现，单侧共循环开挖26环，该步骤主要是为了考虑开挖对对接段应力重分布的影响；③步骤三为对接点两侧盾构机刀盘同时拆除，采用稳态计算方式，通过去除刀盘宽度方向上的壳单位实现，以分析两侧刀盘拆除后围岩稳定后的规律。

2. 压力拱主拱圈范围的确定

选择对接点位置处的纵截面以及横截面作为观测截面，该位置处的地层沉降最大，其大主应力流线成拱效应最显著。通过图 2（a）所示对接点位置纵截面的最大主应力矢量流线偏转图可知，处于隧道外边沿的主应力由水平方向转变为竖直方向，可判定该区域为开挖卸荷导致的围岩松动区，但未见明显的流线成拱，围岩压力在拱顶附近有所下降。从图 2（b）所示对接点位置横截面的最大主应力矢量图可知，在拱顶微风化安山岩和强风化安山岩分界线处的最大主应力发生明显偏转，最大偏转角度约90°。相较于隧道纵截面，横截面成拱效应更为显著。

图 2 大主应力矢量流线偏转图

Figure 2. Streamline deflection of principal stress vector

下载: 全尺寸图片幻灯片

2.1 外边界的确定

从图 3所示刀盘拆卸后拱顶应力曲线可知，在拱顶以上5.5 m处水平应力突然增大，转变为最大主应力垂直应力。当刀盘拆除时，临空面顶部围岩下陷，由此导致隧道外边沿两侧围岩产生挤压作用。此时，上覆围岩压力在拱顶区域发生偏转。由于微风化岩层围岩强度较高，并且临空面两侧盾壳的支撑作用，拱顶围岩有限的沉降使得压力拱仅向上发展至软硬围岩分界面，而非像砂土的活动门试验，滑裂面一直延伸至地表。当超出分界面后，拱顶最大主应力由水平应力转换成垂直应力，表明强风化安山岩基本不承受压力拱水平挤压作用，压力拱外边界位于5.5 m处。

图 3 压力拱主拱圈外边界确定图

Figure 3. Determination of outer boundary main arch ring of pressure arch

下载: 全尺寸图片幻灯片

2.2 内边界的确定

定义水平应力变化指标A的表达式为

。

$A = \frac{{{\sigma_{x1}} - {\sigma_{x0}}}}{{{\sigma_{x0}}}} 。$

(1)

式中： ${\sigma_{x1}}$ 为刀盘拆卸后水平应力（kPa）； ${\sigma_{x0}}$ 为刀盘初始水平应力（kPa）。

从图 4（a）可知，在压力拱外边界，水平应力变化指标A达到最大，两侧围岩挤压程度最高，并随深度减小。具体分界可从图 4（b）水平应力变化指标曲线图得知，在隧道正上方0~0.5 m水平应力变化指标A呈现线性变化，当超出拱顶0.5 m后水平应力指标变化幅度缩小，并且水平应力由下降变成增大。根据上述分析确定主拱圈内边界位于拱顶0.5 m位置处，其并未与隧道拱顶上边沿重合。

图 4 压力拱主拱圈内边界确定图

Figure 4. Determination of inner boundary of main arch ring of pressure arch

下载: 全尺寸图片幻灯片

2.3 侧边界的确定

图 5中h表示水平监测截线与隧道拱顶间距。

，

${B_{\text{3}}} = D(1 + \cot \beta ) \text{，}$

(2)

。

$\beta = {45^ \circ } + \frac{\varphi }{2} 。$

(3)

图 5 隧道上部成拱区域围岩应力分布

Figure 5. Stress distribution of surrounding rock in upper arch area of tunnel

下载: 全尺寸图片幻灯片

式中：B₃为Terzarghi剪切面宽度（m）；D为隧道外直径（m）； $\varphi$ 为土体内摩擦角（°）。

图 5（a）所示剪应力分布表明，相较于按式（2），（3）计算所得Terzarghi地层剪切面^[17]宽度为9.8 m，数值计算得到的剪切面最大宽度仅6.0 m，并且随深度增加而减小，剪切面最小宽度为5.0 m，剪切面倾角与垂直面夹角为5.7°。

图 5（b）所示垂直应力分布图表明，在剪切面内垂直应力均出现下降，剪切面内侧围岩垂直应力通过剪切面转移至拱腰位置。图 5（c）所示水平应力分布图表明水平应力增大区域宽度相对较大，而在h=0.5 m截线处出现了水平应力下降的现象，说明围岩处于松动状态。

2.4 压力拱主拱圈模型

总结上述压力拱的主拱圈范围分析得出如图 6所示的简化模型。相较于传统砂土地层中不完全拱土拱的塔形滑移面，考虑一定黏聚力条件下的不完全压力拱的滑移面呈现出上大下小的盆状，这与叶飞等^[11]和吕伟华等^[18]通过模型试验观测的压力拱形状结果一致。同时，从图 6还可知压力拱外边界与拱顶沉降曲线存在一致性，拱顶沉降曲线的两个转折点分别发生在强风化安山岩和微风化安山岩分界线处。

图 6 压力拱主拱圈简化模型

Figure 6. Simplified model for main arch ring of pressure arch

下载: 全尺寸图片幻灯片

3. 不完全拱条件下的主拱圈围岩垂直压力计算

3.1 平均侧压力系数 ${K_\text{av}}$ 的确定

侧压力系数反映压力拱两侧稳定围岩对滑移面内围岩挤压程度，选取拱顶范围内任一水平微分土条，按下式计算土条上平均侧压力系数 ${K_\text{av}}$ ：

$K_{\mathrm{av}}=\frac{\sigma_{\mathrm{h}}^{\mathrm{A}}}{\sigma_{\mathrm{v}}^{\mathrm{av}}} 。$

(4)

式中： $\sigma_{\mathrm{h}}^{\mathrm{A}}$ 为水平微分土条上端点处的水平应力（kPa）； $\sigma_{\mathrm{v}}^{\mathrm{av}}$ 为微分土条上平均垂直应力（kPa）。

如所示，假设拱顶水平微分土条上各点均处于同一应力莫尔圆，土条上任意一点水平应力 ${\sigma_{\text{h}i}}$ 和垂直应力 ${\sigma_{\text{v}i}}$ 统一表达为关于最大主应力 ${\sigma_1}$ 偏转角 ${\theta_i}$ 的函数，即

${\sigma_{\text{h}i}} = {\sigma_1}{\cos ^2}{\alpha_i} + {\sigma_3}{\sin ^2}{\alpha_i} ={\sigma }_{1}({\mathrm{cos}}^{2}{\alpha }_{i}+{K}_\text{a}{\mathrm{sin}}^{2}{\alpha }_{i})\text{，}$

(5)

${\sigma_{\text{v}i}} = {\sigma_1}{\sin ^2}{\alpha_i} + {\sigma_3}{\cos ^2}{\alpha_i} = {\sigma_1}({\sin ^2}{\alpha_i} + {K_\text{a}}{\cos ^2}{\alpha_i}) \text{，}$

(6)

${K_\text{a}} = \frac{{{\sigma_3}}}{{{\sigma_1}}} = \frac{{1 - \sin \varphi }}{{1 + \sin \varphi }} = \frac{1}{{{K_\text{P}}}} \text{，}$

(7)

${\theta_i} = \frac{\pi }{2} - {\alpha_i} 。$

(8)

图 7 压力拱成拱机理

Figure 7. Forming mechanism of pressure arch

下载: 全尺寸图片幻灯片

式中： ${\sigma_1}$ ， ${\sigma_3}$ 分别为水平微分土条任意一点上的最大主应力和最小主应力（kPa）； ${\alpha_i}$ 为最大主应力与水平方向的夹角(°)，其值域为 $\left[ { - {\alpha_0}, {\alpha_0}} \right]$ ，其中， ${\alpha_i}$ 在水平微分土条中心位置等于0，在极限状态下滑移面位置等于 ${45^ \circ } + \varphi /{\text{2}}$ 。

从可知，随着最大主应力偏转角 ${\theta_i}$ 增大，水平应力逐渐增大，垂直应力逐渐减小，当 ${\theta_i} = {90^ \circ }$ 时，水平应力变为最大主应力，垂直应力变为最小主应力。因此，水平微分土条上端点处的水平应力 $\sigma_{\mathrm{h}}^{\mathrm{A}}$ 和垂直应力 $\sigma_{\mathrm{h}}^{\mathrm{A}}$ 可以表示为

$\sigma_{\mathrm{h}}^{\mathrm{A}} = {\sigma_1}{\cos ^2}{\alpha_0} + {\sigma_3}{\sin ^2}{\alpha_0} \text{，}$

(9)

$\sigma_{\mathrm{v}}^{\mathrm{A}} = {\sigma_1}{\sin ^2}{\alpha_0} + {\sigma_3}{\cos ^2}{\alpha_0} 。$

(10)

为了将平均垂直应力 $\sigma_{\mathrm{v}}^{\mathrm{av}}$ 表达为关于偏转角 ${\theta_i}$ 的函数，需要进一步明确最大主应力 ${\sigma_1}$ 的拱迹线。目前研究认为最大主应力拱迹线主要存在线性、圆弧形和抛物线形三种形式，最大主应力与水平方向夹角 ${\alpha_i}$ 与拱顶宽度B的关系表达式：

${\alpha_i} = \frac{{2y}}{B}\;(线性) ,$

(11)

$\sin {\alpha_i} = \frac{{2y}}{B}\sin {\alpha_0}\;(圆弧),$

(12)

$\tan {\alpha_i} = \frac{{2y}}{B}\tan {\alpha_0}\;(抛物线) 。$

(13)

在同一应力莫尔圆假设下，对水平微分土条上各点垂直应力进行积分可得平均垂直应力 $\sigma_{\mathrm{v}}^{\mathrm{av}}$ ，

$\sigma_{\mathrm{v}}^{\mathrm{av}} = \frac{1}{B}\int_{ - \frac{B}{{\text{2}}}}^{\frac{B}{{\text{2}}}} {({\sigma_1}{{\sin }^2}{\sigma_i} + {\sigma_3}{{\cos }^2}{\sigma_i})} \text{d}x$

$= \frac{{{\sigma_1}}}{B}\int_{ - \frac{B}{{\text{2}}}}^{\frac{B}{{\text{2}}}} {({{\sin }^2}{\sigma_i} + {K_{\text{a}}}{{\cos }^2}{\sigma_i})} \text{d}x 。$

(14)

结合式（11）可得 $\sigma_{\mathrm{v}}^{\mathrm{av}}$ 在线性、圆弧形和抛物线形拱迹线的表达式：

$\sigma_{\mathrm{v}}^{\mathrm{av}} = \frac{1}{2}{\sigma_1}\left[ {1 + {K_\text{a}} - \frac{{(1 - {K_\text{a}})\sin {\alpha_0} \cdot \cos {\alpha_0}}}{{{\alpha_0}}}} \right](线性),$

(15)

$\sigma_{\mathrm{v}}^{\mathrm{av}} = {\sigma_1}\left[ {1 + \frac{1}{3}(1 - {K_\text{a}}){{\sin }^2}{\alpha_0}} \right] \quad (圆弧),$

(16)

$\sigma_{\mathrm{v}}^{\mathrm{av}} = {\sigma_1}\left[ {1 + \frac{{({K_\text{a}} - 1)}}{{\tan {\alpha_0}}} \cdot {\alpha_0}} \right] \quad \quad (抛物线) 。$

(17)

将式（9），（17）代入式（4）即可求得平均侧压力系数 ${K_\text{av}}$ ，线性、圆弧和抛物线拱迹线的平均侧压力系数 ${K_\text{av}}$ 分别为

${K_\text{av}} = \frac{{2 \cdot ({{\cos }^2}{\alpha_0} + {K_\text{a}}{{\sin }^2}{\alpha_0})}}{{1 + {K_\text{a}} - (1 - {K_\text{a}}) \cdot \sin {\alpha_0} \cdot \cos {\alpha_0}/{\alpha_0}}} \;(线性),$

(18)

${K_\text{av}} = \frac{{{{\cos }^2}{\alpha_0} + {K_\text{a}}{{\sin }^2}{\alpha_0}}}{{{K_\text{a}} + 1/3 \cdot (1 - {K_\text{a}}) \cdot {{\sin }^2}{\alpha_0}}} \quad (圆弧)，$

(19)

${K_\text{av}} = \frac{{{{\cos }^2}{\alpha_0} + {K_\text{a}}{{\sin }^2}{\alpha_0}}}{{{K_\text{a}} \cdot {\alpha_0}/\tan {\alpha_0} + 1 - {\alpha_0}/\tan {\alpha_0}}} \quad (抛物线)。$

(20)

应当指出，上述拱迹线均是基于完全拱条件下得出的，而从前面分析可知盾构刀盘拆卸后形成的压力拱为不完全拱，实际最大主应力偏转程度小于上述3种拱迹线。

从所示应力莫尔圆可知，对于土条上任意计算单元，其最大主应力与水平方向夹角 ${\alpha_i}$ 可以表示为关于剪应力 ${\tau_{yz}}$ 、水平应力 ${\sigma_{\text{h}i}}$ 以及最小主应力 ${\sigma_3}$ 的函数，

${\tau_{yz}} = ({\sigma_{\text{h}i}} - {\sigma_3}) \cdot \tan {\alpha_i} 。$

(21)

根据式（21）计算得到如所示数值模拟结果下的最大主应力拱迹线，结合式（8）进一步反映最大主应力偏转角 ${\theta_i}$ 与地层成拱效应的显著程度。为了说明与不完全拱条件下拱迹线的区别，将数值模拟结果下的水平微分土条端点位置处的最大主应力与水平方向夹角 ${\alpha_0}$ 代入式（11）~（13），计算得到其余3种完全拱条件下的拱迹线。从可知，在h=5.5 m时的最大主应力偏转角 ${\theta_i}$ 基本在90°左右，并且由于监测截线中心偏转角和两侧偏转角相差并不大，所以不同形状拱迹线基本重合在一起；从可知，在h=4.5 m时最大主应力偏转角 ${\theta_i}$ 呈现为中间大两边小，最大偏转角 ${\theta_i} =$ 90°，各形状拱迹线差别并不明显；从图 8（c）可知，在h=3.5 m时的圆弧拱迹线条件下最大主应力偏转最显著，数值模拟的拱迹线居于第三位，说明地层压力拱的荷载传递能力出现下降；从图 8（d）~（f）可知，数值模拟拱迹线条件下最大主应力偏转程度随着埋深h减小而不断减小，说明成拱效应不断衰减。

图 8 最大主应力拱迹线

Figure 8. Arch traces of major principal stress

下载: 全尺寸图片幻灯片

根据式（4）计算得到如所示数值模拟结果下的平均侧压力系数 ${K_\text{av}}$ ，并将数值模拟结果下的水平微分土条端点位置处最大主应力与水平方向夹角 ${\alpha_0}$ 代入式（18）~（20），计算得到其余3种拱迹线条件下的平均侧压力系数 ${K_\text{av}}$ 。从可知，现有完全拱理论计算下的 ${K_\text{av}}$ 普遍偏大，高估了主拱圈两侧围岩的挤压作用，在计算上覆围岩压力时偏于不安全，特别是在距离拱顶较远位置，其偏差更加明显。其中，圆弧拱迹线的偏差最大，线性拱迹线的偏差次之，抛物线拱迹线偏差最小，这与前述拱迹线偏转程度的排序一致。

图 9 拱顶平均侧压力系数

Figure 9. Mean lateral pressure coefficients of arch vanlt

下载: 全尺寸图片幻灯片

3.2 调用内摩擦角 $\delta$ 的确定

由于围岩具有较大的黏聚力，压力拱滑移剪切面上的内摩擦角并未完全发挥，调用部分内摩擦角通常就能满足隧道开挖扰动引起的地应力再平衡。从所示的不完全压力拱条件下的应力莫尔圆可知，在无黏聚力情况下（即 ${\tau_\text{c}} - {O_\text{c}} - {\sigma_\text{c}}$ 坐标系），调用内摩擦角 $\delta$ 可表示为剪切面上的剪应力 $\tau_{yz}^\text{c}$ 和水平应力 $\sigma_\text{h}^\text{c}$ 的函数，

$\tau_{y z}^{\mathrm{c}}=\sigma_{\mathrm{h}}^{\mathrm{c}} \cdot \tan \delta 。$

(22)

图 10 不完全拱计算应力莫尔圆

Figure 10. Mohr circle of stress for calculation of incomplete arch

下载: 全尺寸图片幻灯片

在具有黏聚力情况下，通过向右平行移轴 $c \cdot \cot \varphi$ 的距离得到新的坐标系（即 $\tau - O - \sigma$ 坐标系），调用内摩擦角 $\delta$ 可进一步表示为新坐标系下剪切面上的剪应力 $\tau_{yz}^{}$ 和水平应力 $\sigma_\text{h}^{}$ 的函数，

$\tau_{y z}=\tau_{y z}^{\mathrm{c}}=\sigma_{\mathrm{h}}^{\mathrm{c}} \cdot \tan \delta=(\sigma+c \cdot \cot \varphi) \cdot \tan \delta。$

(23)

根据式（23）计算得到如所示距拱顶不同距离下的调用内摩擦角 $\delta$ 。从可知，在不完全拱的外边界（即h=5.5 m），地层调用的内摩擦角 $\delta$ 基本为0，可忽略不计。并且，随着距拱顶距离的减小，调用内摩擦角 $\delta$ 逐渐增大；随着距隧道中心距离的增加，调用内摩擦角 $\delta$ 逐渐增大。在不完全拱的内边界（即h=0.5 m）两侧滑移面位置处，调用内摩擦角 $\delta$ 达到最大，大约为20°，仍然小于地层最大内摩擦角 $\varphi = {36^ \circ }$ 。这说明刀盘拆卸后地层强度仍有较大富余，并未达到极限平衡状态，再平衡后形成的压力拱能够有效保证地层稳定性。

图 11 调用内摩擦角

Figure 11. Call of internal friction angle

下载: 全尺寸图片幻灯片

3.3 拱顶平均垂直应力 $\sigma_{\mathrm{v}}^{\mathrm{av}}$ 的计算

针对不完全压力拱的滑移面呈现倾斜向上的形状，本文通过微分的思想将其分为若干矩形微分土条进行适当简化，每条微分土条的厚度为1 m，其受力分析如图 12所示。

图 12 微分土条受力分析

Figure 12. Force analysis of differential soil strips

下载: 全尺寸图片幻灯片

对于砂土，其垂直受力平衡满足如下方程：

$B \sigma_{\mathrm{av}}^{\mathrm{c}}+\gamma B \mathrm{~d} z=B\left(\sigma_{\mathrm{av}}^{\mathrm{c}}+\mathrm{d} \sigma_{\mathrm{av}}^{\mathrm{c}}\right)+2 \tau_{y z}^{\mathrm{c}} \mathrm{~d} z 。$

(24)

整理得

$\frac{\mathrm{d} \sigma_{\mathrm{av}}^{\mathrm{c}}}{\mathrm{~d} z}+\frac{2 \tau_{y z}^c}{B}=\gamma 。$

(25)

其中：

$\tau_{y z}^{\mathrm{c}}=\sigma_{\mathrm{h}} \tan \delta=K_{\mathrm{av}}^{\mathrm{c}} \sigma_{\mathrm{av}}^{\mathrm{c}} \tan \delta 。$

(26)

解上述微分方程得

$\sigma_{\mathrm{av}}^{\mathrm{c}}=\frac{\gamma B}{2 K_{\mathrm{av}}^{\mathrm{c}} \tan \delta}+C \cdot \mathrm{e}^{-2 K_{\mathrm{av}}^{\mathrm{c}} \cdot \tan \delta \cdot \frac{z}{B}} 。$

(27)

当压力拱外边界存在超载为 $q$ 时，代入初始条件 $\left.\sigma_{\mathrm{av}}^{\mathrm{c}}\right|_{z=0}=q$ 时，得到

$\sigma_{\mathrm{av}}^{\mathrm{c}}=\frac{\gamma B}{2 K_{\mathrm{av}}^{\mathrm{c}} \tan \delta}\left(1-\mathrm{e}^{-2 K_{\mathrm{av}}^{\mathrm{c}} \cdot \tan \delta \cdot \frac{z}{B}}\right)+q \cdot \mathrm{e}^{-2 K_{\mathrm{av}}^{\mathrm{c}} \cdot \tan \delta \cdot \frac{z}{B}} 。$

(28)

对于具有一定黏聚力的岩土体，通过图 10所示平行移轴进行处理，

$\sigma_{\mathrm{av}}+c \cot \varphi=\sigma_{\mathrm{av}}^{\mathrm{c}} 。$

(29)

将式（30）代入式（24）得

$\gamma B \mathrm{~d} z=B \mathrm{~d} \sigma_{\mathrm{av}}+2 K_{\mathrm{av}}^{\mathrm{c}} \cdot\left(\sigma_{\mathrm{av}}+c \cot \varphi\right) \cdot \tan \delta \cdot \mathrm{d} z 。$

(30)

整理得

$\frac{\mathrm{d} \sigma_{\mathrm{av}}}{\mathrm{~d} z}+\frac{2 K_{\mathrm{av}}^{\mathrm{c}} \sigma_{\mathrm{av}} \tan \delta}{B}=\gamma-\frac{2 K_{\mathrm{av}}^{\mathrm{c}} c \cot \varphi \tan \delta}{B} 。$

(31)

解上述微分方程得

$\sigma_{\mathrm{av}}=q \cdot \mathrm{e}^{-2 K_{\mathrm{av}}^{\mathrm{c}} \tan \delta \cdot \frac{Z}{B}}+ \\ \frac{\left(1-\mathrm{e}^{-2 K_{\mathrm{av}}^{\mathrm{c}} \tan \delta \cdot \frac{z}{B}}\right) \cdot\left(\gamma B-2 K_{\mathrm{av}}^{\mathrm{c}} \cdot c \cot \varphi \cdot \tan \delta\right)}{2 K_{\mathrm{av}}^{\mathrm{c}} \tan \delta} 。$

(32)

其中，边界条件 $\sigma_\text{av}^{}\left| {_{z = 0} = } \right.q$ 与压力拱外边界上覆围岩压力有关，按照地层初始垂直应力取值。

所示为拱顶平均垂直应力计算参数。其中，平均侧压力系数 ${K_\text{av}}$ 和调用内摩擦角 $\delta$ 分别根据式（4），（23）确定。从可以看出，拱顶平均侧压力系数 ${K_\text{av}}$ 在压力拱内边界和外边界位置处相对较大，而在拱圈范围内有所减小；拱顶调用内摩擦角 $\delta$ 随着拱顶距离增大而减小，在压力拱的外边界达到0。将上述计算参数代入式（33）可以得到如图 13所示的拱顶各微分土条平均垂直应力。从图 13可知，通过考虑大主应力拱迹线不完全偏转，相较于抛物线、线性、圆弧形拱迹线的17.7%，25.0%，28.6%的最大误差，采用本文理论计算的结果与数值模拟最大误差进一步缩小为10.5%，能够更合理地反映地层不完全成拱状态下上覆围岩压力的计算。总体而言，通过数值模型计算和理论分析可以得出，在微风化安山岩地层中弃壳解体，原则上无需进行加固就可形成自承拱，其压力拱拱圈内垂直应力下降明显，内边界从初始的346 kPa减小为113 kPa（理论计算结果）或119 kPa（数值模型结果），相当于初始围压压力的1/3，能够有效保证刀盘拆卸后临空面的稳定性。

表 2 拱顶平均垂直应力计算参数

Table 2. Parameters of mean vertical stress of arch vault

距拱顶距离/m	土条宽度/m	平均侧压力系数	调用内摩擦角/(°)
5.5	6.0	0.956	0
4.5	5.8	0.895	7.48
3.5	5.6	0.858	11.13
2.5	5.4	0.847	13.59
1.5	5.2	0.861	15.54
0.5	5.0	0.950	17.07

下载: 导出CSV

| 显示表格

图 13 压力拱间平均垂直应力

Figure 13. Mean vertical stresses between arches

下载: 全尺寸图片幻灯片

4. 结论

依托青岛盾构地中对接工程，通过COMSOL Multiphysics建立三维数值模型，确定盾构地中对接过程中压力拱的主拱圈范围，并通过理论推导主拱圈范围内围岩压力计算方法，主要得出4点结论。

（1）由于刀盘厚度相对隧道直径较小，盾构对接段沿隧道轴向的大主应力未见矢量流线成拱，沿隧道径向的成拱效应更为显著。

（2）压力拱外边界位于软硬围岩分界面处，由于拱顶存在松动区域，内边界并未与隧道上边沿重合，而是距离拱顶一定高度，压力拱两侧边界倾斜向上，主拱圈总体上呈现为盆状。

（3）在刀盘拆卸后形成的不完全拱中，采用既有圆弧、线性、抛物线形拱迹线均高估了最大主应力偏转角度，相应所得的平均侧压力系数偏大，计算上覆围岩压力时偏于不安全。并且，随着距拱顶距离的减小，调用内摩擦角 $\delta$ 逐渐增大；随着离隧道中心距离的增加，调用内摩擦角 $\delta$ 逐渐增大。

（4）通过数值模型计算和理论分析得出在微风化安山岩地层中弃壳解体原则上无需进行加固，其压力拱拱圈内垂直应力下降明显，拱圈内边界垂直应力相当于初始围压压力的1/3。

图 1 在不同应力空间的冻土屈服曲线

Figure 1. Yield curves for frozen soils in different stress spaces

下载: 全尺寸图片幻灯片

图 2 潜在破坏应力比M_f的演化规律

Figure 2. Evolution of potential strength M_f

下载: 全尺寸图片幻灯片

图 3 改进Hvorslev面

Figure 3. Improved Hvorslev surface

下载: 全尺寸图片幻灯片

图 4 常规三轴压缩应力路径下屈服面的演化

Figure 4. Evolution of yield surface under triaxial compression stress path

下载: 全尺寸图片幻灯片

图 5 M_f与η的演化规律

Figure 5. Evolution of M_f and $\bar \eta$

下载: 全尺寸图片幻灯片

图 6 塑性流动方向的演化

Figure 6. Evolution of plastic flow direction

下载: 全尺寸图片幻灯片

图 7 剪切过程中的体胀规律

Figure 7. Variation rules of volume during shearing

下载: 全尺寸图片幻灯片

图 8 不同温度下冻土抗剪强度线

Figure 8. Shear strength lines for frozen soils at different temperatures

下载: 全尺寸图片幻灯片

图 9 冻土三向拉伸强度和温度的拟合曲线

Figure 9. Fitting curve of three-dimensional tensile strength and temperature of frozen soils

下载: 全尺寸图片幻灯片

图 10 冻结粉砂试验结果与模型预测对比

Figure 10. Comparison of test results and model predictions for frozen soils

下载: 全尺寸图片幻灯片

表 1 冻土的模型参数

Table 1 Model parameters of frozen soils

λ κ υ M $\varepsilon _{{\text{dr}}}^{\text{p}}$ χ μ

0.042 0.004 0.3 1.35 0.07 -0.392 1.51

下载: 导出CSV

参考文献(33)

[1]	LAI Y M, XU X T, DONG Y H, et al. Present situation and prospect of mechanical research on frozen soils in China[J]. Cold Regions Science and Technology, 2013, 87: 6-18. doi: 10.1016/j.coldregions.2012.12.001
[2]	马巍, 王大雁. 中国冻土力学研究50a回顾与展望[J]. 岩土工程学报, 2012, 34(4): 625-640. http://cge.nhri.cn/cn/article/id/14543 MA Wei, WANG Dayan. Studies on frozen soil mechanics in China in past 50 years and their prospect[J]. Chinese Journal of Geotechnical Engineering, 2012, 34(4): 625-640. (in Chinese) http://cge.nhri.cn/cn/article/id/14543
[3]	ZHAO Y H, ZHANG M Y, GAO J. Research progress of constitutive models of frozen soils: a review[J]. Cold Regions Science and Technology, 2023, 206: 103720. doi: 10.1016/j.coldregions.2022.103720
[4]	XU X T, WANG Y B, BAI R Q, et al. Comparative studies on mechanical behavior of frozen natural saline silty sand and frozen desalted silty sand[J]. Cold Regions Science and Technology, 2016, 132: 81-88. doi: 10.1016/j.coldregions.2016.09.015
[5]	KIM S Y, KIM Y, LEE J S. Effects of frozen water content and silt fraction on unconfined compressive behavior of fill materials[J]. Construction and Building Materials, 2021, 266: 120912. doi: 10.1016/j.conbuildmat.2020.120912
[6]	NIU Y Q, WANG X, LIAO M K, et al. Strength criterion for frozen silty clay considering the effect of initial water content[J]. Cold Regions Science and Technology, 2022, 196: 103521. doi: 10.1016/j.coldregions.2022.103521
[7]	孙晓宇, 齐吉琳, 尹振宇. 冻结饱和标准砂压缩性试验研究[J]. 岩土工程学报, 2018, 40(9): 1723-1728. doi: 10.11779/CJGE201809020 SUN Xiaoyu, QI Jilin, YIN Zhenyu. Experimental study on compressibility of frozen saturated ISO standard sand[J]. Chinese Journal of Geotechnical Engineering, 2018, 40(9): 1723-1728. (in Chinese) doi: 10.11779/CJGE201809020
[8]	高娟, 赖远明, 常丹, 等. 考虑加载速率影响的冻结含盐砂土强度准则研究[J]. 岩土工程学报, 2019, 41(1): 104-110. doi: 10.11779/CJGE201901011 GAO Juan, LAI Yuanming, CHANG Dan, et al. Strength criterion for frozen saline sand considering effects of loading rates[J]. Chinese Journal of Geotechnical Engineering, 2019, 41(1): 104-110. (in Chinese) doi: 10.11779/CJGE201901011
[9]	LI X, YAN Y, JI S Y. Mechanical properties of frozen ballast aggregates with different ice contents and temperatures[J]. Construction and Building Materials, 2022, 317: 125893. doi: 10.1016/j.conbuildmat.2021.125893
[10]	ZHANG D, LIU E L, LIU X Y, et al. A new strength criterion for frozen soils considering the influence of temperature and coarse-grained contents[J]. Cold Regions Science and Technology, 2017, 143: 1-12. doi: 10.1016/j.coldregions.2017.08.006
[11]	LUO F, LIU E L, ZHU Z Y. A strength criterion for frozen moraine soils[J]. Cold Regions Science and Technology, 2019, 164: 102786. doi: 10.1016/j.coldregions.2019.102786
[12]	LAI Y M, LIAO M K, HU K. A constitutive model of frozen saline sandy soil based on energy dissipation theory[J]. International Journal of Plasticity, 2016, 78: 84-113. doi: 10.1016/j.ijplas.2015.10.008
[13]	QI J L, HU W, MA W. Experimental study of a pseudo-preconsolidation pressure in frozen soils[J]. Cold Regions Science and Technology, 2010, 60(3): 230-233. doi: 10.1016/j.coldregions.2009.10.008
[14]	CHANG D, LAI Y M, YU F. An elastoplastic constitutive model for frozen saline coarse sandy soil undergoing particle breakage[J]. Acta Geotechnica, 2019, 14(6): 1757-1783. doi: 10.1007/s11440-019-00775-0
[15]	YAO X L, XU G F, ZHANG M Y, et al. A frozen soil rate dependent model with time related parabolic strength envelope[J]. Cold Regions Science and Technology, 2019, 159: 40-46. doi: 10.1016/j.coldregions.2018.12.006
[16]	ZHANG D, LIU E L. Binary-medium-based constitutive model of frozen soils subjected to triaxial loading[J]. Results in Physics, 2019, 12: 1999-2008. doi: 10.1016/j.rinp.2019.02.029
[17]	LAI Y M, JIN L, CHANG X X. Yield criterion and elasto-plastic damage constitutive model for frozen sandy soil[J]. International Journal of Plasticity, 2009, 25(6): 1177-1205. doi: 10.1016/j.ijplas.2008.06.010
[18]	LAI Y M, YANG Y G, CHANG X X, et al. Strength criterion and elastoplastic constitutive model of frozen silt in generalized plastic mechanics[J]. International Journal of Plasticity, 2010, 26(10): 1461-1484. doi: 10.1016/j.ijplas.2010.01.007
[19]	SUN K, ZHOU A N. A multisurface elastoplastic model for frozen soil[J]. Acta Geotechnica, 2021, 16(11): 3401-3424. doi: 10.1007/s11440-021-01391-7
[20]	YANG Y G, LAI Y M, DONG Y H, et al. The strength criterion and elastoplastic constitutive model of frozen soil under high confining pressures[J]. Cold Regions Science and Technology, 2010, 60(2): 154-160. doi: 10.1016/j.coldregions.2009.09.001
[21]	LU D C, LIANG J Y, DU X L, et al. Fractional elastoplastic constitutive model for soils based on a novel 3D fractional plastic flow rule[J]. Computers and Geotechnics, 2019, 105: 277-290. doi: 10.1016/j.compgeo.2018.10.004
[22]	LU D C, ZHOU X, DU X L, et al. A 3D fractional elastoplastic constitutive model for concrete material[J]. International Journal of Solids and Structures, 2019, 165: 160-175. doi: 10.1016/j.ijsolstr.2019.02.004
[23]	LI H C, TONG C X, CHANG X, et al. Constitutive modelling of temperature-dependent behaviour of soft rocks with fractional-order flow rule[J]. Applied Sciences, 2022, 12(8): 3875. doi: 10.3390/app12083875
[24]	LIANG J Y, LU D C, DU X L, et al. Non-orthogonal elastoplastic constitutive model for sand with dilatancy[J]. Computers and Geotechnics, 2020, 118: 103329. doi: 10.1016/j.compgeo.2019.103329
[25]	LIANG J Y, LU D C, ZHOU X, et al. Non-orthogonal elastoplastic constitutive model with the critical state for clay[J]. Computers and Geotechnics, 2019, 116: 103200. doi: 10.1016/j.compgeo.2019.103200
[26]	QU P F, ZHU Q Z, ZHAO L Y, et al. A micromechanics-based fractional frictional damage model for quasi-brittle rocks[J]. Computers and Geotechnics, 2021, 139: 104391. doi: 10.1016/j.compgeo.2021.104391
[27]	SUN K, TANG L, ZHOU A N, et al. An elastoplastic damage constitutive model for frozen soil based on the super/subloading yield surfaces[J]. Computers and Geotechnics, 2020, 128: 103842. doi: 10.1016/j.compgeo.2020.103842
[28]	姚仰平, 路德春, 周安楠. 岩土类材料的变换应力空间及其应用[J]. 岩土工程学报, 2005, 27(1): 24-29. doi: 10.3321/j.issn:1000-4548.2005.01.003 YAO Yangping, LU Dechun, ZHOU Annan. Transformed stress space for geomaterials and its application[J]. Chinese Journal of Geotechnical Engineering, 2005, 27(1): 24-29. (in Chinese) doi: 10.3321/j.issn:1000-4548.2005.01.003
[29]	姚仰平. UH模型系列研究[J]. 岩土工程学报, 2015, 37(2): 193-217. doi: 10.11779/CJGE201502001 YAO Yangping. Advanced UH models for soils[J]. Chinese Journal of Geotechnical Engineering, 2015, 37(2): 193-217. (in Chinese) doi: 10.11779/CJGE201502001
[30]	YAO Y P, GAO Z W, ZHAO J D, et al. Modified UH model: constitutive modeling of overconsolidated clays based on a parabolic hvorslev envelope[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2012, 138(7): 860-868. doi: 10.1061/(ASCE)GT.1943-5606.0000649
[31]	高娟, 赖远明. 冻结盐渍土三轴剪切试验过程中的损伤及压融分析[J]. 岩土工程学报, 2018, 40(4): 707-715. doi: 10.11779/CJGE201804015 GAO Juan, LAI Yuanming. Damage and pressure melting analysis of frozen saline soilsinprocess of triaxial compression tests[J]. Chinese Journal of Geotechnical Engineering, 2018, 40(4): 707-715. (in Chinese) doi: 10.11779/CJGE201804015
[32]	孙星亮, 汪稔, 胡明鉴. 冻土三轴剪切过程中细观损伤演化CT动态试验[J]. 岩土力学, 2005, 26(8): 1298-1302, 1311. doi: 10.3969/j.issn.1000-7598.2005.08.022 SUN Xingliang, WANG Ren, HU Mingjian. A CT-timely experimental study on meso-scopic structural damage development of frozen soil under triaxial shearing[J]. Rock and Soil Mechanics, 2005, 26(8): 1298-1302, 1311. (in Chinese) doi: 10.3969/j.issn.1000-7598.2005.08.022
[33]	赵淑萍, 马巍, 郑剑锋, 等. 基于CT单向压缩试验的冻结重塑兰州黄土损伤耗散势研究[J]. 岩土工程学报, 2012, 34(11): 2019-2025. http://cge.nhri.cn/cn/article/id/14883 ZHAO Shuping, MA Wei, ZHENG Jianfeng, et al. Damage dissipation potential of frozen remolded Lanzhou loess based on CT uniaxial compression test results[J]. Chinese Journal of Geotechnical Engineering, 2012, 34(11): 2019-2025. (in Chinese) http://cge.nhri.cn/cn/article/id/14883