郭樹起

(石家莊鐵道大學(xué)省部共建交通工程結(jié)構(gòu)力學(xué)行為與系統(tǒng)安全國家重點(diǎn)實(shí)驗(yàn)室,石家莊 050043)

(石家莊鐵道大學(xué)工程力學(xué)系,石家莊 050043)

引言

邊界元方法[1-6]作為一種非常成功的數(shù)值方法,在力、熱、電、磁等多個(gè)物理分支都得到了廣泛的應(yīng)用[7-16].對于難于求得解析解的問題,可以采用邊界元法進(jìn)行數(shù)值求解[17-19].邊界元法通常導(dǎo)致邊界積分方程,也稱邊界積分方程方法.本文則稍稍變化一下邊界元法的求解思路,從Somigliana 等式開始推導(dǎo),利用了線彈性算子零空間的不平凡性,即格林函數(shù)的不唯一性,導(dǎo)出了不同于邊界積分方程方法的另一種方法,使線彈性問題在求解解析解時(shí)提供了一個(gè)新的選擇.首先給出此邊界積分法的基本理論與基本步驟.然后為檢驗(yàn)此法的有效性,求解了一個(gè)典型的彈性問題,圓形夾雜問題.

1 邊界積分法

可以導(dǎo)出邊界元法的邊界積分方程的Somigliana 等式為

其中,ui(P)是在P點(diǎn)的位移,是在任意一點(diǎn)P沿xi方向作用單位集中力時(shí),在任意一點(diǎn)Q處引起的xj方向的位移分量,fj是體力,uj(q)是邊界位移,tj(q)是邊界面力,是相應(yīng)于邊界位移的邊界面力.是在P處沿著任意方向nm的單位集中力作用下,在Q處的引起的位移.的關(guān)系為

將式 (3a)等號兩邊同時(shí)乘以單位方向nm,得即滿足方程

其中Lik是線彈性算子

其中,Cijkl是線彈性常數(shù).線彈性算子Lik與其共軛算子相等.因此式(3a)和式(3b)可以重寫為

對于各向同性彈性體,則進(jìn)一步化為

其中,λ 和G是常數(shù).稱格林函數(shù)[20].

由于線彈性算子Lik的零空間是非平凡的,因此格林函數(shù)不唯一.利用此不唯一性,記兩個(gè)任意的格林函數(shù)為均滿足方程(3a),將其代入到Somigliana 等式(1)中[4,20],將得到的兩個(gè)方程相減得到

或者

各向同性線彈性體則滿足

方程(7)與式(10)所示的Betti 互易定理[21-22]形式上類似.

其中,()(1)和()(2)分別表示同一物體的線彈性邊值問題I,II 的位移uj,面力tj與體力fj.但并不是某個(gè)力學(xué)邊值問題的解,而是需要滿足無體力的齊次方程(8)或(9).Somigliana 等式可以從Betti 互易定理推導(dǎo)而來[4],即把格林函數(shù)定義式(5)代入式(10).式(7)也可以由Betti 互易定理得到.在式(10)中令=0,所得方程即與式(7)含義相同,只是記號不同.將代入應(yīng)力平衡方程中.即得無體力的應(yīng)力平衡方程,在式(8)和式(9)中代入本構(gòu)方程與幾何方程,同樣得無體力的應(yīng)力平衡方程.

無體力時(shí),式(7)變?yōu)槭匠?11)

方程(11)是本文計(jì)算的主要理論依據(jù).為簡單起見,下面本文研究無體力的平面彈性力學(xué)問題.這時(shí)邊界面S進(jìn)一步化為邊界曲線L

利用式(12)求得未知的邊界位移與面力.對于無體力平面問題來說,位移函數(shù)可以采用相應(yīng)于Michell 解的位移函數(shù)[23-25].一般來說,邊界位移與面力假設(shè)中一般包含不止一個(gè)待定常數(shù),因此需要符合題意的多個(gè)位移函數(shù)以及面力函數(shù)這樣的函數(shù)在此稱為試函數(shù).邊界位移與邊界面力全部已知后,邊界元法用Somigliana 等式計(jì)算域內(nèi)位移.本文則繼續(xù)用式(12)求域內(nèi)位移,避免使用格林函數(shù)帶來奇異積分問題.基本步驟是作一條虛擬邊界L?經(jīng)過需要計(jì)算位移與應(yīng)力的點(diǎn),即式(12)變?yōu)?/p>

這類似于矩陣傳遞法.因此方程(13)預(yù)期可以用于多相介質(zhì)問題的分析.表1 簡敘述了邊界元法與此新方法的主要思想與各自特點(diǎn).

表1 邊界元法與本文方法比較Table 1 Comparison between BEM and the proposed method

