彈性地基梁板的徑向基函數(shù)逼近求解方法

林軍志,楊 笛,徐績青,3

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

0 引 言

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

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

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

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

1.1 彈性地基梁

域內(nèi)微分方程的計算如式(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)微分方程的計算如式(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分別為板的長和寬,m。

2 徑向基函數(shù)

筆者根據(jù)文獻[9-10],基于修正再生核逼近思想構(gòu)造了局部緊支撐徑向基函數(shù)(該函數(shù)使離散插值矩陣具有帶疏性)及根據(jù)經(jīng)驗改良后的Gauss公式,如式(7)～式(12)。這樣可以解決因數(shù)值計算中布點過密產(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ù)應滿足式(13)。

(13)

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

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

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

(14)

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

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

(15)

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

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

(16)

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

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

(17)

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

(18)

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

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

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

(19)

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

(20)

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

Aα=U

(21)

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

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

(22)

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

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

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

3.3 方法改進

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

(23)

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

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

(24)

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

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

(25)

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

n個節(jié)點對應n+nQ+nS個未知數(shù),建立了n+nQ+nS個線性方程組,其中Vi,Zi為附加未知量。Vi表示為邊界點在x方向上撓度二階偏導量,即x方向的曲率;Zi表示為邊界點在y方向上撓度二階偏導量,即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 算法實施

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

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

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

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

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

5)對近似解Ui求4階偏導;

6)方法1:將近似解Ui和4階偏導帶入微分方程中,并針對邊界點改變該行的微分方程使其滿足自然邊界條件。顯然具有公因子[a1, …,anb1, …,bnQc1,…,cnS]T;

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

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

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

(29)

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

(30)

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

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

5 算例分析

5.1 彈性地基梁的撓度分析

算例1:假設(shè)一個兩端簡支的彈性地基梁,受均布荷載作用q=2 kN/m;簡支梁長度為L=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ù)配點逼近的方法進行求解,選用文獻[10]提出的一維6階連續(xù)函數(shù)作為主函數(shù)和一維4階連續(xù)函數(shù)作為輔助函數(shù)。一維梁結(jié)構(gòu)不存在y方向的曲率,因此構(gòu)造出近似解Ui和矩陣A。

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

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

5.2 彈性地基板的撓度分析

5.2.1 算例2

若有一塊簡支彈性地基板,板長為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ù)配點逼近的方法進行求解,選用文獻[10]提出的三維4階連續(xù)函數(shù)作為主函數(shù),文獻[9]構(gòu)造的三維4階連續(xù)函數(shù)作為輔助函數(shù)和文獻[10]的三維6階連續(xù)函數(shù)作為輔助函數(shù)構(gòu)造出近似解Ui和矩陣A。

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

得到的各撓度值如圖4。

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

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

5.2.2 算例3

若有一塊四邊簡支的正方形薄板,幾何尺寸為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ù)配點逼近方法進行求解,選用改良后的Gauss函數(shù)作為主函數(shù),輔助函數(shù)為式(8)、 式(9)分別對x、y方向的曲率。在滿足域內(nèi)微分方程和邊界條件的情況下構(gòu)造出近似解Ui,采用方法二求逆矩陣B,并以待求矩陣U作為未知量建立方程組。

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

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

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

(36)

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

(37)

(38)

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

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

6 結(jié) 論

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

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

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

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

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

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