基于瑞利-里茲法的層狀地基應(yīng)變問題的近似解析解

朱桂春，張雯超，陳勝

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

0 引 言

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

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

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

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

1 各向同性地基近似解析解的推導(dǎo)

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

各向同性體的胡克定律為

(1)

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

1.2 推導(dǎo)思路

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

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

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

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

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

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

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

經(jīng)驗(yàn)證，標(biāo)準(zhǔn)瑞利-里茲法中的位移假定函數(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 有限部分推導(dǎo)

為提高假定的多項式精度，可增加待定系數(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個坐標(biāo)，η軸上投影n個坐標(biāo)，從下向上第i行，從左向右第j列的節(jié)點(diǎn)號對應(yīng)的形式數(shù)為Nij=Ni(η)×Nj(ξ)，節(jié)點(diǎn)標(biāo)號順序如圖2所示。

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

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

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

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

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

則有

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

(11)

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

(12)

根據(jù)幾何方程：

(13)

對其插值得

(14)

由于Nu，Nw為關(guān)于ξ，η的函數(shù)，不可直接求關(guān)于x，z的偏導(dǎo)，借助鏈?zhǔn)椒▌t有

(15)

即

(16)

記

(18)

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

(19)

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

(20)

節(jié)點(diǎn)位移矩陣

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

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

1.4 無限部分推導(dǎo)

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

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

滿足兩坐標(biāo)系下的同號節(jié)點(diǎn)一一對應(yīng)，則有

(23)

可見，①②③④圍成的區(qū)域?yàn)閮?nèi)映射，其他區(qū)域?yàn)橥庥成洌蓹z驗(yàn)其正確性。坐標(biāo)變換時，僅需①②③④點(diǎn)的坐標(biāo)，不需要指定無窮遠(yuǎn)處的坐標(biāo)。為建模方便，可以形式上指定無窮遠(yuǎn)節(jié)點(diǎn)坐標(biāo)，方便查看，實(shí)際計算不會影響。

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

2 模型驗(yàn)證

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

2.1 單層各向同性地基的驗(yàn)證

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

圖5 單層地基計算模型圖

由于近似解析解計算程序可實(shí)現(xiàn)對無窮遠(yuǎn)處邊界的高精度描述。為了保證驗(yàn)證的準(zhǔn)確性，將MIADAS有限元地基模型范圍設(shè)置為136 m×136 m，左側(cè)與底部邊界受法向位移約束，右側(cè)邊界上部受x向均布荷載p作用，p=10 kPa，作用深度D=8 m，右側(cè)邊界下部受法向位移約束；而近似解析解計算模型右側(cè)邊界約束條件同有限元模型，左側(cè)及底部邊界為無限邊界，為方便與有限元計算結(jié)果進(jìn)行比對，近似解析解計算程序模型范圍為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é)點(diǎn)，豎向位移取上部邊界節(jié)點(diǎn)，下同。由圖6～7可以看出，兩者結(jié)果吻合度較高，說明本文所用理論的正確性和程序計算的準(zhǔn)確性，可滿足精度要求。

圖6 驗(yàn)證例水平位移對比

圖7 驗(yàn)證例豎向位移對比

2.2 層狀各向同性地基的驗(yàn)證

為驗(yàn)證近似解析解在層狀地基中的適用性，本文對水平向均布荷載作用下三層各向同性地基進(jìn)行求解，計算模型如圖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可以看出，各點(diǎn)吻合度較高，說明本計算程序在求解層狀各向同性地基變形依然具有較高的準(zhǔn)確性和計算精度。

圖9 層狀地基驗(yàn)證例水平位移對比

圖10 層狀地基驗(yàn)證例豎向位移對比

對本文計算程序和MIDAS/GTS軟件響應(yīng)速度進(jìn)行對比，測試操作系統(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ì)地基的位移變化差異較大，因此，工程實(shí)踐中將土層視為各項同性的均質(zhì)體，采用模量與深度的加權(quán)平均近似模擬分層土的特性進(jìn)行設(shè)計計算的方法是不妥當(dāng)?shù)摹?/p>

圖12 各工況豎向位移

4 結(jié) 論

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

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