林軍志,楊 笛,徐績(jī)青,3

(1. 重慶交通大學(xué) 河海學(xué)院,重慶 400074; 2. 山區(qū)公路水運(yùn)交通地質(zhì)減災(zāi)重慶市教委重點(diǎn)實(shí)驗(yàn)室,重慶 400074; 3. 重慶交通大學(xué) 國(guó)家內(nèi)河航道整治工程技術(shù)研究中心 水利水運(yùn)工程教育部重點(diǎn)實(shí)驗(yàn)室,重慶 400074)

0 引 言

在工程中因地基變形而導(dǎo)致建筑物的損壞是不可無視的,因此與土體接觸的結(jié)構(gòu)應(yīng)按彈性地基梁和彈性地基板來考慮。筆者以Winkler模型為基礎(chǔ),假設(shè)地基表面任一點(diǎn)的沉降與該點(diǎn)單位面積上的壓力成正比,建立了撓度微分方程式和邊界條件。

目前針對(duì)彈性地基梁板的撓度求解研究成果頗豐,如經(jīng)典的Navier彈性地基板雙三角級(jí)數(shù)解[1],該方法是將荷載轉(zhuǎn)化為級(jí)數(shù)形式,將滿足邊界條件的試函數(shù)帶回微分方程得到撓度理論解;初參數(shù)法[2]是用撓度ω、轉(zhuǎn)角θ、彎矩M、剪力Q代替撓度函數(shù)中的4個(gè)參數(shù),使積分常數(shù)具有明確的物理意義;處理變截面梁和邊荷載問題的鏈桿法[2],是用鏈桿力代替地基反力;楊維加[3]對(duì)彈性地基梁采用三角級(jí)數(shù)法進(jìn)行求解;楊成永等[4]針對(duì)帶有脫空彈性地基梁?jiǎn)栴}采用傅里葉級(jí)數(shù)法進(jìn)行求解。以上解析方法都是針對(duì)某些特定條件或規(guī)則形狀進(jìn)行的求解,其計(jì)算形式較為復(fù)雜。

有限元方法在樣條函數(shù)空間尋找近似解,但對(duì)于求解高階偏微分問題所構(gòu)造高階連續(xù)樣條函數(shù)基(高階連續(xù)單元)非常困難,且每次計(jì)算都需要剖分網(wǎng)格,計(jì)算工作量大。無網(wǎng)格方法[5]主要有滑動(dòng)最小二乘法[6]、伽遼金法[7]、雜交邊界法[8]等,這些方法均是通過數(shù)值積分得到弱形式解。基于徑向基函數(shù)[9-14]的求解具有如下優(yōu)點(diǎn):表達(dá)與計(jì)算非常簡(jiǎn)單、各向同性、任意多元函數(shù)都可用一元函數(shù)來描述、節(jié)約存儲(chǔ)成本、可以逼近幾乎所有的從各向同性問題中產(chǎn)生的函數(shù)。采用徑向基函數(shù),結(jié)合加權(quán)余量法[15]中的配點(diǎn)法對(duì)線性方程組進(jìn)行離散,通過笛拉克函數(shù)控制殘差可得到計(jì)算簡(jiǎn)單、精度高、通用性強(qiáng)的強(qiáng)形式解。

1 撓度微分方程與邊界條件

1.1 彈性地基梁

域內(nèi)微分方程的計(jì)算如式(1):

(1)

自然邊界條件如式(2):

W|x=0,l=0

(2)

本質(zhì)邊界條件如式(3):

(3)

式中:EI為彈性地基梁的彎曲剛度,kN·m;k為彈性地基系數(shù),kN/m3;b取單寬,m;l為寬度,m;q(x)為施加在彈性地基梁上的荷載,kN/m。

1.2 彈性地基板

域內(nèi)微分方程的計(jì)算如式(4):

D?4W(x,y)+kW(x,y)=q(x,y)

(4)

自然邊界條件如式(5):

W|x=0,2a=0,W|y=0,2b=0

(5)

本質(zhì)邊界條件如式(6):

(6)

