楊釗，許建聰，余俊

（1.同濟大學巖土及地下工程教育部重點實驗室，上海 200092；2.中南大學土木建筑學院，湖南 長沙 410075）

彈性地基Timoshenko梁單元在ABAQUS軟件中的應用

楊釗1，許建聰1，余俊2

（1.同濟大學巖土及地下工程教育部重點實驗室，上海 200092；2.中南大學土木建筑學院，湖南 長沙 410075）

利用ABAQUS軟件中的自定義單元接口，采用Fortran語言開發(fā)彈性地基Timoshenko梁單元程序.通過與經(jīng)典解的比較，驗證所編彈性地基Timoshenko梁單元程序的準確性.該單元不僅考慮了地基彈簧受拉脫開的特點，而且還考慮了曲形梁單元內(nèi)結點不在一條直線上的特點.采用所編寫彈性地基梁單元分析青草沙原水過江輸水工程，計算規(guī)律與已建類似工程實測結果相同.

彈性地基；Timoshenko梁；ABAQUS軟件；單元子程序

在地下結構的計算領域，已有許多用于考慮結構-土體相互作用的計算方法，但是荷載結構法仍是目前使用最廣泛的一種方法.為了在研究周圍地層對結構的約束作用及地層對地下結構的反作用力時考慮地下結構變形，將Winker彈性地基理論引入到荷載結構法中［1］.同時，提出彈性地基上的Timoshenko梁，以考慮橫向剪切變形對厚梁的影響.目前，大量商品化軟件如MARC，SAP等均包含彈性地基梁單元，但均沒有考慮地基彈簧在受拉時的脫開情況.ANSYS，ABAQUS等大型通用有限元法計算分析軟件中未包括彈性地基梁單元，只能通過在梁單元結點上加彈簧單元來近似模擬彈性地基的作用［2-3］.這種做法只能在梁單元足夠小的情況下，才能近似等價于彈性地基作用［4］.本文以ABAQUS軟件為平臺，開發(fā)基于Winker地基理論與Timoshenko梁理論的地基梁單元.

1 Timoshenko梁單元剛度陣與荷載陣

設一Timoshenko梁的彎曲剛度為EI，剪切剛度為S，壓縮剛度為EA，長度為L，梁上有側(cè)向分布荷載q（x），徑向集中力Pi，集中力偶Mj，軸向集中力Fk.因此，其總勢能［5］為

Timoshenko梁單元的基本特點是撓度ω、軸向壓縮量u和截面轉(zhuǎn)動θ各自獨立插值.即

式（2）中：n為單元的結點數(shù)；Ni是Lagrange插值多項式.將式（2）代入式（1）中，由δΠ=0可以得到有限元的求解方程為

2 彈性地基梁

彈性地基上，梁在荷載作用下產(chǎn)生變形的同時，地基土也產(chǎn)生了變形.根據(jù)彈性地基的局部變形理論，地基土對地基梁的反力集度Pd與地基梁的撓度ω間的關系為

式（4）中：kd為地基反力系數(shù).為了考慮地基彈簧受拉脫開，假定梁的撓度值為負值，kd=0.考慮地基土的應變能后，彈性地基梁的總勢能比梁的總勢能多出一項［6］，即

式（5）中：b為梁截面的寬度.由式（2）可得

式（6）中：Bd=［Bd，1，…，Bd，n］；Bd，i=［Ni，0，0］.將式（6）代入式（5），取極值即可得單元剛度陣中地基剛度附加項Kde為

梁單元內(nèi)力計算方程為

式（9）中：N為軸力；Q為剪力；M為彎矩.

3 整體坐標系下的剛度陣

假定單元的方向在整體坐標系下的方向角為β.以單元內(nèi)任意一點i為例，有

式（11）中：Bi，gol=［Bi，gol，…，Bn，gol］，Bi，gol=Bλii.結合式（6），（7）可得

式（12）中：Bd，gol=［Bd，1，gol，…，Bd，n，gol］，Bd，i，gol=Bd，λii.整體坐標系下，彈性地基梁單元剛度矩陣為

當由前處理所得梁單元的內(nèi)部結點位于梁單元兩端結點的連線上，且內(nèi)部結點為梁單元的均分點時，由局部坐標與整體坐標的轉(zhuǎn)換關系可得

在實際工程中，特別是以曲梁形式存在的盾構隧道管片，由前處理軟件所得梁單元內(nèi)部結點位于梁兩端結點連線之外.針對這種情況，式（14）應為

式（15）中：k為單元在整體坐標系下的斜率.當k=∞時，有

結合式（13）～（16），可以得到在整體坐標系下單元剛度矩陣的表達式.

在數(shù)值積分中，若采用精確積分計算剪切變形能項，當梁很薄的情況下，約束條件γ=（dω/dx）-θ不可能在梁單元上處處滿足.采用精確積分計算剪切變形能項，將過分夸大剪切應變能的量級而造成剪切鎖死.為了避免剪切鎖死，對剪切變形項采用縮減積分計算［7］.

