基于有限元法的軸對稱結構極限下限分析

王 崴，王顏輝，孫 萬，李睿智，鞏少鵬

（山西能源學院，山西 晉中 030600）

0 引言

軸對稱結構廣泛應用于航空航天、化工、機械等工程領域中。通過確定各類軸對稱結構的極限載荷，可以為工程設計和安全評定提供準確可靠的理論依據和參考[1-3]。傳統的“許用應力法”對結構極限承載能力的評估往往趨于保守，而彈塑性增量法由于計算難度大、計算成本高等缺點使得其很難應用于實際工程中。相較于前兩種方法，極限分析能夠在考慮材料塑性性質的基礎上得到真實反映材料安全裕度的重要參數，而且計算過程不涉及結構的加載歷史，執行起來簡單易行，是研究軸對稱結構承載能力的一種直接有效的方式。

目前，將數值方法與數學規劃理論相結合已成為分析極限分析問題的主流方法。這種方法將極限分析歸納為一個非線性數學規劃問題，但維數障礙使得該問題的求解異常困難。眾多學者為克服維數障礙問題進行了研究。李春光等[4]對四邊形網格的平衡方程進行積分，減少了自變量和約束方程的數量；錢向東等[5]對路徑跟蹤內點法進行改進，在算法中引入子迭代程序，使大量計算可以在單元一級完成；秦方等[6]將下限分析列式轉化為圓錐二次優化問題進行數值求解。

本文基于塑性極限下限分析理論，結合有限元離散技術，構造了求解軸對稱極限下限問題的相應算法。通過減縮基技術將下限分析列式歸結為一系列數學規劃子問題，并在每個子問題中通過彈塑性增量分析得到的結果形成自平衡應力場搜索子空間。根據數學規劃問題中變量少、約束條件多的特點，采用復合形法進行求解，并編制了相應的計算程序。最后通過對兩個軸對稱結構極限下限載荷系數的計算和對比，驗證了本文方法的可行性和有效性。

1 數學規劃格式

現考慮一個理想彈塑性軸對稱結構V，為簡便計算，取軸對稱結構中的任一對稱面S進行分析。極限下限定理可以表述為[7-9]：結構未達到極限狀態的充要條件是，存在一個自平衡應力場ρij（x），與虛擬彈性應力場（x）相疊加后產生的總應力場處處不違反屈服條件f[·]。而極限下限載荷就是結構在未達到極限狀態時所能承受的最大載荷，由此得到其數學規劃格式為：

其中χ為待求的極限載荷系數，nj表示邊界單位外法線向量，σs是Mises 屈服應力。表達式（2）為結構應滿足的屈服條件，（3）和（4）分別表示自平衡應力場需要滿足的應力邊界條件和平衡方程。

對于一個連續體結構，為減少數學規劃式（1）～（4）的未知數和約束條件的數量，必須對結構進行離散化處理。根據虛位移原理[10]，應力邊界條件（3）和平衡方程（4）的等效積分弱形式表示為：

采用有限元法進行離散，假定將對稱面S劃分為N個單元的的組合，則每個單元中的應變列向量εe的計算a表達式為：

式中，Be為單元應變位移關系矩陣，u表示單元節點位移列向量，根據式（6）可以得到：

將單元節點應變向量εe組合成為總體節點應變向量ε，并代入式（5），由Gauss 積分法得到：

式中，U表示結構的總體節點位移向量，B為擴展之后與U同階的應變位移關系矩陣，N是離散之后結構中的單元總數，G為每個單元內部Gauss 積分點的數目，I表示全部Gauss 積分點的集合，（ρe）i為單元e中第i個Gauss 積分點的自平衡應力，ri和分別為單元e中第i個Gauss 積分點的徑向半徑和Jacobi 矩陣的值。考慮到δU的任意性，式（8）進一步寫為：

令xi表示離散之后第e個單元中第i個Gauss積分點的坐標值，則虛擬彈性應力場。綜合上述分析過程，離散之后的軸對稱極限下限分析數學規劃格式為：

2 減縮基技術

求解非線性規劃問題（12）的最大困難在于未知數和約束條件的數量較為龐大，因此本文采用減縮基技術[11]將自平衡應力場表示為一組自平衡應力基矢量的線性組合：