式中:D為彈性地基板的彎曲剛度,kN·m;2a、2b分別為板的長(zhǎng)和寬,m。

2 徑向基函數(shù)

筆者根據(jù)文獻(xiàn)[9-10],基于修正再生核逼近思想構(gòu)造了局部緊支撐徑向基函數(shù)(該函數(shù)使離散插值矩陣具有帶疏性)及根據(jù)經(jīng)驗(yàn)改良后的Gauss公式,如式(7)～式(12)。這樣可以解決因數(shù)值計(jì)算中布點(diǎn)過密產(chǎn)生的病態(tài)矩陣問題。

φ(r)=(1-r)5×(1+5r+9r2+5r3+r4)

(7)

φ(r)=(1-r)6(6+36r+82r2+72r3+30r4+5r5)

(8)

φ(r)=(1-r)6(3+18r+35r2)

(9)

φ(r)=(1-r)7(5+35r+101r2+147r3+101r4+35r5+5r6)

(10)

φ(r)=(1-r)8(1+8r+25r2+32r3)

(11)

φ(r)=expcr2

(12)

需特別注意,若有r=d/dmax則以上函數(shù)應(yīng)滿足式(13)。

(13)

3 徑向基函數(shù)配點(diǎn)法

3.1 加權(quán)余量配點(diǎn)法

已知基本微分方程式和邊界條件[15]:

(14)

式中:D、G分別為微分算子;d、g分別為不含函數(shù)的項(xiàng);u為待求函數(shù)項(xiàng);m為邊界條件數(shù)。

假設(shè)待定函數(shù)u的近似解為Ui,則有式(15):

(15)

式中:ai為待定系數(shù);fi為形式確定的試函數(shù)項(xiàng);n為試函數(shù)項(xiàng)的數(shù)目。

將近似解代入微分方程,則有內(nèi)部殘量R1和邊界殘量R2[15]:

(16)

選用笛拉克函數(shù)作為權(quán)函數(shù)ωi(x),如式(17):

ωi(x)=δ(x-xi)=0,(x≠xi)

(17)

根據(jù)笛拉克函數(shù)的挑選性,控制殘差在一系列配點(diǎn)上xi=0,則有式(18)。

(18)

求解式(18),即可得到待定系數(shù)ai。

3.2 徑向基函數(shù)配點(diǎn)求解

取徑向基函數(shù)作為試函數(shù),則撓度的近似解Ui[11]可由式(19)表達(dá)。

(19)

所求區(qū)域用n個(gè)節(jié)點(diǎn)離散。設(shè)P為區(qū)域內(nèi)部nP個(gè)節(jié)點(diǎn)的集合;Q為邊界x方向上的nQ個(gè)節(jié)點(diǎn)集合;S為邊界y方向上的nS個(gè)節(jié)點(diǎn)集合;則有n=nP+nQ+nS。n個(gè)節(jié)點(diǎn)對(duì)應(yīng)的未知撓度為n個(gè),試函數(shù)數(shù)量與未知量數(shù)量相同,可建立n個(gè)線性方程組,如式(20)。

(20)

令特征矩陣A=[Φ(r)],系數(shù)矩陣α=[α1,α2,…,αn]T,待求矩陣U=[U1,U2,…,Un]T,則可簡(jiǎn)化如式(21):

Aα=U

(21)

任意點(diǎn)的撓度近似解可表示為式(22):

[φ1(ri)φ2(ri) …φn(ri)]A-1U

(22)

梁、板方程的求解〔式(1)～式(6)〕一般分為兩種方法。

方法1:不求逆矩陣,以待定系數(shù)矩陣α為未知量求解。

方法2:求逆矩陣,以待求矩陣U作為未知量求解。

3.3 方法改進(jìn)

在數(shù)值計(jì)算中,常規(guī)配點(diǎn)法往往會(huì)在邊界處產(chǎn)生嚴(yán)重振蕩[14],針對(duì)這一問題可利用本質(zhì)邊界條件新增未知量,借鑒彈塑性靜力學(xué)的處理方法提出n個(gè)節(jié)點(diǎn)撓度和nQ+nS個(gè)邊界點(diǎn)撓度的二階偏導(dǎo)量聯(lián)合插值的近似解,如式(23),明確附加了未知量的物理意義。

