基于特征線法網格劃分的主動土壓力有限元計算

方宏偉 ，譚卓英 ，方玲玲

(1.北京科技大學 金屬礦山高效開采與安全教育部重點實驗室，北京，100083；2.北京科技大學 土木與環(huán)境工程學院，北京，100083；3.蘇州大學 計算機科學與技術學院，江蘇 蘇州，215006)

有限元法在巖土工程中應用的關鍵是有限元模型的建立和本構關系的選取。有限元建模直接決定計算的精度和規(guī)模，網格劃分處于建模的中心位置[1]。提高建模精度和減少誤差的方法主要是加密網格和增加插值函數的階次，但網格加密到一定程度時，計算精度提高甚小而計算時間急劇增加[2]，插值函數采用高階次到一定程度時，精度增加亦不明顯[3]。另外，有限元解的光滑性是由建模對象的物理本質決定的[4]。通過地質模型建立有限元模型，實質是基于幾何形態(tài)建模[5]。因此，有限元建模應基于其對象的物理力學性質和幾何形態(tài)。目前通用的大型有限元軟件如Ansys等的前處理功能較差，雖開發(fā)有一些專門的網格劃分軟件如 HyperMesh[6]，但劃分網格前仍需作一定的幾何清理工作[7]，計算結果的精度很難確定。特征線法是基于極限平衡理論的一種數值方法，其實質是用差分方法求解滑移線方程組，詳細的公式和理論介紹參見文獻[8]。根據Mohr?Coulomb屈服準則，極限平衡狀態(tài)下土體中任一點有2條滑移線，且兩線與大主應力σ1的跡線夾角為μ=π/4?φ/2(φ為土體內摩擦角)，構成兩族滑移線，并將土體分割成有限單元體。因此，可用特征線法完成極限平衡狀態(tài)下土體的網格劃分。傳統(tǒng)的主動土壓力計算方法主要是Coulomb和Rankine 2個經典理論，為了更加符合工程實際，人們提出了包括對經典理論改進在內的一些新的方法[9]，彭明祥[10]建議采用有限元或滑移線理論求解。在此，本文作者運用特征線法完成土體網格劃分，將土體內部通過墻踵的 1條滑移線(曲線)作為滑裂面，加入邊界條件完成建模后應用位移有限元法計算節(jié)點應力，即將特征線法與有限元法耦合起來計算主動土壓力。采用 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及其近似平行線在內均屬第Ⅱ族滑移線，由特征線法的計算完成土體的網格劃分，得每一個節(jié)點的坐標(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上節(jié)點主應力σ1和σ3分布圖Fig.2 Distribution of principal stress σ1 and σ3 on node of boundary OD

2 有限元法求解主動土壓力

2.1 位移有限元法求解邊界節(jié)點應力

依據特征線法網格劃分結果進行有限元計算時，對于三角形，采用三角形單元，對于四邊形，由于其為非規(guī)則四邊形，故采用四節(jié)點四邊形等參數單元，應用平面應力問題的彈性矩陣進行計算。采用 Gauss數值積分法中的四點法對單元矩陣求解：

整體分析后，可得每個單元的應力向量{σ}e={σx，σy，τxy}e，并對邊界節(jié)點相關的單元應力向量取均值(精確到千分位)，得邊界9個節(jié)點應力向量。

2.2 本構模型參數對邊界節(jié)點應力和位移的影響

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。極限狀態(tài)下的應力水平取 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是反向影響。通過一個參數取極大值，其余參數取極小值的方法對節(jié)點位移和應力進行計算，結果表明：除參數G的變化對45號節(jié)點水平應力稍有正向影響以外，其他參數的變化對節(jié)點應力影響甚?。粎礙和Rf對節(jié)點位移敏感性強，結果見表1，與文獻[18]中結論一致。

由于參數G對υt有正向影響，同時為了綜合評價Et和υt對節(jié)點應力計算結果的影響，可分別按以上6個參數對Et和υt的正、負向影響對其取極大和極小值，考慮到 Ei及 Et對 υt亦有一定反向影響，故只取 Ei和Et極大、υt極小以及 Ei和 Et極小、υt極大兩組數據分析即可，計算結果見表2。由表2可知：本構模型參數對節(jié)點應力影響不大，與文獻[19]應用極限分析有限元法時本構模型參數對強度影響不大的結論一致。

表1 邊界節(jié)點位移水平分量與參數變化的關系Table1 Relationship between horizontal component of boundary nodes displacement and parameter variation cm

表2 邊界節(jié)點應力計算結果Table2 Calculated results of stress on boundary nodes kPa

2.3 主動土壓力的求解與驗證

由特征線法理論可知：在邊界OD上，當墻土摩擦角δ＞0o時，對土體穩(wěn)定性產生影響的是最大主應力σ1，應將其作為主動土壓力來求解。將邊界節(jié)點從下到上重新編號為0，1，2，3，4，5，6，7和8，由于每一點的yi為已知，故可計算邊界節(jié)點主動土壓力pa和總的主動土壓力Ea及其作用點到填土底部的距離y，計算結果見表3。

表3 主動土壓力計算結果Table3 Calculated results of active earth pressure

應注意到：當δ=0o時，η=90o，此時，σ3對土體穩(wěn)定性影響最大，故取

為Rankine土壓力理論表現形式。

為驗證本文耦合方法計算結果的可靠性，可用Coulomb解與其進行對比，Coulomb主動土壓力系數：

式中：ω為墻面與豎直線的夾角，在本例中ω=0o，則

從表3可知：該耦合方法的計算結果與Coulomb解接近，υt=0.49時兩者之間的絕對誤差為 4.101 2 kN/m，相對誤差為 0.72%，因此，其值是可靠的。2種計算方法的節(jié)點主動土壓力pa分布如圖3所示。由前面分析可知：從墻頂到墻踵節(jié)點位移逐漸變大，與文獻[20]中特征線法計算得到的速度場從上到下變大的結論一致，表明擋土墻產生繞頂轉動的模式。在此條件下，已有的試驗研究表明[21?22]：主動土壓力分布為上部大而下部小的拋物線性，且在淺處大于庫侖解，在深處小于庫侖解，更趨近于梯形分布，土壓力合力的作用點大于 0.33H。本文的土壓力合力作用點約為0.64H，由分布圖可見計算結果與實驗結果相符。

圖3 節(jié)點主動土壓力pa分布圖(υt=0.49)Fig.3 Distribution of active earth stress pa on nodes (υt=0.49)

3 結論

(1)運用特征線法進行有限元網格劃分，避免了劃分中人為因素，既可基于土體的幾何形態(tài)，又符合其物理力學本質，且計算精度可以得到保證；根據Coulomb解的對比和已有的實驗研究結果，表明本文的特征線網格劃分的主動土壓力有限元法計算結果可靠，方法可行。 本文的特征線網格劃分和有限元程序具有通用性，除了計算主動土壓力以外，還可用來分析極限平衡條件下均質土體地基和邊坡的穩(wěn) 定性。

(2)特征線網格劃分和有限元程序只適合均質土體極限平衡狀態(tài)下的強度計算，在計算主動土壓力時未考慮地下水的影響。因此，對復雜的邊界條件和非均質土體，在考慮滲流影響的前提下，開發(fā)滑移線法的網格劃分程序，并與現有成熟的大型商業(yè)有限元程序對接是今后的研究方向。

[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].北京: 中國建材工業(yè)出版社, 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]彭明祥.擋土墻主動土壓力的庫侖統(tǒng)一解[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].北京: 中國建筑工業(yè)出版社,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]肖大平, 朱維一, 陳環(huán).滑移線法求解極限承載力問題的一些進展[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]劉發(fā)前.圓形填土土壓力分布模式研究[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].北京: 機械工業(yè)出版社,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].基建優(yōu)化, 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.