基于Mohr-Coulomb准则和二阶锥规划技术的轴对称自适应下限有限元法

孙锐; 张箭; 阳军生; 杨峰

doi:10.11779/CJGE20220781

基于Mohr-Coulomb准则和二阶锥规划技术的轴对称自适应下限有限元法 English Version

孙锐^1,,
张箭^2, ,,
阳军生³,
杨峰³

1.
安徽理工大学土木建筑学院, 安徽淮南 232001
2.
河海大学岩土力学与堤坝工程教育部重点实验室, 江苏南京 210024
3.
中南大学土木工程学院, 湖南长沙 410075

基金项目:

安徽高校自然科学研究重点项目 2022AH050840

安徽理工大学高层次引进人才科研启动基金项目 2022yjrc31

国家自然科学基金面上基金项目 52178386

详细信息

作者简介:
孙锐(1993—)，男，博士，讲师，主要从事隧道与地下工程方面的研究工作。E-mail: sunruilight@163.com

通讯作者:
张箭, E-mail: zhangj0507@163.com

中图分类号: TU43
计量
- 文章访问数: 300
- HTML全文浏览量: 46
- PDF下载量: 80
出版历程
- 收稿日期: 2022-06-19
- 网络出版日期: 2023-11-05
- 刊出日期: 2023-10-31

Axisymmetric adaptive lower bound finite element method based on Mohr-Coulomb yield criterion and second-order cone programming

1.
School of Civil Engineering and Architecture, Anhui University of Science and Technology, Huainan 232001, China
2.
Key Laboratory for Geotechnical Engineering of Ministry of Water Resource, Hohai University, Nanjing 210024, China
3.
School of Civil Engineering, Central South University, Changsha 410075, China

摘要

摘要: 轴对称Mohr-Coulomb准则屈服面的角点问题导致其在数值计算中存在困难，如何高效处理该屈服准则一直是极限分析下限有限元法的重要内容。首先，引入完全塑性假定将轴对称Mohr-Coulomb准则转化为1组不等式约束和3个线性等式约束；然后，将不等式约束直接转化为二阶锥约束，避免了对角点进行光滑近似处理；最后，将基于Mohr-Coulomb准则的轴对称极限分析下限有限元计算模型转化为具有较高计算效率的二阶锥规划数学优化模型。极限分析下限有限元法采用的线性应力单元难以精确模拟破坏区域的应力变化，单元的分布形式对计算精度存在较大影响。因此提出一种基于单元应力的网格自适应加密策略，通过判断单元内节点应力接近屈服的程度，自动识别破坏区域待加密的单元，实现对破坏区域应力分布的精确模拟，进而能够以较少单元获得高精度下限解。通过分析圆形基础承载力及竖向锚板极限抗拔力等典型轴对称岩土工程稳定性问题，表明了所提方法具有较高计算效率及计算精度，具有一定的理论价值和应用前景。
- 下限有限元 /
- 自适应加密 /
- 轴对称 /
- 二阶锥规划 /
- Mohr-Coulomb准则
Abstract: The axisymmetric Mohr Coulomb yield surfaces have edges and corners in the three-dimensional stress space, which leads to difficulties in numerical calculation. Therefore, how to deal with the axisymmetric Mohr-Coulomb criterion efficiently has always been an important research content of the lower bound finite element limit analysis (LB-FELA) method. Firstly, the axisymmetric Mohr-Coulomb criterion is transformed into a set of inequality constraints and three linear equality constraints by introducing the complete plasticity assumption. Then, the inequality constraints are directly transformed into the second-order cone ones, which avoids the smooth approximation of numerical singularities. Finally, the axisymmetric LB-FELA model based on the Mohr Coulomb criterion is transformed into a mathematical optimization one of the second-order cone programming. The linear stress element adopted by the axisymmetric LB-FELA cannot accurately simulate the stress change in the failure region. Therefore, the mesh distribution form has a great impact on the calculation accuracy of the LB-FELA. To solve this problem, an adaptive mesh refinement strategy based on the element yield residual is proposed. By judging the degree of the node stress in the element close to the yield, the elements can be automatically identified and refined in the failure area. By analyzing the stability problems of typical axisymmetric geotechnical projects such as the bearing capacity of circular foundation and the ultimate uplift capacity of vertical anchor plate, it is shown that the proposed method has high calculation efficiency and accuracy, and it has certain theoretical value and application prospect.
- lower bound finite element method /
- adaptive refinement /
- axisymmety /
- second-order cone programming /
- Mohr-Coulomb criterion

HTML全文

0. 引言

非饱和渗透系数是非饱和土水力特性研究的重要参数之一^[1]，它是基质吸力（或体积含水量）的函数，变化范围可达数个数量级。目前，室内试验测定非饱和渗透系数，耗时费力，尤其是在高吸力阶段，测量周期较长。因此，学界建立了不同的模型来间接预测非饱和土的渗透系数。这些模型大致可分为经验模型^[2]、宏观模型^[3]和统计模型^[4-5]等。经验模型和宏观模型一般采用相对简单的数学公式对渗透系数进行预测，需要较多试验数据以获得最佳拟合参数，不便于复杂条件下的非饱和土水力特性理论分析。统计模型往往利用土-水特征曲线（SWCC）确定非饱和土渗透系数，是目前比较流行的方法，常见的基本统计模型有Childs Collis-George（CCG）模型^[4]，Mualem^[5]，Burdine模型^[6]等，许多学者对模型进行了修正，以便更准确地预测非饱和渗透系数。Van Genuchten（VG）^[7]建立了土-水特征曲线模型，并结合Burdine和Mualem提出的传导率模型建立了相对渗透系数方程。Fredlund等^[8]提出了一种新的SWCC模型，并结合CCG模型推导了整个吸力范围内的相对渗透系数函数方程。陶高梁等^[9-10]将SWCC视为反映土体不同尺寸孔隙渗透通道的间接指标，建立了考虑微观渗透通道起始水力梯度的饱和黏性土非线性渗透模型；基于渗透通道分形特性，建立了不同水力梯度下的饱和黏性土非线性渗透分形新模型。

