朱桂春，張雯超，陳勝

(1.蘇州大學 軌道交通學院，江蘇 蘇州 215000；2.揚州工業(yè)職業(yè)技術(shù)學院 建筑工程學院，江蘇 揚州 225000；3.中國地質(zhì)大學(武漢新校區(qū)建設(shè)指揮部)，湖北 武漢 430074)

0 引 言

土體大多自然沉積，形成層狀構(gòu)造[1-2]，但在工程實踐中，一般將土視為均質(zhì)體，采用彈性模量加權(quán)平均方法近似模擬分層土，對其變形特性進行研究，但是這樣處理可能會引起建筑物的不均勻沉降[3]，危害建筑物的安全。

對于任意荷載作用下的層狀地基應力應變狀態(tài)求解，可以簡化為平面應變問題。平面應變問題是三維空間問題的一種特殊情況，即假定y軸方向無位移，研究層狀地基平面應變問題具有一定的現(xiàn)實意義，如在層狀地基中開挖地連墻槽段，槽壁受側(cè)向靜止土壓力與泥漿壓力的共同作用，可簡化為受三角形水平分布荷載作用的平面應變問題求解地層的變形情況。

目前求解層狀結(jié)構(gòu)的方法主要有位移函數(shù)法[4-7]、柔度矩陣法[8-12]、傳遞矩陣法[13-14]、精確的剛度矩陣法[15-17]和解析層元法[18]等。已有層狀結(jié)構(gòu)求解方法中，傳遞矩陣法計算時數(shù)據(jù)溢出，而剛度矩陣法中的矩陣元素表達式復雜，不便于實際應用。

本文利用Mathematica軟件計算程序，基于擴展的瑞利-里茲法分析層狀地基應變問題的近似解答，推導非對稱荷載作用下層狀地基平面應變問題的解答，以及各參數(shù)對變形的敏感性，以期為類似工程的設(shè)計與施工提供參考。

1 各向同性地基近似解析解的推導

1.1 各向同性地基平面問題的基本關(guān)系式

各向同性體的胡克定律為

(1)

式中,λ=E/[(1+v)(1-2v)]。

1.2 推導思路

應變的近似解析解求解思路(瑞利-里茲法)如下：

(1)假定任一層狀區(qū)域內(nèi)任一點的水平位移為u(x,z)，豎直位移為v(x,z)。

(2)構(gòu)造泛函

(3)據(jù)最小勢能原理，對泛函進行變分，求其駐值，同時考慮邊界條件約束，求解待定參數(shù)。

對于地表以下受水平非對稱荷載作用的層狀各向同性土體，由于結(jié)構(gòu)復雜，參數(shù)較多，使用傳統(tǒng)的瑞利-里茲法求解，積分難度較大，難以得出正確的結(jié)果。為求解方便，節(jié)省計算時間，將層狀各向同性土體劃分為多層有限區(qū)域和無限區(qū)域，有限區(qū)域與無限區(qū)域位移采用有理多項式假定，將每層土體作為一個塊體，利用擴展的瑞利-里茲法對其進行近似解析解求解，區(qū)域示意圖見圖1。

圖1 各向同性地基分層示意圖

假定土層為各向同性土層，各層層厚分別為h1，h2，h3，…，hn，HFE為有限區(qū)域，iHFE為無限區(qū)域，模型右側(cè)邊界一定深度內(nèi)受荷載作用，其余深度土層右側(cè)邊界受水平位移約束。左側(cè)及下部邊界為無窮遠處。

經(jīng)驗證，標準瑞利-里茲法中的位移假定函數(shù)

w=C0+C1x+C2x2+C3x3，

(3)

等價于

w=N1(x)wi+N2(x)θi+N3(x)wj+
N4(x)θj，

(4)

兩者待定系數(shù)個數(shù)相同，Ni為插值基或形函數(shù)。

同理可知，二維位移函數(shù)

(5)

等價于

w=N1(ξ,η)w1+N2(ξ,η)w2+N3(ξ,η)w3+
N4(ξ,η)w4，

(6)

式中：N1=(1-ξ)(1-η)/4；N2=(1+ξ)(1-η)/4；N3=(1+ξ)(1+η)/4；N4=(1-ξ)(1+η)/4。

1.3 有限部分推導

為提高假定的多項式精度，可增加待定系數(shù)的個數(shù)，如

