基于自然單元法的軸對稱結(jié)構極限上限分析

陳莘莘，王 崴

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

1 引 言

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

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

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

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

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

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

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

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

(1)

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

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

?j≠i≠k}

(2)

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

圖1 點x的一階和二階Voronoi結(jié)構

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

2.2 Sibson插值

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

(3)

(4)

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

(5)

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

3 無搜索迭代算法

(6)

(7)

ui,i=0, inV

(8)

ui=0， onΓu

(9)

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

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

ε=BU

(10)

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

(11)

(12)

對于軸對稱問題，有

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

(13)

式中

(14)

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

(15)

FTU=1

(16)

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

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

(17)

式中

(18)

(19)

進一歩，式(17)可以表達為

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

(20)

(21)

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

(22)

s.t.FTU=1

(23)

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

(24)

采用拉格朗日乘子法將歸一化條件引入到目標函數(shù)中，則式(22～24)所表示的數(shù)學規(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)成立的前提是對每個積分點都有

UTKiU≠0 (i∈I)

(30)

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

(31)

FTUk=1

(32)

(33)

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

I=Pk∪Rk

(34)

(35)

(36)

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

(37)

s.t.FTUk=1

(38)

(39)

(40)

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

(41)

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

第0步:假設初始狀態(tài)下整個結(jié)構無剛性區(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)

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

(51)

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

4 數(shù)值算例

4.1 厚壁圓筒

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

(52)

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

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

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

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

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

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

4.2 厚壁球殼

為進一歩驗證本文方法的有效性，考慮一受均布內(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)將其作為軸對稱問題進行求解。對于不同的b/a，采用不同的節(jié)點布置方案進行計算。圖5給出了b/a=1.3時，厚壁球殼的節(jié)點布置方案。圖6給出了b/a=1.05～1.3時，厚壁球殼的極限載荷數(shù)值解和解析解的比較。顯然，本文數(shù)值解和解析解吻合很好。

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

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 圓柱-圓錐形組合筒

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

圖7 圓柱-圓錐組合筒

Fig.7 Cylindrical-conical combined cylinder

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

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

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