陳莘莘，王 崴

(華東交通大學(xué) 土木建筑學(xué)院，南昌 330013)

1 引 言

作為塑性力學(xué)最具實(shí)用意義的分支之一，結(jié)構(gòu)極限分析的目的是確定工程結(jié)構(gòu)的極限荷載[1]。結(jié)構(gòu)極限分析一般不需要知道載荷變化的歷史，相對于傳統(tǒng)的彈塑性增量分析往往效率更高更適用。雖然極限分析的數(shù)值方法[2-4]得到了迅速發(fā)展，但發(fā)展高效準(zhǔn)確和切實(shí)可行的極限分析數(shù)值方法仍是當(dāng)前研究的熱點(diǎn)。

作為一種新穎而獨(dú)特的數(shù)值計(jì)算方法，無網(wǎng)格法[5]近似函數(shù)不依賴于網(wǎng)格，且在節(jié)點(diǎn)不規(guī)則分布時(shí)，不會損失多少計(jì)算精度，從而日益得到重視，并得到不斷發(fā)展。目前，無網(wǎng)格法的種類已十分繁多，比較典型的有無單元伽遼金法[6,7]、重構(gòu)核粒子法[8]、無網(wǎng)格局部Petrov-Galerkin法[9]和自然單元法[10,11]等。其中，發(fā)展較晚的自然單元法是一種以自然鄰近插值作為試函數(shù)，以Galerkin法作為加權(quán)殘量法而得到的新型無網(wǎng)格方法。與大多數(shù)無網(wǎng)格法不同，自然單元法的形函數(shù)具有插值性且在邊界節(jié)點(diǎn)間是線性變化的，從而可以直接施加本質(zhì)邊界條件。此外，自然單元法的形函數(shù)構(gòu)造簡單，不涉及復(fù)雜的矩陣求逆運(yùn)算，更不需要任何人為的參數(shù)。蓋珊珊等[12]提出了線彈性斷裂問題的自然單元法。董軼等[13]將應(yīng)力雜交的思想引入自然單元法中，提出了彈性力學(xué)問題的雜交自然單元法。丁道紅等[14]提出了材料非線性問題分析的自然單元法。周書濤等[15]提出了平面問題極限上限分析的自然單元法。但目前尚未見到自然單元法針對軸對稱結(jié)構(gòu)極限上限分析的研究成果。

本文根據(jù)極限分析上限定理，詳細(xì)推導(dǎo)了軸對稱結(jié)構(gòu)極限上限分析的自然單元法求解格式，并編制了相應(yīng)的計(jì)算程序。通過典型算例的計(jì)算和對比分析，驗(yàn)證了所提軸對稱結(jié)構(gòu)極限上限分析方法的有效性和合理性。

2 無網(wǎng)格自然單元法

無網(wǎng)格自然單元法是采用自然鄰近插值構(gòu)造函數(shù)近似空間的一種新型數(shù)值方法，具有不依賴網(wǎng)格的優(yōu)勢，同時(shí)形函數(shù)構(gòu)造簡單，且在邊界上滿足插值性質(zhì)。自然鄰近插值按插值基函數(shù)的不同可分為Sibson插值[16]和Laplace插值[17]。

2.1 Voronoi結(jié)構(gòu)和Delaunay三角形

設(shè)平面中有一點(diǎn)集N={x1,x2,x3,…,xN}，則任意點(diǎn)xi的一階Voronoi 結(jié)構(gòu)實(shí)際上是以點(diǎn)xi為最近離散點(diǎn)的空間位置的集合，其數(shù)學(xué)表達(dá)式為

Ti={x∈R2d(x,xi)≤d(x,xj),?j≠i}

(1)

二階Voronoi結(jié)構(gòu)Ti j表示以xi為最近點(diǎn)，以xj為次近點(diǎn)的空間位置的集合，其數(shù)學(xué)表達(dá)式為

Ti j={x∈R2d(x,xi)

?j≠i≠k}

(2)

圖1所示為平面7個(gè)節(jié)點(diǎn)的一階Voronoi結(jié)構(gòu)和待插點(diǎn)x的二階Voronoi結(jié)構(gòu)。作為Voronoi結(jié)構(gòu)的對偶，Delaunay三角化是把有公共Voronoi結(jié)構(gòu)邊界的節(jié)點(diǎn)連接起來將整體域劃分成三角形區(qū)域。Delaunay三角形的一個(gè)重要特性是每個(gè)Delaunay三角形的外接圓中無任何其他節(jié)點(diǎn)，即空圓準(zhǔn)則，常用于確定計(jì)算點(diǎn)的自然鄰近點(diǎn)。

