基于特征线法网格划分的主动土压力有限元计算

方宏伟^{1, 2}，谭卓英^{1, 2}，方玲玲³

(1. 北京科技大学 金属矿山高效开采与安全教育部重点实验室，北京，100083；

2. 北京科技大学 土木与环境工程学院，北京，100083；

3. 苏州大学 计算机科学与技术学院，江苏 苏州，215006)

摘要：针对有限元法网格划分存在人为因素且不能保证精度的不足，提出特征线法划分网格的新方法。假定通过墙踵的滑移线为滑裂面，采用Duncan-Zhang(E-υ)模型，将其运用于主动土压力的有限元分析计算。研究结果表明：本构模型参数对边界节点应力影响不大，但实验常数K和破坏比R_f对其位移有影响；当墙土摩擦角δ＞0°时，最大主应力σ₁对边界的土体稳定性影响最大，宜将其作为主动土压力p_a进行计算；当δ=0°时，取最小主应力σ₁=p_a，此时为Rankine理论的表现形式；绘制节点p_a分布图，求得主动土压力合力E_a；以底部计算点为矩心，用合力矩定理求得E_a作用点的位置y。经过与Coulomb解的对比和已有实验结论的验证，表明计算结果可靠，方法可行。

关键词：特征线法；网格划分；主动土压力；有限元法

中图分类号：TU43         文献标志码：A         文章编号：1672-7207(2011)09-2825-06

Characteristic line mesh division for active earth pressure calculation in finite-element method

FANG Hong-wei^{1, 2}, TAN Zhuo-ying^{1, 2}, FANG Ling-ling³

(1. State Key Laboratory of High-Efficient Mining and Safety of Metal Mines of Ministry of Education,

University of Science and Technology Beijing, Beijing 100083, China;

2. School of Civil and Environmental Engineering, University of Science and Technology Beijing, Beijing 100083, China;

3. College of Computer Science and Technology, Soochow University, Suzhou 215006, China)

Abstract: In order to overcome the disadvantages of human factors and less precision in mesh division, a new approach of characteristic curve mesh generation was established. Using the hypothesis of glide lines crossing at the wall-heel as a slip crack surface, the Duncan-Chang (E-υ) nonlinear elastic model was used for calculation of active earth pressure in finite-element method. The results show that constitutive model parameters play a minor role in boundary node stresses but experimental constant K and break ratio R_f have a larger impact on displacement. The maximum principal stress σ₁ plays the biggest influence in stability of earth mass boundary when the internal friction angle of earth mass δ is more than 0° and it is appropriate to take σ₁ as active earth pressure p_a. The minimum principal stress σ₃ equals p_a and has an expression of Rankine’s in case of δ=0°. From the distribution of p_a on nodes, the total active earth pressure E_a can be decided and the corresponding position y of application point of E_a can be determined by resultant moment theorem at the bottom of earth mass. The results are compared with those from Coulomb’s earth pressure theory and verified by the existing ones from experiments, indicating that the calculated precision is realizable and the new method is feasible.

Key words: method of characteristic curve; mesh division; active earth pressure; finite-element method