除了基于SWCC模型的方法外，陈正汉等^[11]、姚志华等^[12]结合自主设计的原状黄土取样设备对原状和重塑黄土的渗水特性进行了一维水平土柱入渗试验，形成了成套的直接取样测量方法，该方法过程简明，原理清晰，所得结果较为可靠。

值得一提的是，直接测量方法一般试样尺寸较大，试验过程耗时相对较长，工作量较大；对于基于SWCC的间接预测方法，由于需要获得SWCC，同样需要较长的测量时间与较大的工作量。因此，提出一种可快速准确预测非饱和渗透系数的方法具有重要的实际意义和工程应用价值。

核磁共振（NMR）技术是一种灵敏度高、测试速度快、无破坏性的先进水分探测技术，核磁共振技术可应用于岩土工程领域中土料水分和孔隙变化的准确测量，具有测量时间快的优点，需要指出的是，NMR方法设备昂贵，相关成本较高。

目前利用NMR技术探究土壤中孔隙水分布和迁移的微观机理的研究已有一定基础^[13]。在多孔介质渗透特性预测方面，Seevers^[14]推导了基于弛豫时间原理的砂岩饱和渗透率方程。在此基础上，Kleinberg等^[15]提出了Schlumberger-Doll Research（SDR）方程，并采用平均弛豫时间作为预测岩土饱和渗透系数的参数。然而，现有的NMR试验经典模型主要通过连通孔隙的宏观总体积来描述岩石的渗透率，并未充分考虑核磁共振中横向弛豫时间分布几何形态等信息。综上所述，上述模型大多依赖于经验模型来推导渗透率函数，适用范围较小，并且预测对象都是饱和土，对于非饱和土渗透系数的预测较少，而本文的研究重点是通过NMR技术提出一种能够快速预测非饱和渗透系数的新方法。

结合NMR理论和渗流理论，本文提出了一种基于NMR曲线快速预测饱和/非饱和渗透系数的方法。为验证该方法的准确性，以湖南黏土为例，制备不同初始孔隙比土样，采用瞬时剖面法得到了非饱和渗透系数与体积含水量的关系；利用核磁共振技术测量了不同初始孔隙比试样脱湿及加湿过程的NMR曲线，分析了基于试样不同状态下NMR曲线的非饱和渗透系数预测结果，并与已有基于SWCC的预测模型预测结果进行了比较。

1. 基于NMR曲线的饱和/非饱和渗透系数新模型

假设土样由不同孔径的连通孔隙通道组成，土料中水的渗流运动皆发生于这些连通的孔隙通道中。饱和土的渗透系数即为海量连通孔隙通道饱和渗透系数的叠加值，可表示为^[16]

${k}_{\text{s}}=\frac{Q}{AJ}=\frac{{\displaystyle \sum _{i=1}^{n}\frac{\gamma {d}_{i}{}^{2}}{32\mu }\times {A}_{i}J}}{AJ}=\frac{{\displaystyle \sum _{i=1}^{{}_{n}}\frac{\gamma {d}_{i}{}^{2}}{32\mu }\times {A}_{i}}}{A}\text{ }。$

(1)

式中：Q为单位时间流量（cm³/s）；J为水力梯度（cm）；A为分析土样的横截面总面积（cm²）；μ为黏度（Pa·s）； $\gamma$ 为水的重度（N/m³）；d_i（cm）为第i级孔隙通道的直径；A_i（cm²）为第i级孔隙通道的横截面积。

假设第i级孔隙通道的实际长度与土样长l（cm）之比为p_i，则孔隙通道实际长度可记为p_il。若第i级孔隙通道的所含水分对应的体积含水量为θ_i，则相应的水分总体积为θ_iV_T，其中V_T（m³）为土样的总体积。则第i级孔隙通道截面积为

${A_i} = \frac{{{\theta _i}{V_{\text{T}}}}}{{{p_i}l}} 。$

(2)

将式（2）代入式（1）中可得

${k_{\text{s}}} = \sum\limits_{i = 1}^{_n} {\frac{{\gamma {d_i}^2}}{{32\mu }} \times \frac{{{\theta _i}{V_{\text{T}}}}}{{{p_i}lA}}} = \sum\limits_{i = 1}^{_n} {\frac{{\gamma {d_i}^2}}{{32\mu }} \times \frac{{{\theta _i}}}{{{p_i}}}} 。$

(3)

通过核磁共振（NMR）技术可获得土样的弛豫时间（T₂）分布曲线，该曲线反映了信号幅度和弛豫时间的对应关系。已有研究^[17]表明，弛豫时间T₂与土料孔隙特性直接相关：

$\frac{1}{{{T_2}}} \approx {\rho _2}\frac{S}{V} = \frac{{\lambda {\rho _2}}}{d} 。$

(4)

式中：S，V分别为水分所处孔隙的表面积（cm²）与体积（cm³）； ${\rho _2}$ 为横向弛豫率（Hz），其与土颗粒表面的物理化学性质有关；d为孔隙直径（cm），与式（3）中d_i都表示孔径， $\lambda$ 为与孔隙形状相关的因子，柱状孔隙为4，球形孔隙为6。

将式（4）简化后可得

$d = \lambda {\rho _2}{T_2} 。$

(5)

由上式可知，横向弛豫时间T₂与孔隙直径d成正比，因此弛豫时间T₂分布曲线实际上反应了孔隙分布特性。将T₂分布曲线的横坐标从右至左依次划分为n，n-1，n-2， $\cdots$ ，2级，1级，假设大小相近的第i级孔隙组成连续孔隙通道，其等效孔径为d_i（n级孔隙最大、1级孔隙最小），通过标定T₂分布曲线总信号幅度为土样总体积含水量，可得到该级通道对应的体积含水量θ_i，如图 1所示，图中编号表示了孔隙通道的分级方式。

