文穎，曾慶元

(中南大學 土木工程學院，湖南 長沙，410075)

傳統(tǒng)的結構彈塑性極限強度分析一般采用經典的塑性鉸理論，但該方法只能計算極限荷載，無法考慮結構構件彈塑性發(fā)展的全過程[1]，特別是用于大型復雜結構的彈塑性分析時，需事先知道塑性鉸出現的先后次序及位置，分析過程繁瑣。隨著矩陣結構分析理論應用的深入和高速電子計算機的快速發(fā)展，有限元方法已廣泛用于求解一般結構彈塑性響應分析問題[2]。目前，國內外研究者采用的彈塑性有限元分析模型主要分為2類：塑性鉸法和塑性區(qū)域法。塑性鉸法假定塑性變形集中在構件端部而構件中間部分保持彈性工作，該法分析過程簡單，但是，考慮構件塑性發(fā)展過程非常粗糙和近似。為了保證塑性鉸法分析模型與構件實際彈塑性分布盡量接近，提高分析精度，Kim等[3]提出三鉸梁模型，在構件中點增設1個塑性鉸考慮塑性沿構件長度方向的發(fā)展，采用纖維單元模型考慮塑性沿構件截面擴展；Attalla等[4]提出偽塑性鉸法，通過引入基于內力的初始屈服面和繼生屈服面等概念以及構件整體受力平衡條件考慮塑性在構件內部發(fā)展全過程；Ren等[5]提出剛體-彈簧模型，建立平面框架結構內力塑性屈服面方程以取代基于應力的塑性分析模型，提高分析效率；Krishnan[6]提出改進的塑性鉸法，認為構件塑化部分不再局限于節(jié)點，而是集中分布在靠近節(jié)點的有限區(qū)域，采用靜力凝聚消除構件內部自由度；顏全勝[7]提出內力塑性系數法，假定構件塑性系數沿長度方向線性分布；Liew等[8]提出精細化塑性鉸法，引入截面塑性發(fā)展系數考慮塑性沿截面的發(fā)展，假定該系數沿桿長拋物線變化近似考慮塑性沿長度方向的分布。塑性區(qū)域法[9-10]將構件內部分成完全彈性區(qū)、塑性區(qū)和彈塑性區(qū)，考慮塑性沿整個構件的逐步開展，由于沒有引入簡化近似假定，一般被視為結構彈塑性分析的精確方法。其缺陷是由于將構件沿截面以及長度方向劃分為若干有限段，增加數據處理量和計算機時而使得分析過程繁瑣。在此，本文作者根據有限元分析的基本概念，提出桿系結構彈塑性分析的單元節(jié)點截面增量內力塑性系數法，通過引入單元平衡條件，簡便地保證了單元的整體平衡，也自動考慮塑性沿單元長度方向的發(fā)展。該法顯式推導了單元彈塑性增量割線剛度矩陣，避免了數值積分，提高了分析效率?；谠摲椒?，提出求解結構非線性增量平衡方程的基于增量割線剛度的直接迭代算法。該算法在平衡迭代中不需計算結構切線剛度矩陣，也不需計算結構整體不平衡力，簡化了迭代過程。

1 框架結構彈塑性分析的增量內力塑性系數法

1.1 彈塑性增量平衡方程

假定任意離散化結構在節(jié)點荷載P作用下節(jié)點變位為u，結構內部應力為σ，應變?yōu)棣?。假定結構產生虛位移 δu和虛應變 δε，根據虛功原理，結構的虛應變能等于外力虛功，可以得到：

式中：ΔP為施加的荷載增量；Δu為節(jié)點位移增量；Δσ和Δε分別為結構應力和應變增量。由虛功原理，忽略幾何位形改變對結構平衡的影響，在荷載 P+ΔP作用下，結構平衡方程如下：

將(2)式減(1)式，可得：

