求解地基极限承载力的上限有限元滑动面搜索法

王宇航 杨峰

摘 要：提出一种基于上限有限元的滑动面搜索法用于求解地基极限承载力。计算思路为：通过将破坏区域按照假定可能的破坏模式的方式进行网格划分并且网格参数化，然后利用上限有限元求解对应参数条件下的上限解，再从中选取最优上限解及其对应的滑动面，即达到利用上限有限元搜索最优滑动面的目的。该法综合了上限有限元计算灵活、刚体滑块上限法可获得直观滑动面的特点，可高效快速的获得一系列上限解并从中优选，具有计算成本低，且能保证计算精度等优势。通过条形基础地基极限承载力的计算分析，验证了该法的有效性，并可望应用于类似问题。

关键词：极限分析上限有限元上限法刚体滑动面地基极限承载力

1 引言

在进行岩土工程稳定性分析时，采用极限分析上限法能直接快速获得上限解。传统的上限法（刚体滑块上限法）需预先假定破坏模式[1]，而上限有限元法通过将破坏区域离散成三角形单元，通过系统耗能最小原理直接搜索获得上限解[2,3]。

刚体滑块上限法对于特定问题求解简单，并可得到具体的滑动面；上限有限元法计算结果与网格单元存在依赖性，计算获得的滑动面形状不明显。

本文试图结合两者的优势，即采用假定破坏模式的方式划分上限有限元网格，将破坏区域的滑动面表示成带参数的多段线，调整滑动面参数并利用上限有限元进行计算获得一系列上限解，最后从中选取最优解。

本文以经典的条形基础地基承载力问题为例，说明采用上述思路的实现过程，计算结果与现有文献进行对比分析，论证基于上限有限元的滑动面搜索法对特定问题的可行性和有效性。

2 极限分析上限有限元

极限分析上、下限理论应用于岩土工程稳定性分析的原理和方法在Chen W F的著作中有详细的论述[1]。由于计算简便，极限分析上限法长期以来均采用与极限平衡法类似的假定刚体滑面模型的形式。

然而，当破坏模式不易假定时，采用Sloan等提出的上限有限元模型就变得更加有效[2]。本文所采用的上限有限元原理和流程参考了Sloan等所作的工作[2,3]，区别在于网格划分按照假定破坏模式的方法进行。

当采用线性规划模型时，上限有限元需对摩尔-库伦屈服准则进行线性化，之后模型即可转化为如下线性规划问题：

Minimize(1)

Subject to(2)

式(1)、(2)中为目标函数(1)的系数向量；为决策变量，由单元节点速度（实域）、速度间断线辅助速度参数（非负值）和单元内部塑性乘子（非负值）组成；为等式约束系数矩阵，为等式约束右侧向量。

由虚功率平衡方程获得目标函数(1)；等式 (2)由单元内部和速度间断线塑性流动约束条件、速度和应力边界约束条件组成。上限有限元基本原理及具体实现过程参见文献[3]。

3 条形基础地基承载力

3.1 上限有限元模型建立

地基极限承载力可方便的转化为求解承载力系数,和，其中与自重相关的承载力系数的解答需通过数值计算获得。为减少篇幅，本文仅进行的计算。

地基承载力上限有限元模型网格划分如图1所示。利用对称性，只考虑模型右侧的一半。可以看出，三角形单元的布置与现有的破坏模式近似。破坏模式由5个三角形组成，其中过渡区三角形个数；破坏模式可由图中的角度,,,唯一确定。为使刚体运动更加灵活，过渡区中每个三角形的角度,可取不同数值。以角度为决策变量，耗散能最小化为目标函数，对应的速度矢量闭合图为几何约束条件，即构成刚体滑块上限法的求解方法。

刚体滑块上限法认为耗散能仅在间断线上发生[4]。采用上限有限元时，三角形单元也允许发生塑性变形，而不再是刚体。如图1，上限有限元求解承载力问题时，建立图示的坐标系，模型边界条件为：基础下方三角形单元左侧边节点向速度分量，上边向速度分量，向速度分量；下边界滑动面对应的边界条件（刚性边界）为，此约束并未将临近的三角形单元节点速度置零，而是在刚性边界添加虚拟节点，并将这些节点的速度分量均置零，然后在滑动面上施加间断线约束条件。

