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

方宏伟 ，谭卓英 ，方玲玲

(1.北京科技大学 金属矿山高效开采与安全教育部重点实验室，北京，100083；2.北京科技大学 土木与环境工程学院，北京，100083；3.苏州大学 计算机科学与技术学院，江苏 苏州，215006)

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

1 特征线法划分土体网格

1.1 边界条件

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

1.2 网格划分

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

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

图1 有限元网格划分分布图Fig.1 Distribution of finite element mesh

图2 边界OD上节点主应力σ1和σ3分布图Fig.2 Distribution of principal stress σ1 and σ3 on node of boundary OD

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

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

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

整体分析后，可得每个单元的应力向量{σ}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，Rf=0.6～0.85，G=0.2～0.6，F=0.1～0.2，D=1～20。极限状态下的应力水平取 s=0.95[17]，并由最大主应力σ1=γH+q=253 kPa ， 求 得 最 小 主 应 力 为 σ3=σ1(1−0.9sin φ)/(1+sin φ)=93 kPa，按公式可计算弹性模量和泊松比初始值 Ei，υi和切线值 Et，υt(υt＞0.49 时取0.49)，并可知K，G，F和D对Et和υt是正向影响，n和Rf对Et是反向影响。通过一个参数取极大值，其余参数取极小值的方法对节点位移和应力进行计算，结果表明：除参数G的变化对45号节点水平应力稍有正向影响以外，其他参数的变化对节点应力影响甚小；参数K和Rf对节点位移敏感性强，结果见表1，与文献[18]中结论一致。

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

表1 边界节点位移水平分量与参数变化的关系Table1 Relationship between horizontal component of boundary nodes displacement and parameter variation cm

表2 边界节点应力计算结果Table2 Calculated results of stress on boundary nodes kPa

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

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

表3 主动土压力计算结果Table3 Calculated results of active earth pressure

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

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

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

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

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

图3 节点主动土压力pa分布图(υt=0.49)Fig.3 Distribution of active earth stress pa 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 andSoil 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.