基于 MATLAB 的网格划分和下限有限元法实现
黄 琳
1
肖洪天 2
1. 长沙理工大学 土木建筑学院，长沙 417100 2. 山东省土木工程防灾减灾重点实验室，青岛 266590
摘要：针对极限分析下限法中静力许可应力场建立困难这一问题，本文在下限定理的基 础上，引入有限元的思想来构造静力许可应力场。本文通过自编的 MATLAB 程序自动 进行映射网格划分，并利用 MATLAB 软件优化工具箱内置的内点算法求解线性规划模 型。通过对经典算例进行求解，计算结果表明该方法是一种合理有效的方法。 关键词：极限分析；下限定理；有限元法；线性规划；网格划分
作 者：黄 琳，岩土工程本科生. E-mail: franklyn0601@
1
1 引言
20 世纪 50 年代，D.C.Drucker 和 W.Prager 基于极值原理建立了土体极限分析理论； W.F.Chen[1]进一步阐明了极限分析理论在土工问题中的应用。目前研究土工问题的常用 方法主要有极限平衡法与有限元法。前者以概念简单明了，实际工程经验丰富而在工程 中得到广泛应用。但极限平衡法需要预先假设滑动面，且并不一定满足精力平衡条件， 因此它的解即缺乏理论基础同时又不是严格意义上的上下限解； 后者虽能弥补前者的不 足，但却无法回避本构关系的盲目性以及相应力学参数的不确定性。 S.W.Sloan[2] 提出的在二维条件下基于线性规划的下限分析有限单元法在理论上克 服了上述两种方法的不足。该方法引入有限单元法构造静力许可的应力场，可直接得到 材料到达塑性流动时所对应的极限荷载或者安全系数。换言之，应力场在满足静力许可 的约束条件下，将对应的极限荷载或安全系数作为目标函数进行优化求解，这样下限问 题可以转化为一个等价的约束非线性规划问题。 又可以通过用内接正 P 边形去拟合屈服 面，将非线性规划问题线性化。另外国内学者在下限分析中也做了很多工作，如：陈祖 煜证明了边坡稳定极限分析的垂直条分法和斜条分法的基础分别是塑性力学下限和上 限原理；李国英等在采用下限原理有限元解决结构面问题，三维非线性问题等方面做了 一定的探索。 上述学者大多采用 Fortran 语言编写塑性极限下限法程序并形成结构土体的约束规 划模型，调用商业的或自编的优化程序求解抽象得到的数学规划模型。本文则直接基于 MATLAB 平台求解整个计算过程，针对边坡和条形基础承载力问题进行计算。

e
Байду номын сангаасe x1
e e e e e y xy x y xy 1 1 2 2 2
这样每条特定的存在边界应力条件的边界 l 都会产生最多四个等式约束条件。
q2 t2
q1
t1
⑴
l
⑵
y
side l
x
4
2.4 应力连续性条件 静力许可的间断线允许应力连续性条件与经典的下限法一样， 在相邻单元的公共边 界上的正应力和剪应力的分量相等，可得：
2 下限原理有限元法
2.1 下限原理 下限定理认为一个满足平衡条件， 应力边界条件和屈服准则的静力场所对应的荷载 将不大于真实的破坏荷载。因此选定一个合适的静力场是解决工程实际问题的关键。根 据下限定理，设定求解区域 V 及其应力边界 S，满足以下三个条件的应力场ij 称为静力 许可应力场：
ij , j Fi 0
sin 2 d cos 2 d

