方宏偉 ，譚卓英 ，方玲玲

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

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

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

1.1 邊界條件

某擋土墻后的填土主要為砂土。為了說(shuō)明特征線劃分網(wǎng)格的方法，將邊界條件進(jìn)行簡(jiǎn)化，如圖1所示。由圖1可知：墻面垂直，填土面水平，即兩者夾角β=90o，填土深H==13.36 m，γ=17 kN/m3，φ=30o，c=0o(c的取值不影響兩族滑移線的形狀和夾角[11])，按Coulomb理論取土體與擋土墻摩擦角 δ=φ/3=10o，μ=π/4?φ/2=30o，填土表面均布條形荷載q=15 kPa，且全部作用于土楔體上，主動(dòng)滑裂面傾角λ的計(jì)算如下[12]：

1.2 網(wǎng)格劃分

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

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

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

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

2 有限元法求解主動(dòng)土壓力

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

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

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

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

Duncan?Zhang (E?υ)非線性彈性模型共有8個(gè)參數(shù)，其中φ和c為已知，參考文獻(xiàn)[15?16]給出其余6個(gè)參數(shù)工程經(jīng)驗(yàn)取值范圍為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，按公式可計(jì)算彈性模量和泊松比初始值 Ei，υi和切線值 Et，υt(υt＞0.49 時(shí)取0.49)，并可知K，G，F(xiàn)和D對(duì)Et和υt是正向影響，n和Rf對(duì)Et是反向影響。通過(guò)一個(gè)參數(shù)取極大值，其余參數(shù)取極小值的方法對(duì)節(jié)點(diǎn)位移和應(yīng)力進(jìn)行計(jì)算，結(jié)果表明：除參數(shù)G的變化對(duì)45號(hào)節(jié)點(diǎn)水平應(yīng)力稍有正向影響以外，其他參數(shù)的變化對(duì)節(jié)點(diǎn)應(yīng)力影響甚?。粎?shù)K和Rf對(duì)節(jié)點(diǎn)位移敏感性強(qiáng)，結(jié)果見(jiàn)表1，與文獻(xiàn)[18]中結(jié)論一致。

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

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

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

2.3 主動(dòng)土壓力的求解與驗(yàn)證

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

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

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

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

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

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

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

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

3 結(jié)論

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

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

[1]于亞婷, 杜平安, 王振偉.有限元法的應(yīng)用現(xiàn)狀研究[J].機(jī)械設(shè)計(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]杜平安.有限元網(wǎng)格劃分的基本原則[J].機(jī)械設(shè)計(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].機(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].長(zhǎng)沙: 湖南科學(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ì)模型的巖土工程有限元自動(dòng)建模方法[J].巖土工程學(xué)報(bào), 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]張萍, 鄭東健, 張獻(xiàn)法.基于Hyper Mesh軟件的復(fù)雜地質(zhì)工程有限元建模[J].長(zhǎng)江科學(xué)院院報(bào), 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]劉榮軍, 吳新躍, 鄭建華.有限元建模中的幾何清理問(wèn)題[J].機(jī)械設(shè)計(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]陳震.散體極限平衡理論基礎(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]顧慰慈.擋土墻土壓力計(jì)算手冊(cè)[M].北京: 中國(guó)建材工業(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]彭明祥.擋土墻主動(dòng)土壓力的庫(kù)侖統(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]龔曉南.土工計(jì)算機(jī)分析[M].北京: 中國(guó)建筑工業(yè)出版社,2000: 198.GONG Xiao-nan.Geotechnical computer analysis[M].Beijing:China Architecture and Building Press, 2000: 198.

[12]應(yīng)宏偉, 蔣波, 謝康和.條形荷載下?lián)跬翂χ鲃?dòng)土壓力計(jì)算[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).滑移線法求解極限承載力問(wèn)題的一些進(jìn)展[J].巖土工程學(xué)報(bào), 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].北京: 中國(guó)水利水電出版社, 2007:228.YIN Zong-ze.Earthwork principle[M].Beijing: China Water Conservancy and Hydropower Press, 2007: 228.

[16]廖紅建.巖土工程數(shù)值分析[M].北京: 機(jī)械工業(yè)出版社,2006: 47, 90?140.LIAO Hong-jian.Numerical analysis of geotechnical engineering[M].Beijing: China Machine Press, 2006: 47,90?140.

[17]潘樹(shù)來(lái), 王全鳳, 涂帆.土體破壞時(shí) 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é)報(bào), 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]姜朋明, 尚羽, 陸長(zhǎng)鋒.基于滑移線法擋土墻土壓力問(wèn)題的討論[J].江蘇科技大學(xué)學(xué)報(bào): 自然科學(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)英, 任美龍.剛性擋土墻主動(dòng)土壓力的試驗(yàn)研究[J].巖土工程學(xué)報(bào), 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.