w=C0+C1ξ+C2η+C3ξη+C4ξ2+C5η2+C6ξ3+
C7ξ2η+C8ξη2+C9η3……+Cnξmηm， (7)

可等價假定為

(8)

在[-1,1]×[-1,1]區(qū)域內(nèi),ξ軸上投影n個坐標，η軸上投影n個坐標，從下向上第i行，從左向右第j列的節(jié)點號對應的形式數(shù)為Nij=Ni(η)×Nj(ξ)，節(jié)點標號順序如圖2所示。

圖2 節(jié)點投影示意圖

Ni為ξ軸[-1,1]內(nèi)，n個點構(gòu)成一維拉格朗日插值基，

式中：-1=ξ1<ξ2<ξ3…<ξn=1；ξi的分布符合擴展切比雪夫多項式零點分布。同理Nj將ξ替換成η即可,式(8)假定完畢。它的多項式成分可以從帕斯卡三角形中去數(shù)n項。

由于以上假定都是基于[-1,1]×[-1,1]范圍的單位坐標系，而地層結(jié)構(gòu)建立在x-z坐標系內(nèi)，因此可引入坐標變換，建立(x-z)與(ξ-η)之間的一一對應關(guān)系?？烧J為x與z都是(ξ，η)的函數(shù)，即x=x(ξ，η)，z=z(ξ，η)，則該函數(shù)不唯一，只要連續(xù)、一一對應即可，最簡單的關(guān)系為多項式關(guān)系，有限區(qū)域的坐標變換見圖3。

圖3 有限區(qū)域坐標變換

則有

基于以上分析，一切均可用(ξ-η)坐標描述，方便后續(xù)計算。

(11)

若(ξ-η)坐標可對應得(x-z)坐標及u,v；根據(jù)上述位移模式，可利用最小勢能原理：

(12)

根據(jù)幾何方程：

(13)

對其插值得

(14)

由于Nu，Nw為關(guān)于ξ，η的函數(shù)，不可直接求關(guān)于x，z的偏導，借助鏈式法則有

(15)

即

(16)

記

(18)

于是，式(14)中的B矩陣便可解出，從而根據(jù)[σ]=[D]·[ε]=[D][B]·[δe]，與1.1中，D含義是否相同

(19)

將應變[ε]根據(jù)最小勢能原理代回,可得

(20)

節(jié)點位移矩陣

由于δδe的任意性，式(12)等價為

經(jīng)多次調(diào)試，盡量平衡計算時間與計算精度，文中式(7)采用8次多項式假定。

1.4 無限部分推導

無限部分形函數(shù)與有限部分中一致，節(jié)點編號也一致，但其坐標變換有所不同，如圖4所示。

圖4 無限區(qū)域坐標變換

滿足兩坐標系下的同號節(jié)點一一對應，則有

(23)

可見，①②③④圍成的區(qū)域為內(nèi)映射，其他區(qū)域為外映射，可檢驗其正確性。坐標變換時，僅需①②③④點的坐標，不需要指定無窮遠處的坐標。為建模方便，可以形式上指定無窮遠節(jié)點坐標，方便查看，實際計算不會影響。

由于無限區(qū)域與有限區(qū)域采用相同的節(jié)點插值和自由度，所以在單元組裝時，將它們視為相同的單元,即K1=K有限，K2=K無限，這樣就可以考慮兩種單元共節(jié)點的情況，自動滿足兩者交界處變形協(xié)調(diào)。

2 模型驗證

為驗證近似解析解計算結(jié)果的準確性，本文采用MIDAS/GTS有限元模型進行驗證。

2.1 單層各向同性地基的驗證

模型設(shè)置如圖5所示。

圖5 單層地基計算模型圖

由于近似解析解計算程序可實現(xiàn)對無窮遠處邊界的高精度描述。為了保證驗證的準確性，將MIADAS有限元地基模型范圍設(shè)置為136 m×136 m，左側(cè)與底部邊界受法向位移約束，右側(cè)邊界上部受x向均布荷載p作用，p=10 kPa，作用深度D=8 m，右側(cè)邊界下部受法向位移約束；而近似解析解計算模型右側(cè)邊界約束條件同有限元模型，左側(cè)及底部邊界為無限邊界，為方便與有限元計算結(jié)果進行比對，近似解析解計算程序模型范圍為136 m×12 m，有限及無限區(qū)域的劃分如圖1所示,其中有限區(qū)域范圍為100 m×8 m。橫觀各向同性地基彈性參數(shù)分別為E=20 MPa，v=0.25，G=8 MPa。