有限元法在岩土工程中应用的关键是有限元模型的建立和本构关系的选取。有限元建模直接决定计算的精度和规模，网格划分处于建模的中心位置^[1]。提高建模精度和减少误差的方法主要是加密网格和增加插值函数的阶次，但网格加密到一定程度时，计算精度提高甚小而计算时间急剧增加^[2]，插值函数采用高阶次到一定程度时，精度增加亦不明显^[3]。另外，有限元解的光滑性是由建模对象的物理本质决定的^[4]。通过地质模型建立有限元模型，实质是基于几何形态建模^[5]。因此，有限元建模应基于其对象的物理力学性质和几何形态。目前通用的大型有限元软件如Ansys等的前处理功能较差，虽开发有一些专门的网格划分软件如HyperMesh^[6]，但划分网格前仍需作一定的几何清理工作^[7]，计算结果的精度很难确定。特征线法是基于极限平衡理论的一种数值方法，其实质是用差分方法求解滑移线方程组，详细的公式和理论介绍参见文献[8]。根据Mohr-Coulomb屈服准则，极限平衡状态下土体中任一点有2条滑移线，且两线与大主应力σ₁的迹线夹角为μ=π/4-φ/2(φ为土体内摩擦角)，构成两族滑移线，并将土体分割成有限单元   体。因此，可用特征线法完成极限平衡状态下土体的网格划分。传统的主动土压力计算方法主要是Coulomb和Rankine 2个经典理论，为了更加符合工程实际，人们提出了包括对经典理论改进在内的一些新的方法^[9]，彭明祥^[10]建议采用有限元或滑移线理论求解。在此，本文作者运用特征线法完成土体网格划分，将土体内部通过墙踵的1条滑移线(曲线)作为滑裂面，加入边界条件完成建模后应用位移有限元法计算节点应力，即将特征线法与有限元法耦合起来计算主动土压力。采用Duncan-Zhang(E-υ)非线性弹性本构模型及总应力分析法，通过VisualC++6.0编制计算程序分析求解，依据Coulomb解和已有的试验结论对计算结果进行验证。

1  特征线法划分土体网格

1.1  边界条件

某挡土墙后的填土主要为砂土。为了说明特征线划分网格的方法，将边界条件进行简化，如图1所示。由图1可知：墙面垂直，填土面水平，即两者夹角β=90^o，填土深H==13.36 m，γ=17 kN/m³，φ=30^o，c=0^o(c的取值不影响两族滑移线的形状和夹角^[11])，按Coulomb理论取土体与挡土墙摩擦角δ=φ/3=10^o，μ= π/4-φ/2=30^o，填土表面均布条形荷载q=15 kPa，且全部作用于土楔体上，主动滑裂面倾角λ的计算如下^[12]：

     (1)

式中：为墙面与填土底水平面的夹角；为填土上表面与水平面的夹角，本例中为=90^o，=0^o。可得λ=30^o，由正切函数得水平宽度=   8 m。

1.2  网格划分

设填土顶部为坐标原点，x轴为水平方向，y轴为竖直方向，如图1所示，对OA取8个等分点，设A点、OB线、OC线及其近似平行线和D点属于第Ⅰ族滑移线，包括O点、曲线ABPCD及其近似平行线在内均属第Ⅱ族滑移线，由特征线法的计算完成土体的网格划分，得每一个节点的坐标(x, y)。由于误差累积，计算得=13.94 m。同时可知边界OD上最大主应力σ₁与正x轴夹角为θ=β+(Δ-δ)/2=95^o，式中Δ=arcsin(sin δ/sin φ)=20^o，则最大主应力σ₁与负x轴(竖直边界外法线)夹角为η=π-θ=85^o，如图2所示。

为了便于有限元分析，将原点重新设置在底部D点，即将滑裂土体放在第一象限内，进行坐标变换得到新的坐标(x，13.94-y)，同时从下到上、从左到右按逆时针进行编号，共88个网格，可知特征线法的计算误差在2%~5%之间^[13]。对于非均质土也可采用滑移线法^[14]划分网格。

图1  有限元网格划分分布图

Fig.1  Distribution of finite element mesh

图2  边界OD上节点主应力σ₁和σ₃分布图

Fig.2  Distribution of principal stress σ₁ and σ₃ on node of boundary OD

2  有限元法求解主动土压力

2.1  位移有限元法求解边界节点应力

依据特征线法网格划分结果进行有限元计算时，对于三角形，采用三角形单元，对于四边形，由于其为非规则四边形，故采用四节点四边形等参数单元，应用平面应力问题的弹性矩阵进行计算。采用Gauss数值积分法中的四点法对单元矩阵求解：

     (2)

式中：ξ₁= η₁=0.577 35；ξ₂= η₂=-0.577 35；α₁=α₂=1.0。

整体分析后，可得每个单元的应力向量{σ}^e=  {σ_x，σ_y，τ_xy}^e，并对边界节点相关的单元应力向量取均值(精确到千分位)，得边界9个节点应力向量。