图 1 孔隙通道分级示意图

Figure 1. Schematic diagram of classification of pore channels

下载: 全尺寸图片幻灯片

将式（5）代入式（3），得到

${k_{\text{s}}} = \sum\limits_{i = 1}^n {\frac{{\gamma {\lambda ^2}\rho _2^2{\theta _i}{T_{2i}}^2}}{{32\mu {p_i}}}} 。$

(6)

陶高梁等^[16]认为对于同一种土而言，p_i可假设为常数，于是上式可简化为

${k_{\text{s}}} = \sum\limits_{i = 1}^{_n} {{k_{\text{c}}}{\theta _i}T_{2i}^2} 。$

(7)

式中： ${k_{\text{c}}}$ = $\gamma {\lambda ^2}\rho _2^2$ /(32μ ${p_i}$ )，对于同一种土为常数。

假设总孔隙通道分为n级，现只有1~m级通道充满水（m＜n），此时试样总体积含水量为θ（θ＜θ_s，θ_s为试样饱和体积含水量），根据式（7）可得到此时非饱和渗透系数k(θ)：

$k(\theta ) = \sum\limits_{i = 1}^m {{k_{\text{c}}}{\theta _i}{T_{2i}}^2} 。$

(8)

结合式（7），（8），非饱和相对渗透系数k_r可以表示为

${k_{\text{r}}} = {{\sum\limits_{i = 1}^m {{\theta _i}{T_{2i}}^2} } \mathord{\left/ {\vphantom {{\sum\limits_{i = 1}^m {{\theta _i}{T_{2i}}^2} } {\sum\limits_{i = 1}^n {{\theta _i}} }}} \right. } {\sum\limits_{i = 1}^n {{\theta _i}} }}{T_{2i}}^2 。$

(9)

综上所述，基于NMR曲线，本文建立的饱和土渗透系数预测模型如（7）式所示，建立的非饱和土绝对渗透系数和相对渗透系数模型分别如式（8），（9）所示。

2. 湖南黏土饱和/非饱和渗透试验与NMR试验

2.1 饱和土渗透试验

本试验研究对象为湖南邵阳某地非饱和黏性土，属于原生黏土经过再次搬运堆积形成的次生黏土，土粒相对质量密度为2.76，液限w_L为46.3%，塑限w_p为27.8%。

将土料风干后碾碎，过2 mm圆孔筛，配置含水率一定的土料，选用直径为45 mm，高度为2 mm的环刀制样，制备初始孔隙比为1.12，1.04，0.97，0.90，0.84的试样，抽真空饱和，采用变水头法测量其饱和渗透系数。为提高试验精度，经多次重复试验取得平均值，进行温度修正后，得到20℃条件下不同初始孔隙比土样的饱和渗透系数^[18-19]，如表 1所示。

表 1 不同初始孔隙比试样的饱和渗透系数（20℃）

Table 1. Saturated permeability coefficients of soil samples with different initial void ratios (20℃)

初始孔隙比	1.12	1.04	0.97	0.90	0.84
渗透系数k_s/(cm·s^-1)	7.72×10⁻⁴	4.15×10⁻⁴	2.49×10⁻⁴	1.73×10⁻⁴	7.63×10⁻⁵

下载: 导出CSV

| 显示表格

2.2 非饱和土渗透试验

采用瞬时剖面法测量非饱和渗透系数。试验装置为课题组自主设计的有机玻璃桶，直径23 cm，高1 m，沿桶周均匀钻打五列竖向排列的圆形孔，孔与孔之间的间距为5 cm，直径为1 cm，并用A、B、C、D、E进行标记。试验开始前将这些小孔封闭以防漏水，本试验采用土仍为2.1节所述湖南黏土。在该试验的准备阶段，配置含水率为19%的土样待用，待水分迁移均匀后复测含水率，采用特制击实器分层击实，并在土柱顶部铺上8 cm厚的细砂层，随后从土柱顶部持续加水20 min，共计1500 mL，使水分在整个界面下渗均匀，观察湿润锋动态变化特点；分不同时间间隔于打好的孔中取土烘干以实测A、B、C、D、E列不同深度土样的含水率。结合所测试样的土-水特征曲线^[18]，计算得到非饱和渗透系数和体积含水量的关系^[18-19]。

将图 2中非饱和渗透系数实测值除以表 1中饱和渗透系数实测值，得到非饱和相对渗透系数实测值并用于后续验证。值得指出，瞬态剖面法试验中，上部土料以脱湿为主，下部土料以吸湿为主，研究表明土料脱湿与吸湿过程非饱和渗透系数存在一定差异，但相对于含水率的影响，该因素可以忽略不计^[8]。故一般情况下本文不考虑非饱和土渗透系数的滞回效应。

图 2 湖南红黏土非饱和渗透系数实测数据

Figure 2. Measured data of unsaturated permeability coefficient of Hunan clay

下载: 全尺寸图片幻灯片

2.3 核磁共振试验

（1）试验设备及原理

本文核磁共振试验于中科院武汉岩土所进行，所用仪器为苏州纽迈公司研发的PQ-001低磁场核磁共振仪。该试验过程中，试样置于试管中，其有效测试范围是60 mm×Φ60 mm。

核磁共振（NMR）的信号来源主要为氢质子，氢质子越多，表明含水率越大，反之则越低。因此通过信号量定标的方法，该技术可以用来测量物质中水的质量。NMR通过测定土样弛豫时间与信号强度的关系，计算峰面积，进而获得土样的水分尺度分布及孔隙分布特性。为了消除磁体磁场的非均匀性对T₂测量的影响，运用NMR测量分析软件及CPMG序列采集样品T₂衰减曲线，运用反演软件得到NMR试验数据。

（2）试验方案