式(3)實際上就是小變形條件下增量形式的虛功原理，反映了結構增量平衡條件。不考慮材料非線性的影響，由式(3)可以直接得出線彈性結構平衡方程：

考慮彈塑性分析,由式(3)可得到結構彈塑性增量平衡方程：

式中：Kep為結構彈塑性增量割線剛度矩陣。

1.2 梁單元彈塑性剛度列式的增量內力塑性系數法

從有限元法的物理概念來說，式(5)表示結構節(jié)點截面的抗力增量和外荷載增量的平衡。隨著外荷載的增加，結構材料彈塑性行為集中反映在節(jié)點截面抗力增量降低，這是相對于單元為完全彈性時的抗力增量而言的。為此，可以定義節(jié)點截面彈塑性抗力增量和彈性抗力增量的比值(稱為單元增量內力塑性系數)來反映構件彈塑性發(fā)展情況，并利用該系數來確定單元彈塑性剛度矩陣，最后采用“對號入座”法則[11]得到結構彈塑性剛度矩陣。

首先，考慮單元為完全彈性時的情況，單元的平衡條件可寫為：

式中，等號右端各項表示節(jié)點截面彈性抗力增量，下標“e”表示單元處于彈性狀態(tài)，上標“e”表示任意結構單元。若考慮單元已經屈服，則單元的平衡條件可改寫為：

式中，等號右端各項表示節(jié)點截面彈塑性抗力增量，下標“ep”表示單元處于彈塑性狀態(tài)。需要說明的是：與式(6)類似，式(7)的表現形式可由經典有限單元法得到。

式(6)直接給出了彈性節(jié)點內力的計算表達式。一般假定單元的彈塑性狀態(tài)只由正應力決定，彎曲剪應力不影響屈服，對于平面結構，多個實例已證實該假定合理性[2,9]。為了考慮塑性區(qū)沿單元截面高度的擴展，這里采用纖維單元法[2,9,12-14]即將單元截面細分成很多小塊 (如圖1所示)，用小塊中點的應變、應力代表小塊的應變、應力水平。累加所有小塊對于截面內力的貢獻，則有：

其中：下標i為節(jié)點截面編號；r為截面上小塊編號；ΔAr為第r小塊的面積；xr和yr為第r小塊幾何中心的截面形心坐標；Δσir為i節(jié)點截面上第r小塊彈塑性正應力增量，該應力增量一般根據材料單向應力狀態(tài)下本構關系模型和應變增量得到。

圖1 梁截面分塊示意圖Fig.1 Discretized cross-section of beam

需要說明的是：目前比較流行的考慮塑性沿截面開展的途徑是通過截面數值積分[10]，盡管數值積分與纖維單元法相比減少了數據存儲量和計算工作量，但應用范圍較窄。例如，對于組合材料(鋼筋混凝土)截面的塑性分析就不適用，因此，本文考慮到分析方法的通用性，選擇纖維單元法?；谏鲜龇治鼋Y果，定義下列單元增量塑性系數：

式中：αi和 βi分別為i節(jié)點截面的軸力和彎矩增量塑性系數；αj和 βj與 j節(jié)點截面相關。這里“增量”一詞的用意是強調計算這些塑性系數是基于內力增量而不是全量。

與正應力的計算不同，彈塑性彎曲剪應力無法通過材料本構關系得到。為了簡化，一些研究者通過計算彈性彎曲剪力增量作為替代[2]，這顯然破壞了單元整體平衡條件，導致出現計算誤差。為了克服這種不足，可以利用處于彈塑性狀態(tài)的單元平衡條件計算節(jié)點截面剪力增量，考慮式(7)等號右端的彈塑性抗力增量在單元內部自相平衡，有：

式中：各節(jié)點內力的正方向和傳統(tǒng)有限元的規(guī)定一致；Le為單元的長度。將式(10)代入式(7)有：

由式(6),(9)和(11)，可得彈塑性剛度矩陣和線彈性剛度矩陣有如下關系：