圖1 點(diǎn)x的一階和二階Voronoi結(jié)構(gòu)

Fig.1 First-order and second-order Voronoi cells aboutx

2.2 Sibson插值

設(shè)插值點(diǎn)x的一階Voronoi結(jié)構(gòu)Tx和二階Voronoi結(jié)構(gòu)Tx i的面積分別為A(x)和Ai(x)，則插值點(diǎn)x的形函數(shù)及其導(dǎo)數(shù)為[16]

(3)

(4)

定義了各節(jié)點(diǎn)的插值函數(shù)后，位移場函數(shù)u(x)可寫為

(5)

式中ui表示點(diǎn)x周圍n個(gè)自然鄰近點(diǎn)的節(jié)點(diǎn)位移。

3 無搜索迭代算法

(6)

(7)

ui,i=0, inV

(8)

ui=0， onΓu

(9)

式中σs表示材料的屈服應(yīng)力,εi j表示應(yīng)變率。為了簡便起見，通常將位移速度場和應(yīng)變率分別稱為位移和應(yīng)變。

將軸對稱面Ω離散成N個(gè)任意分布的節(jié)點(diǎn)，則應(yīng)變矩陣ε可以寫為

ε=BU

(10)

式中U為節(jié)點(diǎn)位移列向量，且

(11)

(12)

對于軸對稱問題，有

εi jεi j=εTDε=UTBTDBU=UTKU

(13)

式中

(14)

根據(jù)式(13)，目標(biāo)函數(shù)(6)可離散為

(15)

FTU=1

(16)

式中F為等效載荷列向量。不可壓條件(8)可表示為

ui,i=εr+εz+εθ=BcU=0

(17)

式中

(18)

(19)

進(jìn)一歩，式(17)可以表達(dá)為

(εr+εz+εθ)2=UTKcU=0

(20)

(21)

經(jīng)過上述推導(dǎo)，離散化的軸對稱結(jié)構(gòu)極限上限分析的數(shù)學(xué)規(guī)劃格式為

(22)

s.t.FTU=1

(23)

UT(Kc)iU=0 (i∈I)

(24)

采用拉格朗日乘子法將歸一化條件引入到目標(biāo)函數(shù)中，則式(22～24)所表示的數(shù)學(xué)規(guī)劃問題可寫為

(25)

s.t.UT(Kc)iU=0 (i∈I)

(26)

式中q為拉格朗日乘子。根據(jù)式(25,26)的極小化條件，可得

(27)

FTU=1

(28)

UT(Kc)iU=0 (i∈I)

(29)

式(27)成立的前提是對每個(gè)積分點(diǎn)都有

UTKiU≠0 (i∈I)

(30)

為便于求解式(27～29)所示的非線性方程組，構(gòu)造如下的迭代格式，

(31)

FTUk=1

(32)

(33)

式中Uk和qk分別表示第k迭代步的節(jié)點(diǎn)位移向量和拉格朗日乘子，Uk - 1是第k-1迭代步的節(jié)點(diǎn)位移向量。然而，當(dāng)某些積分點(diǎn)上的應(yīng)變?yōu)?時(shí)，式(31)是沒有意義的。因此，在進(jìn)行第k歩迭代前，需將全部積分點(diǎn)分為剛性區(qū)子集Rk和塑性區(qū)子集Pk：

I=Pk∪Rk

(34)

(35)

(36)

綜上所述，極限上限分析的離散數(shù)學(xué)規(guī)劃格式可表達(dá)為

(37)

s.t.FTUk=1

(38)

(39)

(40)

采用拉格朗日乘子法和罰函數(shù)法分別將式(38～40)引入到目標(biāo)函數(shù)后，式(37～40)的數(shù)學(xué)規(guī)劃格式可以等效成易求解的極小化問題：

(41)

式中α和β為罰因子，通常取104～108可得到較高的計(jì)算精度[15]。根據(jù)以上構(gòu)想，極限上限分析的迭代算法可歸納為

第0步:假設(shè)初始狀態(tài)下整個(gè)結(jié)構(gòu)無剛性區(qū)，求解如下的極小問題，

(42)

s.t.FTU0=1

(43)

(44)

通過求解以上的極小化問題，可得到初始的極限載荷乘子為

(45)

第k步：根據(jù)第k-1步的位移向量Uk - 1確定剛性區(qū)子集Rk和塑性區(qū)子集Pk，并求解如下極小化問題，

(46)

s.t.FTUk=1