本试验同样采用2.1节所述湖南黏土，制作两组平行试样，每组试样分别包含初始孔隙比为1.12，1.04，0.97，0.90，0.84的5个试样，用于脱湿与吸湿试验，分别模拟上层土料脱湿与下层土料吸湿过程。脱湿过程中，对饱和试样进行称重并记录，将不同初始孔隙比试样放在小铝盒中，自然风干。开始时，需每小时一次对试样进行称重，避免结果差异过大。土样脱水速度逐渐变慢后，可逐渐增大测量间隔。吸湿过程中，将初始与饱和含水率的间隔分为十级，每级之间的差相同以确定每级加湿目标含水率。由初始含水率开始，用滴水球分级向试样加水，密封72 h使水分均匀分布，最终使每个试样达到该初始含水率下的饱和含水率。利用核磁共振测试系统分别测定不同初始孔隙比土样在吸湿和脱湿过程中不同含水率试验的NMR曲线，NMR试验对应试样的含水率如表 2所示，NMR试验共计95次。

表 2 脱湿及吸湿过程NMR试验方案及其对应的试样含水率

Table 2. Moisture contents of samples corresponding to NMR tests of desorption and absorption process

孔隙比	脱湿过程					吸湿过程
孔隙比	1.12	1.04	0.97	0.90	0.84	1.12	1.04	0.97	0.90	0.84
体积含水量/%	53.16	52.35	50.39	49.18	48.36	19.85	20.61	21.38	22.14	22.91
	46.49	45.62	43.69	42.62	41.91	22.87	23.45	24.36	24.48	25.62
	40.89	40.16	38.44	36.74	36.53	29.17	26.43	27.15	27.49	28.07
	36.49	35.59	34.38	32.58	32.33	29.77	29.85	30.30	30.01	30.84
	36.15	35.21	34.08	32.23	31.98	32.75	32.85	33.28	32.70	33.65
	30.67	29.09	28.10	26.74	26.75	36.50	36.38	36.34	35.44	36.36
	29.69	28.05	27.13	25.69	25.62	40.57	39.62	39.40	38.21	39.23
	25.62	23.92	23.25	22.07	21.90	43.33	42.63	42.30	40.86	41.78
	14.35	13.64	14.36	13.17	14.13	46.72	45.56	45.26	43.50	44.52
	11.91	11.64	12.22	11.49	12.00

下载: 导出CSV

| 显示表格

（3）试验结果

a）不同初始孔隙比试样脱湿过程NMR曲线

按上节试验方案进行脱湿过程的核磁共振试验，获得不同初始孔隙比核磁共振试验反演后获得的核磁共振曲线，由于篇幅限制，本节选用e=1.12与e=0.84两种初始孔隙比结果作为代表，如图 3所示。

图 3 脱湿过程中核磁共振分布曲线试验数据

Figure 3. Experimental data of NMR distribution curves during desorption process

下载: 全尺寸图片幻灯片

由图 3可知，从饱和含水率到残余含水率过程中，曲线所围成的峰面积逐渐减小，表明含水率逐渐减少。不同初始孔隙比试样的NMR曲线变化规律存在相似之处，即均存在两个峰值，第一个峰值在0.5~0.8 ms处，第二个峰值在30~80 ms处。在脱湿过程中，随着含水率的降低，第二峰逐渐消失，说明在脱湿过程前期，大孔隙优先排水；而后第一峰值逐渐降低，表明在脱湿过程后期以小孔隙排水为主。

b）不同初始孔隙比试样吸湿过程NMR曲线

按上节试验方案进行吸湿过程的核磁共振试验，同样以e=1.12与e=0.84两种初始孔隙比反演后获得的核磁共振曲线为代表，如图 4所示。

图 4 吸湿过程中核磁共振分布曲线试验数据

Figure 4. Experimental data of NMR distribution curves in absorption process

下载: 全尺寸图片幻灯片

可以看出，与脱湿过程类似，不同初始孔隙比的NMR曲线所包围的峰面积随着含水率的增加而增加，不同初始孔隙比试样的NMR曲线同样有两个峰值，且第一个峰值的位置与脱湿过程大致相同，第二峰值点在20~30 ms。随着含水率的增加，第一峰值逐渐升高；然而不同初始孔隙比试样的NMR曲线在后五级吸湿过程中，其曲线峰值几乎不变。第二峰在前期吸湿过程中并未出现，直到含水率增加至较大含水率时，峰值略有起势，之后随含水率上升逐渐升高。这一试验结果表明，在吸湿过程中，水分优先进入小孔隙，由于小孔隙数量远多于大孔隙，因此小孔隙含水率总和大于大孔隙。NMR曲线表现为：第一峰面积明显大于第二峰面积，随着初始孔隙比的减小，大孔隙体积压缩量高于小孔隙，使得二者间差值增大。

c）不同初始孔隙比试样饱和NMR曲线

图 5为不同初始孔隙比下饱和试样的核磁共振反演曲线。

图 5 饱和状态试样NMR曲线

Figure 5. NMR curves under saturated state

下载: 全尺寸图片幻灯片

随着初始孔隙比的减小，NMR曲线所围成的总面积逐渐减小，即饱和土样的含水率减少。饱和试样的NMR曲线同样有两个峰值，大约分别位于0.8 ms和30~50 ms。随着初始孔隙比的减小，第二峰峰面积明显减小，而第一峰峰面积略有增加。上述现象表明，随着土样初始孔隙比的降低，大孔隙体积急剧减小，导致孔隙总体积减小，孔隙水质量随之减少。

3. NMR非饱和渗透系数模型验证

利用NMR试验预测渗透系数，关键是通过弛豫时间T₂分布曲线获得各孔隙通道对应的含水率。将不同初始孔隙比试样的总含水率除以各自的总信号幅度便可确定不同初始孔隙比试样中单位信号幅度对应的含水率, 随后便可将各弛豫时间 $T_{2i}^{}$ 对应的信号幅度转变为相应的体积含水量 ${\theta _i}$ 。