2.2  本构模型参数对边界节点应力和位移的影响

Duncan-Zhang (E-υ)非线性弹性模型共有8个参数，其中φ和c为已知，参考文献[15-16]给出其余6个参数工程经验取值范围为K=300~1 000，n=0.3~0.6，R_f =0.6~0.85，G=0.2~0.6，F=0.1~0.2，D=1~20。极限状态下的应力水平取s=0.95^[17]，并由最大主应力σ₁=γH+q=253 kPa，求得最小主应力为σ₃=  σ₁(1-0.9sin φ)/(1+sin φ)=93 kPa，按公式可计算弹性模量和泊松比初始值E_i，υ_i和切线值E_t，υ_t(υ_t＞0.49时取0.49)，并可知K，G，F和D对E_t和υ_t是正向影响，n和R_f对E_t是反向影响。通过一个参数取极大值，其余参数取极小值的方法对节点位移和应力进行计算，结果表明：除参数G的变化对45号节点水平应力稍有正向影响以外，其他参数的变化对节点应力影响甚小；参数K和R_f对节点位移敏感性强，结果见表1，与文献[18]中结论一致。

由于参数G对υ_t有正向影响，同时为了综合评价E_t和υ_t对节点应力计算结果的影响，可分别按以上6个参数对E_t和υ_t的正、负向影响对其取极大和极小值，考虑到E_i及E_t对υ_t亦有一定反向影响，故只取E_i和E_t极大、υ_t极小以及E_i和E_t极小、υ_t极大两组数据分析即可，计算结果见表2。由表2可知：本构模型参数对节点应力影响不大，与文献[19]应用极限分析有限元法时本构模型参数对强度影响不大的结论一致。

表1  边界节点位移水平分量与参数变化的关系

     Table 1  Relationship between horizontal component of boundary nodes displacement and parameter variation    cm

表2  边界节点应力计算结果

                           Table 2  Calculated results of stress on boundary nodes                         kPa

2.3  主动土压力的求解与验证

由特征线法理论可知：在边界OD上，当墙土摩擦角δ＞0^o时，对土体稳定性产生影响的是最大主应力σ₁，应将其作为主动土压力来求解。将边界节点从下到上重新编号为0，1，2，3，4，5，6，7和8，由于每一点的y_i为已知，故可计算边界节点主动土压力p_a和总的主动土压力E_a及其作用点到填土底部的距离y，计算结果见表3。

    (3)

         (4)

         (5)

表3  主动土压力计算结果

Table 3  Calculated results of active earth pressure

应注意到：当δ=0^o时，η=90^o，此时，σ₃对土体稳定性影响最大，故取

     (6)

为Rankine土压力理论表现形式。

为验证本文耦合方法计算结果的可靠性，可用Coulomb解与其进行对比，Coulomb主动土压力系数：

式中：ω为墙面与竖直线的夹角，在本例中ω=0^o，则

。

从表3可知：该耦合方法的计算结果与Coulomb解接近，υ_t=0.49时两者之间的绝对误差为4.101 2 kN/m，相对误差为0.72%，因此，其值是可靠的。2种计算方法的节点主动土压力p_a分布如图3所示。由前面分析可知：从墙顶到墙踵节点位移逐渐变大，与文献[20]中特征线法计算得到的速度场从上到下变大的结论一致，表明挡土墙产生绕顶转动的模式。在此条件下，已有的试验研究表明^[21-22]：主动土压力分布为上部大而下部小的抛物线性，且在浅处大于库仑解，在深处小于库仑解，更趋近于梯形分布，土压力合力的作用点大于0.33H。本文的土压力合力作用点约为0.64H，由分布图可见计算结果与实验结果相符。

图3  节点主动土压力p_a分布图(υ_t=0.49)

Fig.3  Distribution of active earth stress p_a on nodes (υ_t=0.49)

3  结论

