方宏偉 ，譚卓英 ，方玲玲

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

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

1 特征線法劃分土體網(wǎng)格

1.1 邊界條件

某擋土墻后的填土主要為砂土。為了說明特征線劃分網(wǎng)格的方法，將邊界條件進行簡化，如圖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 網(wǎng)格劃分

設(shè)填土頂部為坐標原點，x軸為水平方向，y軸為豎直方向，如圖1所示，對OA取8個等分點，設(shè)A點、OB線、OC線及其近似平行線和D點屬于第Ⅰ族滑移線，包括O點、曲線ABPCD及其近似平行線在內(nèi)均屬第Ⅱ族滑移線，由特征線法的計算完成土體的網(wǎng)格劃分，得每一個節(jié)點的坐標(x, y)。由于誤差累積，計算得=13.94 m。同時可知邊界OD上最大主應(yīng)力 σ1與 正 x 軸 夾 角 為 θ=β+(Δ?δ)/2=95o， 式 中Δ=arcsin(sin δ/sin φ)=20o，則最大主應(yīng)力 σ1與負 x軸(豎直邊界外法線)夾角為η=π?θ=85o，如圖2所示。

為了便于有限元分析，將原點重新設(shè)置在底部D點，即將滑裂土體放在第一象限內(nèi)，進行坐標變換得到新的坐標(x，13.94?y)，同時從下到上、從左到右按逆時針進行編號，共88個網(wǎng)格，可知特征線法的計算誤差在 2%～5%之間[13]。對于非均質(zhì)土也可采用滑移線法[14]劃分網(wǎng)格。

圖1 有限元網(wǎng)格劃分分布圖Fig.1 Distribution of finite element mesh

圖2 邊界OD上節(jié)點主應(yīng)力σ1和σ3分布圖Fig.2 Distribution of principal stress σ1 and σ3 on node of boundary OD

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

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

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

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

2.2 本構(gòu)模型參數(shù)對邊界節(jié)點應(yīng)力和位移的影響

Duncan?Zhang (E?υ)非線性彈性模型共有8個參數(shù)，其中φ和c為已知，參考文獻[15?16]給出其余6個參數(shù)工程經(jīng)驗取值范圍為K=300～1 000，n=0.3～0.6，Rf=0.6～0.85，G=0.2～0.6，F(xiàn)=0.1～0.2，D=1～20。極限狀態(tài)下的應(yīng)力水平取 s=0.95[17]，并由最大主應(yīng)力σ1=γH+q=253 kPa ， 求 得 最 小 主 應(yīng) 力 為 σ3=σ1(1?0.9sin φ)/(1+sin φ)=93 kPa，按公式可計算彈性模量和泊松比初始值 Ei，υi和切線值 Et，υt(υt＞0.49 時取0.49)，并可知K，G，F(xiàn)和D對Et和υt是正向影響，n和Rf對Et是反向影響。通過一個參數(shù)取極大值，其余參數(shù)取極小值的方法對節(jié)點位移和應(yīng)力進行計算，結(jié)果表明：除參數(shù)G的變化對45號節(jié)點水平應(yīng)力稍有正向影響以外，其他參數(shù)的變化對節(jié)點應(yīng)力影響甚??；參數(shù)K和Rf對節(jié)點位移敏感性強，結(jié)果見表1，與文獻[18]中結(jié)論一致。

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

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

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

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

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

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

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

為Rankine土壓力理論表現(xiàn)形式。

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

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

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

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

3 結(jié)論

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

(2)特征線網(wǎng)格劃分和有限元程序只適合均質(zhì)土體極限平衡狀態(tài)下的強度計算，在計算主動土壓力時未考慮地下水的影響。因此，對復(fù)雜的邊界條件和非均質(zhì)土體，在考慮滲流影響的前提下，開發(fā)滑移線法的網(wǎng)格劃分程序，并與現(xiàn)有成熟的大型商業(yè)有限元程序?qū)邮墙窈蟮难芯糠较颉?/p>

[1]于亞婷, 杜平安, 王振偉.有限元法的應(yīng)用現(xiàn)狀研究[J].機械設(shè)計, 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]杜平安.有限元網(wǎng)格劃分的基本原則[J].機械設(shè)計與制造,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]陳傳淼.有限元超收斂構(gòu)造理論[M].長沙: 湖南科學(xué)技術(shù)出版社, 2001: 205.CHEN Chuan-miao.Structure theory of super convergence of finite elements[M].Changsha: Hunan Science and Technology Press, 2001: 205.

[5]李新星, 朱合華, 蔡永昌, 等.基于三維地質(zhì)模型的巖土工程有限元自動建模方法[J].巖土工程學(xué)報, 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軟件的復(fù)雜地質(zhì)工程有限元建模[J].長江科學(xué)院院報, 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].機械設(shè)計與制造, 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]陳震.散體極限平衡理論基礎(chǔ)[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].巖土力學(xué), 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]應(yīng)宏偉, 蔣波, 謝康和.條形荷載下?lián)跬翂χ鲃油翂毫τ嬎鉡J].巖土力學(xué), 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].巖土工程學(xué)報, 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].上海: 上海交通大學(xué)土木工程系, 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]廖紅建.巖土工程數(shù)值分析[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模型應(yīng)用的若干關(guān)鍵技術(shù)[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]尹蓉蓉, 朱合華.鄧肯—張模型參數(shù)敏感性分析[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]鄭穎人, 趙尚毅.巖土工程極限分析有限元法及其應(yīng)用[J].土木工程學(xué)報, 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].江蘇科技大學(xué)學(xué)報: 自然科學(xué)版, 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]周應(yīng)英, 任美龍.剛性擋土墻主動土壓力的試驗研究[J].巖土工程學(xué)報, 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.