2 圓形夾雜問題

圓形夾雜問題是彈性力學(xué)的一個(gè)經(jīng)典問題,通常采用Airy 應(yīng)力函數(shù)法與復(fù)變函數(shù)法求解[26-30].理論結(jié)果被實(shí)驗(yàn)所驗(yàn)證[31].在此用邊界積分法求其解析解.如圖1 所示,無限大平板上有半徑r=a的圓形夾雜,夾雜內(nèi)的彈性常數(shù)與基體不同.基體的彈性常數(shù)用E,ν,μ分別來表示楊氏模量、泊松比、剪切模量,夾雜內(nèi)的彈性常數(shù)用下標(biāo)1 以示區(qū)別,如E1,ν1,μ1.夾雜所在區(qū)域?yàn)?f,其他不含夾雜的區(qū)域?yàn)?.

圖1 圓形夾雜示意圖Fig.1 Schematic diagram of a circular inclusion

2.1 界面位移與面力

界面處認(rèn)為是完美連接的,位移、應(yīng)力連續(xù)

其中,上標(biāo)()f表示與夾雜有關(guān)的量.平板在無窮遠(yuǎn)處有均布應(yīng)力σx=q,其他應(yīng)力分量為零.相應(yīng)的極坐標(biāo)下位移如式(15a)所示;徑向應(yīng)力周向應(yīng)力與切應(yīng)力如式(15b)所示

先研究基體部分,即夾雜之外的區(qū)域?.把位移與應(yīng)力分解為均勻場部分與擾動部分[33],ui=ui+并考慮式(15a)和式(15b),有

這時(shí),無窮遠(yuǎn)處單向均勻應(yīng)力邊界條件化為了位移ui與應(yīng)力τij在無窮遠(yuǎn)處為零,即

采用極坐標(biāo)系,式(12)在本節(jié)中改為

考慮邊界L即為r=a,以及ds=rdθ=adθ,并將式(16)代入式(18)得到

其中,a0,a2,b2,c0,c2,d2是待定常數(shù).

表2 位移與應(yīng)力函數(shù)φTable 2 Displacementsand stress function φ

表2 位移與應(yīng)力函數(shù)φTable 2 Displacementsand stress function φ

這樣在夾雜與基體的交界處r=a處的位移u與面力t為

將表達(dá)式(22)代入式(18),得到式(23).然后將階數(shù)n=0 與n=2 階的位移試函數(shù)與面力代入,得到關(guān)于待定系數(shù)方程式(24).

從式(24)解出待定常數(shù)a0,a2,b2

若c0=c2=d2=0,得a0=qa/E,a2=2qa/E,b2=?2qa/E,將其代入式(20),即得夾雜退化為圓孔時(shí),圓孔邊界的位移、面力

式(24)中只有3 個(gè)方程,少于式(20)中待定系數(shù)的個(gè)數(shù).為此需要單獨(dú)考慮圓形夾雜的情況.研究ra的夾雜區(qū)域?f,由于此時(shí)所考慮區(qū)域不含無窮遠(yuǎn)處,因此不再將位移、面力分解為均勻場部分與擾動部分,此時(shí)邊界積分方程(18)更改為

邊界位移與面力仍如式(20)所示.圓形夾雜的彈性常數(shù)與基體不同,所取位移的彈性常數(shù)需要相應(yīng)改變.另外由于ra區(qū)域包含原點(diǎn),所取位移試函數(shù)在原點(diǎn)的位移值以及相應(yīng)的應(yīng)力值都應(yīng)是有限的.所取位移試函數(shù)如表3 所示,其中μ1=E1/(2+2ν1).平面應(yīng)力,κ1=(3-ν1)/(1+ν1);平面應(yīng)變κ1=3-4ν1.所得方程為式(28),解出待定常數(shù)a0,a2與b2,如式(29)所示.

表3 位移與應(yīng)力函數(shù)φTable 3 Displacementsand stress function φ

表3 位移與應(yīng)力函數(shù)φTable 3 Displacementsand stress function φ

根據(jù)位移與應(yīng)力連續(xù)性條件,兩式(25)和式(29)中的a0,a2,b2,c0,c2,d2不僅記號相同、數(shù)值也相等.聯(lián)立式(24)和式(27)得

由式(30)知道d2=?c2,代入式(25)和式(29)中可以進(jìn)一步簡化為式(31)和式(32)

由式(30)知道,當(dāng)E1=0,c0=c2=d2=0,界面面力為零,簡化為圓孔解答.當(dāng)E1=+∞,表示剛性夾雜,式(30)化為c0=q/(1+ν),c2=2q/(3-ν),d2=?2q/(3-ν),將其代入式(25)中得到a0=a2=b2=0,這表示界面位移等于0;將其代入式(20)得剛性夾雜的r=a處的面力