4 程序編制

ABAQUS提供用戶單元接口子程序UEL，用戶通過自定義UEL接口與求解器Standard的接口實現(xiàn)數(shù)據(jù)傳遞.UEL有其固定的書寫格式與規(guī)范，與主程序共享的變量必須在子程序開關予以定義，而主程序通過ABAQUS輸入文件（.inp）中的關鍵字“element，type=Un”來判斷是否使用自定義單元.

依據(jù)ABAQUS軟件二次開發(fā)的約定，用戶單元子程序UEL應至少包括5部分［8］，分別為：ABAQUS約定的子程序題名說明、ABAQUS定義的參數(shù)聲明表、用戶自定義的局部變量聲明表、用戶編寫的程序代碼段和子程序返回與結束語句等.在UEL中，用戶需要給出單元的結點數(shù)目、結點自由度、材料參數(shù)，通過主程序傳送給UEL的結點位移及結點位移增量更新單元應力，并最終將單元剛度矩陣（AMA TRX）及單元不平衡力矩陣（RHS）提供給ABAQUS主程序進行迭代求解.

4.1 程序流程

ABAQUS主程序進行迭代求解有如下8個步驟.

（1）計算坐標轉(zhuǎn)換矩陣λi.

（2）計算高斯積分點數(shù)與高斯積分點坐標.對于n結點Timensheno梁、軸力項、彎矩項剛度，可由n-1個高斯積分點精確求解；剪切項剛度需要n個高斯積分點才能精確求解.但為了避免剪切鎖死的發(fā)生，剪切項剛度在計算中也采用n-1個高斯積分點積分.

（3）計算形函數(shù)矩陣，并由式（14）～（16）計算積分系數(shù).

（4）由式（11）計算整體坐標系下梁單元的單元剛度矩陣.

（5）由式（9）計算積分點處的軸力、剪力、彎矩.

（6）將主程序傳入的位移值代入式（2），（10），計算地基梁單元在積分點處的撓度值.如計算撓度值為負值，則取此積分點處的地基剛度為零，再由式（12）計算地基對梁單元的附加剛度矩陣.

（7）由式（13）計算地基梁單元的整體剛度矩陣.

（8）計算單元的殘余力，并判斷收斂.如果不收斂，返回ABAQUS主程序進行第i+1次迭代.主程序?qū)⒏鶕?jù)UEL子程序第i次迭代所得到的單元剛度矩陣與殘余力項，計算位移增量與總位移量，然后跳到第（1）步進行UEL的第i+1次迭代計算.

4.2 程序驗證

梁荷載和彈性地基梁豎向位移圖，分別如圖1，2所示.兩端自由的彈性地基梁參數(shù)：長度l為10m，寬度b為0.5m，高度h為0.5m，梁身的彈性模量為10.0GPa，剪切模量為5.0GPa，地基的剛度系數(shù)K為4.0GN·m-3.求梁截面A，B和C的彎矩與撓度［9］.

采用自編彈性地基梁單元計算，將梁劃分為20個3結點地基梁單元.計算過程中，可得到每個單元內(nèi)積分點處的彎矩、剪力、軸力.將積分點處的內(nèi)力值外推，可以得到單元結點處的內(nèi)力值.

彎矩和撓度的理論解與數(shù)值解的對比，如表1所示.表1中：ωth，Mth分別為彎矩，撓度的理論解；ωc，Mc分別為彎矩，撓度的數(shù)值解；σ為相對誤差.由表1可知，理論解與數(shù)值解相差較小，該誤差可能源于有限元數(shù)值計算中網(wǎng)格的劃分、迭代收斂判斷準則，以及其他諸多綜合因素的影響.由此可見，編制的有限元程序是可靠的.

圖1 梁荷載示意圖Fig.1 Loads on the beam

圖2 彈性地基梁豎向位移圖Fig.2 Vertical displacement of the elastic foundation beam

表1 彎矩和撓度的理論解與數(shù)值解的對比Tab.1 Comparison of bending moment and deflection between analytical solution and numerical solution

5 工程實例計算

青草沙原水過江隧洞工程位于上海長江隧道下游約80m處，浦東側(cè)越江點在五號溝，長興島越江點在該島新開河附近，全長7.23km.越江輸水管道采用全斷面隧道掘進機（TBM）施工，有壓輸水，設計為圓形斷面，襯砌結構外直徑為6.8m，管片厚為480mm.考慮沖刷后上部垂直水土壓力為439.1kPa，上部水平土壓力為313.0kPa，下部水平土壓力為375.0kPa，隧道內(nèi)水壓力為404.1kPa.隧道周圍地層的地基剛度系數(shù)為10.0MN·m-3.采用3根彈簧分別模擬接頭的抗彎、抗壓與抗剪性能，其剛度系數(shù)分別為500（MN·m）·rad-1，5.0TN·m-1，0.5TN·m-1.