图1 地基承载力上限有限元模型网格划分及边界条件

由于基础下方的三角形向速度置零，意味着基础与地基之间完全粗糙。

上限有限元计算前需确定三角形单元的数目()和摩尔-库伦屈服准则线性化对应的塑性乘子数目，之后将模型的几何和力学参数导入已编制的上限有限元程序，通过求解线性规划问题即可获得承载力系数的上限解。

3.2 地基破坏时滑动面搜索策略

按照图1模型设置的三角形单元较少，因此对于特定网格所获上限解精度不高。于是，将上述上限有限元的网格参数化，即由三个角度参数,,确定一系列的网格，计算对应的上限解，再从中选取最小值，与其对应的,,所确定的破坏模式和滑动面即为搜索过程获得的最终结果。

从图1可知，参数,,的取值存在合理范围：

(3)

因此，参数,,的取值可在式(3)所示的范围内均匀选取。同时，为了减少无效的计算量，在试算的基础上也可将式(3)中的参数取值范围进一步缩小。

4 计算结果讨论

按照上述方法，以下采用上限有限元进行地基承载力的计算以及最优滑动面的搜索。

4.1 计算参数的设置

选取破坏模式过渡区三角形单元数目为，塑性乘子数目。于是模型单元总数为102，间断线数目202。以内摩擦角为 例，参数的取值范围选取为；参数的取值范围选取为；参数的取值范围为。

计算步骤为：①按照2°间隔依次选取参数；②对于每个值，按1°间隔依次选取参数；③对于每组和值，按1°间隔依次选取参数；④对于每组,和值，进行一次上限有限元计算并记录获得的值；⑤进行循环计算，获得取值范围内每个,和值组合对应的值；⑥得到所有值的最小值。

4.2 计算结果分析与讨论

当时，采用上限有限元法获得的承载力系数计算结果如表1。

表1 地基承载力系数Nγ计算参数和结果

/° /° /° Nγ

128 100 100 40 53 21.00

102 40 51 20.63

104 40 48 20.43

106 40 46 20.40

108 40 43 20.51

110 41 41 20.78

112 41 39 21.20

114 42 37 21.74

116 41 36 22.34

118 42 32 22.84

120 43 30 23.54

表1列出了不同参数取值时计算结果。其中最优解为，对应的,和值分别为106°, 40°和46°。

最优解对应的三角形单元网格划分和破坏模式见图2。如图2(a)，其中地基破坏模式的过渡区内划分了密集的三角形单元，以形成速度间断线并减小滑动面范围。图2(b)为破坏模式，即表示基础向下移动单位长度1时，地基内部破坏区域的运动形态。可以看出，单元之间错动，说明速度间断线的作用明显。被动区的三角形单元向上挤出，说明单元发生了显著的变形，不再为刚体。

(a) 三角形单元网格划分

(b) 破坏模式

圖2 计算承载力系数时的网格划分和所获破坏模式

按4.1节计算步骤，表1列出了时不同取值对应的较优上限解和参数,值。将表中每个取值对应的破坏区域滑动面形状绘制如图3所示。其中值越大，对应的滑动面范围越大。最优解对应于从上至下第四条滑动面。可以想象，破坏范围越大，对应的值越大；然而破坏范围越小，其破坏区域的速度场变化越剧烈，自重功率也将增加。于是，最优解对应的滑动面包含在图3所示的最大与最小破坏区域之间。

图3 不同取值对应的破坏区域滑动面

为说明不同参数取值时计算结果的差异，以下将时，不同和取值对应的值示意如图4。

图4 和不同取值对应的承载力系数计算结果

