彈性地基Timoshenko梁?jiǎn)卧贏BAQUS軟件中的應(yīng)用

楊釗，許建聰，余俊

（1.同濟(jì)大學(xué)巖土及地下工程教育部重點(diǎn)實(shí)驗(yàn)室，上海 200092；2.中南大學(xué)土木建筑學(xué)院，湖南 長(zhǎng)沙 410075）

彈性地基Timoshenko梁?jiǎn)卧贏BAQUS軟件中的應(yīng)用

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

（1.同濟(jì)大學(xué)巖土及地下工程教育部重點(diǎn)實(shí)驗(yàn)室，上海 200092；2.中南大學(xué)土木建筑學(xué)院，湖南 長(zhǎng)沙 410075）

利用ABAQUS軟件中的自定義單元接口，采用Fortran語(yǔ)言開(kāi)發(fā)彈性地基Timoshenko梁?jiǎn)卧绦?通過(guò)與經(jīng)典解的比較，驗(yàn)證所編彈性地基Timoshenko梁?jiǎn)卧绦虻臏?zhǔn)確性.該單元不僅考慮了地基彈簧受拉脫開(kāi)的特點(diǎn)，而且還考慮了曲形梁?jiǎn)卧獌?nèi)結(jié)點(diǎn)不在一條直線上的特點(diǎn).采用所編寫(xiě)彈性地基梁?jiǎn)卧治銮嗖萆吃^(guò)江輸水工程，計(jì)算規(guī)律與已建類似工程實(shí)測(cè)結(jié)果相同.

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

在地下結(jié)構(gòu)的計(jì)算領(lǐng)域，已有許多用于考慮結(jié)構(gòu)-土體相互作用的計(jì)算方法，但是荷載結(jié)構(gòu)法仍是目前使用最廣泛的一種方法.為了在研究周圍地層對(duì)結(jié)構(gòu)的約束作用及地層對(duì)地下結(jié)構(gòu)的反作用力時(shí)考慮地下結(jié)構(gòu)變形，將Winker彈性地基理論引入到荷載結(jié)構(gòu)法中［1］.同時(shí)，提出彈性地基上的Timoshenko梁，以考慮橫向剪切變形對(duì)厚梁的影響.目前，大量商品化軟件如MARC，SAP等均包含彈性地基梁?jiǎn)卧?，但均沒(méi)有考慮地基彈簧在受拉時(shí)的脫開(kāi)情況.ANSYS，ABAQUS等大型通用有限元法計(jì)算分析軟件中未包括彈性地基梁?jiǎn)卧?，只能通過(guò)在梁?jiǎn)卧Y(jié)點(diǎn)上加彈簧單元來(lái)近似模擬彈性地基的作用［2-3］.這種做法只能在梁?jiǎn)卧銐蛐〉那闆r下，才能近似等價(jià)于彈性地基作用［4］.本文以ABAQUS軟件為平臺(tái)，開(kāi)發(fā)基于Winker地基理論與Timoshenko梁理論的地基梁?jiǎn)卧?

1 Timoshenko梁?jiǎn)卧獎(jiǎng)偠汝嚺c荷載陣

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

Timoshenko梁?jiǎn)卧幕咎攸c(diǎn)是撓度ω、軸向壓縮量u和截面轉(zhuǎn)動(dòng)θ各自獨(dú)立插值.即

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

2 彈性地基梁

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

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

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

式（6）中：Bd=［Bd，1，…，Bd，n］；Bd，i=［Ni，0，0］.將式（6）代入式（5），取極值即可得單元?jiǎng)偠汝囍械鼗鶆偠雀郊禹?xiàng)Kde為

梁?jiǎn)卧獌?nèi)力計(jì)算方程為

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

3 整體坐標(biāo)系下的剛度陣

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

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

式（12）中：Bd，gol=［Bd，1，gol，…，Bd，n，gol］，Bd，i，gol=Bd，λii.整體坐標(biāo)系下，彈性地基梁?jiǎn)卧獎(jiǎng)偠染仃嚍?/p>

當(dāng)由前處理所得梁?jiǎn)卧膬?nèi)部結(jié)點(diǎn)位于梁?jiǎn)卧獌啥私Y(jié)點(diǎn)的連線上，且內(nèi)部結(jié)點(diǎn)為梁?jiǎn)卧木贮c(diǎn)時(shí)，由局部坐標(biāo)與整體坐標(biāo)的轉(zhuǎn)換關(guān)系可得

在實(shí)際工程中，特別是以曲梁形式存在的盾構(gòu)隧道管片，由前處理軟件所得梁?jiǎn)卧獌?nèi)部結(jié)點(diǎn)位于梁兩端結(jié)點(diǎn)連線之外.針對(duì)這種情況，式（14）應(yīng)為

式（15）中：k為單元在整體坐標(biāo)系下的斜率.當(dāng)k=∞時(shí)，有

結(jié)合式（13）～（16），可以得到在整體坐標(biāo)系下單元?jiǎng)偠染仃嚨谋磉_(dá)式.

在數(shù)值積分中，若采用精確積分計(jì)算剪切變形能項(xiàng)，當(dāng)梁很薄的情況下，約束條件γ=（dω/dx）-θ不可能在梁?jiǎn)卧咸幪帩M足.采用精確積分計(jì)算剪切變形能項(xiàng)，將過(guò)分夸大剪切應(yīng)變能的量級(jí)而造成剪切鎖死.為了避免剪切鎖死，對(duì)剪切變形項(xiàng)采用縮減積分計(jì)算［7］.

