基于共旋坐標法的空間桁架幾何非線性分析

朱麗君, 顏海泉, 董正方*

(1.河南大學土木建筑學院, 開封 475004; 2.上海市政工程設(shè)計研究總院(集團)有限公司, 上海 200082)

傳統(tǒng)的幾何非線性分析方法完全拉格朗日列式法(total Lagrange formulation, TL)與更新拉格朗日列式法(updated Lagrange formulation, UL)均未考慮變形后的坐標系的變化[1],單元發(fā)生較大剛體轉(zhuǎn)動時的計算誤差較大。共旋坐標法基于UL法,通過建立隨動坐標系[2-4]來解決這一問題,單元剛體轉(zhuǎn)角任意大時計算都是精確的;共旋坐標法迭代過程中,整體坐標系與局部坐標系的轉(zhuǎn)換矩陣不斷更新,幾何非線性通過坐標轉(zhuǎn)換計算而被考慮進來;目前諸多工程實例中均需考慮幾何非線性與材料非線性的耦合作用[5-6],相較于TL法和UL法,共旋坐標法可使耦合非線性問題的剛度矩陣脫耦、計算更加簡潔[7],在大位移小應(yīng)變幾何非線性分析中更具有優(yōu)勢。

共旋坐標法[8]的基本原理是將單元剛體位移扣除,從而得到單元純變形與桿端力的關(guān)系,因其高效精確的優(yōu)點而被廣泛應(yīng)用于梁桿單元的幾何非線性分析中。常用的共旋坐標法有隨動坐標系法[7,9]與隨轉(zhuǎn)坐標系法[10]兩種:前者提出以單元i節(jié)點為原點建立隨動坐標系,將梁桿的大位移分解為隨動坐標的剛性位移與彈性位移;后者將隨動坐標系的原點設(shè)在初始局部坐標系原點處,此坐標系隨單元轉(zhuǎn)動而不平動,進而得到幾何非線性與線性平衡方程的統(tǒng)一表達式。蔡松柏等[11]、鄧繼華等[12]、馮曉東等[13]利用桁架單元的幾何與物理關(guān)系,通過變分法直接推導出桁架單元的剛度矩陣,進一步簡化了桁架幾何非線性的計算,并編制基于面向過程技術(shù)的程序或?qū)ζ涠伍_發(fā)。馬玉全[14]對傳統(tǒng)共旋坐標法中的單元長度、荷載及抗力進行改進,使其更好地適用于梁桿單元,提高了計算精度,并編制相關(guān)程序。李東升等[15]使用逆向運動法將梁單元的剛體旋轉(zhuǎn)分解,從而計算出精確的旋轉(zhuǎn)矩陣,并編制了相應(yīng)有限元程序。

但共旋坐標法的隨動與整體坐標系的轉(zhuǎn)換矩陣需以整體與初始局部坐標系的轉(zhuǎn)換矩陣為中介矩陣[9-10,14,16],若能進一步簡化轉(zhuǎn)換矩陣的計算,將會減少軟件中矩陣的存儲步驟,尤其對于大型算例,計算效率將顯著提高。且共旋坐標法在程序中的實現(xiàn)諸多基于面向過程技術(shù)[11,14,17],然而面向?qū)ο蠹夹g(shù)具有更高的數(shù)據(jù)抽象和封裝能力,更易實現(xiàn)程序功能的擴充,也是有限元程序設(shè)計的主流工具[18]。此外,桁架幾何非線性的判斷條件還缺乏研究,經(jīng)驗理論認為應(yīng)變達到5%時考慮幾何非線性,但此理論在桁架中的適用性有待驗證。

現(xiàn)通過取整體坐標系原點為隨動坐標系的原點,得到更簡單高效的轉(zhuǎn)換矩陣計算方法,使用面向?qū)ο驝++語言在程序Aduts中實現(xiàn)此法桁架單元。在程序中進行算例分析并與理論解對比來驗證其正確性,并對其進行參數(shù)分析,總結(jié)桁架幾何非線性的判斷依據(jù)。