(23)

式(23)對(duì)于徑向基函數(shù)的高階連續(xù)性有一定要求,為了避免徑向基函數(shù)高階求導(dǎo)后形成的特征矩陣條件數(shù)增大,應(yīng)對(duì)式(23)進(jìn)行簡(jiǎn)化。

分別采用線性無關(guān)的輔助徑向基函數(shù)βi(ri),γi(ri)來代替高階偏導(dǎo)項(xiàng),如式(24)。

(24)

式(24)中:右邊第1項(xiàng)代表n個(gè)節(jié)點(diǎn)的撓度;第2項(xiàng)代表在邊界點(diǎn)上對(duì)x方向的曲率;第3項(xiàng)代表在邊界點(diǎn)上對(duì)y方向的曲率。

此時(shí)式(24)可寫作式(25):

(25)

n+nQ+nS個(gè)待定系數(shù)的n個(gè)方程組為超定解,因此為滿足一一對(duì)應(yīng)關(guān)系,構(gòu)造了附加未知量、矩陣A的附加行。數(shù)值計(jì)算中為降低方程組求解難度,從減小誤差思想出發(fā),構(gòu)建附加未知量與本質(zhì)邊界條件等價(jià)可以使方程組的附加行列成為類單位陣。改造后的矩陣A及待求矩陣U如式(26)、 式(27):

n個(gè)節(jié)點(diǎn)對(duì)應(yīng)n+nQ+nS個(gè)未知數(shù),建立了n+nQ+nS個(gè)線性方程組,其中Vi,Zi為附加未知量。Vi表示為邊界點(diǎn)在x方向上撓度二階偏導(dǎo)量,即x方向的曲率;Zi表示為邊界點(diǎn)在y方向上撓度二階偏導(dǎo)量,即y方向的曲率。得到的近似解如式(28):

(26)

U=[U1…UnV1…VnQZ1…ZnS]T

(27)

[φ1(ri) …φn(ri)β1(ri) …βnQ(ri)γ1(ri) …γnS(ri)]A-1U

(28)

4 算法實(shí)施

1)建立4階偏微分方程和邊界條件方程;

2)在規(guī)定區(qū)域合理均勻配點(diǎn),并求得個(gè)節(jié)點(diǎn)之間的距離及其影響半徑;

3)選取各線性無關(guān)的徑向基函數(shù)作為主函數(shù)和輔助函數(shù)構(gòu)造近似解Ui;

4)方法1:根據(jù)各節(jié)點(diǎn)位置結(jié)合邊界條件構(gòu)造特征矩陣A使得附加未知量具有物理意義;

方法2:根據(jù)各節(jié)點(diǎn)位置結(jié)合邊界條件構(gòu)造特征矩陣A,使得附加未知量具有物理意義并求得逆矩陣B=A-1;

5)對(duì)近似解Ui求4階偏導(dǎo);

6)方法1:將近似解Ui和4階偏導(dǎo)帶入微分方程中,并針對(duì)邊界點(diǎn)改變?cè)撔械奈⒎址匠淌蛊錆M足自然邊界條件。顯然具有公因子[a1, …,anb1, …,bnQc1,…,cnS]T;

方法2:將近似解和4階偏導(dǎo)函數(shù)帶入微分方程中,并通過化零置一的方法使其滿足自然邊界條件。顯然具有公因子[U1, …,UnV1, …,VnQZ1,…,ZnS]T;

7)方法1:利用A附加行的意義構(gòu)造方程組,滿足微分方程和本質(zhì)邊界條件;

方法2:由于附加未知量與本質(zhì)邊界條件等價(jià),因此在聯(lián)立方程時(shí),式(3)、式(5)無需再次計(jì)算2階偏導(dǎo)式,而是構(gòu)建類單位陣滿足邊界條件。若以第一個(gè)邊界點(diǎn)x方向曲率為例,如式(29):

(29)

轉(zhuǎn)換為式(30):

