余能原理基面力元法在桁架彈塑性分析中的應(yīng)用

李瑞雪，彭一江

(北京工業(yè)大學(xué)建筑工程學(xué)院，北京 100124)

余能原理基面力元法在桁架彈塑性分析中的應(yīng)用

李瑞雪，彭一江

(北京工業(yè)大學(xué)建筑工程學(xué)院，北京 100124)

基面力元法(BFEM)是一種新型的有限單元法，是以"基面力"來(lái)描述一點(diǎn)的應(yīng)力狀態(tài)。為了編制適用于桁架彈塑性分析的MATLAB基面力元程序，以余能原理基面力元法為基礎(chǔ)，將基面力元法的平面4節(jié)點(diǎn)的實(shí)體單元模型向桁架單元進(jìn)行退化，列出了退化桁架單元的柔度矩陣、單元軸力和節(jié)點(diǎn)位移的具體表達(dá)式，并對(duì)其控制方程進(jìn)行了必要的處理。根據(jù)退化的桁架單元模型，采用增量法的思想和線(xiàn)性強(qiáng)化彈塑性模型，以三桿超靜定桁架和復(fù)雜靜定桁架為算例進(jìn)行了計(jì)算。將計(jì)算結(jié)果與大型有限元分析軟件ABAQUS的T2D2計(jì)算結(jié)果進(jìn)行了對(duì)比分析，發(fā)現(xiàn)所用方法與ABAQUS的計(jì)算結(jié)果高度吻合，進(jìn)一步驗(yàn)證了基面力理論的正確性和程序的適用性。

塑性力學(xué)；有限單元法；基面力元法；余能原理；桁架；彈塑性分析

有限單元法[1]從正式提出至今已經(jīng)歷了半個(gè)多世紀(jì)的發(fā)展。在過(guò)去的幾十年中，以節(jié)點(diǎn)位移作為基本未知量的“位移型有限元”，因其清晰易懂的物理概念率先被發(fā)展并得到廣泛應(yīng)用[2]，它是以最小勢(shì)能原理作為理論基礎(chǔ)，而且要求它所采用的離散化結(jié)構(gòu)必須滿(mǎn)足位移協(xié)調(diào)性。但是隨著日益發(fā)展的需要，待解決的問(wèn)題也越來(lái)越多樣化和復(fù)雜化。不同的有限元模型，例如基于最小余能原理的平衡元[3]、雜交元[4]、非協(xié)調(diào)元[5]以及廣義協(xié)調(diào)元[6]等被一大批學(xué)者提出并被研究。因此，針對(duì)有限元問(wèn)題探索具有更簡(jiǎn)潔的應(yīng)力和變形描述體系，建立不依賴(lài)網(wǎng)格且收斂性好的新型有限元模型是一些學(xué)者所關(guān)心的。

2003年，GAO[7]提出了一個(gè)描述應(yīng)力狀態(tài)的新概念——“基面力(base force，BF)”,該基面力為一階張量，可以像其他傳統(tǒng)的二階張量一樣來(lái)描述一點(diǎn)的應(yīng)力狀態(tài)，其表達(dá)公式簡(jiǎn)潔，推導(dǎo)方程方便。在高玉臣[8-9]工作的基礎(chǔ)上，彭一江[10]開(kāi)始了基于基面力概念的新型有限元研究，分別對(duì)平面問(wèn)題、空間問(wèn)題和彈性大變形問(wèn)題建立了相應(yīng)的基面力元模型，并于2009年將該方法簡(jiǎn)稱(chēng)為“基面力元法(BFEM)”[11]。經(jīng)過(guò)近幾年的研究，基面力元法理論已經(jīng)不斷地應(yīng)用在許多領(lǐng)域中[11-13]，其相關(guān)研究工作也引起了國(guó)際同行的關(guān)注與研究[14-16]。

本文以余能原理基面力元法為基礎(chǔ)，對(duì)平面4節(jié)點(diǎn)的實(shí)體單元模型[17]進(jìn)行退化，得到了平面2節(jié)點(diǎn)的桁架基面力元模型，推導(dǎo)了其柔度矩陣的具體表達(dá)式，并進(jìn)一步對(duì)控制方程進(jìn)行了處理。利用增量法的思想和MATLAB編程語(yǔ)言，編制了桁架的彈塑性基面力元程序，結(jié)合典型的平面桁架算例進(jìn)行了數(shù)值計(jì)算，驗(yàn)證了該模型的正確性。