1 桁架單元共旋坐標法理論

常用桁架單元的共旋坐標法主要有以單元端點、初始局部坐標系原點為隨動坐標系的原點兩種[7-8],且多基于平面問題。

1.1 隨動坐標系法

此理論需用到三個坐標系:整體坐標系、初始局部坐標系、隨動坐標系,隨動坐標系的原點為單元的i節(jié)點,此坐標系隨單元平動和轉(zhuǎn)動,以平面問題為例。

首先需計算整體坐標系與初始局部坐標系的轉(zhuǎn)換矩陣R0,計算式為

(1)

式(1)中:lxX、lxY為局部坐標系的x軸對整體坐標系X、Y軸的方向余弦;lyX、lyY為局部坐標系的y軸對整體坐標系X、Y軸的方向余弦。

再計算初始時刻局部坐標系與t時刻局部坐標系的轉(zhuǎn)換矩陣Rt,計算式為

(2)

式(2)中:α為初始局部坐標系x軸與t時刻局部坐標系的夾角(圖1),計算式為

(3)

最后計算整體坐標系與t時刻局部坐標系的轉(zhuǎn)換矩陣R=Rt×R0。

XOY為整體坐標系;xoy為初始局部坐標系;xtotyt為隨動坐標系; i0、j0表示初始時刻單元的兩個端點位置;it、jt表示t時刻后單元的 兩個端點位置圖1 單元位形及隨動坐標系Fig.1 Configuration of element and concomitant coordinate system

1.2 隨轉(zhuǎn)坐標系法

如圖2所示,XOY為整體坐標系,xoy為初始局部坐標系,xtotyt為隨動坐標系,以平面問題為例,此理論與隨動坐標系法的計算步驟相同,只是初始局部坐標系與隨動坐標系的原點重合,隨動坐標系相對于初始局部坐標系只轉(zhuǎn)動,但各轉(zhuǎn)換矩陣的計算列式與1.1節(jié)中一致。

由此可見,這兩種建立隨動坐標系的方法都會使隨動與整體坐標的轉(zhuǎn)換矩陣R的計算分成三步:①計算整體與初始局部坐標系的轉(zhuǎn)換矩陣R0;②計算隨動與初始局部坐標系的轉(zhuǎn)換矩陣Rt;③整體與隨動坐標系的轉(zhuǎn)換矩陣R=Rt×R0。

桁架共旋坐標法的核心為轉(zhuǎn)換矩陣的更新,而轉(zhuǎn)換矩陣在每個迭代步都需根據(jù)位移值來同步更新,若是能對轉(zhuǎn)換矩陣的計算做進一步的簡化,幾何非線性的計算將更加高效,尤其是對于大型工程,將會有效節(jié)省計算時間。

圖2 單元位形及隨轉(zhuǎn)坐標系Fig.2 Configuration of element and co-rotational coordinate system

1.3 基于整體坐標系的隨轉(zhuǎn)坐標系法

在對空間桁架進行大變形幾何非線性分析時,將隨動坐標系的原點設(shè)在整體坐標系原點處,此隨動坐標系只隨單元發(fā)生轉(zhuǎn)動而不平動。

如圖3所示,將t時刻的單元平移,使單元的i節(jié)點與整體坐標系的原點重合,單元兩節(jié)點的連線方向為隨動坐標系xt軸方向,且假設(shè)整體坐標系的Z軸始終在隨動坐標系的xtotzt面內(nèi),xt軸的方向向量由單元與整體坐標系三個坐標軸的夾角即可計算出,yt軸的方向向量可由xt的方向向量與xtotzt面

圖3 空間桁架坐標系與位形Fig.3 Coordinate system and configuration of space truss