(30)

8)方法1:聯(lián)立方程組求得未知系數(shù)并帶回式(24),求得近似解Ui;

方法2:聯(lián)立方程組直接求得未知量Ui。

5 算例分析

5.1 彈性地基梁的撓度分析

算例1:假設(shè)一個(gè)兩端簡(jiǎn)支的彈性地基梁,受均布荷載作用q=2 kN/m;簡(jiǎn)支梁長(zhǎng)度為L(zhǎng)=4 m;EI=2.5×106kN·m2;地基彈性系數(shù)k=4×104kN/m3。如圖1。

圖1 彈性地基梁受均布荷載Fig. 1 Elastic foundation beam subjected to uniformly distributed load

已知彈性地基梁的撓度微分方程與邊界條件為式(1)～式(3),可得到理論解,如式(31):

(31)

(32)

(33)

則無量綱撓度理論解如式(34):

(34)

利用徑向基函數(shù)配點(diǎn)逼近的方法進(jìn)行求解,選用文獻(xiàn)[10]提出的一維6階連續(xù)函數(shù)作為主函數(shù)和一維4階連續(xù)函數(shù)作為輔助函數(shù)。一維梁結(jié)構(gòu)不存在y方向的曲率,因此構(gòu)造出近似解Ui和矩陣A。

取節(jié)點(diǎn)間距為0.05均勻布點(diǎn),則一共有n=21個(gè)節(jié)點(diǎn),其中兩個(gè)邊界點(diǎn)nQ=2。采用方法一直接聯(lián)立式(1)～式(3),求解得到未知系數(shù)代入近似解的撓度值如圖2。從數(shù)學(xué)精度方面考慮全局相對(duì)誤差1.38%;工程上為滿足安全性原則,一般考慮最大撓度變形點(diǎn)。由圖2可見:最大撓度變形發(fā)生在跨中0.012 5 m處,與理論解相比,其相對(duì)誤差僅為0.084%,優(yōu)于最小滑動(dòng)二乘解。

圖2 彈性地基梁受均布荷載撓度分析Fig. 2 Deflection analysis of elastic foundation beam subjected to the uniformly distributed load

5.2 彈性地基板的撓度分析

5.2.1 算例2

若有一塊簡(jiǎn)支彈性地基板,板長(zhǎng)為10 m(即a=5),板寬為10 m(即b=5);板的彎曲剛度為D=2.5×104kN·m;地基彈性系數(shù)為k=104kN/m3;板上受到均布荷載為q=2 kN/m2,如圖3。

圖3 彈性地基板受均布荷載Fig. 3 Elastic foundation beam and plate subjected to the uniformly distributed load

已知彈性地基板的撓度微分方程與邊界條件為式(4)～式(6),則可得到理論解,如式(35):

(35)

利用徑向基函數(shù)配點(diǎn)逼近的方法進(jìn)行求解,選用文獻(xiàn)[10]提出的三維4階連續(xù)函數(shù)作為主函數(shù),文獻(xiàn)[9]構(gòu)造的三維4階連續(xù)函數(shù)作為輔助函數(shù)和文獻(xiàn)[10]的三維6階連續(xù)函數(shù)作為輔助函數(shù)構(gòu)造出近似解Ui和矩陣A。

板面上4個(gè)角點(diǎn)同時(shí)考慮對(duì)x、y兩個(gè)方向求曲率,利用方法2求出逆矩陣B,利用附加未知量的意義構(gòu)造一個(gè)類似單位矩陣的大小為2m行n+2m列的矩陣M(M=[RT]),R為2m行n列的零矩陣,T為2m行2m列的單位矩陣。使其乘以n+2m個(gè)未知量,并通過化零置一轉(zhuǎn)化成未知系數(shù)與本質(zhì)邊界條件等價(jià)的未知量。取節(jié)點(diǎn)間距為0.25,均勻布點(diǎn),則一共有n=1 681個(gè)節(jié)點(diǎn),其中邊界點(diǎn)160個(gè)。

得到的各撓度值如圖4。

