基于D-P条件等效的黏性土隧道开挖地层力学行为有限元分析

摘要

基于三轴试验应力‒应变结果，获取了依托工程黏性土强度参数。通过德鲁克‒普拉格条件（D-P条件）等效，实现了摩尔‒库伦（M-C）黏性土材料的D-P准则模拟方法。采用“间隙约束法”模拟衬砌支护时机对地层力学行为的影响。通过数值计算讨论了黏性土盾构隧道地层参数以及施工特性对地层力学行为及支护反力的影响。研究发现：黏聚力、内摩擦角以及开挖洞径对围岩塑性区半径及支护反力的影响基本呈线性关系，随黏聚力和内摩擦角的增大，塑性区半径和支护反力明显减小，内摩擦角的影响显著于黏聚力。开挖洞径由6 m增大到14 m，塑性区半径增大约210 %，支护反力增大约230 %。土侧压力系数将显著影响开挖后塑性区形状、范围以及支护反力的分布。随着支护时间的延后，塑性区逐渐增大，支护反力逐渐减小，当释放的位移超过总位移的40 %时，两者增减幅度均有明显降低。就依托工程而言，当将地层含水率控制在8 %~15 %间，当地层收敛约为总位移的20 %时，可充分发挥地层自承能力从而降低支护反力。

关键词

隧道; 黏性土; 塑性区; 支护反力; 德鲁克‒普拉格条件等效

隧道工程的力学问题是土体与隧道结构相互作用、互相影响下的复杂问题。土体特性、隧道开挖及支护时机对地层和隧道结构的稳定性与受力有重要影响，在长期的研究中也受到了广泛的关注。其中数值模拟方法以其相对低廉的成本^［1］以及对复杂模型的适用性，成为了大量隧道力学问题研究中的重要手段。

有限单元法（FEM）自Clough ^［2］于1960年正式提出以来，在工程计算领域已有大量的应用。在隧道工程领域Atkinson和 Cairncross^［3］较早地使用有限单元法计算了浅埋隧道的稳定性。类似的问题也得到了Davis 等 ^［4］的关注。1983年时，FEM的提出者Clough ^［5］也针对土压平衡盾构隧道的土体响应进行了分析。上述早期的研究多采用线弹性的本构，在目前国内大量的工程研究报告乃至学位论文中也仍多有采用。后续，杨林德等^［6］给出了洞室开挖平面应变问题非线性有限元计算流程。Lyamin 和 Sloan ^［7］， Yamamoto等^［8］采用M-C非线性本构，得到了更为理想的隧道开挖围岩受力解析。郑颖人等^［9］通过多款有限元商业软件采用理想弹塑性本构，提出了使用有限元强度折减法求解安全系数评价隧道的稳定性。

再进一步发展，几何复杂的三维模型^［10］，更加精确的弹塑性、黏弹塑性本构^［11- 12］，以及基于节点的有限元法（NS-FEM） ^［13-15］、基于稳定节点的光滑有限元法（SNS-FEM） ^［16］等新的拓展方法的应用更加巩固了有限元数值模拟在隧道工程力学分析中的重要地位。

总结目前的研究可以发现，尽管复杂本构及几何模型的隧道力学有限元分析方法已经较为成熟，但是这类方法难度大且针对性强，往往需要自行编制相应计算程序。因此大量计算依旧通过商业软件，选用简单的弹性本构模拟隧道围岩。此类模型获取的所谓结果仅能从探寻规律的角度获得极为有限的结论。随着商业软件内置本构的丰富多元，采用M-C、D-P等模型获取隧道开挖力学行为的研究也越来越多，但是本构模型及屈服准则选用的准确性对数值模拟结果的可靠性有重要影响，但是在大量研究中多被忽略，缺少讨论^［17-19］。

综上所述，为了获取相对可靠的黏性土隧道土性参数及施工对地层受力的影响规律，本研究拟首先基于三轴试验获取依托工程原状土样的强度参数及应力‒应变发展规律。其后，探寻适用黏性土的屈服准则及数值模拟方法，通过对比三轴试验结果，证明本研究提出的屈服准则选用及数值模拟方法的适用性。最后，基于所提出的数值模拟方法，通过数值计算试验讨论土性参数（黏聚力、内摩擦角、侧压力系数）以及施工参数（开挖洞径、支护时机）对地层受力的影响，以期探明合理的规律，对黏性土地层的隧道结构设计及施工参数设置提供可靠参考。