(47)

(48)

(49)

通過求解以上的極小化問題，可得到第k歩的極限載荷乘子為

(50)

對于預(yù)先給定的誤差容限ε，如果滿足

(51)

則迭代終止。本文誤差容限取ε=2.0×10-4。

4 數(shù)值算例

4.1 厚壁圓筒

考慮一受均布內(nèi)壓P作用的厚壁圓筒，其內(nèi)半徑為a，外半徑為b，其極限載荷解析解為[1]

(52)

對于不同的b/a，采用不同的節(jié)點(diǎn)布置方案進(jìn)行計(jì)算，圖2給出了b/a=2時(shí)的節(jié)點(diǎn)布置方案。

圖2b/a=2時(shí)厚壁圓筒的節(jié)點(diǎn)布置

Fig.2 Nodal distribution for a thick-walled cylinder whenb/a=2

圖3給出了厚壁圓筒極限載荷的本文數(shù)值解和解析解的比較。顯然，本文數(shù)值解與解析解吻合很好，說明本文方法的計(jì)算精度較高。

圖3 厚壁圓筒極限載荷數(shù)值解與解析解的比較

Fig.3 Numerical limit loads compared with analytical solutions for the thick-walled cylinder

4.2 厚壁球殼

為進(jìn)一歩驗(yàn)證本文方法的有效性，考慮一受均布內(nèi)壓P作用的厚壁球殼，如圖4所示，其內(nèi)外半徑分別為a和b，其極限載荷的解析解為[1]

圖4 厚壁球殼

Fig.4 A thick-walled spherical shell

P=2σsln (b/a)

(53)

此問題是球?qū)ΨQ問題，現(xiàn)將其作為軸對稱問題進(jìn)行求解。對于不同的b/a，采用不同的節(jié)點(diǎn)布置方案進(jìn)行計(jì)算。圖5給出了b/a=1.3時(shí)，厚壁球殼的節(jié)點(diǎn)布置方案。圖6給出了b/a=1.05～1.3時(shí)，厚壁球殼的極限載荷數(shù)值解和解析解的比較。顯然，本文數(shù)值解和解析解吻合很好。

圖5b/a=1.3時(shí)厚壁球殼的節(jié)點(diǎn)布置

Fig.5 Nodal distribution for a thick-walled spherical shell whenb/a=1.3

圖6 厚壁球殼極限載荷數(shù)值解與解析解的比較

Fig.6 Numerical limit loads compared with analytical solutions for the thick-walled spherical shell

4.3 圓柱-圓錐形組合筒

本算例考慮一個(gè)圓柱-圓錐形組合筒，結(jié)構(gòu)的尺寸和受載情況如圖7所示，這也是一個(gè)軸對稱問題。

圖7 圓柱-圓錐組合筒

Fig.7 Cylindrical-conical combined cylinder

采用圖8所示的節(jié)點(diǎn)布置方案計(jì)算得到了該組合筒頂部均布荷載P的極限載荷為1.065σs。圖9 給出了極限上限分析的迭代收斂示意圖，可以看出，本文算法具有較好的收斂性。

圖8 圓柱-圓錐組合筒的節(jié)點(diǎn)布置

Fig.8 Nodal distribution for a cylindrical-conical combined cylinder

圖9 圓柱-圓錐形組合筒的極限上限分析迭代過程

Fig.9 Iterative process for upper bound limit analysis of the cylindrical-conical combined cylinder

5 結(jié) 論

無網(wǎng)格自然單元法不需要對求解域進(jìn)行網(wǎng)格劃分，同時(shí)本質(zhì)邊界條件的施加十分方便，是一種兼具了有限元法和無網(wǎng)格法優(yōu)點(diǎn)的數(shù)值方法。本文根據(jù)極限上限分析定理，采用自然鄰近插值構(gòu)造軸對稱結(jié)構(gòu)的位移場，并通過將積分點(diǎn)劃分為剛性區(qū)和塑性區(qū)子集的辦法克服目標(biāo)函數(shù)非線性且不可微的困難，建立了軸對稱結(jié)構(gòu)極限上限分析的自然單元法。本文的數(shù)值計(jì)算結(jié)果表明，采用自然單元法進(jìn)行軸對稱結(jié)構(gòu)的極限上限分析是可行的，具有前處理方便、穩(wěn)定性好和精度高等優(yōu)點(diǎn)。本文方法不僅推廣了自然單元法的應(yīng)用范圍，而且對拓展軸對稱結(jié)構(gòu)極限上限分析的數(shù)值方法有著積極的意義。