将中不同初始孔隙比试样NMR曲线的纵坐标分别求和，得到总信号幅度，随后将不同初始孔隙比试样的饱和含水率分别除以相对应的总信号幅度，得到单位信号幅度对应的含水率，信号幅度与单位信号幅度对应的含水率的乘积即为每个纵坐标对应的含水率，如此便将每级横向弛豫时间T_2i对应的的信号幅度转换为对应 ${\theta _i}$ 。基于式（7），将每级横向弛豫时间 $T_{2i}^{}$ 平方后乘以对应 ${\theta _i}$ 求和，便得到 $\sum\limits_{i = 1}^n {{\theta _i}} T_{2i}^2$ ，将其作为横坐标，表 1的饱和渗透系数实测值作为纵坐标进行线性拟合，设置截距为0，结果如图 6所示。

图 6 不同初始孔隙比湖南黏性土k_c线性拟合

Figure 6. Linear fitting of k_c of Hunan clay with different initial void ratios

下载: 全尺寸图片幻灯片

对于湖南黏土，其计算结果k_c=0.1046 cm/s³，均方根误差为7.151×10^-5，相关系数为0.9311，说明拟合结果良好，k_c对于同种类土样来说可视为定值。图 6证明了式（7）的合理性，式（9）是基于式（7）推导而来，因而证明了式（9）的理论基础是科学的。

利用脱湿、吸湿过程及饱和试样的NMR曲线，基于式（9）分别预测非饱和相对渗透系数。值得说明的是，本文为了准确刻画土料脱、吸湿路径不同含水率状态水分分布特性，共进行了95次核磁共振试验，在利用式（9）预测非饱和相对渗透系数时，计算方法有3种，即：吸湿NMR曲线预测方法、脱湿NMR曲线预测方法和饱和NMR曲线预测方法。3种方法的具体预测计算过程分别如下：

（a）对于脱湿、吸湿核磁共振方法而言，分母部分 $\sum\limits_{i = 1}^n {{\theta _i}} T_{2i}^2$ 按照土料饱和状态NMR曲线（）进行计算，采取的方法与中横坐标的获取方式一致，分子部分 $\sum\limits_{i = 1}^m {{\theta _i}} T_{2i}^2$ 的计算方法与分母相同，只不过需要分别选取非饱和相对渗透系数对应含水率的NMR曲线（图 3与）计算得到的 ${\theta _i}$ 和 $T_{2i}^{}$ 。

（b）对于饱和NMR曲线预测方法，同样是利用式（9）进行计算，分母 $\sum\limits_{i = 1}^n {{\theta _i}} T_{2i}^2$ 仍然按照饱和状态的NMR曲线进行同样的计算，不同的是，在计算分子 $\sum\limits_{i = 1}^m {{\theta _i}} T_{2i}^2$ 时，假设土料排水为大孔先排出，再排小孔的水（图 3脱湿NMR曲线也证明了这一现象），因此在计算某一特征含水率θ对应的非饱和相对渗透系数时，以饱和状态NMR曲线为基础，如图 1所示，根据NMR数据所得的T₂弛豫时间分布将孔隙通道分为1~n级，若θ值与第1~i级孔隙的累计体积含水量相等（或处于第1~（i-1）级孔隙的累计体积含水量与第1~i级孔隙的累计体积含水量之间），则分子 $\sum\limits_{i = 1}^m {{\theta _i}} T_{2i}^2$ 计算时只计算饱和NMR曲线1~i级孔隙的对应数值，如此便得到利用饱和NMR曲线的不同特征含水率下的非饱和相对渗透系数预测值。

将脱湿、吸湿及饱和状态NMR曲线的非饱和相对渗透系数式（9）预测值和实测值进行比较，见图 7。

图 7 式（9）非饱和渗透系数预测值与实测值对比

Figure 7. Comparison between predicted values from three kinds of NMR curves and experimental results

下载: 全尺寸图片幻灯片

由图 7可知，基于脱湿、吸湿过程及饱和状态NMR预测值整体上与实测值基本吻合，吸湿和脱湿过程NMR曲线预测值整体低于实测值。其中，当含水率大于35%时，吸湿过程NMR曲线预测误差最大，其预测值低于实测值、脱湿过程NMR曲线预测值和饱和NMR曲线预测值；当含水率达到35%左右时，脱湿过程NMR曲线预测值要低于实测值和饱和状态NMR曲线预测值，并与吸湿过程NMR曲线预测值相近，与实测值相差较远，此时饱和状态NMR曲线预测值更接近实测值。整体上看，饱和状态NMR曲线预测值误差最小，其预测值与实测值总体较为吻合，尤其在初始孔隙比为1.04时基本重合。

在含水率大于35%时，瞬态剖面法中土料一般处于上部土体脱湿过程中，因此脱湿过程NMR预测方法所得误差较小。在吸湿过程中，滴水球滴入试样中的水分不能完全被试样吸收，部分水分会在空气中挥发；而在密封迁移水分的过程中，部分水分会残留这密封膜上，导致试样的实际含水率偏小。在脱湿过程中，随着含水率的降低，试样慢慢收缩，逐渐与环刀脱离，试验过程中需要对试样进行频繁称重，导致部分土颗粒脱落，形成质量误差。对于饱和状态NMR曲线预测方法而言，该方法并没有上述试验操作过程，因此避免了试验操作带来的各项误差，对于其深层次的机理原因，有待进一步深入研究与分析。

4. 分析与讨论

4.1 核磁共振曲线缺级的影响

由2.3节（3）中的T₂分布曲线可知，较为密实的饱和试样及低含水率的非饱和试样容易出现孔隙水含量缺级，如图 8所示，连续NMR曲线出现中断的现象被称为缺级。

图 8 核磁共振曲线缺级示意图

Figure 8. Schematic diagram of missing grade of NMR curves

下载: 全尺寸图片幻灯片