1 三轴试验及强度参数的规划求解

西安地铁9号线疗养院‒华清池区间，盾构埋深处地层以黏性土（粉质黏土）为主。通过现场探井取土，制备土样，通过三轴剪切试验（图1）得到的土体的应力‒应变曲线如图2所示。

图1  土样制备及三轴试验

Fig. 1  Soil specimens preparation and triaxial tests

图2  不同含水率黏性土应力-应变关系

Fig. 2  Relationship of stress-strain of cohesive soil at different moisture contents

土体破坏点取偏应力峰值点，该土体在测试中，除8 %含水率下100 kPa和200 kPa围压测试曲线，其余均没有明显的峰值，按15 %应变对应主应力差值取值。得到土体的强度指标如表1所示。

表1  土体强度参数

Tab. 1  Soil strength parameters

序号	含水率	粘聚力/kPa	内摩擦角/º
1	8%	25.3	22.54
2	15%	26.4	21.70
3	20%	30.3	18.68

2 黏性土力学行为的数值模拟方法

强度准则的选用对黏性土地层相关工程力学问题的求解结果有重要影响。Mohr-Coulomb（M-C）强度理论对岩土工程材料的使用性被大量的工程实践证实，因此也得到了极为广泛的应用。大量土性参数的基本力学试验，如直接剪切试验、三轴试验等，最终获取的均为M-C准则下的强度参数。随着计算机的发展，有限单元法成为了岩土工程问题解析的重要分支。然而，M-C准则在π平面上的屈服曲线为存在角点的六边形，给数值计算的便捷性和准确性带来了影响。Drucker-Prager（D-P）屈服准则在π平面上的屈服面为圆形，在众多有限元软件中均有应用。陈鹏等^［20］的研究表明，将M-C准则中的强度参数进行一定的换算，从而使用D-P准则对岩土体求解，仍然能够获得合理的解答。前人学者为了降低使用D-P准则分析岩土工程问题的误差，发展形成了D-P系列准则（图3），包括了：M-C外角点外接圆准则（DP1）；M-C内角点外接圆准则（DP2）；M-C内切圆准则（DP3）；M-C等面积圆准则（DP4）；M-C匹配圆准则（DP5）。目前大量有限元软件如ANSYS中使用的D-P准则为DP1。

图3  各准则在π平面上的屈服曲线

Fig. 3  Yield of each criterion in deviator plane

对于M-C准则，用常见的主应力形式表达为

$$\frac{\mathrm{1}}{\mathrm{2}}\left({\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{3}}\right)=c\mathrm{c}\mathrm{o}\mathrm{s}\phi -\frac{\mathrm{1}}{\mathrm{2}}\left({\sigma }_{\mathrm{1}}+{\sigma }_{\mathrm{3}}\right)\mathrm{s}\mathrm{i}\mathrm{n}\phi$$

（1）

式中：σ₁为最大主应力；σ₃为最小主应力；c为黏聚力；φ为内摩擦角。

通过单元体平衡可以得到主应力和应力偏量及球应力的关系：

$$\left(\begin{array}{c}{\sigma }_{\mathrm{1}}\\ {\sigma }_{\mathrm{2}}\\ {\sigma }_{\mathrm{3}}\end{array}\right)=\frac{\mathrm{2}}{\sqrt[]{\mathrm{3}}}\sqrt[]{J_{\mathrm{2}}}\left[\begin{array}{c}\mathrm{s}\mathrm{i}\mathrm{n}\left({\theta }_{\sigma }+\frac{\mathrm{2}}{\mathrm{3}}\pi \right)\\ \mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\sigma }\\ \mathrm{s}\mathrm{i}\mathrm{n}\left({\theta }_{\sigma }-\frac{\mathrm{2}}{\mathrm{3}}\pi \right)\end{array}\right]+\left(\begin{array}{c}{\sigma }_{\mathrm{m}}\\ {\sigma }_{\mathrm{m}}\\ {\sigma }_{\mathrm{m}}\end{array}\right)$$

（2）

式中：σ₂为中间主应力；J₂为偏应力第二不变量；θ_σ为应力洛德角；σ_m为三个正应力的平均值。

将（2）代入（1）式有：