1 桁架單元的基面力元模型

1.1 柔度矩陣

對(duì)于平面問(wèn)題，考慮如圖1所示的平面4節(jié)點(diǎn)單元。I，J，K和L分別表示單元4個(gè)邊的邊中節(jié)點(diǎn)，TI，TJ，TK和TL分別為作用在相應(yīng)節(jié)點(diǎn)上的力，也稱(chēng)單元面力。

圖1 平面4節(jié)點(diǎn)單元Fig.1 Four-node plane element

單元的柔度矩陣可表示為

I,J=1,2,3,4。

(1)

退化的平面桁架單元如圖2所示，假設(shè)應(yīng)力均勻的分布在整個(gè)桿單元中，α，β分別表示單元的第α節(jié)點(diǎn)和第β節(jié)點(diǎn)，且α，β=1，2；Tα,Tβ為作用在第α節(jié)點(diǎn)和第β節(jié)點(diǎn)上的力。

圖2 平面桁架單元Fig.2 Planar truss element

從桁架結(jié)構(gòu)的基本假設(shè)上說(shuō)，不考慮其剪切能力，而泊松比的本質(zhì)是對(duì)剪切效果的研究，所以單純的桁架結(jié)構(gòu)分析不需要泊松比也可以求解出其軸向的力學(xué)響應(yīng)，所以本文假設(shè)泊松比v=0，于是桁架單元的余能可寫(xiě)為

(2)

平面桁架單元的柔度矩陣也可由平面4節(jié)點(diǎn)單元退化為

(3)

式中：A為單元的面積；U為單位張量。Pαβ為單元的第α節(jié)點(diǎn)和第β節(jié)點(diǎn)的徑矢rα和rβ的點(diǎn)積，有

U=ex?ex+ey?ey,

(4)

rαβ=rα·rβ=rαxrβx+rαyrβy。

(5)

最終可以得到Cαβ的具體形式

(6)

1.2 單元軸力表達(dá)式

根據(jù)文獻(xiàn)[17]，基線(xiàn)力Tα與Cauchy應(yīng)力張量σ的關(guān)系式為

(7)

由于桁架單元中假設(shè)只有軸力，且其應(yīng)力在整個(gè)桿中都是均勻分布的，單元的應(yīng)力可寫(xiě)為

(8)

本文計(jì)算軸力的思路一開(kāi)始是在求出桿端面力后沿著軸線(xiàn)方向合成得到，后來(lái)考慮到桁架不存在剪力，推斷其軸力必是沿軸線(xiàn)方向的主應(yīng)力與橫截面積的乘積，并且驗(yàn)證了推斷的正確性，所以有

(9)

軸力取不為零的主應(yīng)力與橫截面積的乘積，有

(10)

式中Acs為桿單元的橫截面面積。

1.3 控制方程及其處理

余能基面力元的控制方程可以寫(xiě)成

(11)

式(11)的第1式的物理意義是位移協(xié)調(diào)條件，主要包括：1)不同單元在連接節(jié)點(diǎn)處的相對(duì)位移為零；2)單元的位移要滿(mǎn)足邊界條件。

式(11)的第2式的物理意義是各單元力的平衡，即每個(gè)單元2個(gè)節(jié)點(diǎn)的力加起來(lái)為零。

式(11)的第3式的物理意義是各單元轉(zhuǎn)動(dòng)的平衡，即每個(gè)單元2個(gè)節(jié)點(diǎn)的力繞原點(diǎn)的力矩之和為零。

對(duì)于之前研究的平面實(shí)體單元，每個(gè)公共邊界由于只包含2個(gè)單元，可由計(jì)算機(jī)自動(dòng)判別并賦予兩單元接觸邊上大小相等、方向相反的單元面力T和-T，從而滿(mǎn)足這一約束條件。而對(duì)于桁架結(jié)構(gòu)來(lái)說(shuō)，每個(gè)公共節(jié)點(diǎn)至少包含2個(gè)及以上的多個(gè)單元，此時(shí)未知力的個(gè)數(shù)也會(huì)相應(yīng)增加，各個(gè)未知力不能再簡(jiǎn)單設(shè)置為等大反向。假如有n個(gè)桿交于一點(diǎn)，則需要假設(shè)的未知力數(shù)是n-1個(gè)。