图4右侧坐标轴表示角度，取值范围为；左侧坐标轴为角度，取值范围为；竖向坐标轴表示承载力系数。从图中可知，当取较小值、取较大值时，计算得到的值较大。值组成的曲面呈凹形，尽管曲面底部较平缓，但仍存在最小值为20.40，说明通过设定参数取值范围并逐次计算确实可以搜索到最优解。

表2列出了内摩擦角取不同值时所得的承载力系数计算结果以及对应的角度参数,和值，计算选取, 。其中右侧文献[5]为Soubra采用多刚体块上限法获得的值；而文献[6]为Martin采用滑移线法经高精度数值计算获得的值，可认为是精确解。

表2 地基承载力系数Nγ计算参数和结果

/° /° /° /° 本文

Nγ 文献[5]

Nγ 文献[6]

Nγ

10 86 3 60 0.71 0.85 0.43

20 92 26 56 4.24 4.67 2.84

30 106 40 46 20.40 21.88 14.75

40 118 52 38 113.77 120.96 85.57

从表2可看出，本文所得上限解优于采用刚体滑块上限法的文献[5]，这与单元允许发生塑性变形有关；本文上限解较之文献[6]的精确解仍存在误差，可通过进一步细化网格，特别是设置多层单元的方式进行弥补，这也是后续需展开的工作。

为了说明塑性乘子数目对计算结果的影响，将值取4~128时对应结果列表如表3。可以看出，值较小时，所得值变大，计算结果精度降低；而当值大于64时，计算结果的变化不再明显。

表3 地基承载力系数Nγ计算参数和结果

/° /° /° Nγ

4 100 106 40 46 21.43

8 20.96

16 20.55

32 20.44

64 20.40

128 20.40

地基破坏模式中的过渡区的三角形数目亦对计算结果有影响，将取不同值时对应的值列表如表4。可以看出，随着值增加，计算结果变小，计算精度增加。这主要得益于三角形数目越多，破坏范围越小且过渡区速度间断线也越多。

表4 地基承载力系数Nγ计算参数和结果

/° /° /° Nγ

128 10 106 40 46 22.32

20 21.72

40 20.99

60 20.68

80 20.50

100 20.40

为了揭示值对破坏模式以及滑动面的影响，将时的三角形单元网格划分和计算得到的破坏模式示意如图5。与图2对比可发现，取值较小时，破坏范围相应变大，破坏区域的变形特别是基础角点附近的变形剧烈程度降低，因此所得上限解精度降低。

(a) 三角形单元网格划分

(b) 破坏模式

图5 计算承载力系数时的网格划分和所获破坏模式

5 结论

以条形基础地基承载力问题为例，提出采用假定破坏模式的方式建立上限有限元计算网格，并将破坏模式的滑动面参数化。通过上限有限元求解对应参数条件下的一系列承载力系数的上限解，再从中选取最优上限解和对应滑动面，即达到采用上限有限元搜索最优滑动面的目的。

基于上限有限元的滑动面搜索法具有计算成本低，且计算精度有保证，能直观获得破坏区域滑动面的优势；该法可望通过进一步的改进应用于特定的岩土工程稳定性课题。

参考文献：

[1]Chen W F, Limit analysis and soil mechanics [M]. Elsevier Scientific Publishing Company, New York, 1975.

[2]Sloan S W, Kleeman P W. Upper bound limit analysis using discontinuous velocity fields [J]. Computer Methods in Applied Mechanics and Engineering, 1995, 127: 293–314.

[3]楊峰, 阳军生, 张学民. 基于线性规划模型的极限分析上限有限元的实现 [J]. 岩土力学, 2011, 32(3): 914–921.

[4]陈祖煜. 土力学经典问题的极限分析上、下限解 [J]. 岩土工程学报, 2002, 24(1): 1–11.

[5]Soubra A H. Upper-bound solutions for bearing capacity of foundations [J]. Journal of Geotechnical and Geoenvironmental Engineering, 1999, 125: 59–68.

[6]Martin C M, Exact bearing capacity calculations using the method of characteristics [C]// Proceedings of the 11th International Conference of IACMAG, Turin, 2005, 4: 441–450.

注：文章内所有公式及图表请用PDF形式查看。