4 程序編制

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

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

4.1 程序流程

ABAQUS主程序進(jìn)行迭代求解有如下8個(gè)步驟.

（1）計(jì)算坐標(biāo)轉(zhuǎn)換矩陣λi.

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

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

（4）由式（11）計(jì)算整體坐標(biāo)系下梁?jiǎn)卧膯卧獎(jiǎng)偠染仃?

（5）由式（9）計(jì)算積分點(diǎn)處的軸力、剪力、彎矩.

（6）將主程序傳入的位移值代入式（2），（10），計(jì)算地基梁?jiǎn)卧诜e分點(diǎn)處的撓度值.如計(jì)算撓度值為負(fù)值，則取此積分點(diǎn)處的地基剛度為零，再由式（12）計(jì)算地基對(duì)梁?jiǎn)卧母郊觿偠染仃?

（7）由式（13）計(jì)算地基梁?jiǎn)卧恼w剛度矩陣.

（8）計(jì)算單元的殘余力，并判斷收斂.如果不收斂，返回ABAQUS主程序進(jìn)行第i+1次迭代.主程序?qū)⒏鶕?jù)UEL子程序第i次迭代所得到的單元?jiǎng)偠染仃嚺c殘余力項(xiàng)，計(jì)算位移增量與總位移量，然后跳到第（1）步進(jìn)行UEL的第i+1次迭代計(jì)算.

4.2 程序驗(yàn)證

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

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

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

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

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

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

5 工程實(shí)例計(jì)算

青草沙原水過(guò)江隧洞工程位于上海長(zhǎng)江隧道下游約80m處，浦東側(cè)越江點(diǎn)在五號(hào)溝，長(zhǎng)興島越江點(diǎn)在該島新開(kāi)河附近，全長(zhǎng)7.23km.越江輸水管道采用全斷面隧道掘進(jìn)機(jī)（TBM）施工，有壓輸水，設(shè)計(jì)為圓形斷面，襯砌結(jié)構(gòu)外直徑為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.

計(jì)算采用三結(jié)點(diǎn)Timeshenko地基梁.為簡(jiǎn)化計(jì)算，在建立地基梁模型的同時(shí)也建立三結(jié)點(diǎn)Timeshenko梁模型.此梁模型與地基梁模型共結(jié)點(diǎn)且結(jié)點(diǎn)編號(hào)一致，但單元編號(hào)不同.將地基梁上的荷載施加到共結(jié)點(diǎn)的梁上，取共結(jié)點(diǎn)梁?jiǎn)卧獜椥阅Ａ繛橐粯O小數(shù).此時(shí)，梁?jiǎn)卧拇嬖趯⒑?jiǎn)化荷載的施加且對(duì)計(jì)算結(jié)果無(wú)影響.由于梁?jiǎn)卧Y(jié)點(diǎn)與地基梁?jiǎn)卧步Y(jié)點(diǎn)，因此，ABAQUS后處理中梁?jiǎn)卧奈灰茍?chǎng)即為地基梁?jiǎn)卧奈灰茍?chǎng).

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

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

襯砌截面的彎矩、軸力和剪力圖，如圖4～6所示.從圖4～6可知，管片彎矩的峰值出現(xiàn)在管頂、管底和兩腰，其管頂、管底為正值，兩腰彎矩為負(fù)值.最大正彎矩位于管頂處，其值為102.7kN·m；最大負(fù)彎矩位于管腰處，其值為-95.3kN·m.設(shè)軸力以受拉為正，受壓為負(fù).管片大部分截面受壓，軸壓值管頂、管底小，而管腰大.軸壓最大值位于管腰處，其最大軸壓值為150.8kN；管頂部分截面受拉，最大軸拉值為-34.6kN.設(shè)剪力以截面呈順時(shí)針轉(zhuǎn)動(dòng)為正，反之為負(fù)，剪力最大值為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 結(jié)束語(yǔ)

針對(duì)絕大多數(shù)商品化軟件不能考慮彈性地基梁在受拉時(shí)地基彈簧脫開(kāi)的不足，基于ABAQUS軟件平臺(tái)開(kāi)發(fā)了彈性地基Timeshenko梁?jiǎn)卧?該單元不僅考慮了地基彈簧受拉脫開(kāi)的特點(diǎn)，而且還考慮了曲形梁?jiǎn)卧獌?nèi)結(jié)點(diǎn)不在一條直線上的特點(diǎn).采用所編寫(xiě)彈性地基梁?jiǎn)卧治銮嗖萆吃^(guò)江輸水工程，計(jì)算規(guī)律與已建類似工程實(shí)測(cè)結(jié)果相同.研究結(jié)果表明，所編寫(xiě)的彈性地基梁?jiǎn)卧雀?，可供?shí)際工程計(jì)算應(yīng)用.

［1］孫鈞.地下結(jié)構(gòu)［M］.北京：科技出版社，1987.

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

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

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

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

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

［7］馮紫良.桿系結(jié)構(gòu)的計(jì)算機(jī)分析［M］.上海：同濟(jì)大學(xué)出版社，1991.

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

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

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

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

2009-09-19

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

國(guó)家自然科學(xué)基金資助項(xiàng)目（40872179）；中國(guó)博士后科研基金資助項(xiàng)目（20080440652）