桁架結(jié)構(gòu)的上述特點(diǎn)，導(dǎo)致控制方程與基面力元實(shí)體單元主要有兩點(diǎn)不同之處：1)必須考慮各個(gè)桿在交點(diǎn)處力的平衡；2)各個(gè)桿在交點(diǎn)處的位移協(xié)調(diào)，即各個(gè)桿的相對(duì)位移為零。

考慮一個(gè)由m個(gè)公共節(jié)點(diǎn)組成的桿系結(jié)構(gòu)，其部分系統(tǒng)與公共節(jié)點(diǎn)i相連的n個(gè)桿件及其節(jié)點(diǎn)荷載如圖3所示。

圖3 部分桿系系統(tǒng)Fig.3 Partial rods system

根據(jù)交點(diǎn)處力的平衡條件，有

(12)

對(duì)于i節(jié)點(diǎn)來(lái)說(shuō)，須設(shè)有n-1個(gè)未知力，第n個(gè)桿的力Tn可由其他桿的力來(lái)表示，即

Tn=-P-T1-T2-…-Tn-1。

(13)

各桿在交點(diǎn)處的位移協(xié)調(diào)條件由式(11)的第1式表示，對(duì)上述的部分桿系系統(tǒng)可展開(kāi)為

(14)

式(14)表達(dá)的物理意義可以理解為第k根桿與第n根桿在相對(duì)位移為零，即n根桿是連在一起的。

1.4 節(jié)點(diǎn)位移表達(dá)式

桁架結(jié)構(gòu)節(jié)點(diǎn)I的位移δI可由該點(diǎn)所包含的任意一個(gè)單元的支配方程求得，其表達(dá)式為

(15)

可得到節(jié)點(diǎn)位移δI的顯式表達(dá)式

δI=CIJ·TJ+λ+λ3ε·rI。

(16)

其矩陣形式可寫(xiě)為

(17)

根據(jù)推導(dǎo)的基面力元公式和控制方程，并且采用增量法的思想和線(xiàn)性強(qiáng)化彈塑性模型，編制出了適用于平面桁架彈塑性分析的MATLAB程序。在程序的每一個(gè)增量步的計(jì)算中，假定方程組是線(xiàn)性的，柔度矩陣也是常數(shù)，那么，每步施加一個(gè)荷載增量ΔP，得到一個(gè)位移增量Δδ，累計(jì)后即得到最終的位移。

2 算 例

本文算例采用的材料本構(gòu)是線(xiàn)性強(qiáng)化彈塑性模型，為研究方便，取各桿的楊氏模型為210GPa，屈服極限為250MPa，并服從Mises屈服條件，屈服后的彈性模量取為42MPa，各桿的橫截面積取單位10cm2。

2.1 三桿超靜定桁架

平面三桿超靜定桁架的單元尺寸和所受荷載如圖4所示。將應(yīng)用本文編制的MATLAB基面力元程序計(jì)算該算例，求得的單元軸力解和節(jié)點(diǎn)位移解與ABAQUS的計(jì)算結(jié)果進(jìn)行對(duì)比，并分別列于表1和表2。

圖4 三桿超靜定桁架的彈塑性分析Fig.4 Elastic-plastic analysis of three bar statically indeterminate truss

kN

表2 節(jié)點(diǎn)位移

從表1中可以看出，(1)單元和(2)單元均率先進(jìn)入了塑性階段，表2中也可以看出，A點(diǎn)產(chǎn)生了較大的塑性位移。本文對(duì)加載后的三桿桁架進(jìn)行卸載，并將卸載后A節(jié)點(diǎn)的位移值與ABAQUS的計(jì)算結(jié)果對(duì)比列于表3。

表3 卸載后的節(jié)點(diǎn)位移

最后給出了A點(diǎn)荷載與位移的關(guān)系，如圖5所示。