$$\begin{array}{l}f=\frac{\mathrm{1}}{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}\phi I_{\mathrm{1}}+\\ \left(\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\sigma }-\frac{\mathrm{1}}{\sqrt[]{\mathrm{3}}}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\sigma }\mathrm{s}\mathrm{i}\mathrm{n}\phi \sqrt[]{J_{\mathrm{2}}}\right)-\\ c\mathrm{c}\mathrm{o}\mathrm{s}\phi =\mathrm{0}\end{array}$$

（3）

式中：f为塑性势函数；I₁为应力张量第一不变量。

在平面应变问题中D-P准则可以描述为

$$f=\alpha I_{\mathrm{1}}+\sqrt[]{J_{\mathrm{2}}}-k$$

（4）

式中：α是与岩土材料内摩擦角φ有关的常数；k是与c和φ有关的常数；α与k可统称为D-P条件。

将式（4）变形有：

$$\begin{array}{l}f=\frac{\mathrm{s}\mathrm{i}\mathrm{n}\phi }{\sqrt[]{\mathrm{3}}\left(\sqrt[]{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{\sigma }}-\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{\sigma }}\mathrm{s}\mathrm{i}\mathrm{n}\phi \right)}{I}_{\mathrm{1}}+\\ \sqrt[]{J_{\mathrm{2}}}-\frac{\sqrt[]{\mathrm{3}}c\mathrm{c}\mathrm{o}\mathrm{s}\phi }{\sqrt[]{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{\sigma }}-\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{\sigma }}\mathrm{s}\mathrm{i}\mathrm{n}\phi }\end{array}$$

（5）

可得：

$$\begin{array}{l}\alpha =\frac{\mathrm{s}\mathrm{i}\mathrm{n}\phi }{\sqrt[]{\mathrm{3}}\left(\sqrt[]{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{\sigma }}-\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{\sigma }}\mathrm{s}\mathrm{i}\mathrm{n}\phi \right)}\\ k=\frac{\sqrt[]{\mathrm{3}}c\mathrm{c}\mathrm{o}\mathrm{s}\phi }{\sqrt[]{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{\sigma }}-\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{\sigma }}\mathrm{s}\mathrm{i}\mathrm{n}\phi }\end{array}$$

（6）

当材料处于单纯剪切状态时，D-P准则为M-C准则的匹配圆，θ_σ=0º，对应D-P条件α、k为

$$\begin{array}{l}\alpha =\frac{\mathrm{s}\mathrm{i}\mathrm{n}\phi }{\mathrm{3}}\\ k=c\mathrm{c}\mathrm{o}\mathrm{s}\phi \end{array}$$

（7）

ANSYS内置的D-P准则的屈服面为M-C屈服准则六角形外角点的外接圆（DP1），对应的θ_σ=30º，常数α、k为

$$\begin{array}{l}\alpha =\frac{\mathrm{2}\mathrm{s}\mathrm{i}\mathrm{n}\phi }{\sqrt[]{\mathrm{3}}\left(\mathrm{3}-\mathrm{s}\mathrm{i}\mathrm{n}\phi \right)}\\ k=\frac{\mathrm{6}c\mathrm{c}\mathrm{o}\mathrm{s}\phi }{\sqrt[]{\mathrm{3}}\left(\mathrm{3}-\mathrm{s}\mathrm{i}\mathrm{n}\phi \right)}\end{array}$$

（8）

考虑实际参与数值计算的D-P条件与M-C匹配圆等效，即取式（8）α、k与式（7）相等，可以得到用于数值计算的D-P模型参数与实际的强度参数的换算关系为

$$\begin{array}{l}\phi =\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}\left(\frac{\mathrm{3}\sqrt[]{\mathrm{3}}\alpha }{\mathrm{2}+\sqrt[]{\mathrm{3}}\alpha }\right)\\ c=\frac{\sqrt[]{\mathrm{3}}\left(\mathrm{3}-\mathrm{s}\mathrm{i}\mathrm{n}\phi \right)k}{\mathrm{6}\mathrm{c}\mathrm{o}\mathrm{s}\phi }\end{array}$$

（9）

为了验证上述理论过程的合理性，通过FEM模拟含水率为15 %的三轴试验。试验测试值用线表示，通过上述方法换算得到的计算结果用实心符号表示，未经过上述换算得到的计算结果用空心符号表示。显然采用D-P条件等效换算后，数值模拟得到的应力应变关系与试验值拟合性良好（图4）。不通过换算直接计算，数值模拟的结果，在应变超过5 %以后，偏应力明显高于三轴试验结果，20 %应变下的计算误差超过200 %。

