基于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
计量
- 文章访问数: 301
- 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-9]。其中大部分模型描述的都是土水特征曲线的边界曲线,即主干燥线和主湿化线,但实际工程中土体经历干湿循环的过程大多位于扫描线,由于理论上存在无数条扫描线,很难得知实际条件下扫描线的变化,目前针对扫描线的预测模型主要有域模型^[10]、边界面模型^[11-14]和接触角模型^[14-17]。以上模型可以很好地预测土水特征曲线的扫描线并反映其滞后效应,但是这些模型都是在平衡条件下建立,并没有考虑非饱和土的渗流过程中由于流速变化引起的动态效应。动态效应与滞后效应均与水力历史相关,区别在于：滞后效应与饱和度改变的方向有关,而动态效应与饱和度改变的速率有关。

受到气候环境影响,实际工程中土的含水状态将不断改变,在孔隙水瞬态流动的过程中,土的持水特征会表现出明显的流速相关的率效应,土水特征曲线和水力参数都与流速有关^[18-20]。大量试验已经证实^[21-22],土水特征曲线存在动态效应,与静态土水特征曲线相比,相同饱和度时,动态土水特征曲线的吸力在脱湿过程中要更大,而吸湿过程中则更低。对于产生动态效应的原因有很多不同的解释,如气液接触面的动态运动^[23]、接触角的动态变化^[24]、土体的不均性^[25]等。不论何种原因,这种动态效应必将对土体的行为产生影响。以降雨入渗为例,不同降雨强度下的土水特征曲线并不相同,此时如果仍采用平衡条件下土水特征曲线来建立渗流方程,就会产生误差。研究发现基于平衡条件的土水特征曲线求解的Richard方程无法准确预测动态多步排水过程中的水头变化^[26]。

为了更精确地描述非平衡过程中土的力学行为,需要建立动态条件下的土水特征曲线模型。目前的模型主要有两类^[19]：一是孔隙尺度的网格模型^[27],需要考虑微观的交界面,并借助统计学理论和电镜扫描等试验手段来进行计算,较为复杂；二是宏观尺度的动态模型,主要基于试验现象和经验方程来建立,便于应用,本文将主要基于宏观模型来建立滞回模型。这类模型主要有两种不同的形式,一种是给出动态毛细力的方程^[23],另一种是给出动态含水率的方程^[28]。尽管目前已有关于动态土水特征曲线的试验研究和理论模型,但这些模型只能描述主干燥和主湿化线,无法预测任意扫描线的变化。

建立动态土水特征曲线滞后模型具有重要的现实意义,实际过程中非饱和土力学的含水率大多处于动态变化的非平衡状态,此时采用平衡条件下得到的土水特征曲线很难准确地描述其力学行为。本文将利用边界面模型的思想,利用动态边界曲线方程来建立动态扫描线模型。首先对动态土水特征曲线的理论基础进行简单的论述,然后建立动态滞回模型,最后将利用已有试验结果对模型进行验证。

1. 动态土水特征曲线的理论基础

赵成刚等ADDIN NE.Ref.{BC7CB835-6C0E-453C-BC0E-60F1A5D71D7A}^[29]以热力学基本平衡方程出发,在不考虑温度影响和质量交换,忽略固体颗粒和流体压缩性的假设下,推导得到了单位体积的非饱和土变形功的基本表达式。在此基础上进一步考虑交界面影响,刘艳等^[30]给出了非饱和土流体的熵增不等式：

$[n (p_{g} - p_{w}) + γ^{wg} \frac{\partial a^{wg}}{\partial S_{r}}] d S_{r} - d φ_{w} \geq 0$ 。

(1)

式中 n为孔隙率； $p_{g}$ 和 $p_{w}$ 分别为孔隙气压和孔隙水压； $S_{r}$ 为饱和度； $γ^{wg}$ 为表面张力； $a^{wg}$ 为气液交界面面积； $φ_{w}$ 为水的自由能。式中第一项为孔隙流体的压力差,第二项为宏观毛细力,用 $p_{c}$ 表示。