圖5 A點(diǎn)荷載與位移關(guān)系曲線(xiàn)Fig.5 Curves between displacement and load at node A

由表1—表3可知，應(yīng)用彈塑性基面力元程序計(jì)算的軸力和位移解與ABAQUS的計(jì)算結(jié)果高度吻合。荷載和位移的關(guān)系曲線(xiàn)圖也比較符合理論情況。

2.2 復(fù)雜靜定桁架

平面復(fù)雜靜定桁架的單元尺寸和所受荷載如圖6所示。將應(yīng)用本文編制的MATLAB基面力元程序計(jì)算該算例，求得的單元軸力解和節(jié)點(diǎn)位移解與ABAQUS的計(jì)算結(jié)果進(jìn)行對(duì)比，并分別列于表4和表5。

圖6 復(fù)雜靜定桁架的彈塑性分析Fig.6 Elastic-plastic analysis of complex statically determinate truss

kN

表5 節(jié)點(diǎn)位移

對(duì)加載后的復(fù)雜桁架算例進(jìn)行卸載，將卸載后各節(jié)點(diǎn)的位移值與ABAQUS的計(jì)算結(jié)果列于表6。

表6 卸載后的節(jié)點(diǎn)位移

最后給出了F點(diǎn)荷載與位移的關(guān)系，如圖7所示。

圖7 F點(diǎn)荷載與位移關(guān)系曲線(xiàn)Fig.7 Curves between displacement and load at node F

由表4—表6可知，應(yīng)用彈塑性基面力元程序計(jì)算的軸力和位移解與ABAQUS的計(jì)算結(jié)果高度吻合。荷載和位移的關(guān)系曲線(xiàn)圖也比較符合理論情況。

3 結(jié) 論

本文把增量法的思想應(yīng)用于平面桁架單元的余能原理基面力元模型中，編制了能夠?qū)ζ矫骅旒苓M(jìn)行彈塑性分析的MATLAB程序。與傳統(tǒng)的位移型有限元模型相比，該計(jì)算程序誤差較小，具有較高的精度。同時(shí)發(fā)現(xiàn)，在本文柔度矩陣、軸力以及位移的求解中均不用進(jìn)行局部坐標(biāo)系與整體坐標(biāo)系的轉(zhuǎn)換。表明以基面力作為基本未知量的桁架基面力元模型的簡(jiǎn)潔性和相應(yīng)程序的適用性，進(jìn)一步豐富和完善了余能原理在桁架中的應(yīng)用，同時(shí)也表明了基面力元理論的正確性和廣泛的應(yīng)用前景。

/References:

[1] CLOUGH R W. The finite element method in plane stress analysis[J]. Proceedings of American Society of Civil Engineers, 1960, 23: 345-378.

[2] WASHIZU K. Variational Methods in Elasticity and Plasticity[M]. 2nd ed. New York:Pergamon Press, 1975.

[3] VEUBEKE B F D. Displacement and equilibrium models in the finite element method[J]. International Journal for Numerical Methods in Engineering， 1965, 9(3):145-197.

[4] PIAN T H H,CHEN D P, KANG D. A new formulationof hybrid/mixed finite element[J]. Computers & Structures， 1983, 16(1/2/3/4): 81-87.

[5] WILSON E L, TAYLOR R L, DOHERTY W P, et al. Incompatible displacement models[J]. Numerical and Computer Methods in Structural Mechanics, 1973, 43-57.

[6] 龍馭球.新型有限元論[M].北京：清華大學(xué)出版社，2004.

[7] GAO Y C. A new description of the stress state at a point with applications[J].Archive of Applied Mechanics， 2003, 73(3): 171-183.

[8] 高玉臣.固體力學(xué)基礎(chǔ)[M].北京：中國(guó)鐵道出版社，1999.