图4  D-P条件等效换算前后结果与试验对比

Fig. 4  Comparison of results before and after equivalent conversion of D-P condition and experiments

3 土性参数对地层屈服特性影响规律

在第2节提出的数值模拟方法下，通过计算不同的黏聚力、内摩擦角、侧压力系数等下隧道开挖后的塑性区半径及支护反力，研究土性参数对地层屈服特性的影响情况。

3.1　黏聚力影响

数值计算模型内摩擦角φ取20º，侧压力系数取1，隧道开挖直径取6 m，塑性区计算时不考虑支护效应，计算支护反力是考虑盾尾间隙为10 cm，且盾壳相较土体为刚性支护，分别计算c为10、15、18、20、25、30、35、40 kPa时的围岩的塑性区发展及支护反力。计算结果如图5所示。

图5  黏聚力对塑性区半径及支护反力的影响

Fig.5  Effect of cohesion on radius of plastic zone and supporting reaction

由图5可知，随着黏聚力的增大，塑性区半径和支护反力均基本呈线性减小。黏聚力从10 kPa增大到40 kPa，塑性区半径由4.59 m减小为2.59 m，斜率为－0.067；黏聚力从10 kPa增大到27 kPa，支护反力由6.6 kPa减小为0.81 kPa，斜率为－0.34。

3.2　内摩擦角影响

数值计算模型黏聚力c取20 kPa，侧压力系数取1，隧道开挖直径取6 m，塑性区计算时不考虑支护效应，计算支护反力时考虑盾尾间隙为10 cm，且盾壳相较土体为刚性支护，分别计算内摩擦角为10 º、15º、19 º、20 º、25 º、30 º、35 º、40 º时，围岩的塑性区发展及支护反力。计算结果如图6所示。

图6  内摩擦角对塑性区半径及支护反力的影响

Fig. 6  Effect of internal friction angle on radius of plastic zone and supporting reaction

由图6可知，随着内摩擦角的增大塑性区半径和支护反力均基本呈线性减小。内摩擦角从10 º增大到40 º，塑性区半径由5.68 m减小为0.48 m，斜率为－0.17；内摩擦角从10 º增大到22 º，支护反力由23.8 kPa减小为0.72 kPa，斜率为－1.92。

3.3　侧压力系数

数值计算模型黏聚力c取30 kPa，内摩擦角φ取30 º，隧道开挖直径取6 m，塑性区计算时不考虑支护效应，计算支护反力时考虑盾尾间隙为10 cm，且盾壳相较土体为刚性支护，分别计算侧压力系数为0.2、0.5、1、1.5、2时的围岩的塑性区发展及支护反力。计算结果如图7所示。图中，MN和MX分别表示塑性应力的最小值点和最大值点。

图7  侧压力系数对塑性区的影响

Fig. 7  Effect of lateral pressure coefficient on radius of plastic zone

由图7可知，侧压力系数对塑性区的形态及塑性区范围有显著影响。当竖向压力与水平向压力的差异增大时，土体开挖后更容易发生破坏。当水平应力明显小于竖向应力时，塑性区呈“X”带状分布，从拱顶至边墙约60 º位置初始产生，与竖直方向呈小角度延伸。当竖向压力与水平压力相等时，隧道开挖引发的塑性破坏呈环形。当竖向压力小于水平压力时，塑性破坏区域由拱墙外侧转移至拱顶和拱底外侧。

支护反力随土侧压力系数的变化情况如图8所示，可知侧压力系数由0.2向2变化时，支护反力由竖“8”形变为横“8”形，与塑性区形状的发展恰好相反。侧压力较小时，拱顶和拱底支护反力较大，侧压力较大时，两侧边墙部位支护反力较大。

图8  侧压力系数对支护反力的影响

Fig. 8  Effect of lateral pressure coefficient on supporting reaction

4 盾构施工对地层屈服特性影响规律

4.1　开挖直径的影响

数值计算模型黏聚力c取20 kPa，内摩擦角φ取20 º，隧道开挖直径取6 m，塑性区计算时不考虑支护效应，计算支护反力时不考虑盾尾间隙且盾壳相较土体为刚性支护，分别计算开挖直径为6、8、10、12、14 m时的围岩的塑性区发展及支护反力。计算结果如图9所示。

图9  开挖洞径对塑性区半径及支护反力的影响