Hassanizadeh等^[31]从热力学角度探讨了非饱和土的毛细力,指出毛细力本质上是气液交界面上的一种作用力,并不完全等同于界面两侧流体的压力差,并可将宏观毛细力表示为

$p_{c} = - \frac{γ^{wg}}{n} \frac{\partial a^{wg}}{\partial S_{r}}$ 。

(2)

在线性假设条件下,根据熵不等式（1）可以得到液相广义力与广义流之间本构关系如下：

$n [(p_{g} - p_{w}) - p_{c}] = - τ^{'} {\dot{S}}_{r}$ 。

(3)

式（3）与Hassanizadeh等^[31]给出的表达式是一致的, $τ^{'}$ 称为毛细阻尼系数（capillary damping coefficient）,它控制了饱和度改变的速率。如果 $τ^{'}$ 取值很小,说明系统受扰动后会迅速地恢复新的平衡。

将式（3）改写一下,可以得到动态土水特征曲线的方程：

$p_{c}^{d} = p_{c} - τ {\dot{S}}_{r}$ 。

(4)

式中, $τ = τ^{'} / n$ , $p_{c}^{d} = p_{g} - p_{w}$ 表示动态毛细力,等于流体的压力差, $p_{c}$ 为平衡时的毛细力,可以用传统的土水特征曲线方程来计算。等式右边第二项的存在说明此时动态毛细力与饱和度的变化历史有关,已有试验证明这种水力历史会对土的变形和强度产生影响。由于流体压力差( $p_{g} - p_{w}$ )与毛细力 $p_{c}$ 的不平衡,导致了系统饱和度的改变,使系统向着新的平衡状态发展,最终恢复流体压力差( $p_{g} - p_{w}$ )与毛细力 $p_{c}$ 的平衡。只有在平衡条件下,饱和度不随时间变化,右边第二项等于0,式（4）退化为传统的土水特征曲线方程。

2. 动态土水特征曲线模型

利用热力学理论推导得到的式（4）,可用于计算动态的土水特征曲线,由于平衡条件下土水特征曲线的边界线可以利用经验方程给出（如VG模型等）,因此式（4）可直接计算出动态边界线,而扫描线的计算则需要借助于边界面模型。

2.1 边界面模型的理论基础

实际工程中,土体受气候和环境影响经历干湿循环,导致其持水特征曲线往往处于扫描线上,一般通过室内试验比较容易确定边界线的方程,而扫描线与水力历史密切相关,难以确定。按照边界面塑性理论,加载面上的塑性反应取决于加载面上的应力点与其在边界面上的映射点之间的距离。只要知道了边界面的塑性模量,根据塑性边界面理论就可得出其中任意一条扫描线的塑性模量。

如图1所示,两条湿化扫描线BC和B'C',尽管它们起点和终点的饱和度相同,但吸力的变化是不一样的。原因就在于由于B'所在的动态干燥曲线与B所在的静态干燥曲线不重合,导致湿化扫描线也不重合。对于同一个点A,在静态和动态条件下可能处于不同的扫描线上。根据Wei等^[12],处于静态的土水特征曲线,其扫描线的斜率与O,A,S三点的距离的比值有关,基于这个思想,刘艳等ADDIN NE.Ref.{30D53DDE-64E0-477C-8C1C-91BCA8348D6B}^[13]给出了静态扫描线斜率的计算方法如下：

图 1 土水特征曲线边界面模型示意图

Figure 1. Schematic bounding surface model for SWCC

下载: 全尺寸图片幻灯片

$K (P_{c}, S_{r}) = \bar{K} (S_{r}) + c \frac{ϕ}{r}$ ,

(5)

式中,c为材料参数, $\bar{K} (S_{r})$ 为土水特征曲线边界线的斜率, $K$ 为扫描线的斜率。r为当前饱和度下两边界曲线的吸力之差（如图1中的OS）,可以表示为

$r = κ_{D} (S_{r}) - κ_{W} (S_{r})$ 。

(6)

式中, $κ_{α} (S_{r})$ 表示静态条件下边界土水特征曲线方程,下标 $α$ =D代表干燥过程,下标 $α$ =W代表湿化过程。

φ表示当前吸力与所对应的边界曲线的吸力之差,如果是从B至C湿化时为图1中的AS段,如果是从C到B的干燥则是图1中OA段,因此可以表示为

$ϕ = | p_{c} - κ_{α} (S_{r}) |$ 。

(7)

式（7）中 $p_{c}$ 代表当前扫描线上的吸力。

边界线斜率可对边界曲线方程求导得到,即

${\bar{K}}_{α} (S_{r}) = - d κ_{α} (S_{r}) / d S_{r}$ 。

(8)

把式（6）～（8）代入到式（5）中则可得扫描线斜率为

$K_{α} = {\bar{K}}_{α} (S_{r}) + c \frac{| p_{c} - κ_{α} (S_{r}) |}{κ_{D} (S_{r}) - κ_{W} (S_{r})}$ 。

(9)

利用式（9）,给定初始条件,就可以利用

${\dot{p}}_{c} = - K_{α} {\dot{S}}_{r}$

(10)

给出扫描线。式（9）预测扫描线的合理性已在文献[13]中得到了验证。以这个模型为基础,可以建立动态条件土水特征曲线的滞回模型。

2.2 动态条件下边界面模型的建立

在瞬态渗流过程中,土体的土水特征曲线与流速相关,此时可以利用动态的边界曲线来预测动态的扫描线。因此计算原理与静态曲线相同,仍然可以采用式（5）的方法来计算,但此时 $\bar{K} (S_{r})$ 需采用动态土水特征曲线边界线的斜率,r和φ的取值也将发生变化。

根据式（4）可知,动态边界曲线与静态条件下的边界曲线不再相同,分别为

$p_{c}^{d} = κ_{W}^{} (S_{r}) - τ_{W} {\dot{S}}_{rW}^{} (边界湿化曲线)$ ,

(11)

$p_{c}^{d} = κ_{D} (S_{r}) - τ_{D} {\dot{S}}_{rD}^{} (边界干燥曲线)$ 。

(12)

预测扫描线时,一般要采用与其相同速率的边界线,如后文图2和图3所示。相同饱和度时,两边界线的吸力之差 $r^{d}$ （即图1中O'S'）可以表示为

$r^{d} = κ_{D} (S_{r}) - κ_{W} (S_{r}) - (τ_{D} {\dot{S}}_{rD}^{} - τ_{W} {\dot{S}}_{rW}^{})$ 。

(13)

图 2 干燥条件下动态扫描线预测

Figure 2. Prediction of dynamic scanning curve during drying process

下载: 全尺寸图片幻灯片

图 3 湿化条件下动态扫描线预测

Figure 3. Prediction of dynamic scanning curve during wetting process

下载: 全尺寸图片幻灯片

对于当前吸力与所对应的边界曲线的吸力之差 $ϕ^{d}$ ,如果是从B'至C'湿化时为图1中的AS'段,如果是从C'到B'的干燥则是图1中O'A段,可以统一表示为

$ϕ^{d} = | p_{c} - [κ_{α}^{} (S_{r}) - τ_{α} {\dot{S}}_{r α}^{}] |$ 。

(14)

把式（13）,（14）代入到式（5）中则可得动态扫描线斜率为

$K_{α}^{d} = {\bar{K}}_{α}^{d} (S_{r}) + c \frac{ϕ^{d}}{r^{d}}$ 。

(15)

对式（11）或（12）求导,假设饱和度变化率不随时间变化,可得到边界线的斜率变为

${\bar{K}}_{α}^{d} (S_{r}, {\dot{S}}_{r}) = - \frac{d κ_{α} (S_{r})}{d S_{r}} + \frac{d τ_{α}}{d S_{r}} {\dot{S}}_{r}$ 。

(16)

将式（8）,（16）代入式（15）,可得

$K_{α}^{d} = {\bar{K}}_{α} (S_{r}) + c \frac{ϕ^{d}}{r^{d}} + \frac{d τ_{α}}{d S_{r}} {\dot{S}}_{r α}$ 。

(17)

与静态条件下式（9）相比,式（17）中r和φ的取值有所改变,同时还多了一个第三项,该项与饱和度变化率有关,同时也与参数τ有关。现有研究发现τ并不是一个常数,而是会随着饱和度变化而改变^[32],因此第三项不能忽略。通常可将τ表示为饱和度的函数,目前主要做法是根据试验结果给出一个经验函数,函数形式以线性和指数型为主（如下文的式（20）和（21））。

结合式（17）给定初始条件,利用

${\dot{p}}_{c}^{d} = - K_{α}^{d} {\dot{S}}_{r}$

(18)

可以预测动态条件下的扫描线。

通过以上步骤,建立了动态扫描线的模型,计算时需首先给出边界面曲线的方程,然后再标定出扫描线参数c和细阻尼系数τ,即可对动态扫描线进行预测。

3. 模型验证

为了验证模型,下面将利用上述模型跟已有的实验结果进行比较。首先要确定出边界面,这里将选用Feng等^[33]提出的方程作为边界面的方程：

$P_{c} = b {(\frac{S_{r}^{s a t} - S_{r}}{S_{r} - S_{r}^{irr}})}^{1 / d}$ ,

(19)

式中,b和d为材料参数,对于干燥和浸湿曲线分别取不同值, $S_{r}^{s a t}$ 为饱和时饱和度, $S_{r}^{irr}$ 为残余饱和度。

模型的验证主要包含以下步骤：①根据静态边界线试验结果,利用式（19）拟合出干燥和湿化边界的参数b和d；②根据任意一条静态扫描线试验结果,利用式（9）标定出参数c；③根据动态边界曲线和静态边界曲线试验结果,利用式（4）标定出参数τ；④将以上到的参数,代入式（17）,（18）可预测动态扫描线。

按照上述步骤,接下来将利用Zhuang等^[34]给出的干燥曲线结果和Milatz等^[35]给出的湿化曲线的结果对模型进行验证,模型的部分参数如表1所示。

表 1 土样材料参数列表

Table 1. Summary of material parameters

文献	$S_{r}^{sat}$	$S_{r}^{irr}$	干燥		浸湿		c
文献	$S_{r}^{sat}$	$S_{r}^{irr}$	b/kPa	d	b/kPa	d	c
Zhuang等^[34]	0.86	0.21	0.35	9.06	0.23	10.04	200
Milatz等^[35]	1.00	0.18	5.59	8.73	2.95	3.28	18

下载: 导出CSV

| 显示表格

3.1 干燥扫描线预测

首先利用模型对干燥条件下的扫描线进行预测,选用试验数据来源于Zhuang等^[34]。Zhuang等^[34]分别测量了砂土在静态和动态条件下的土水特征曲线,如图2中的试验点所示。利用静态边界线试验点标定土水特征曲线参数b和d,再利用静态扫描线试验点可以标定参数c,所得参数如表1所示。

根据试验结果,动态主干燥曲线的饱和度变化率为 ${\dot{S}}_{rD} = - 0.05 s^{- 1}$ 。利用动态干燥线和静态干燥线可以标定参数τ(kPa•s),Zhuang等^[34]将其表示为

$τ = 200 \times 10^{9.4 (0.85 - S_{r})}$ 。

(20)

边界线的确定对于预测扫描线非常关键,但Zhuang等^[34]只给出了动态干燥曲线,而动态湿化曲线由于试验失败,没有相关数据,可利用静态湿化边界曲线计算得到。假设参数τ在湿化与干燥条件下取值相同,即可将式（20）得到的参数τ代入式（4）中计算出动态湿化边界线,如图2中的虚线所示。然后将表1中的参数代入式（17）,（18）,可得到动态扫描线的预测结果,如图2中的虚线所示,所得结果与试验结果相符,表明模型可以较好地反映干燥过程中动态试验的结果。

从图2可以看到动态干燥线位于静态干燥线的右侧,表明相同饱和度条件下,动态曲线的吸力要大于静态曲线的吸力。流速越大,饱和度变化速率越快,动态吸力会越大。这对于模拟地基降水过程中非饱和土体的变形问题有重要的意义,快速降水将导致土体中的吸力更大,由吸力更大导致变形也将更大,说明快速降水可能增加非饱和土地基的沉降。

3.2 湿化扫描线预测

为了说明模型在湿化过程中的合理性,接下来采用Milatz等^[35]的试验结果对模型进行验证。Milatz等^[35]采用改进的直剪仪对3种不同密实程度的砂土开展了测试,这里选择动态效应最明显的细砂来进行验证。试验给出了水的流速为5 mm³/s,相应的饱和度变化率为 ${\dot{S}}_{r} = 2.9 \times 10^{- 4} s^{- 1}$ 。

由于试验只给了干燥边界线,并没有湿化边界线,这里只能假定给出一个静态的湿化边界面参数,动态湿化曲线采用与动态干燥线相同的参数τ。计算所需的参数如表1所示。参数τ(kPa·s)采用Milatz等^[35]给出的表达式：

$τ = 5.2812 \times 10^{3} - 3.795 \times 10^{2} \times S_{r}$ 。

(21)

利用所得参数,对动态扫描线进行预测,所得结果如图3所示。虚线为动态扫描线的预测曲线,与动态扫描线试验点比较吻合,表明模型可以较好地反映湿化过程中动态试验的结果。

动态湿化扫描线对于模拟降雨滑坡问题具有重要意义。目前模拟降雨滑坡问题采用的是静态曲线,但是通过以上计算可以看到,不同降雨强度使得饱和度的变化率不同,此时扫描线不重合,从图3可以看到动态湿化曲线位于静态湿化曲线下方。降雨强度越大,相同的饱和度对应的吸力将更小,根据非饱和土的强度理论可知,吸力对强度是有贡献的,吸力的减小可能降低土体的强度,从而可以解释为什么强降雨时更容易出现滑坡。

4. 结语

本文对非饱和土的动态土水特征曲线进行了探讨,建立了动态土水特征曲线的边界面模型,利用模型可以较好地预测动态土水特征曲线的扫描线变化。

建立动态土水特征曲线模型对于实际应用具有重要的意义,目前的土水特征曲线大多是基于室内试验获得,这些试验都是在平衡条件下给出的,但实际问题可能一直处在一个动态过程,用平衡条件下得到的土水特征曲线和水力参数去预测一个动态或非平衡的实际问题,就会存在误差。

数值预测结果表明,相同饱和度,动态条件下所对应的吸力与静态条件不同,动态干燥过程吸力要更大,而动态湿化过程的吸力要更小,这对于分析实际工程问题具有重要意义,应用动态的土水特征曲线来分析实际问题是下一步要开展的工作。

图 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/.

施引文献(14)

期刊类型引用(6)

1.	张高翔，刘艳，刘志强. 土-水特征曲线的动态效应试验及模型研究. 岩土力学. 2025(01): 178-186+198 . 百度学术
2.	郑健龙，刘绍平，胡惠仁. 公路路基湿度计算理论研究进展. 中外公路. 2023(01): 1-10 . 百度学术
3.	王欣. 基于孔隙水分布的非饱和土SWCC及相对渗透系数预测模型. 长江科学院院报. 2022(02): 89-93+101 . 百度学术
4.	伏映鹏，廖红建，吕龙龙，柴小庆. 考虑接触角及粒径级配影响的土水特征曲线滞回模型. 岩土工程学报. 2022(03): 502-513 . 本站查看
5.	罗婷婷，王林，裴鹏，杨斌，邹行. 换热孔回填材料含水率特征变化规律及其对热物性影响研究. 太阳能学报. 2022(07): 485-492 . 百度学术
6.	许韬，白冰. 考虑循环温度驱动的土-水特征曲线滞后模型. 岩土力学. 2022(12): 3393-3402 . 百度学术

其他类型引用(8)

资源附件(1)

其他相关附件
- PDF格式
  19-202311-G22-0781一文审稿意见与作者答复点击下载(1700KB)

图(8) / 表(2)

计量

文章访问数: 301
HTML全文浏览量: 46
PDF下载量: 80
被引次数: 14

0. 引言
1. 动态土水特征曲线的理论基础
2. 动态土水特征曲线模型
2.1 边界面模型的理论基础
2.2 动态条件下边界面模型的建立
3. 模型验证
3.1 干燥扫描线预测
3.2 湿化扫描线预测
4. 结语

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

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

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

计量

出版历程