式中：α為單元增量內力塑性系數矩陣，它只與軸力和彎矩的增量塑性系數有關。

1.3 非線性增量平衡方程求解步驟

用增量內力塑性系數法分析結構彈塑性響應的具體步驟如下。

(1) 假定結構彈性極限節(jié)點荷載為P0，節(jié)點位移、內部應變和應力分別為u0，ε0和σ0。在結構開始屈服時就采取增量加載模式，施加荷載增量ΔPi，i表示荷載增量步。選取初始剛度

(5) 繼續(xù)加載，i=i+1，求解流程回到步驟(1)，直至結構達到承載極限狀態(tài)為止。結構極限狀態(tài)由本步剛度參數法[13]確定。

增量-迭代求解過程的圖解形式如圖2所示。

圖2 基于增量割線剛度的增量-迭代過程Fig.2 Incremental secant stiffness based on incremental-iterative procedure

2 數值算例

算例1 文獻[1]計算了簡支梁在跨中集中荷載作用下的彈塑性響應，為了與文獻[1]中的數值結果進行比較，本算例采用相同的幾何和材料參數(見文獻[1])。文獻[1]將梁橫截面等分成8層，沿梁軸方向劃分為20個梁段，而這里也將截面分成上下8層，沿梁跨方向劃分為 2個梁段，公共節(jié)點位于跨中。簡支梁荷載-撓度曲線如圖3所示。文獻[1]采用分段分塊變剛度法計算所得到的結構極限荷載P≈2.8×106N，由增量內力塑性系數法得P≈2.689 3×106N，而結構彈塑性極限強度的理論值為P≈2.688×106N。文獻[1]中的極限荷載值偏大的原因可能是計算剪切恢復力時用采用彈性剪力而不是彈塑性剪力，導致結構平衡條件只在節(jié)點處滿足，結構單元內部平衡條件不滿足。

算例2 如圖4所示，一單層兩榀平面框架承受面內若干集中荷載作用。至于其彈塑性響應，Franchi等[15]采用數學規(guī)劃法計算彈塑性極限荷載，Wong等[16]采用分段線性化內力屈服面法分析結構彈塑性響應全過程。該框架幾何與材料見在文獻[15]。這里將所有柱子離散為1個梁單元，橫梁用2個單元模擬。所有桿件節(jié)點截面沿其厚度方向統(tǒng)一劃分成10個條帶。圖5所示為荷載-框架頂部側移曲線，并與現有其他文獻結果進行對比。

圖3 簡支梁的荷載-位移曲線Fig.3 Load deflection curves of simple supported beam

圖4 平面框架Fig.4 One storey portal frame of simple supported beam

圖5 平面框架的荷載-側移曲線Fig.5 Load-top horizontal displacement curve of portal frame

按本文方法算得懸臂梁的極限荷載因子λ=0.768 7，而文獻[15]和文獻[16]給出的極限荷載因子分別為0.768 5和0.768 9，它們都非常接近。值得注意的是：當荷載因子λ=0.745時，3和4號節(jié)點處塑性鉸出現彈性卸載。這也是結構發(fā)生塑性變形后內力重新調整的結果，這一現象在文獻[15]和[16]中得到證明。

3 結論

(1) 提出了平面框架彈塑性分析的增量內力塑性系數法。通過引入單元平衡條件，自動地考慮塑性沿單元長度方向的發(fā)展，得到的單元彈塑性增量割線剛度矩陣可以顯式給出，無需數值積分，提高了計算效率。

(2) 增量內力塑性系數法的計算步驟與經典彈塑性有限元分析法的相比無需計算結構切線剛度矩陣，也不需在每次迭代中計算不平衡力，簡化了迭代計算過程。

(3) 增量內力塑性系數法不僅物理概念明確，而且計算結果正確、可靠。在分析一般框架結構時，每個構件用至多2個梁單元離散便可以獲得滿意結果。