內(nèi)且與xt軸不重合的向量(0,0,1)叉乘求得,zt軸的方向向量可由yt軸與xt軸的方向向量叉乘得到。

圖3中,α、β、γ分別為變形后的單元與整體坐標系X、Y、Z軸的夾角,也即隨動坐標系x軸與整體坐標系夾角,平面P1與xt軸垂直,平面P2與整體坐標系Z軸垂直,平面P1與P2的交線即為隨動坐標系yt軸的方向,zt軸方向由右手定則確定。

現(xiàn)做以下假設(shè):(Xi,Yi,Zi)、(Xj,Yj,Zj)分別為變形前單元i、j端的坐標;(Xi+Ui,Yi+Vi,Zi+Wi)、(Xj+Uj,Yj+Vj,Zj+Wj)分別為變形后單元i、j端的坐標;Uλ,Vλ,Wλ(λ=i,j)為單元i、j端沿X、Y、Z軸方向的位移值;L0、Lt分別為變形前后的單元長度。則xt軸的三個方向余弦值為

cosα=(Xj-Xi+Uj-Ui)/Lt

(4)

cosβ=(Yj-Yi+Vj-Vi)/Lt

(5)

cosγ=(Zj-Zi+Wj-Wi)/Lt

(6)

(7)

Lt=[(Uj-Ui+Xj-Xi)2+

(Vj-Vi+Yj-Yi)2+

(8)

則隨動坐標系xt軸的方向向量為(cosα,cosβ,cosγ),yt軸的方向向量為(-cosβ,cosα,0),zt軸的方向向量為(-cosαcosγ,-cosβcosγ,cos2α+cos2β)由yt、xt兩坐標軸的方向向量叉乘得來。

則可得到隨動坐標系與整體坐標系的轉(zhuǎn)換矩陣R[19]。即

(9)

由局部坐標系下桿端軸力計算式(10)及式(4)～式(6),可得軸力在整體坐標系六個方向上的分力,再經(jīng)變分運算可得單元剛度矩陣K=KL+KN,其中KL為剛度矩陣的線性部分,KN為剛度矩陣的非線性部分,也稱幾何剛度矩陣。表達式為

(10)

式(10)中:E為桿件的彈性模量;A為其截面面積。

[KL]6×6=RTTlbTK3×3TlbR

(11)

(12)

(13)

(14)

(15)

整體坐標系下剛度矩陣的具體表達式為

(16)

(17)

式中:Δ為桿件單元長度變化,Δ=Lt-L0; symmetry表示關(guān)于對角線對稱的矩陣的另一部分。

(18)

將隨動坐標系原點取為整體坐標系原點的優(yōu)點是可直接得出轉(zhuǎn)換矩陣,不需整體坐標系與初始局部坐標系的轉(zhuǎn)換矩陣做中介,減少了計算步驟,提高了計算效率。

2 共旋坐標法桁架單元的程序?qū)崿F(xiàn)

2.1 單元接口基類

以桁架的共旋坐標法為原理,使用面向?qū)ο驝++語言在程序Aduts中編寫了桁架單元。面向?qū)ο笳Z言由于具有動態(tài)內(nèi)存分配的功能,程序的數(shù)據(jù)組織相互獨立,修改某一處的數(shù)據(jù)不會對其他部分數(shù)據(jù)造成影響,使程序更易擴充;且面向?qū)ο蠹夹g(shù)中的數(shù)據(jù)類型可自行定義,不會像面向過程技術(shù)一樣受數(shù)據(jù)類型的限制,故更易編制類型豐富的單元及本構(gòu)材料;面向?qū)ο蟪绦蛟O(shè)計方法將程序框架劃分為不同的類,每個類只進行抽象描述,具體的定義在各子類中實現(xiàn),使用時不需知道子類的實現(xiàn)細節(jié),具有很好的數(shù)據(jù)封裝與代碼可重用性,易于程序開發(fā)。