2.2 域內(nèi)位移與應(yīng)力

2.2.1 基體位移與應(yīng)力(ra)

有了位移與應(yīng)力在界面處的值,就可以計(jì)算域內(nèi)位移場與應(yīng)力場.先研究夾雜之外的區(qū)域?,虛設(shè)一條邊界r=ρ >a.研究由r=a與r=ρ 圍成的環(huán)形區(qū)域,所考慮區(qū)域不含無窮遠(yuǎn)處,因此不再將位移、面力分解為均勻場部分與擾動部分,此時(shí)邊界積分方程式(18)變?yōu)?/p>

在r=a處,位移ur,uθ與面力tr,tθ如式(20)所示,其中待定系數(shù)如式(25)和式(30)所示.在r=ρ 處,假設(shè)位移ur,uθ與面力tr,tθ設(shè)為式(33),即類似式(20)的形式

將式(20)和式(35)代入式(30),并選取適當(dāng)?shù)奈灰圃嚭瘮?shù)與面力由于此環(huán)形區(qū)域既不包含原點(diǎn)也沒有無窮遠(yuǎn)處做邊界,因此只滿足位移單值條件即可.所選取的應(yīng)力函數(shù)φ 如表2 與表3 所示,此時(shí)表3 中的μ1,κ1應(yīng)替換為μ,κ.面力由式(21)求出或者直接查表得到[33].并利用式(25)、c2+d2=0 與a2+b2=0,所得方程如下

聯(lián)立此6 個(gè)方程,得待定系數(shù)A0,A2,B2,C0,C2與D2,考慮式(30),并用彈性常數(shù)κ 與μ來替換E與ν,用κ1與μ1來替換E1與ν1,結(jié)果如式(36)所示

其中,α,β 是彈性常數(shù)的組合,如式(37)所示

將式(36)代入式(35)得到位移場ur,uθ與除σθ之外的應(yīng)力場.σθ可以由位移ur,uθ代入幾何方程,再代入本構(gòu)方程得出.并將ρ 仍寫為r,結(jié)果如式(38)所示

2.2.2 夾雜中位移與應(yīng)力(r < a)

研究夾雜區(qū)域?f,虛設(shè)一條邊界r=ρ

式(39)形式上與式(34)完全相同,除了r=ρ

將式(20)和式(40)代入式(39),并選取適當(dāng)?shù)奈灰婆c面力由于此環(huán)形區(qū)域既不包含原點(diǎn)也沒有無窮遠(yuǎn)處做邊界,因此只滿足位移單值條件即可.因此預(yù)期可以得到6 個(gè)方程,解出6 個(gè)未知系數(shù).所選取的位移試函數(shù)可以分為兩類,一類如表3 所示,一類如表2 所示.先考慮表3 中的這時(shí)有式(27)成立,將其代入式(39)得到

式(41)與式(27)形式完全類似,從數(shù)學(xué)上看只是做了符號替換,r=a換成了r=ρ、式(35)換成了式(40),彈性常數(shù)也相同.結(jié)果上也完全類似,式(29)替換相應(yīng)的符號得式(40)

從式(44)得知,夾雜內(nèi)部的應(yīng)力與界面處相同,不隨極徑變化.利用式(44)簡化式(42),并用r替換ρ,得到

式(45)代入式(40),并利用式(30),得到夾雜內(nèi)的位移與應(yīng)力

結(jié)果與文獻(xiàn)[17]相同.同邊界元法一樣,此邊界積分法也可以用于求解振動問題[35].

3 結(jié)論

本文參考了邊界元的思路,從Somigliana 等式出發(fā),利用格林函數(shù)的不唯一性,得到了適用于本文的邊界積分方程.可以用此方程求解彈性力學(xué)問題的解析解.求解步驟也與邊界元類似,先求出未知邊界位移與面力,再求域內(nèi)位移與面力.對于含有有限邊界的平面彈性問題來說,第一步可以將邊界上位移、面力進(jìn)行傅里葉級數(shù)展開,利用三角函數(shù)的正交性,求解邊界上未知位移與面力.第二步類似矩陣傳遞法.知道邊界上的位移與面力后,虛設(shè)邊界通過域內(nèi)的點(diǎn),再次利用第一步的方法求解域內(nèi)位移與應(yīng)力.

為了驗(yàn)證本方法,具體說明本方法求解思路.求解了單向均布應(yīng)力下平板含圓形夾雜問題,得到了精確解.在特殊情況下可以化為含圓孔或剛性夾雜的解析解.結(jié)果說明了此方法的有效性.