d T T d equil
a x1
a a b b b a a a b b b y xy x y xy y xy x y xy 1 1 2 2 2 x1 1 1 2 2 2
b 0
2.5 屈服条件
MESH GENERATION AND LOWER BOUND ANALYSIS USING FINITE ELEMENT METHOD ON MATLAB
HUANG Lin 1 XIAO Hong-tian 2
1. School of Civil Engineering and Achitecture, Changsha University of Science and Technology, Qingdao 266590; 2. Shandong Provincial Key Laboratory of Disaster Prevention and Mitigation, Qingdao 200092 Abstract: Given the difficulty of finding a statically admissible field, this paper describes a technique for computing lower bound limit loads using finite element method in soil mechanics under conditions of plane strain. Integrated with the optimization toolbox of Matlab and automatically mapping mesh method, a new calculating method is proposed to solve linear programming problems. Several examples of the ultimate bearing capacity of a rigid strip footing and slope in homogeneous soil was presented to illustrate the validation and effectiveness of the present method with the comparison between numerical lower bound analysis and analytical solutions. Key words: limit analysis; lower bound analysis; finite elements; linear programming; mesh generation
A2 b1
c A1 b1
T
对于工程问题，其静力许可场的构造非常困难，更不用说计算极限荷载。不过随着 有限元的引入和计算机的应用，该问题已经可以求解。在对计算区域进行离散后，取节 点应力为优化变量，这里我们借助于有限元中三角形三节点单元的插值函数构造应力 场，并施加平衡条件，应力间断条件，边界条件和屈服条件。这样就构成了有限自由度 的规划问题来求解极限荷载。 由于使用的是三角形三节点单元进行离散，应采纳对应单元的插值函数如下：
A b
d equil d d equil
b

n
⑵ ⑴
⑷ ⑶
d
a
y
x
式中
T A 0
d equil
T 0
0 T

0 T
sin 2 d T 1 sin 2 d 2
cos 2 d 1 sin 2 d 2
2k Ck 2 sin p ， D 2c cos cos p
y 2 xy
k 6
k 5
Mohr-Coulomb 屈服函数

p
k 1 k 4
x x y
k=2
k=3
线性 Mohr-Coulomb 屈服准则 （p=6）
3 基于 MATLAB 的网格划分和求解策略
本文基于 MATLAB 7.10.0 编写程序。特点如下： 1）网格的自动划分； 其基本思想为：将待划分的整个区域分成若干个超单元，要求每个超单元为单连通 区域，而整个区域可为单连通的或复连通的。每个超单元由 4 个节点描述之。然后对每 个超单元进行再划分。 2）数据的可视化； 数据可视化是 MATLAB 的一项重要功能。在后处理中，作者为直观分析样本数据 的分布和趋势特性，将单元间各个节点的应力值绘制成二维平面等值线图。 3）矩阵的稀疏存储； 矩阵的稀疏存储对于塑性有限元极限分析具有重要意义。 因为约束条件通常是由一 系列高度稀疏的矩阵组合而成，利用稀疏存储可以显著提高计算速度，减少所需存储空 间。结果比对证明，MATLAB 可以有效准确的处理塑性极限分析中的线性规划求解问 题。
6
4 程序编制与算例
4.1 程序编制 根据上述推导的数学模型，本文采用 MATLAB 7.10.0 编写了基于线性规划的有限 元塑性分析下限法主体计算程序 CALORS.m 和网格划分程序 Meshgenerator.m，主体计 算程序框图如下图所示。程序可根据给定边界自动进行网格划分，并手动控制网格的精 细程度，并计算外力超载系数以及生成应力场。
式中： c 为材料的黏聚力； 为材料的内摩擦角。 用内接正 p 边形去拟合屈服面，如图所示，那么每个节点的屈服条件均可用 p 个线 性方程代替，则正多边形第 k 条边表示的约束条件为
Ak x Bk y Ck xy D, k 1,2, p
5
2k 2k 其中： Ak cos p sin cos p ， Bk sin cos p cos p ，
A b
l bound l l bound
其中
T A 0
l bound

0 T sin 2 l cos 2 l
sin 2 l cos 2 l [T ] 1 1 sin 2 l sin 2 l 2 2 l bbound q1 t1 q2 t2
0 0 0
因此每条相邻单元的公共边界的应力间断线都会产生四个等式约束条件。
对于平面应变问题，假设拉应力为正，则 Mohr-Coulomb 屈服准则可表示为
F x y 2 xy 2c cos x y sin 0