(1) 运用特征线法进行有限元网格划分，避免了划分中人为因素，既可基于土体的几何形态，又符合其物理力学本质，且计算精度可以得到保证；根据Coulomb解的对比和已有的实验研究结果，表明本文的特征线网格划分的主动土压力有限元法计算结果可靠，方法可行。 本文的特征线网格划分和有限元程序具有通用性，除了计算主动土压力以外，还可用来分析极限平衡条件下均质土体地基和边坡的稳 定性。

(2) 特征线网格划分和有限元程序只适合均质土体极限平衡状态下的强度计算，在计算主动土压力时未考虑地下水的影响。因此，对复杂的边界条件和非均质土体，在考虑渗流影响的前提下，开发滑移线法的网格划分程序，并与现有成熟的大型商业有限元程序对接是今后的研究方向。

参考文献：

[1] 于亚婷, 杜平安, 王振伟. 有限元法的应用现状研究[J]. 机械设计, 2005, 22(3): 6-8.
YU Ya-ting, DU Ping-an, WANG Zheng-wei. Research on the current application status of finite element method[J]. Journal of Machine Design, 2005, 22(3): 6-8.

[2] 杜平安. 有限元网格划分的基本原则[J]. 机械设计与制造, 2000(1): 34-36.
DU Ping-an. The basic principle for finite element mesh generation[J]. Machinery Design and Manufacture, 2000(1): 34-36.

[3] 杜平安. 建立有限元模型的基本原则[J]. 机械与电子, 2001(4): 40-42.
DU Ping-an. The basic principle for creation of finite element model[J]. Machinery and Electronics, 2001(4): 40-42.

[4] 陈传淼. 有限元超收敛构造理论[M]. 长沙: 湖南科学技术出版社, 2001: 205.
CHEN Chuan-miao. Structure theory of super convergence of finite elements[M]. Changsha: Hunan Science and Technology Press, 2001: 205.

[5] 李新星, 朱合华, 蔡永昌, 等. 基于三维地质模型的岩土工程有限元自动建模方法[J]. 岩土工程学报, 2008, 30(6): 855-862.
LI Xin-xing, ZHU He-hua, CAI Yong-chang, et al. Automatic modeling method of numerical analysis in geotechnical engineering based on 3D geologic model[J]. Chinese Journal of Geotechnical Engineering, 2008, 30(6): 855-862.

[6] 张萍, 郑东健, 张献法. 基于Hyper Mesh软件的复杂地质工程有限元建模[J]. 长江科学院院报, 2008, 25(2): 84-86.
ZHANG Ping, ZHENG Dong-jian, ZHANG Xian-fa. Finite element model of complex geologic engineering based on Hyper Mesh software[J]. Journal of Yangtze River Scientific Research Institute, 2008, 25(2): 84-86.

[7] 刘荣军, 吴新跃, 郑建华. 有限元建模中的几何清理问题[J]. 机械设计与制造, 2005(9): 145-147.
LIU Rong-jun, WU Xin-yue, ZHENG Jian-hua. The problem of geometry clean in constructing finite element model[J]. Machinery Design and Manufacture, 2005(9): 145-147.

[8] 陈震. 散体极限平衡理论基础[M]. 北京: 水利电力出版社, 1987: 11-36, 104-123.
CHEN Zhen. Granular materials limit equilibrium theory foundation[M]. Beijing: Hydraulic Press, 1987: 11-36, 104-123.

[9] 顾慰慈. 挡土墙土压力计算手册[M]. 北京: 中国建材工业出版社, 2004: 263-290, 473-490.
GU Wei-ci. Calculation of earth pressure on retaining wall[M]. Beijing: China Material Industry Press, 2004: 263-290, 473-490.

[10] 彭明祥. 挡土墙主动土压力的库仑统一解[J]. 岩土力学, 2009, 30(2): 379-386.
PENG Ming-xiang. Coulomb’s unified solution of active earth pressure on retaining wall[J]. Rock and Soil Mechanics, 2009, 30(2): 379-386.

[11] 龚晓南. 土工计算机分析[M]. 北京: 中国建筑工业出版社, 2000: 198.
GONG Xiao-nan. Geotechnical computer analysis[M]. Beijing: China Architecture and Building Press, 2000: 198.