此时，大孔隙水分总含量较少，大孔隙间很难形成直接的联通孔隙，需借道小孔隙形成联通孔隙通道，因此连通孔隙通道的等效孔径由小孔隙孔径控制。本文将这部分水分所对应的孔隙通道称为“非连续孔隙通道”，其对土样的渗透系数不起控制作用。若按照式（7）~（9）计算，会产生较大误差，应舍弃这些数据。图 9分别分析了在上述e=1.12试样中，脱湿NMR曲线和吸湿NMR曲线预测非饱和相对渗透系数时，NMR曲线缺级的影响。

图 9 NMR曲线缺级对脱湿和吸湿过程非饱和渗透系数预测的影响

Figure 9. Effects of lack of level in NMR curves on predicted results during desorption and absorption processes

下载: 全尺寸图片幻灯片

可以明显看出，若不舍弃非连续孔隙通道（NMR曲线中缺级部分）的数据，预测值将非常不稳定，并且从10^-5至10²数量级之间变化，总体上远远偏离实测值；当舍弃非连续孔隙通道数据时，预测结果均处于10^-4~1数量级区间，与实测值吻合较好。因此在模型验证部分，非连续孔隙通道数据，本文均予以舍弃。

4.2 脱湿、吸湿及饱和状态NMR曲线预测方法的比较

非饱和渗透系数可通过脱湿、吸湿和饱和试样的NMR曲线来预测。虽然这3种情况下的预测值都与试验值吻合较好，但整体上来看基于饱和状态下的NMR曲线总体预测效果更好。另一方面，采用脱湿或吸湿过程的NMR曲线对某一含水率条件下的非饱和渗透系数进行预测时，需要先对相应含水率试样进行NMR试验测量相应NMR曲线，预测不同含水率下的非饱和渗透系数则需要进行多次试验，试验周期长，工作量大，成本高；而饱和NMR曲线预测方法仅需要一次核磁共振试验便可预测全含水率区间的所有非饱和相对渗透系数，其工作量和成本远低于其他两种方法。因此，本文建议采用饱和NMR曲线直接预测饱和/非饱和渗透系数。

4.3 NMR曲线预测方法与已有方法对比

目前为止，基于土-水特征曲线（SWCC）预测非饱和渗透系数是最常用的方法，然而，SWCC本身测量周期较长，获取全吸力范围的试验数据较为困难，不利于实际工程应用。值得一提的是，核磁共振（NMR）方法具有快速方便的优点，值得进一步研究。前文已对NMR方法的有效性进行了论证，其中饱和NMR曲线预测方法是最值得推广的方法。为进一步比较SWCC预测方法与NMR预测方法，本文选取了两种SWCC预测方法：CCG（Childs Collis-George）模型和文献[16]模型，对湖南黏土非饱和相对渗透系数进行预测。图 10给出了饱和NMR预测法、CCG模型预测法和文献[16]模型预测法结果与非饱和相对渗透系数实测值。

图 10 NMR预测方法与已有SWCC预测方法对比

Figure 10. Comparison of prediction methods based on NMR curves and SWCCs

下载: 全尺寸图片幻灯片

由图 10可知，CCG模型预测值均高于实测值，且预测误差较大，土样含水率越小，预测误差越大，当含水率接近30%时，预测值偏离较大，可见CCG模型对湖南黏土适用性不佳；在含水率高于35%时，文献[16]与NMR预测方法预测结果接近，且与实测值吻合较好，而在体积含水量低于35%时，文献[16]模型预测值远低于实测值，而NMR预测方法仍具有较好预测效果。

本文提出的饱和NMR曲线预测方法能够快速与较为准确地预测湖南黏土的非饱和相对渗透系数，但对于其他土类的适用性仍有待验证。同时，对存在缺级现象的土料，直接应用时误差可能较大，应舍弃缺级部分。值得一提的是，在应用式（9）预测非饱和相对渗透系数时不需要获得k_c，若需要应用式（7）或式（8）获得非饱和/饱和（绝对）渗透系数，需要对不同初始孔隙比的试样进行多次饱和渗透试验与核磁共振试验以确定k_c，增加了试验的工作量。此外，由图 10可以看出，SWCC预测方法（CCG、文献[16]模型）无法给出低含水率区间的非饱和相对渗透系数预测结果，其原因在于，SWCC通常由压力板仪测得，一般无法给出高基质吸力（即低含水率）区间的实测结果，而饱和NMR预测方法可以同时给出全含水率范围内的渗透系数预测结果。尽管如此，饱和NMR预测方

法在低含水率区间的预测效果仍需进一步验证。

5. 结论

结合核磁共振（NMR）理论和渗流理论建立了基于核磁共振曲线的饱和/非饱和渗透系数预测模型。为验证该方法的合理性，以湖南黏土为研究对象，进行了渗透试验及核磁共振试验，利用试验数据验证了该模型的合理性，得到以下4点结论。

（1）首次提出了通过核磁共振曲线直接预测饱和/非饱和渗透系数的模型（式（7）~（9）），并结合湖南黏土的饱和/非饱和渗透系数实测值进行了验证，证明该模型可靠有效。

（2）脱湿、吸湿及饱和NMR曲线均能较好预测湖南黏性土非饱和相对渗透系数，其中饱和NMR曲线预测效果最佳。此外，饱和NMR曲线预测方法具有成本低、工作量小、速度更快等优点，本文建议采用饱和NMR曲线预测方法预测非饱和渗透系数。

（3）对于缺级的NMR曲线，非连续孔隙通道数据对预测结果具有较大影响，计算时应舍弃NMR曲线缺级之后的大孔隙试验数据，即“非连续孔隙通道”数据。

（4）采用CCG模型、文献[16]模型两种基于土-水特征曲线的预测方法对湖南黏性土非饱和渗透系数进行了预测，并与本文提出的饱和NMR曲线预测方法的预测结果、湖南黏土非饱和渗透系数实测结果进行对比，发现本文提出的饱和NMR法的预测结果与实测值吻合更好。