[9] 高玉臣.彈性大變形的余能原理[J].中國(guó)科學(xué)G輯:物理學(xué)、力學(xué)、天文學(xué)， 2006，36(3)：298-311. GAO Yuchen. Complementary energy principle for large elastic deformation[J].Science in China(Series G. Physics, Mechanics & Astronomy, 2006, 36(3): 298-311.

[10]彭一江.基于基面力概念的新型有限元方法[D].北京:北京交通大學(xué),2006. PENG Yijiang. A New Finite Element Method Based on the Concept of Base Force[D]. Beijing: Beijing Jiaotong of Technology, 2006.

[11]PENG Yijiang, LIU Yinghua. Base force element method of complementary energy principle for large rotation problems[J]. Acta Mechanica Sinica, 2009, 25(4): 507-515.

[12]張昭鋒.基于余能原理的基面力元法在巖體力學(xué)中的應(yīng)用[D]. 北京：北京工業(yè)大學(xué),2014. ZHANG Zhaofeng. The Application of Base Force Element Method on Complementary Energy Principle in Rock Mechanics[D]. Beijing:Beijing University of Technology, 2014.

[13]王耀.三維基面力元法及其在再生混凝土細(xì)觀(guān)損傷分析中的應(yīng)用[D].北京：北京工業(yè)大學(xué),2015. WANG Yao. Research of Three-Dimensional Base Force Element Method and Application on Meso-Demage Analysis for Recycled Concrete[D]. Beijing: Beijing University of Technology, 2015.

[14]SANTOS H A F A, PIMENTA P M, ALMEIDA J P M. A hybrid-mixed finite element formulation for the geometrically exact analysis of three-dimensional framed structures[J]. Computational Mechanics,2011, 48: 591-613.

[15]SANTOS H A F A. Variationally consistent force-based finite element for the geometrically non-linear analysis of Euler-Bernoulli framed structures[J]. Finite Elements in Analysis and Design, 2012, 53: 24-36.

[16]SANTOS H A F A. A novel updated Lagrangian complementary energy-based formulation for the elastica problem: force-based finite element model[J]. Acta Mechanica,2015, 226(4): 1133-1151.

[17]彭一江,雷文賢,彭紅濤.基于基線(xiàn)力概念的平面4節(jié)點(diǎn)余能有限元模型[J].北京工業(yè)大學(xué)學(xué)報(bào),2009,35(11)：1484-1490. PENG Yijiang, LEI Wenxian, PENG Hongtao. Four-node plane element model on complementary energy based on the concept of base line forces[J]. Journal of Beijing University of Technology, 2009,35(11)：1484-1490.

Application of base force element method with complementary energy principle for elasto-plastic truss analysis

LI Ruixue, PENG Yijiang

(College of Architecture and Civil Engineering, Beijing University of Technology, Beijing 100124, China)

Base force element method(BFEM) is a new type of finite element method, which uses the "base force" as fundamental variables to describe the stress state at a point. In order to compile a base force element program for the elasto-plastic analysis of truss using MATLAB language, based on the BFEM with complementary energy principle, the solid element of the plane 4 nodes in BFEM is degenerated into the truss element. The concrete expressions of the flexibility matrix, axial force and displacement of the degenerated truss elements are derived, and necessary modification of the governing equations is conducted. According to the degenerated truss element model, the idea of incremental method and linear strengthened elasto-plastic model are adopted for calculation by using 3-pole indeterminate truss and complicated static determinacy as samples. It is found that through comparison of the calculated results and the results from T2D2 of ABAQUS, the calculated results using the method are in highly agreement with those obtained with ABAQUS, which further verifies the correctness of the base force element method and the applicability of the program.

plastic mechanics; finite element method; base force element method(BFEM); complementary energy principle; truss; elasto-plastic analysis

1008-1534(2017)04-0276-06

2017-04-07;

2017-04-27；責(zé)任編輯：陳書(shū)欣

國(guó)家自然科學(xué)基金(10972015)

李瑞雪(1991—)，男，河北邯鄲人，碩士研究生，主要從事余能原理基面力元法方面的研究。

彭一江教授。E-mail：pengyijiang@bjut.edu.cn

O344

A

10.7535/hbgykj.2017yx04008

李瑞雪，彭一江.余能原理基面力元法在桁架彈塑性分析中的應(yīng)用[J].河北工業(yè)科技,2017,34(4):276-281. LI Ruixue, PENG Yijiang.Application of base force element method with complementary energy principle for elasto-plastic truss analysis[J].Hebei Journal of Industrial Science and Technology,2017,34(4):277-281.