用于上限有限元的非结构化网格重划加密方法研究
杨峰，阳军生
(中南大学 土木工程学院，湖南 长沙，410075)
摘要：提出一种用于二维上限有限元的非结构化网格数据转化方法并编制数据处理和上限有限元程序，借助MESH 2D非结构化网格划分工具建立基于耗散能密度为评判参数的网格重划加密方法，可达到不明显增加单元数量而获得较优上限解的目的。同时，分别实现用于上限有限元的六节点三角形单元无间断线和三节点三角形单元+速度间断线这2种网格模式。研究结果表明：六节点三角形单元无间断线模式结合非结构化网格重划加密方法能在同等单元数量的前提下显著提高上限解精度，显示出明显的优势和实用性。
关键词：极限分析；上限有限元；三角形单元；速度间断线；破坏模式；非结构化网格；网格加密
中图分类号：TU43 文献标志码：A 文章编号：1672-7207(2014)10-3571-07
Investigation of unstructured mesh regeneration and refinement method for finite element upper bound solution
YANG Feng, YANG Junsheng
(School of Civil Engineering, Central South University, Changsha 410075, China)
Abstract: A mesh data transforming method with general unstructured mesh for finite element upper bound solution was proposed and corresponding program was compiled. By means of unstructured meshing tool MESH 2D,a mesh regeneration and refinement method was built based on evaluation of dissipated energy to obtain favorable upper bound solution without the increase of triangular elements. Meanwhile, the two modes of six nodal triangular elements without velocity discontinuities and three nodal triangular elements with velocity discontinuities were implemented for finite element upper bound solution. The results show that the precisions of results are improved significantly by using six nodal triangular elements mode with unstructured mesh regeneration and refinement method.
Key words: limit analysis; finite element upper bound solution; triangular element; velocity discontinuity; failure mechanism; unstructured mesh; mesh refinement
目前，与有限元技术相结合的极限分析上限法即上限有限元法用于岩土工程稳定性分析已成为研究热点[1-5]。该法假定岩土体为理想弹塑性材料，利用上限原理将岩土稳定性问题转化为线性或非线性规划模型进行求解[6-8]。上限有限元法具有数值求解稳定、快速、无需假定破坏模式、适应复杂模型等诸多优势。然而，即便对于二维问题，目前该法还存在如下困难：1) 上限有限元形成的线性或非线性规划模型含有大量的决策变量和约束条件；2) 常采用三节点或六节点三角形单元，其他类型单元应用难度大；3) 常需在相邻单元间设速度间断线，网格数据格式与共节点的有限元不同；4) 岩土材料失稳破坏往往形成局部化的剪切带和塑性区，获得较优上限解和破坏形态常需划分密集的网格；5) 目前尚无成熟的上限有限元软件可用。对于剪切带和塑性区的精确模拟问题，李大钟等[9]通过衡量屈服准则残余和等效应变，对网络进行剖分加密以实现极限分析有限元网络自适应。为改进上限有限元研究，本文作者提出1种二维上限有限元非结构化网格数据转化方法并编制数据处理和计算程序，将MESH 2D[10]网格划分工具用于上限有限元分析；提出基于耗散能密度评判参数的网格重划加密技术，通过3次网格重划加密和上限有限元计算，达到不显著增加计算规模并获得较优上限解的目的；实现用于上限有限元的三节点三角形单元+速度间断线和六节点三角形单元无间断线这2种网格模式。最后，通过对比已有文献中关于边坡稳定性算例，验证上述方法的有效性。
1 上限有限元与六节点三角形单元
三节点三角形单元+速度间断线模式用于上限有限元的相关理论见文献[6，8]。对于六节点三角形单元用于上限有限元，Yu等[11]以及Makrodimopoulos等[12-13]均进行了研究，前者研究了内摩擦角
为0°且模型为线性规划的情况，后者采用二阶锥规划模型进行研究。本文旨在建立基于线性规划模型和六节点三角形单元无间断线模式的上限有限元并用于
型土体，其基本理论[14]如下。
1.1 六节点三角形单元
六节点三角形单元如图1所示，六节点三角形单元每边中点增设1个节点。每个节点含2个速度分量，单元顶点对应的节点1，2和3还分别包含p个塑性乘子
。
单元内部任一点的速度可表示为
(1)
式中：ui和vi分别为节点x和y向速度分量；Ni为单元形函数，若用面积坐标L1，L2和L3表示，则
(2)
1.2 单元内部约束条件
在六节点三角形单元中，x向应变率
，剪应变率
和塑性乘子
在单元内部呈线性变化。
图1 六节点三角形单元
Fig. 1 Six nodal triangular element
对摩尔-库仑屈服函数线性化时，单元内部应服从相关联流动法则，于是，节点i(i=1，2，3)处塑性应变率(
，
，
)需满足
(3)
其中：Fk为线性化后的屈服函数；Ak，Bk和Ck为系数；σx，σy和τxy为应力分量；
为第i节点在屈服函数外接多边形第k条边对应的塑性乘子；k=1，2，…，p。于是，
(4)
其中：
式(4)为节点i上的塑性流动约束条件，单元顶点对应的节点(i=1，2，3)均需施加上述约束。
1.3 速度边界条件
设节点i位于与x轴夹角为θ的边界上，其切向和法向速度分别为
和
，则ui和vi须满足
(5)
其中：
；
；
1.4 自重功率
在如图1所示坐标系中，每个单元自重功率Pt为
(6)
其中：
γ为土体容重；A为三角形单元面积。
1.5 单元内部耗散功率
每个单元内部耗散功率Pc为
(7)
其中：
c为土体黏聚力；f为内摩擦角。
1.6 上限有限元的线性规划模型
根据上限定理，由耗散能、外力功率和虚功率平衡方程即可建立目标函数。其决策变量为节点速度分量ui和vi，以及单元塑性乘子。对于模型包含的所有单元、节点和边界条件，上限有限元线性规划模型如下。
最小化：
服从于：
(8)
式(8)中各量均为全局变量，与式(4)~(7)中的局部变量对应，如X1为所有节点速度向量，X2为所有单元塑性乘子向量。
2 非结构化网格重划加密方法
2.1 六节点三角形单元无间断线模式
本文非结构化网格重划加密方法系指由前次计算结果重新划分网格且变形较大的区域依网格尺寸权重加密，其中网格尺寸权重h与单元边长直接关联。以下以六节点三角形单元为例，将耗散能密度(这里与塑性乘子一致)作为网格尺寸权重的评判参数进行网格重划与计算。
1) 初次计算。初次计算可获得模型破坏区域的大致信息。以均匀网格划分模型，网格尺寸权重h为定值。然后转换网格数据，再进行上限有限元计算。三角形单元数目合理选取，若过大则将使计算变缓甚至不可行，若过小则计算结果精度降低。
2) 二次计算。二次计算前的网格重划需设置网格权重h，以特征点塑性乘子作为权重h评判参数。特征点为三角形单元形心，该点塑性乘子
为
(9)
为便于设置评判阀值
，将每个单元的归一化，于是，评判参数η为
(10)
式中：
为模型所有单元中的最大值。
网格重划时，对
的区域，划分较密的网格，即权重h取较小值；而对η＜△的区域，划分稀疏的网格，即h取较大值。之后利用该网格进行上限有限元二次计算。
3) 最终计算。最终计算前，网格重划所需尺寸权重h依据二次计算所得塑性乘子值，权重计算点仍为三角形单元形心。评判参数η的计算与前面的相同。所不同的是：此次网格重划时，由0≤η≤1取值范围设多级阀值△1，△2，…，△n，不同的阀值区间递减设置权重(h1，h2，…，hn)，以此划分网格并进行上限有限元最终计算。
2.2 三节点三角形单元+间断线模式
对于三节点三角形单元+间断线模式，其网格重划加密方法和步骤与六节点三角形单元的一致，不同的是评判参数η需考虑速度间断线耗散能。
仍选取耗散能密度作为网格权重h的评判参数，计算点为三角形单元形心，将单元相邻间断线的耗散能的一半叠加到该三角形单元耗散能上。
每个三角形单元耗散能密度
为[8]
(11)
该单元相邻间断线耗散能密度为[8]
(12)
式中：
，
，
和
均为速度间断线上非负优化变量。
设
，
和
分别为该单元三条边所在间断线耗散能密度，将其一半叠加到该单元可得
(13)
将模型耗散能密度ρ归一化，网格权重h的评判参数η为
(14)
2.3 MESH 2D网格划分与数据转换
MESH 2D是Engwirda[10]编制的一个二维非结构网格划分工具，可嵌入MATLAB。网格划分所得数据为节点编号与坐标、三角形单元编号及内部节点编号。对于六节点三角形单元无间断线模式，上限有限元所需的网格数据转化思路如下：
1) 原型每个节点含x和y向速度分量。
2) 获取单元的边编号，每条边中点附加1个节点，该节点含x和y向速度分量。
3) 塑性乘子归属于每个单元的3个顶点对应的节点，分属不同单元的同位置顶点节点处塑性乘子不共有。
4) 单元节点及附加量顺序按图1所示排序。
5) 依模型边界标识搜寻并记录边界节点、单元及局部坐标法向等信息，以便施加边界条件。
对于三节点三角形单元+间断线模式，上限有限元所需的网格数据转化思路为：
1) 速度分量归属于单元节点，相邻单元同位置处的节点不共有，塑性乘子归属于单元。
2) 所有单元的边依边界标识分为速度间断线和边界边2种。速度间断线附加辅助变量，边界边记录所属单元、所含节点及局部坐标法向等信息。
3) 搜索、判断并建立速度间断线两侧单元编号、间断线节点在两侧单元节点的位置等信息，以便施加速度间断线塑性流动约束。
2.4 上限有限元计算及网格重划加密流程
为便于编程，将上限有限元线性规划模型决策变量之一(节点速度分量)转化为非负量。上限有限元主程序采用MATLAB编制，网格划分利用MESH 2D[10]，线性规划模型采用SeDuMi[15]求解。
当采用网格重划加密和六节点三角形单元无间断线模式时，上限有限元的计算流程如图2所示。
图2 上限有限元计算流程图
Fig. 2 Flow chart for finite element upper bound solution
3 算例验证与讨论
选取文献[12-13]中的边坡稳定性算例进行上限有限元计算验证。该算例求解不同内摩擦角对应的边坡稳定性系数Ns。模型的初次网格划分、边界条件、几何尺寸如图3所示。边坡坡度为70°，左右和下边界均约束2个方向的速度，即u=0 m/s，v=0 m/s。边坡稳定性系数Ns(无量纲)定义为
(15)
图3 边坡稳定性计算模型网格划分及边界条件
Fig. 3 Mesh and boundary condition of calculating model with slope stability
本例边坡坡高H=25 m，黏结力c=1 kPa，内摩擦角=35°。以Ns为目标函数，令单位容重(γ=1)时自重功率为1，利用式(6)施加自重约束条件。摩尔-库仑屈服准则线性化参数p取32。
采用网格重划加密的上限有限元法，按六节点三角形单元无间断线和三节点三角形单元+间断线2种模式所得边坡稳定性系数Ns分别如表1和表2所示。表1和表2中也列出了文献[12]和[16]中的计算结果，并以文献[16]中的计算结果作为评判相对误差的参考值。
从表1可看出：采用六节点三角形单元进行上限有限元计算，初次和二次计算结果Ns的相对误差分别为14.90%和5.50%，Ns值均与文献[12]中的相应值14.97和14.45接近，不过文献[12]中的第2种网格单元数约为本文二次计算网格单元数的2.5倍。此外，本文最终计算的单元数为4 006，Ns相对误差为1.86%，而文献[12]中采用的网格单元数为28 864，Ns相对误差为2.38%。本文结果相对误差小且模型网格单元数少，这是由于采用了网格重划加密方法。
表1 边坡稳定性系数Ns上限有限元计算结果(六节点三角形单元无间断线)
Table 1 Results of slope stability factor Ns by finite element upper bound solution with six nodal triangular elements
表2 边坡稳定性系数Ns上限有限元计算结果(三节点三角形单元+间断线)
Table 2 Results of slope stability factor Ns by finite element upper bound solution with three nodal triangular elements and velocity discontinuities
从表2可以看出：采用三节点三角形单元+间断线进行上限有限元计算，二次和最终计算所得Ns的相对误差分别为11.75%和8.74%，均小于文献[12]中的13.13%和11.54%，而本文采用的网格单元数少，特别是最终计算时采用的网格单元数更少。这说明基于三节点三角形单元+间断线模式的网格重划加密方法对降低计算规模也具有优势。
对比表1和表2最终计算的网格参数可知：表1中单元、节点和约束条件数目均比表2中的少，仅变量数目比表2中的大，这说明二者计算规模差异不大。然而，所得Ns相对误差分别为1.86%和8.74%，说明采用网格重划加密方法时，六节点三角形单元无间断线模式比三节点三角形单元+间断线模式精度更高。
此外，本文认为对该边坡算例，文献[16]采用基于假定刚性滑体破坏模式所得上限解更优，这是由于其模型和破坏模式均简单明确，且形成明显剪切带。对于难以预先找到合理破坏模式的复杂模型，假定刚体滑块破坏模式的上限法精度优势将丧失。
以下对初次、二次和最终三阶段的网格塑性变形及最终所得耗散能密度分布进行讨论。图4所示为采用六节点三角形单元无间断模式所得网格塑性变形图(仅显示了模型局部范围)。初次计算采用均匀网格(见图4(a))。二次计算对耗散能较大的区域网格加密，加密区域网格相对均匀，其网格变形在坡脚处较剧烈(见图4(b))。最终计算采用的网格在坡脚附近相对密集，向上至坡顶形成1条密集的带状区域，稍远处网格稀疏(见图4(c))。图5所示为最终计算所得耗散能密度分布图。耗散能密度以塑性乘子替换并进行归一化，使之成为无量纲参数，在坡脚处其值最大，向上呈圆滑、狭窄的带状直至地表。耗散能密度所反映的带状区域即为边坡滑动面。
图4 不同计算阶段网格塑性变形图(六节点三角形单元)
Fig. 4 Plastic deformation of mesh for various calculating stages with six nodal triangular elements
图6所示为采用三节点三角形单元+间断线模式所得网格塑性变形图。其初次、二次和最终网格塑性变形图与图4所示的基本类似。不同之处在于由于存在间断线，坡脚处单元之间错动明显，显示了间断线模拟剧烈变形区的优势。然而，三节点单元模拟剧烈变形区效果差，最终计算前网格重划加密形成的密集区较宽，计算结果精度也较低。图7所示为采用该模式最终计算所得耗散能密度分布图，耗散能密度同样进行了无量纲处理，该图显示的耗散能密度分布形状与图5所示的类似，不过其表示的剪切带更宽泛。
图5 最终计算所得耗散能密度分布图(六节点三角形单元)
Fig. 5 Distribution of dissipated energy density for last calculating stage with six nodal triangular elements
图6 不同计算阶段网格塑性变形图 (三节点三角形单元+间断线)
Fig. 6 Plastic deformation of mesh for various calculating stage with three nodal triangular elements and velocity discontinuities
图7 最终计算所得耗散能密度分布图 (三节点三角形单元+间断线)
Fig. 7 Distribution of dissipated energy density for last calculating stage with three nodal triangular elements and velocity discontinuities
4 结论
1) 通过网格单元数据转换，可将MESH 2D作为上限有限元的网格非结构化网格划分工具，从而减小复杂模型建模的工作量。
2) 基于耗散能密度为评判参数的网格重划加密方法，能依据前次结果信息重划网格对应变率较大的区域网格加密，达到降低计算规模且不影响计算精度的目的。
3) 配合网格重划加密方法，采用六节点三角形单元无间断线模式比三节点三角形单元+间断线模式更能有效提高计算精度并降低计算规模。
参考文献：
[1] Augarde C E, Lyamin A V, Sloan S W. Stability of an undrained plane strain heading revisited[J]. Computers and Geotechnics, 2003, 30: 419-430.
[2] 王均星, 王汉辉, 吴雅峰. 土坡稳定的有限元塑性极限分析上限法研究[J]. 岩石力学与工程学报, 2004, 23(11): 1867-1873.
WANG Junxing, WANG Hanhui, WU Yafeng. Stability analysis of soil slope by finite element method with plastic limit upper bound[J]. Chinese Journal of Rock Mechanics and Engineering, 2004, 23(11): 1867-1873.
[3] 殷建华, 陈健, 李焯芬. 岩土边坡稳定性的刚体有限元上限分析法[J]. 岩石力学与工程学报, 2004, 23(6): 898-905.
YIN Jianhua, CHEN Jian, LI Zhuofen. Upper bound limit analysis of stability of rock and soil slopes using rigid finite elements[J]. Chinese Journal of Rock Mechanics and Engineering, 2004, 23(6): 898-905.
[4] 杨峰, 阳军生, 张学民, 等. 黏土不排水条件下浅埋隧道稳定性上限有限元分析[J]. 岩石力学与工程学报, 2010, 29(S2): 3952-3959.
YANG Feng, YANG Junsheng, ZHANG Xuemin, et al. Finite element analysis of upper bound solution of shallow-buried tunnel stability in undrained clay[J]. Chinese Journal of Rock Mechanics and Engineering, 2010, 29(S2): 3952-3959.
[5] YANG Feng, YANG Junsheng. Stability of shallow tunnel using rigid blocks and finite element upper bound solutions[J]. International Journal of Geomechanics, 2010, 10(6): 242-247.
[6] Sloan S W. Upper bound limit analysis using finite elements and linear programming[J]. International Journal for Nummerical and Analytical Methods in Geomechanics, 1989, 13(3): 263-282.
[7] Sloan S W, Kleeman P W. Upper bound limit analysis using discontinuous velocity fields[J]. Computer Methods in Applied Mechanics and Engineering, 1995, 127(1): 293-314.
[8] 杨峰, 阳军生, 张学民. 基于线性规划模型的极限分析上限有限元的实现[J]. 岩土力学, 2011, 32(3): 914-921.
YANG Feng, YANG Junsheng, ZHANG Xuemin. Implementation of finite element upper bound solution of limit analysis based on linear programming model[J]. Rock and Soil Mechanics, 2011, 32(3): 914-921.
[9] 李大钟, 郑榕明, 王金安, 等. 自适应有限元极限分析及岩土工程中的应用[J]. 岩土工程学报, 2013, 35(5): 922-929.
LI Dazhong, ZHENG Rongming, WANG Jinan, et al. Application of finite-element-based limit analysis with mesh adaptation in geotechnical engineering[J]. Chinese Journal of Geotechnical Engineering, 2013, 25(5): 922-929.
[10] Engwirda D. MESHZD-automatic mesh generation[EB/OL]. http://www. mathworks. com/ matlabcentral/ fileexchange/ 25555 -mesh2d –automatic –mesh –generation, 2009-10-12.
[11] Yu H S, Sloan S W, Kleeman P W. A quadratic element for upper bound limit analysis[J]. Engineering Computations, 1994, 11(3): 195-212.
[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.
[13] Makrodimopoulos A, Martin C M. Upper bound limit analysis using discontinuous quadratic displacement fields[J]. Communications in Numerical Methods in Engineering, 2008, 24(11): 911-927.
[14] 杨峰, 阳军生, 李昌友, 等. 基于六节点三角形单元和线性规划模型的上限有限元研究[J]. 岩石力学与工程学报, 2012, 31(12): 2556-2563.
YANG Feng, YANG Junsheng, LI Changyou, et al. Investigation of finite element upper bound solution based on six nodal triangular elements and linear programming modal[J]. Chinese Journal of Rock Mechanics and Engineering, 2012, 31(12): 2556-2563.
[15] Sturm J F. Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones[J]. Optimization Methods and Software, 1999, 11(12): 625-653.
[16] Chen W F. Limit analysis and soil mechanics[M]. New York: Elsevier Scientific Publishing Company, 1975: 254-255.
(编辑 陈灿华)
收稿日期：2013-09-15；修回日期：2013-11-23
基金项目(Foundation item)：国家自然科学基金资助项目(51008309)；“十二五”国家科技支撑计划项目(2012BAK24B02)(Project (51008309) supported by the National Natural Science Foundation of China; Project (2012BAK24B02) supported by the National Science & Technology Pillar Program During the 12th Five-year Plan Period)
通信作者：阳军生(1969-)，男，湖南永兴人，博士，教授，博士生导师，从事隧道与地下工程的研究；电话：13973182363；E-mail：jsyang@mail.csu.edu.cn