其中M1，M2，…，MK表示基矢量對應的待定系數。通過減縮基技術，整個極限問題的求解將在維數很低的自平衡應力場子空間中進行迭代。具體過程如下：首先根據4.1 節的計算得到結構的彈性極限載荷，之后在此基礎上施加一個載荷增量△χ1，結合4.2 節的步驟得到一組自平衡應力場基矢量，并通過復合形法[12]對該數學規劃問題進行求解，得到近似解χ1。由于此時規劃式中的未知數只有K+1 個（即待求的極限載荷系數和基矢量對應的K個待定系數），因此計算規模得到很大程度的減小。χ1的精度一般來說并不高，還需進行第2 次數學規劃問題的計算。以第1 次數學規劃問題得到的平衡應力場為基礎，施加第2 個載荷增量△χ2，同時將求解第1 次數學規劃問題得到的自平衡應力場作為第K個基矢量，由此得到一組新的自平衡應力場基矢量，并再次通過復合形法計算得到近似解χ2。如此重復進行，直到滿足收斂條件：

則計算終止，式中vol 表示誤差容限。

3 平衡應力場構造

3.1 虛擬彈性應力場

假設在軸對稱面S上分布著體力b= [br bz]T，面力邊界Γt作用有給定的面力，位移邊界Γu有已知位移，則通過加權余量法可得到軸對稱彈性力學問題控制方程的等效弱形式為：

式中，r代表軸對稱半徑，σ表示應力向量：

根據有限元法插值方法，將求解域內的位移場函數代入式（15）中，并注意到位移變分δu具有任意性，由此進一步推導出軸對稱結構線彈性問題的控制方程：

式中u代表求解域中所有節點的位移向量，K和f分別表示總體剛度矩陣和總體載荷向量。

3.2 自平衡應力場

約束條件式（10）表明，求解極限下限問題的關鍵在于找到一個能讓載荷系數χ取最大值的自平衡應力場。為了讓迭代得到良好的收斂效果，本文利用彈塑性增量分析中平衡迭代的結果構造自平衡應力基矢量，具體過程如下：

以第m次迭代為例，根據上一次迭代得到的載荷系數χ（m-1）和自平衡應力場ρm-1i，則各Gauss 積分點總應力為：

現在此基礎上增加一個載荷增量△χm使結構進一步屈服。為提高計算效率，本文選用修正的Newton-Raphson 迭代法[11]進行計算。則第1 次和第p迭代的平衡表達式為：

其中F代表基準載荷作用時對應的節點載荷向量，De為軸對稱問題的彈性矩陣。用（20）減去（19），得到：

式中ρp是一個與外載荷有關的自平衡應力場。將每一次迭代得到的ρp作為一個基矢量，就形成了自平衡應力場的搜索子空間。

4 算例驗證

4.1 厚壁圓筒

一個理想彈塑性厚壁圓筒結構，受到均布內壓P的作用，分別用r和R表示厚壁圓筒的內徑和外徑，其計算模型如圖1 所示。根據塑性力學理論，厚壁圓筒極限載荷的解析解計算表達式可寫為[13]：

圖1 受內壓厚壁圓筒

本文對R/r在2.0 ～4.0 的一組厚壁圓筒進行極限分析。對于不同的內外徑之比，設置不同的單元數量進行計算。圖2 所示為當R= 2r時有限元網格的劃分情況。為驗證本文算法的有效性，將采用本文算法得到的數值解與解析解進行對比，具體情況如圖3所示。從圖中可以看出，數值解和解析解形成的曲線吻合很好，說明了所提算法的有效性。

圖2 R/r=2 時厚壁圓筒單元布置

圖3 本文數值解與解析解的比較

4.2 厚壁球殼

對一個內徑為r，外徑為R的厚壁球殼進行分析，球殼內部受到均布載荷P的作用。圖4 為R=1.3r時厚壁球殼的計算模型，則球殼結構的極限載荷解析解計算表達式為[13]：

圖4 受內壓厚壁球殼

為進一步驗證所提方法的有效性，算例中計算了外徑與內徑之比為1.1～1.3 之間的一組厚壁球殼，計算結果如表1 所示。表1 結果顯示，采用本文算法得到的數值解和理論計算得到的解析解之間的誤差比較小，都在1%以內。為了探究基矢量數目對計算精度的影響，算例中以R/r= 1.3 時的厚壁球殼為研究對象，分別設置不同數目的基矢量進行計算，結果如表2 所示。表2 結果顯示，不同基矢量數目得到的計算結果均穩定在0.52 附近，這充分說明了本文算法的合理性和穩定性。

表1 厚壁球殼極限下限載荷系數的計算結果

表2 R=1.3r 時厚壁球殼不同基矢量數目的計算結果

5 結語

本文基于塑性極限分析理論，將有限元法應用于求解軸對稱結構的極限下限問題中，并建立了相應算法。通過將整個下限分析轉化為一系列未知數較少的非線性數學規劃子問題，解決了維數障礙問題。每個子問題中通過彈塑性增量法迭代的結果構造出相應的自平衡應力場，并采用復合形法直接進行求解。算例的計算結果表明，本文提出的算法具有數值穩定性好、計算精度高的優點。