計算采用三結點Timeshenko地基梁.為簡化計算，在建立地基梁模型的同時也建立三結點Timeshenko梁模型.此梁模型與地基梁模型共結點且結點編號一致，但單元編號不同.將地基梁上的荷載施加到共結點的梁上，取共結點梁單元彈性模量為一極小數(shù).此時，梁單元的存在將簡化荷載的施加且對計算結果無影響.由于梁單元結點與地基梁單元共結點，因此，ABAQUS后處理中梁單元的位移場即為地基梁單元的位移場.

圖3 襯砌變位矢量圖Fig.3 Displacement vector diagram of lining

襯砌變位矢量圖，如圖3所示.在外周水土壓力與內(nèi)水壓力的聯(lián)合作用下，襯砌結構豎直方向內(nèi)縮，最大壓縮量為2.4mm；水平方向伸長，最大伸長量為2.2mm.襯砌結構的形狀由原先的圓形變成扁平的橢圓形.

襯砌截面的彎矩、軸力和剪力圖，如圖4～6所示.從圖4～6可知，管片彎矩的峰值出現(xiàn)在管頂、管底和兩腰，其管頂、管底為正值，兩腰彎矩為負值.最大正彎矩位于管頂處，其值為102.7kN·m；最大負彎矩位于管腰處，其值為-95.3kN·m.設軸力以受拉為正，受壓為負.管片大部分截面受壓，軸壓值管頂、管底小，而管腰大.軸壓最大值位于管腰處，其最大軸壓值為150.8kN；管頂部分截面受拉，最大軸拉值為-34.6kN.設剪力以截面呈順時針轉(zhuǎn)動為正，反之為負，剪力最大值為68.3kN，剪力最小值為-59.7kN.

圖4 襯砌截面彎矩圖Fig.4 Moment of lining section

圖5 襯砌截面軸力圖Fig.5 Axial force of lining section

圖6 襯砌截面剪力圖Fig.6 Shear force of lining section

6 結束語

針對絕大多數(shù)商品化軟件不能考慮彈性地基梁在受拉時地基彈簧脫開的不足，基于ABAQUS軟件平臺開發(fā)了彈性地基Timeshenko梁單元.該單元不僅考慮了地基彈簧受拉脫開的特點，而且還考慮了曲形梁單元內(nèi)結點不在一條直線上的特點.采用所編寫彈性地基梁單元分析青草沙原水過江輸水工程，計算規(guī)律與已建類似工程實測結果相同.研究結果表明，所編寫的彈性地基梁單元精度高，可供實際工程計算應用.

［1］孫鈞.地下結構［M］.北京：科技出版社，1987.

［2］黃群賢，林建華.液化側(cè)擴地基中樁基的有限元分析［J］.華僑大學學報：自然科學版，2004，25（3）：328-330.

［3］賈瑞華，陽軍生，馬濤，等.既有管線下盾構施工地層沉降監(jiān)測和位移加載數(shù)值分析［J］.巖土工程學報，2009，31（3）：425-430.

［4］楊釗，潘曉明，余俊.盾構輸水隧洞復合襯砌計算模型研究［C］//2009全國土木工程博士生學術論壇優(yōu)秀論文集.長沙：中南大學出版社，2009.

［5］王勖成.有限單元法［M］.北京：清華大學出版社，2006.

［6］朱伯芳.有限單元法原理與應用［M］.北京：中國水利水電出版社，2004.

［7］馮紫良.桿系結構的計算機分析［M］.上海：同濟大學出版社，1991.

［8］葉志才，徐磊，王超.基于ABAQUS的三維錨桿單元的開發(fā)［J］.三峽大學學報，2008，30（5）：29-32.

［9］黃義，何芳杜.彈性地基上的梁板殼［M］.北京：科學出版社.

Application of Elastic Foundation Timoshenko Beam Element in ABAQUS

YAN G Zhao1，XU Jian-cong1，YU Jun2
（1.Key Laboratory of Geotechnical and Underground Engineering of Ministry of Education，Tongji University，Shanghai 200092，China；2.School of Civil Engineering and Architecture，Central South University，Changsha 410075，China）

Using the user defined element interface in ABAQUS，the elastic foundation Timoshenko beam element is developed with Fortran language.Comparing with classic analytical solution，the result indicates that the accuracy of the element is high enough.That element not only considers the separation of the foundation spring under tension，but also considers the internal node of curve beam element is not on the line of two end nodes.Using the elastic foundation beam element，the Qingcaosha river-cross water diversion project is analyzed，the results agrees with the measured data of similar projects.

elastic foundation；Timoshenko beam；ABAQUS；element subroutine

TU 471.2

A

1000-5013（2010）04-0448-05

（責任編輯：黃仲一 英文審校：方德平）

2009-09-19

楊釗（1984-），男，博士研究生，主要從事盾構隧道數(shù)值計算和模型實驗的研究.E-mail：yangzhaolp@126.com.

國家自然科學基金資助項目（40872179）；中國博士后科研基金資助項目（20080440652）