本文程序中抽象出了單元接口基類ElementInterface.h,開發(fā)人員通過重寫其內(nèi)部的虛函數(shù)來添加所需要的單元子類。接口定義了繼承該接口的子類需要實現(xiàn)的方法,其虛函數(shù)如表1所示。

本程序的單元接口基類Topology2()函數(shù)中設(shè)置了附加參數(shù),各單元子類可根據(jù)功能需求自行添加與修改附加參數(shù),而不會對其他單元子類產(chǎn)生影響,這也是本程序的一大優(yōu)點,附加參數(shù)的設(shè)置簡化了單元子類的開發(fā)、維護了程序的封裝性與多態(tài)性。

表1 ElementInterface類中的虛函數(shù)Table 1 Virtual functions of ElementInterface

2.2 桁架單元的擴展

通過繼承單元接口基類的虛函數(shù)實現(xiàn)桁架單元的擴展,關(guān)鍵步驟為每個迭代步對單元轉(zhuǎn)換矩陣、剛度矩陣、抗力向量的更新,這就需處理好程序框架對單元子類中各函數(shù)的調(diào)用工作。故Aduts程序中桁架單元的編寫工作分為兩部分,一是對單元接口基類中函數(shù)的具體定義,二是在不修改程序框架的前提下,完成對單元子類中各函數(shù)的調(diào)用,從而實現(xiàn)各數(shù)據(jù)的更新與計算流程。

在重寫基類函數(shù)Topology1()時,因所加桁架單元是考慮幾何非線性的,故幾何非線性標簽改為true。重寫函數(shù)Topology2()時,附加參數(shù)有四個:“l(fā)ength0”“Initial”“l(fā)engthT”與“CorotTag”,分別表示桁架單元的初始單元長度、調(diào)用轉(zhuǎn)換矩陣標簽(用于判斷是進行迭代前轉(zhuǎn)換矩陣的調(diào)用,還是進入迭代步后轉(zhuǎn)換矩陣的調(diào)用)、t時刻的單元長度與共旋坐標法的標簽(用于判斷是否使用共旋坐標法),比普通單元(不考慮幾何非線性的單元)多了后三個附加參數(shù)。

重寫轉(zhuǎn)換矩陣ComputeTransformMatrixData()時,與普通單元不同的是,共旋坐標法桁架單元需隨著單元位移的更新而更新局部坐標系與整體坐標系的轉(zhuǎn)換矩陣,普通單元只需根據(jù)初始位形計算一次轉(zhuǎn)換矩陣。

重寫剛度矩陣ComputeElementStiffMatrix()時,由于幾何非線性單元的剛度矩陣分為線性剛度與非線性幾何剛度兩部分,故比普通單元多了幾何剛度矩陣的定義部分。

其余單元接口基類函數(shù)的重寫與普通單元基本相同。

實現(xiàn)單元子類中各函數(shù)的準確調(diào)用是單元擴展的難點之一,本文中對共旋坐標法桁架單元計算流程的處理如圖4所示。

圖4 桁架的計算流程Fig.4 Calculational flow of truss

圖4中計算流程與普通單元的區(qū)別在于轉(zhuǎn)換矩陣的更新,這需通過在單元子類內(nèi)部定義附加參數(shù)來實現(xiàn)。不使用共旋坐標法的單元只在迭代之前計算一次整體坐標與局部坐標系的轉(zhuǎn)換矩陣,使用此法時則在每個迭代步都需更新轉(zhuǎn)換矩陣。由此可見,通過在單元子類中設(shè)置不同附加參數(shù),滿足了單元子類特殊的功能需求,對單元子類的擴展有重要作用。

3 算例

3.1 單元正確性的驗證

如圖5所示桁架,各桿件的彈性模量與橫截面積大小一致,圖5中桁架單元桁架跨長a=100 cm,b=-2 cm,彈性模量E=200 GPa,橫截面積A=4 cm2,桁架頂端作用豎直向下的集中荷載,劃分10個荷載步,位移收斂條件為|ΔU/U|≤1×10-3。