[1] 潘家英, 張國政, 程慶國. 大跨度橋梁極限承載力的幾何與材料非線性耦合分析[J]. 土木工程學報, 2000, 33(1): 5-9.PAN Jia-yin, ZHANG Guo-zheng, CHENG Qing-guo.Geometrical and material nonlinear analysis for determining ultimate load capacity of long-pan bridges[J]. China Civil Engineering Journal, 2000, 33(1): 5-9.

[2] Hinton E, Owen D.R.J. Finite elements in plasticity: Theory and practice[M]. Swansea: Pineridge Press, 1979: 3-7.

[3] Kim S E, Kim Y, Choi S H. Nonlinear analysis of 3-D steel frames[J]. Thin Walled Structures, 2001, 39(6): 445-461.

[4] Attalla M R, Deierlein G G, McGuire W. Spread of plasticity:Quasi-plastic hinge approach[J].Journal of Structural Engineering, 1994, 120(8): 2451-2473.

[5] Ren W X, Tan X G, Zheng Z C. Nonlinear analysis of plane frames using rigid body-spring discrete element method[J].Computers & Structures, 1999, 71(1): 105-119.

[6] Krishnan S. Modified elastofiber element for steel slender column and brace modeling[J]. Journal of Structural Engineering,ASCE, 2010, 136(11): 1350-1366.

[7] 顏全勝. 大跨度鋼斜拉橋極限承載力分析[D]. 長沙: 長沙鐵道學院土木工程系, 1994: 33-37.YAN Quan-sheng. Analysis of ultimate load capacity of the steel cable-stayed bridge[D]. Changsha: Changsha Railway University.Department of Civil Engineering, 1994: 33-37.

[8] Liew J Y R, White D W, Chen W F. Second-order refined plasti-hinge analysis for frame design: Part I and Part II [J].Journal of Structural Engineering, ASCE, 1993, 119(11):3196-3237.

[9] Toma S. European calibration frames for second-order inelastic analysis[J]. Engineering Structures, 1992, 14(1): 7-14.

[10] Teh L H, Clarke M J. Plastic-zone analysis of 3D steel frames using beam elements[J]. Journal of Structural Engineering,ASCE, 1999, 125(11): 1328-1337.

[11] 曾慶元, 楊平. 形成矩陣的“對號入座”法則與桁梁空間分析的桁段有限元法[J]. 鐵道學報, 1986, 8(2): 48-59.ZENG Qing-yuan, YANG Ping. The “set-in-right-position” rule for formulating matrix and the truss finite element method for spatial truss analysis[J]. Journal of the Railway Society, 1986,8(2): 48-59.

[12] Zhang Y X, Bradford M A, Gilbert R I. A new triangular layered plate element for the non-linear analysis of reinforced concrete slabs[J]. Communications in Numerical Methods in Engineering,2005, 22(7): 699-709.

[13] 韋成龍. 大跨度板桁結合主梁斜拉橋極限承載力分析[D]. 長沙: 中南大學土木工程學院, 2003: 80-83.WEI Cheng-long. Analysis of the ultimate load for long span cable-stayed bridge with a plate-truss composite beam[D].Changsha: Central South University. School of Civil Engineering,2003: 80-83.

[14] 呂烈武, 沈世釗, 沈祖炎, 等. 鋼結構構件穩(wěn)定理論[M]. 北京: 中國建筑工業(yè)出版社, 1983: 246-252.Lü Lie-wu, SHEN Shi-zhao, SHEN Zu-yan, et al. The stability theory of steel members[M]. Beijing: China Architecture &Building Press, 1983: 246-252.

[15] Franchi A, Cohn M Z. Computer analysis of elastic-plastic structures[J]. Computer Methods in Applied Mechanics and Engineering, 1980, 21(3): 271-294.

[16] Wong M B, Loi F T. Yield surface linearization in elasto-plastic analysis[J]. Computers and Structures, 1987, 26(6): 951-956.