Fig. 9  Effect of tunnel diameter on radius of plastic zone and supporting reaction

由图9可知，随着开挖洞径的增大塑性区半径和支护反力均基本呈线性增大。洞径从6 m增大到14 m，塑性区半径由4.04 m增大到8.63 m，斜率为0.57，支护反力由10.97 kPa增大到25.45 kPa，斜率为1.81。

4.2　支护时机的影响

数值计算模型黏聚力c取18.906 kPa，内摩擦角φ取18.2 º，隧道开挖直径取6 m。采用“间隙约束法”模拟支护时机。首先计算无支护状态开挖轮廓线的最大位移，此时对开挖洞径一周的节点不施加约束。开挖后即刻支护对应围岩位移释放为0，即约束开挖轮廓一周节点的位移为0。开挖后位移释放20 %时，对应围岩位移释放为0.2，即约束开挖轮廓线上各节点沿径向发生最大位移的20 %，同理计算位移释放为总位移的0.4、0.6、0.8倍时围岩的塑性区及支护反力的发展情况。

由图10可知，随着支护时间的增长，塑性区逐渐增大，但是作用在支护结构上的反力随之减小。当释放的位移超过总位移的40 %时，塑性区半径的增加和支护反力的减小幅度均有明显降低。

图10  支护时机对塑性区半径及支护反力的影响

Fig. 10  Effect of supporting time on radius of plastic zone and supporting reaction

5 工程应用

对本工程原状土在8 %、15 %、20 %含水率下的隧道开挖及不同支护时机下地层稳定情况及支护反力进行计算。发现无支护状态下，8%含水率在第0.95计算子步时不收敛，对应拱顶下沉值为1.84 m；20 %含水率从在第0.83计算子部时不收敛，对应拱顶下沉值为1.7 m；1 5%含水率计算收敛，拱顶下沉值为2.11 m。表明本工程地层在8 %及20 %含水率下的自稳能力差，无支护状态下不能保障地层稳定。15 %含水率下，围岩自稳能力较好。不同支护时机下的支护反力如图11所示。8 %和15 %含水率支护特性曲线较为接近，在变形释放初期（20 %总位移内），支护力有明显的下降，释放后期支护力下降速率较为缓慢。20 %含水率地层自稳性差，其支护反力明显高于8 %及15 %含水率。因此，就本工程而言，将地层含水率控制在15 %以内，隧道开挖对地层稳定性的影响将相对较小。同时，如果将盾尾间隙控制在隧道总收敛的20%左右，可较为充分的发挥地层的自稳能力，保障地层稳定的同时，可较大幅度地降低管片后期承受的地层荷载。

图11  不同含水率土体支护特性曲线

Fig.11  Support characteristic of soil at different moisture contents

6 结论

（1）通过M-C准则匹配圆与M-C准则外角点的外接圆的D-P条件等效，可换算得到适用数值计算的土体黏聚力及内摩擦角值，获得较为准确的数值解析。

（2）黏聚力、内摩擦角以及开挖洞径对围岩塑性区半径及支护反力的影响基本呈线性关系，随黏聚力和内摩擦角的增大塑性区半径和支护反力明显减小，内摩擦角的影响显著于黏聚力。

（3）土侧压力系数将显著影响开挖后塑性区形状、范围以及支护反力的分布。水平应力明显小于竖向应力时，塑性区呈“X”带状分布，当竖向压力与水平压力相等时，塑性破坏呈环形。当竖向压力小于水平压力时，塑性破坏区域由拱墙外侧转移至拱顶和拱底外侧。支护反力呈现出与塑性区分布相反的特性。

（4）随着开挖洞径的增大，塑性区半径和支护反力均线性增大。开挖洞径由6 m增大到14 m，塑性区半径增大约210 %，支护反力增大约230 %。

（5）随着支护时间的延后，塑性区半径逐渐增大，支护反力逐渐减小，当释放的位移超过总位移的40%时，两者增减幅度均有明显降低。

（6）就依托工程而言，当将地层含水率控制在8 %~15 %间，当地层收敛约为总位移的20 %时，可充分发挥地层自承能力从而降低支护反力。

作者贡献声明

何乔： 论文初稿，算法和程序，以及绘图和修改论文。

韩兴博： 指导论文思路和框图，以及绘图和校核论文。

邓念兵： 调研土工试验的相关资料。

杨海挺： 补充土工试验的相关资料。

朱建朝： 调研及参与论文的修改。