图 1 单元形式及应力间断

Figure 1. Element type and stress discontinuity

下载: 全尺寸图片幻灯片

图 2 基础初始网格及边界条件

Figure 2. Initial mesh and boundary conditions for footing

下载: 全尺寸图片幻灯片

图 3 圆形基础地基承载力下限有限元加密网格(接触光滑)

Figure 3. Adaptive refined meshes of bearing capacity using lower bound finite element method with smooth interface

下载: 全尺寸图片幻灯片

图 4 圆形基础地基承载力下限有限元加密网格(接触粗糙)

Figure 4. Adaptive refined meshes of bearing capacity using lower bound finite element method with rough interface

下载: 全尺寸图片幻灯片

图 5 网格自适应加密次数对地基承载力系数的影响

Figure 5. Effects of mesh refinement times on bearing capacity coefficient

下载: 全尺寸图片幻灯片

图 6 圆形锚板竖向抗拔力问题初始网格及边界条件

Figure 6. Initial mesh and boundary conditions for vertical uplift of circular anchor plate

下载: 全尺寸图片幻灯片

图 7 圆形锚板竖向极限抗拔力问题加密网格（局部破坏区域）

Figure 7. Refinement meshes for vertical uplift of circular anchor plates

下载: 全尺寸图片幻灯片

${N_q}$ 计算结果与文献结果对比

Figure 8. Comparison between present results and those available in Reference[1]

下载: 全尺寸图片幻灯片

表 1 计算效率对比(地基与基础直径接触粗糙)

Table 1 Comparison of calculation efficiencies between present method and those available in Reference [18]

方法	单元数	节点数	优化变量数	锥约束数量	线性约束数量	$\phi {\text{ = }}5^\circ$			$\phi {\text{ = }}20^\circ$
方法	单元数	节点数	优化变量数	锥约束数量	线性约束数量	承载力系数	计算误差/%	计算时长/s	承载力系数	计算误差/%	计算时长/s
本文方法	22846	11592	479766	68538	593794	7.73	3.37	31	22.35	3.74	37
文献[18]	22846	11592	890994	205614	799408	7.74	3.25	68	22.37	3.67	87

下载: 导出CSV

表 2 N_c计算结果与文献结果对比

Table 2 Comparison between present results and those available in references

$\phi$	本文方法	Kumar等^[5]	Zhang等^[25]	Martin等^[8]
5	8.03(7.41)	8.00(7.31)	8.00(7.47)	8.06(7.43)
10	11.04(9.94)	10.99(9.78)	11.16(10.05)	11.09(9.99)
15	15.76(13.86)	15.66(13.51)	16.03(13.99)	15.84(13.87)
20	23.53(19.94)	23.22(19.38)	24.02(20.32)	23.67(20.07)
25	36.98(30.25)	36.17(29.06)	37.96(30.97)	37.31(30.52)
30	61.92(48.78)	61.48(47.10)	64.27(50.48)	62.70(49.29)
35	112.07(84.59)	112.47(81.47)	117.99(89.01)	113.99(85.88)
注：括号外为接触粗糙时地基承载力系数，括号内为接触光滑时的承载力系数；自适应加密次数均为15次。

下载: 导出CSV

参考文献(28)