圖4 彈性地基板受到均布荷載的撓度分析Fig. 4 Deflection analysis of elastic foundation beam and plate subjected to the uniformly distributed load

板最大撓度變形發(fā)生在板中央,近似解為2.42×10-4m,略大于理論解,撓度最大值相對(duì)誤差為0.012%。此時(shí)逆矩陣的求逆精度為3.05×10-5,從工程實(shí)際而言效果較好;節(jié)點(diǎn)全局平均相對(duì)誤差為0.46%,從數(shù)學(xué)層面考慮精度較好。

5.2.2 算例3

若有一塊四邊簡(jiǎn)支的正方形薄板,幾何尺寸為1 m×1 m(即a=0.5,b=0.5)。板的彈性模量為E=2.1×1010Pa,泊松比為μ=0.3,地基彈性系數(shù)為k=4.9×104kN/m3,板上受到均布荷載為q=1 kN/m2。

利用徑向基函數(shù)配點(diǎn)逼近方法進(jìn)行求解,選用改良后的Gauss函數(shù)作為主函數(shù),輔助函數(shù)為式(8)、 式(9)分別對(duì)x、y方向的曲率。在滿足域內(nèi)微分方程和邊界條件的情況下構(gòu)造出近似解Ui,采用方法二求逆矩陣B,并以待求矩陣U作為未知量建立方程組。

取節(jié)點(diǎn)間距為0.025,共有n=1 681個(gè)節(jié)點(diǎn),其中160個(gè)邊界點(diǎn)。由于地基反力變大,撓度最大值不會(huì)發(fā)生在板中央,由此得到撓度結(jié)果[7],見圖5。與理論解對(duì)比,撓度最大值相對(duì)誤差為0.178%,與整體相對(duì)誤差嚴(yán)重不匹配,這是因?yàn)镚auss函數(shù)求逆精度低。

圖5 馬鞍形撓度值(調(diào)整前)Fig. 5 Saddle-type deflection value (before the adjustment)

故從代數(shù)方面考慮,讓邊界點(diǎn)除了滿足邊界條件外同時(shí)滿足邊界微分方程,如式(36),通過改變矩陣大小及排列來提高精度。

(36)

主函數(shù)與輔助函數(shù)不變,此時(shí)x、y兩個(gè)方向的曲率都用式(8)表示,4階偏導(dǎo)用式(9)表示,構(gòu)造出近似解Ui和矩陣A如式(37)、式(38),其余步驟與方法2相同。

(37)

(38)

調(diào)整后的撓度值如圖6。與理論解對(duì)比,其整體相對(duì)誤差明顯減小,撓度最大值的相對(duì)誤差為0.013%,比調(diào)整前提高10倍,精度滿足工程需要;板中心撓度相對(duì)誤差為0.013 5%,與雜交邊界法相當(dāng)。

圖6 馬鞍形撓度值(調(diào)整后)Fig. 6 Saddle-type deflection value (after the adjustment)

6 結(jié) 論

筆者為求解彈性地基梁和彈性地基板的撓度值,采用了徑向基函數(shù)配點(diǎn)逼近的方法,得出如下結(jié)論:

1)徑向基函數(shù)配點(diǎn)法形式簡(jiǎn)單,無網(wǎng)格化,計(jì)算方便,求解精度高;

2)在構(gòu)建矩陣A時(shí),可根據(jù)邊界條件編寫附加行,使得附加未知量具有一定的物理意義,從而通過簡(jiǎn)化方程求解來減小誤差;

3)針對(duì)方程求解提出兩種算法:不求逆矩陣,以待定系數(shù)作為未知數(shù)進(jìn)行求解;求逆矩陣,以撓度值作為未知數(shù)進(jìn)行求解,根據(jù)矩陣求逆精度選擇算法;

4)參考緊支柱思想,改造了與全局支撐域相關(guān)聯(lián)的Gauss函數(shù),使其無量綱化,根據(jù)經(jīng)驗(yàn)確定其參數(shù)c=-23.026;

5)為解決改造后的Gauss函數(shù)求逆性差的問題,可利用邊界微分方程建立輔助函數(shù),以提高求解精度。