本文程序的拱頂位移計算結(jié)果與理論結(jié)果[20]的對比如圖6所示,由圖6可知,本文程序解與理論值吻合度很高,驗證了本文程序中桁架單元的正確性。

圖5 計算結(jié)構(gòu)模型Fig.5 Structural model of calculation

圖6 拱頂豎向位移對比Fig.6 Vertical displacement comparison of vault

3.2 參數(shù)分析

由應(yīng)變與軸力的公式可知,應(yīng)變大小與拱跨比b/a、P/EA的值緊密相關(guān),無論單元長度、截面面積如何變化,只要b/a與P/EA的值不變,應(yīng)變值不會改變。

因此,對桁架進行幾何非線性分析時,通過改變圖5中桁架的拱跨比b/a與P/EA的值,使應(yīng)變保持在指定值,得到圖7中幾何非線性與線性位移誤差與拱跨比關(guān)系圖。位移相對誤差為線性與幾何非線性位移之差與非線性位移的比值,其中,幾何非線性計算的第一次結(jié)果即為線性結(jié)果,此參數(shù)分析得出的結(jié)論對任意長桁架均適用。

位移誤差是線性計算未考慮單元的位形變化、將剛體轉(zhuǎn)動計入到線性應(yīng)變中,使得線性位移的計算結(jié)果出錯而引起的。

如圖7所示,當應(yīng)變分別為1%、2%、3%、4%、5%、8%,拱跨比分別為0.073、0.103、0.125、0.144、0.161、0.201時,位移相對誤差等于工程精度誤差要求5%;各應(yīng)變值對應(yīng)的拱跨比分別為0.081、0.115、0.140、0.161、0.179、0.220時幾何非線性與線性的相對誤差為0;各應(yīng)變值對應(yīng)的拱跨比分別為0.093、0.130、0.159、0.182、0.202、0.251時幾何非線性與線性的相對誤差為-5%。

由圖7可知,應(yīng)變不變,拱跨比越小位移的相對誤差越大,故即使應(yīng)變未超過5%,位移誤差可能會超過工程精度誤差要求,也需考慮幾何非線性。因此,判斷桁架是否考慮幾何非線性時,應(yīng)將桁架的應(yīng)變與拱跨比的具體數(shù)值結(jié)合來討論,不能單純以應(yīng)變的大小作為判斷標準。

圖7 幾何非線性與線性位移的相對誤差Fig.7 Relative error of displacement of geometric non-linearity and linearity

4 結(jié)論

本文給出一種更簡單的空間桁架共旋坐標法的轉(zhuǎn)換矩陣及抗力向量,并用數(shù)值算例分析了桁架的幾何非線性,得到如下結(jié)論。

(1)基于整體坐標系建立共旋坐標法的隨動坐標系,能夠直接給出轉(zhuǎn)換矩陣表達式,此法轉(zhuǎn)換矩陣的列式更簡單、計算步驟更少,適用于任意大小剛體轉(zhuǎn)角的梁桿單元。

(2)本文中的共旋坐標法比以往方法少了一步中介轉(zhuǎn)換矩陣的計算,在有限元程序中實現(xiàn)時會節(jié)省數(shù)據(jù)的存儲空間、提高程序計算效率,對于大型工程結(jié)構(gòu)更加明顯。

(3)面向?qū)ο蠹夹g(shù)實現(xiàn)的共旋桁架單元比面向過程技術(shù)具有更好的數(shù)據(jù)管理、更易開發(fā)與維護的特點,且單元子類中的附加參數(shù)使此程序有更好的封裝性與多態(tài)性。

(4)桁架不能單純以應(yīng)變達到5%為幾何非線性的判斷依據(jù),需將應(yīng)變與拱跨比結(jié)合來判斷是否考慮幾何非線性。