[1]	KHATRI V N, KUMAR J. Vertical uplift resistance of circular plate anchors in clays under undrained condition[J]. Computers and Geotechnics, 2009, 36(8): 1352-1359. doi: 10.1016/j.compgeo.2009.06.008
[2]	KHATRI V N, KUMAR J. Bearing capacity factor Nc under ϕ=0 condition for piles in clays[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2009, 33(9): 1203-1225. doi: 10.1002/nag.763
[3]	PASTOR J, TURGEMAN S. Limit analysis in axisymmetrical problems: numerical determination of complete statical solutions[J]. International Journal of Mechanical Sciences, 1982, 24(2): 95-117. doi: 10.1016/0020-7403(82)90041-8
[4]	COX A D. Axially-symmetric plastic deformation in soils—Ⅱ. Indentation of ponderable soils[J]. International Journal of Mechanical Sciences, 1962, 4(5): 371-380. doi: 10.1016/S0020-7403(62)80024-1
[5]	KUMAR J, KHATRI V N. Bearing capacity factors of circular foundations for a general c-ϕ soil using lower bound finite elements limit analysis[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2011, 35(3): 393-405. doi: 10.1002/nag.900
[6]	KUMAR J, CHAKRABORTY D. Stability numbers for an unsupported vertical circular excavation in c-ϕ soil[J]. Computers and Geotechnics, 2012, 39: 79-84. doi: 10.1016/j.compgeo.2011.08.002
[7]	KUMAR J, CHAKRABORTY M. Upper-bound axisymmetric limit analysis using the mohr-coulomb yield criterion, finite elements, and linear optimization[J]. Journal of Engineering Mechanics, 2014, 140(12): 06014012. doi: 10.1061/(ASCE)EM.1943-7889.0000820
[8]	MARTIN C. M. User guide for ABC: Analysis of bearing capacity, Version 1.0[R]. Oxford: Department of Engineering Science, University of Oxford, 2004.
[9]	孙锐, 阳军生, 赵乙丁, 等. 基于Drucker-Prager准则的高阶单元自适应上限有限元研究[J]. 岩土工程学报, 2020, 42(2): 398-404. doi: 10.11779/CJGE202002022 SUN Rui, YANG Junsheng, ZHAO Yiding, et al. Upper bound adaptive finite element method with higher-order element based on Drucker-Prager yield criterion[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(2): 398-404. (in Chinese) doi: 10.11779/CJGE202002022
[10]	孙锐, 杨峰, 阳军生, 等. 基于二阶锥规划与高阶单元的自适应上限有限元研究[J]. 岩土力学, 2020, 41(2): 687-694. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX202002040.htm SUN Rui, YANG Feng, YANG Junsheng, et al. Investigation of upper bound adaptive finite element method based on second-order cone programming and higher-order element[J]. Rock and Soil Mechanics, 2020, 41(2): 687-694. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX202002040.htm
[11]	MAKRODIMOPOULOS A, MARTIN C M. Lower bound limit analysis of cohesive-frictional materials using second-order cone programming[J]. International Journal for Numerical Methods in Engineering, 2006, 66(4): 604-634. doi: 10.1002/nme.1567
[12]	MAKRODIMOPOULOS A, MARTIN C M. Upper bound limit analysis using simplex strain elements and second-order cone programming[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2007, 31(6): 835-865. doi: 10.1002/nag.567
[13]	杨昕光, 周密, 张伟, 等. 基于二阶锥规划的边坡稳定上限有限元分析[J]. 长江科学院院报, 2016, 33(12): 61-67. https://www.cnki.com.cn/Article/CJFDTOTAL-CJKB201612013.htm YANG Xinguang, ZHOU Mi, ZHANG Wei, et al. Upper bound finite element limit analysis of slope stability using second-order cone programming[J]. Journal of Yangtze River Scientific Research Institute, 2016, 33(12): 61-67. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-CJKB201612013.htm
[14]	杨昕光, 迟世春. 土石坝坡极限抗震能力的下限有限元法[J]. 岩土工程学报, 2013, 35(7): 1202-1209. http://www.cgejournal.com/cn/article/id/15099 YANG Xinguang, CHI Shichun. Lower bound FEM for limit aseismic capability of earth-rockfill dams[J]. Chinese Journal of Geotechnical Engineering, 2013, 35(7): 1202-1209. (in Chinese) http://www.cgejournal.com/cn/article/id/15099
[15]	刘锋涛, 张绍发, 戴北冰, 等. 边坡稳定分析刚体有限元上限法的锥规划模型[J]. 岩土力学, 2019, 40(10): 4084-4091, 4100. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201910044.htm LIU Fengtao, ZHANG Shaofa, DAI Beibing, et al. Upper bound limit analysis of soil slopes based on rigid finite element method and second-order cone programming[J]. Rock and Soil Mechanics, 2019, 40(10): 4084-4091, 4100. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201910044.htm
[16]	周锡文, 刘锋涛, 戴北冰, 等. 基于混合常应力-光滑应变单元的极限分析方法[J]. 岩土力学, 2022, 43(6): 1660-1670. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX202206019.htm ZHOU Xiwen, LIU Fengtao, DAI Beibing, et al. Limit analysis method based on mixed constant stress-smoothed strain element[J]. Rock and Soil Mechanics, 2022, 43(6): 1660-1670. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX202206019.htm
[17]	杨昕光, 迟世春. 基于非线性破坏准则的土坡稳定有限元上限分析[J]. 岩土工程学报, 2013, 35(9): 1759-1764. http://www.cgejournal.com/cn/article/id/15293 YANG Xinguang, CHI Shichun. Upper bound FEM analysis of slope stability using a nonlinear failure criterion[J]. Chinese Journal of Geotechnical Engineering, 2013, 35(9): 1759-1764. (in Chinese) http://www.cgejournal.com/cn/article/id/15293
[18]	TANG C, TOH K C, PHOON K K. Axisymmetric lower-bound limit analysis using finite elements and second-order cone programming[J]. Journal of Engineering Mechanics, 2014, 140(2): 268-278.
[19]	LYAMIN A V, SLOAN S W, KRABBENHØFT K, et al. Lower bound limit analysis with adaptive remeshing[J]. International Journal for Numerical Methods in Engineering, 2005, 63(14): 1961-1974.
[20]	ZHANG R, CHEN G H, ZOU J F, et al. Study on roof collapse of deep circular cavities in jointed rock masses using adaptive finite element limit analysis[J]. Computers and Geotechnics, 2019, 111: 42-55.
[21]	李大钟, 郑榕明, 王金安, 等. 自适应有限元极限分析及岩土工程中的应用[J]. 岩土工程学报, 2013, 35(5): 922-929. http://www.cgejournal.com/cn/article/id/15048 LI Dazhong, CHENG Yungming, WANG Jinan, et al. Application of finite-element-based limit analysis with mesh adaptation in geotechnical engineering[J]. Chinese Journal of Geotechnical Engineering, 2013, 35(5): 922-929. (in Chinese) http://www.cgejournal.com/cn/article/id/15048
[22]	ZHANG R, LI L A, ZHAO L H, et al. An adaptive remeshing procedure for discontinuous finite element limit analysis[J]. International Journal for Numerical Methods in Engineering, 2018: 287-307.
[23]	MUÑOZ J J, BONET J, HUERTA A, et al. Upper and lower bounds in limit analysis: adaptive meshing strategies and discontinuous loading[J]. International Journal for Numerical Methods in Engineering, 2009, 77(4): 471-501.
[24]	SUN R, YANG J S, ZHAO Y D, et al. Upper bound finite element limit analysis method with discontinuous quadratic displacement fields and remeshing in non-homogeneous clays[J]. Archive of Applied Mechanics, 2021, 91(3): 1007-1020.
[25]	ZHANG J, GAO Y F, FENG T G, et al. Upper-bound finite-element analysis of axisymmetric problems using a mesh adaptive strategy[J]. Computers and Geotechnics, 2018, 102: 148-154.
[26]	LIU F Q, WANG J H, ZHANG L L. Analytical solution of general axisymmetric active earth pressure[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2009, 33(4): 551-565.
[27]	LIU F Q, WANG J H. A generalized slip line solution to the active earth pressure on circular retaining walls[J]. Computers and Geotechnics, 2008, 35(2): 155-164.
[28]	MOSEK ApS. The MOSEK C optimizer API manual, version10.1[EB/OL]. 2023-10-01. https://www.mosek.com/.