近似解析解計算程序與MIDAS/GTS有限元軟件計算結(jié)果如圖6～7所示。其中土體水平位移取右側(cè)邊界節(jié)點，豎向位移取上部邊界節(jié)點，下同。由圖6～7可以看出，兩者結(jié)果吻合度較高，說明本文所用理論的正確性和程序計算的準確性，可滿足精度要求。

圖6 驗證例水平位移對比

圖7 驗證例豎向位移對比

2.2 層狀各向同性地基的驗證

為驗證近似解析解在層狀地基中的適用性，本文對水平向均布荷載作用下三層各向同性地基進行求解，計算模型如圖8所示。

圖8 多層地基計算模型圖

MIADAS有限元地基模型范圍設(shè)置為136 m×136 m，左側(cè)與底部邊界受法向位移約束，右側(cè)邊界上部受x向均布荷載p作用，p=10 kPa，作用深度D=8 m，右側(cè)邊界下部受法向位移約束，而近似解析結(jié)計算模型右側(cè)邊界約束條件同有限元模型，左側(cè)及底部邊界為無限邊界，近似解析解計算模型如2.1節(jié)。有限元模型各層地基彈性參數(shù)如表1所示。

表1 各地層參數(shù)

將求解結(jié)果與MIDAS/GTS計算的結(jié)果對比，如圖9～10所示。從圖9～10可以看出，各點吻合度較高，說明本計算程序在求解層狀各向同性地基變形依然具有較高的準確性和計算精度。

圖9 層狀地基驗證例水平位移對比

圖10 層狀地基驗證例豎向位移對比

對本文計算程序和MIDAS/GTS軟件響應速度進行對比，測試操作系統(tǒng)版本：Win7旗艦版Service Pack1 64位；CPU版本：AMD A8-7650K Radeon R7 3.3 Ghz；內(nèi)存：8.00 GB。經(jīng)測試，本文計算程序計算時間約17.5 s，MIDAS/GTS軟件計算時間為11.68 s，但本文計算程序在工況發(fā)生變化時，參數(shù)輸入簡便，可自動重新劃分網(wǎng)格，設(shè)置邊界條件，較傳統(tǒng)的有限元軟件建模時間大幅縮短，操作更為便捷。

3 土的分層特性對位移的影響

由以上分析可知，對地層變形敏感度較高的參數(shù)為彈性模量，所以下述算例考慮層狀均質(zhì)地基在水平荷載作用下每層彈性參數(shù)不同對地層變形的影響，主要控制變量為各層的彈性模量E。

算例模型尺寸及邊界條件同2.1節(jié)，土體分為5層，有限區(qū)域各層地基E,v,h如表2所示。各工況地層分布如表3所示。

表2 各地層參數(shù)

表3 各工況地層分布

另外增加工況6，工況6為均質(zhì)地基，采用模量與深度的加權(quán)平均近似模擬分層土的特性，其彈性參數(shù)E由式(24)計算。各工況位移如圖11～12所示。

(24)

由圖11～12可以看出，不同工況下的地基水平位移和豎向位移明顯不同，土的分層特性對地基的位移有顯著影響：上層土的物理性質(zhì)對整體的變形影響較大，上層土的彈性模量越大，其水平及豎向的最大位移越小。對比工況1～5和工況6可以看到，在相同的外部荷載及邊界條件下，分層地基與均質(zhì)地基的位移變化差異較大，因此，工程實踐中將土層視為各項同性的均質(zhì)體，采用模量與深度的加權(quán)平均近似模擬分層土的特性進行設(shè)計計算的方法是不妥當?shù)摹?/p>

圖12 各工況豎向位移

4 結(jié) 論

(1)本文基于擴展的瑞利-里茲法分析層狀地基應變問題的近似析解，求解過程采用高階多項式和無窮坐標變換,實現(xiàn)了對無窮遠處邊界的高精度描述，克服了傳統(tǒng)有限元方法對無窮遠處結(jié)構(gòu)盲目截斷的缺點。

(2)層狀結(jié)構(gòu)對地基變形有較大影響，其中上層土的物理性質(zhì)對整體的變形影響較大，工程實踐中常規(guī)的將土層視為各項同性的均質(zhì)體，采用模量與深度的加權(quán)平均近似模擬分層土的特性進行設(shè)計計算的方法是不妥當?shù)?，所以在工程實際應用中考慮土的分層特性是十分必要的。