[12] 应宏伟, 蒋波, 谢康和. 条形荷载下挡土墙主动土压力计算[J]. 岩土力学, 2007, 28(增刊): 183-186.
YING Hong-wei, JIANG-Bo, XIE Kang-he. Active earth pressure on retaining wall due to strip surcharge[J]. Rock and Soil Mechanics, 2007, 28(S): 183-186.

[13] 肖大平, 朱维一, 陈环. 滑移线法求解极限承载力问题的一些进展[J]. 岩土工程学报, 1998, 20(4): 25-29.
XIAO Da-ping, ZHU Wei-yi, CHEN Huan. Progress in slip lines method to solve the bearing capacity problem[J]. Chinese Journal of Geotechnical Engineering, 1998, 20(4): 25-29.

[14] 刘发前. 圆形填土土压力分布模式研究[D]. 上海: 上海交通大学土木工程系, 2008: 14-15, 76-92.
LIU Fa-qian. Study on the distribution of the earth pressure of circular pit[D]. Shanghai: Shanghai Jiaotong University. Department of Civil Engineering, 2008: 14-15, 76-92.

[15] 殷宗泽. 土工原理[M]. 北京: 中国水利水电出版社, 2007: 228.
YIN Zong-ze. Earthwork principle[M]. Beijing: China Water Conservancy and Hydropower Press, 2007: 228.

[16] 廖红建. 岩土工程数值分析[M]. 北京: 机械工业出版社, 2006: 47, 90-140.
LIAO Hong-jian. Numerical analysis of geotechnical engineering[M]. Beijing: China Machine Press, 2006: 47, 90-140.

[17] 潘树来, 王全凤, 涂帆. 土体破坏时Duncan-Zhang模型应用的若干关键技术[J]. 基建优化, 2007, 28(6): 157-160.
PAN Shu-lai, WANG Quan-feng, TU Fan. Some key techniques in applying Duncan-Chang model to analyze soil failure[J]. Optimization of Capital Construction, 2007, 28(6): 157-160.

[18] 尹蓉蓉, 朱合华. 邓肯—张模型参数敏感性分析[J]. 地下空间, 2004, 24(4): 434-437.
YIN Rong-rong, ZHU He-hua. Analysis of parameters sensitivity of Duncan-Chang model[J]. Underground Space, 2004, 24(4): 434-437.

[19] 郑颖人, 赵尚毅. 岩土工程极限分析有限元法及其应用[J]. 土木工程学报, 2005, 38(1): 91-98.
ZHENG Ying-ren, ZHAO Shang-yi. Limit state finite element method for geotechnical engineering analysis and its application[J]. China Civil Engineering Journal, 2005, 38(1): 91-98.

[20] 姜朋明, 尚羽, 陆长锋. 基于滑移线法挡土墙土压力问题的讨论[J]. 江苏科技大学学报: 自然科学版, 2008, 22(6): 9-12.
JIANG Peng-ming, SHANG Yu, LU Chang-feng. Discussion on soil pressure against retaining wall based on slip line method[J]. Journal of Jiangsu University of Science and Technology: Natural Science Edition, 2008, 22(6): 9-12.

[21] 周应英, 任美龙. 刚性挡土墙主动土压力的试验研究[J]. 岩土工程学报, 1990, 12(2): 19-26.
ZHOU Ying-ying, REN Mei-long. An experiment study on active earth pressure behind rigid retaining wall[J]. Chinese Journal of Geotechnical Engineering, 1990, 12(2): 19-26.

[22] Fang Y S, Ishibashi I. Static earth pressure with various wall movements[J]. Journal of Geotechnical Engineering, 1986, 112(3): 317-333.

(编辑 杨幼平)

收稿日期：2010-09-13；修回日期：2010-12-19

基金项目：国家重点基础研究发展规划(“973”计划)项目(2010CB731501)；国家自然科学基金资助项目(51174013)

通信作者：谭卓英(1965-)，男，湖南涟源人，教授，博士生导师，从事岩土工程、采矿工程及防灾减灾研究；电话：010-62332939；E-mail: markzhy_tan@163.com