應(yīng)用ANSYS用戶可編程特性實現(xiàn)梁單元開發(fā)研究*

黃 斌 鄧洪洲

(同濟(jì)大學(xué)建筑工程系 上海 200092)

應(yīng)用ANSYS用戶可編程特性實現(xiàn)梁單元開發(fā)研究*

黃 斌 鄧洪洲

(同濟(jì)大學(xué)建筑工程系 上海 200092)

應(yīng)用ANSYS平臺提供的UPFs二次開發(fā)工具，對考慮幾何非線性梁元的開發(fā)過程進(jìn)行探索.對于用到的子程序文件uec102，uel102，以及uep102進(jìn)行了研究，編寫代碼實現(xiàn)了梁單元開發(fā)，并通過經(jīng)典算例校驗了程序的可靠性.同時，在一般梁元開發(fā)研究基礎(chǔ)上，進(jìn)一步創(chuàng)建了能夠考慮節(jié)點剛度影響的新型梁元，應(yīng)用于一單層球面網(wǎng)殼結(jié)構(gòu)中，與普通梁元加彈簧模擬方法進(jìn)行了比較，兩者計算結(jié)果吻合較好，證明了該方法的可行性.

ANSYS二次開發(fā)；幾何非線性；梁單元開發(fā)；節(jié)點剛度

0 引 言

對于網(wǎng)殼結(jié)構(gòu)中構(gòu)件通常按照兩端鉸接的軸心受力模型進(jìn)行驗算.但實際結(jié)構(gòu)節(jié)點并不完全鉸接，存在一定的抗彎剛度.在用ANSYS軟件進(jìn)行有限元分析時，對于考慮節(jié)點剛度影響的構(gòu)件常使用梁單元加彈簧方式進(jìn)行模擬，但就目前ANSYS單元庫中的彈簧元而言，采用梁元兩端連接彈簧，每個節(jié)點至少要連接2個彈簧，其中1個用于控制平動位移，另1個用于控制轉(zhuǎn)角位移，如此，在進(jìn)行建模時，每個節(jié)點要額外增加2個節(jié)點用于構(gòu)建彈簧元，這樣節(jié)點數(shù)量較多，且出現(xiàn)節(jié)點位置重合，增加建模難度.

文中應(yīng)用ANSYS平臺提供的用戶可編程特性UPFs(user programmable features)從一般梁元創(chuàng)建出發(fā)，探索梁元開發(fā)過程.對于考慮節(jié)點影響構(gòu)件，可將構(gòu)件與半剛性節(jié)點作為一體，研究其等效單元剛度矩陣，作為一新型梁元單元剛度矩陣，在建模時不用再增加節(jié)點模擬彈簧元，仍可按梁元看待，降低了建模難度，提高了建模效率.

UPFs是對ANSYS程序提供的可供用戶修改子程序文件(uex100～uex105)進(jìn)行調(diào)整或修改，用于建立新型材料的本構(gòu)關(guān)系，定義復(fù)合材料的失效準(zhǔn)則，或根據(jù)用戶要求開發(fā)具有特定力學(xué)性能的單元等[1].

1 梁單元的開發(fā)過程探索

1.1 單元形函數(shù)

文中一般梁單元為3節(jié)點線型單元，其中僅節(jié)點i，j對單元剛度有貢獻(xiàn)，k節(jié)點為輔助節(jié)點，用于截面朝向的控制.

梁元形函數(shù)采用2節(jié)點(除去輔助節(jié)點k)3次hermite插值函數(shù)，單元場函數(shù)在節(jié)點處具有一階導(dǎo)數(shù)連續(xù)性，為C1型單元，形函數(shù)如下所示[2].

采用量綱一的量坐標(biāo)

1.2 幾何非線性分析

1.2.1 空間大轉(zhuǎn)動處理

文獻(xiàn)[3-5]給出了考慮幾何非線性時空間大轉(zhuǎn)動的處理方法.三維空間旋轉(zhuǎn)變換除了指定旋轉(zhuǎn)角外，還需指定旋轉(zhuǎn)軸.

給定具有單位長度的旋轉(zhuǎn)軸A=(ax,ay,az)和旋轉(zhuǎn)角θ，則繞OA軸旋轉(zhuǎn)變換的旋轉(zhuǎn)矩陣Tθ表示為

(1)

式中：ax，ay,az為旋轉(zhuǎn)軸的方向余弦.

1.2.2 單元坐標(biāo)轉(zhuǎn)換矩陣

坐標(biāo)轉(zhuǎn)換矩陣取決于單元坐標(biāo)系與整體坐標(biāo)系之間的相對位置關(guān)系.考慮幾何非線性時，單元兩端截面的位置不斷變化見圖1，因此需要更新每次迭代的坐標(biāo)轉(zhuǎn)換矩陣.

圖1 單元坐標(biāo)系示意圖

梁單元坐標(biāo)系與整體坐標(biāo)系的轉(zhuǎn)換矩陣為

式中：r為單元定向矩陣，可表示為

(2)

式中:rij(i，j=1,2,3)為單元局部坐標(biāo)軸j(j=x,y,z)與整體坐標(biāo)軸i(i=X,Y,Z)的夾角的方向余弦，r的第一列為單元局部坐標(biāo)x軸與整體坐標(biāo)系中各坐標(biāo)軸的方向余弦，由現(xiàn)時構(gòu)形下單元兩端節(jié)點的位置即可得到，此時單元的長度為

單元矩陣的y軸與z軸在整體坐標(biāo)系中相對位置的確定比較復(fù)雜.設(shè)節(jié)點坐標(biāo)系的任意時刻的坐標(biāo)矩陣為

式中：α為節(jié)點坐標(biāo)系的坐標(biāo)軸任意時刻與整體坐標(biāo)軸夾角的方向余弦矩陣.

初始構(gòu)形下節(jié)點坐標(biāo)系分別平行于3根整體坐標(biāo)系的坐標(biāo)軸，則初始節(jié)點坐標(biāo)系矩陣為

假設(shè)在t時刻αt已經(jīng)求得t+Δt時刻的轉(zhuǎn)角增量為Δθ(Δθx，Δθy，Δθz)，則根據(jù)式(1)可以得到t+Δt時刻的節(jié)點定向矩陣αt+Δt.

αt+Δt=TΔθαt

(3)

式中：αt+Δt為節(jié)點坐標(biāo)系相對于整體坐標(biāo)系的位置.定義端截面定向矩陣S.S的第一列為截面法線的方向余弦，其余兩列為截面主軸的方向余弦.則有

Sit+Δt=αit+Δtr0

Sjt+Δt=αjt+Δtr0

(4)

式中：r0為初始單元坐標(biāo)的定向矩陣.初始時刻，截面坐標(biāo)系平行于單元坐標(biāo)系，也即為截面坐標(biāo)系的矩陣，r0的確定可參閱文獻(xiàn)[6].

截面坐標(biāo)系S確定之后，在上文單元坐標(biāo)系的x軸的方向已經(jīng)確定，因此，便可獲得單元坐標(biāo)系的x軸與截面坐標(biāo)系的x軸之間的夾角，假定單元的相對變形是微小的，如下式.

(5)

為求r的第二、三列，令：

由此可根據(jù)兩端截面的定向矩陣得到兩單元坐標(biāo)系的定向矩陣

Pit+Δt=Sit+Δteit+Δt

Pjt+Δt=Sjt+Δtejt+Δt

(6)

取單元截面主軸方向為兩端截面主軸方向的平均值，則可得

(7)

1.3 UPFs子程序文件說明

主要用到的子程序文件有：uec102， uel102，uep102.

uec102文件用于描述所創(chuàng)建單元的基本特征，該子程序文件中主要是對一維數(shù)組ielc()進(jìn)行賦值.uel102子程序文件是單元計算的核心程序.在該文件中單元剛度矩陣計算，內(nèi)力計算，單元坐標(biāo)轉(zhuǎn)換矩陣的更新，非線性計算時Newton-Raphson恢復(fù)力的計算，非線性分析迭代的邏輯都需要在uel102中編寫相應(yīng)代碼.uep102主要是用于控制線單元輸出結(jié)果，用戶可在此子程序中指定輸出內(nèi)容，如單元桿端力，單元應(yīng)力等.

2 經(jīng)典算例檢驗

為檢驗程序的可靠性及有效性，與經(jīng)典算例進(jìn)行了對比分析.

2.1 William平面剛架幾何非線性分析

William平面剛架是結(jié)構(gòu)幾何非線性分析中的經(jīng)典算例，其結(jié)構(gòu)布置簡圖見圖2.文中對2種不同高度(h=9.804，8.128 mm)的框架進(jìn)行了計算分析.

圖2 William框架結(jié)構(gòu)簡圖(尺寸單位：mm)

圖3～4分別對應(yīng)h=9.804 mm時點A處荷載P-豎向位移δA曲線及荷載P-水平支座反力H曲線.

圖3 點A處P-δA曲線

圖4 P-H曲線

圖5～6分別對應(yīng)h=8.128 mm時點A處荷載p-豎向位移δA曲線及荷載p-水平支座反力H曲線.單元劃分?jǐn)?shù)量均為4單元.

圖5 點A處P-δA曲線

圖6 P-H曲線

2.2 六角星形穹頂平衡路徑追蹤

六角星形穹頂結(jié)構(gòu)平衡路徑分析也是幾何非線性分析的經(jīng)典算例，結(jié)構(gòu)布置見圖7.

圖7 六角形猩猩穹頂結(jié)構(gòu)布置簡圖(尺寸單位：cm)

文中在進(jìn)行分析時，考慮到文中使用的形函數(shù)與ANSYS單元庫中BEAM4的形函數(shù)相同，因此與使用了ANSYS的BEAM4梁元模擬的結(jié)果進(jìn)行了對比.計算結(jié)果見圖9～12.

圖8～9為結(jié)構(gòu)點1處的荷載位移曲線.圖8單元劃分?jǐn)?shù)量均為1，圖9單元劃分?jǐn)?shù)量為4.圖10為4單元時，荷載與點2處的x向水平位移ux關(guān)系曲線，圖11為4單元時點2處的z向位移uz關(guān)系曲線.

圖8 點1處P-δ1曲線

圖9 點1處P-δ1曲線

圖10 點2處P-ux曲線

圖11 點2處P-uz曲線

通過比較可知：2種方法的計算所得的荷載位移曲線基本重合，說明文中使用UPFs編寫的程序在邏輯上是合理的.

2.3 六邊形12構(gòu)件空間框架幾何非線性分析

文中選取了文獻(xiàn)[7]中的六邊形框架(見圖12)算例進(jìn)行了計算分析.，并與ANSYS單元庫中的BEAM4單元進(jìn)行對比.

圖12 六邊形空間框架(尺寸單位：mm)

計算結(jié)果見圖13～14，圖13為單元劃分?jǐn)?shù)量均為1單元，圖14為單元劃分?jǐn)?shù)量均為4單元.兩者的荷載位移曲線均較好的吻合.

圖13 點A處P-δA曲線

圖14 點A處P-δA曲線

3 新型梁元的網(wǎng)殼分析應(yīng)用

對于網(wǎng)殼結(jié)構(gòu)，許多學(xué)者為了考慮節(jié)點對結(jié)構(gòu)計算結(jié)果的影響，提出了多種計算模型.在一般梁元開發(fā)研究的基礎(chǔ)上，文中基于剛臂-彈簧模型，研究了考慮節(jié)點尺寸以及節(jié)點半剛性影響的構(gòu)件等效剛度矩陣，并作為應(yīng)用ANSYS的UPFs開發(fā)的新型梁元單元剛度矩陣，應(yīng)用于一單層球面網(wǎng)殼結(jié)構(gòu)計算分析中.

當(dāng)不考慮節(jié)點尺寸的影響，僅考慮節(jié)點半剛性時，計算簡圖見圖15.

圖15 考慮節(jié)點半剛性單元計算簡圖

彈簧的力矩-轉(zhuǎn)角關(guān)系可以表示為

Ms=Ksθs

(8)

Ms=[MyaMzaMyaiMzaiMyb

MzbMybiMzbi]T

(9)

θs=[θyaθzaθyaiθzaiθybθzbθybiθzbi]T

(10)

式中：θyai，θzai，θybi，θzbi為增加的內(nèi)部節(jié)點自由度；Myai，Mzai，Mybi，Mzbi為與其對應(yīng)的力矩.

式中：Rki為彈簧的轉(zhuǎn)動剛度，與圓管相連的焊接球的轉(zhuǎn)動剛度可根據(jù)文獻(xiàn)[8-9]提供的公式計算.

構(gòu)件采用普通梁元模擬每個單元有12個自由度，考慮彈簧之后，增加了4個內(nèi)部自由度，此時，構(gòu)件的自由度變?yōu)?6個.將原構(gòu)件的12×12剛度矩陣與彈簧的8×8剛度矩陣均擴(kuò)展為16×16階矩陣，并且進(jìn)行合并，則有：

fa=Kaua

(12)

式中：

fa=[FxaFyaFzaMxaMyaMza

FxbFybFzbMxbMybMzb

MyaiMzaiMybiMzbi]T

(13)

ua=[uavawaθxaθyaθzaubvbwb

θxbθybθzbθyaiθzaiθybiθzbi]T

(14)

式(12)可寫為

(15)

通過凝聚內(nèi)部自由度，將U4消掉，由式(15)可得

(16)

將式(16)代入式(15)可得

(17)

式(17)中

(18)

(19)

(20)

進(jìn)一步再考慮焊接球節(jié)點的尺寸影響，采用剛臂彈簧梁元模型，計算見圖16.

圖16 剛臂彈簧模型的計算簡圖

根據(jù)向量運算法則，a點的位移向量和荷載向量與A點的位移向量和荷載向量有如下關(guān)系

(21)

(22)

(23)

(24)

(25)

Kab與KAB之間的關(guān)系為

(26)

式(26)中，轉(zhuǎn)換矩陣的表達(dá)式為

(27)

考慮焊接球尺寸影響時，可取exA為焊接球的半徑長度，eyA=0，ezA=0.

文獻(xiàn)[8]采用梁元連接彈簧元再加剛臂模擬網(wǎng)殼中的構(gòu)件，考慮焊接球影響，對一單層球面網(wǎng)殼結(jié)構(gòu)進(jìn)行了分析計算.網(wǎng)殼結(jié)構(gòu)布置圖見圖17～18.結(jié)構(gòu)跨度30 m，矢跨比為1/4，桿件外徑為0.077 5 m、內(nèi)徑為0.069 5 m，焊接球的直徑為0.45 m、厚度為0.014 m.

圖17 單層球面網(wǎng)殼結(jié)構(gòu)布置圖

圖18 單層球面網(wǎng)殼1/8結(jié)構(gòu)

文中將基于剛臂彈簧模型的等效剛度矩陣，作為一新型梁元的單元剛度矩陣，應(yīng)用在單層球面網(wǎng)殼結(jié)構(gòu)分析中，并分別計算對比了文獻(xiàn)[8]中的2種工況下構(gòu)件應(yīng)力.工況1為網(wǎng)殼各節(jié)點均受豎向力10 kN，工況2為僅網(wǎng)殼的頂點受豎向力10 kN.

計算結(jié)果見表1～2，表中的“組合模型”是指文獻(xiàn)[8]中構(gòu)件采用梁元加彈簧元再加剛臂方法建立的模型，由表1～2可知，2種模型計算結(jié)果的軸應(yīng)力比較接近，說明使用新型梁元模擬考慮節(jié)點剛度的構(gòu)件的力學(xué)模型是可行的.

表1 工況1時的構(gòu)件應(yīng)力對比 ×108 N/m2

表2 工況2下的構(gòu)件應(yīng)力對比 ×108 N/m2

4 結(jié) 論

1) 在創(chuàng)建單元時，用到了3個子程序文件，分別為：uec102，uel102，uep102.uec102用于描述單元的基本特征，主要是給ielc數(shù)組賦值；uel102用于編寫單元計算的核心程序；uep102主要用于線單元在后處理中的打印輸出，需要輸出的結(jié)果必須已經(jīng)寫入了結(jié)果文件中.

2) 通過經(jīng)典算例檢驗了文中使用UPFs二次開發(fā)工具進(jìn)行梁元開發(fā)所編寫程序的可靠性，同時說明文中在開發(fā)過程中關(guān)鍵步驟處理的合理性.

3) 在網(wǎng)殼結(jié)構(gòu)中常采用彈簧與梁元組合模擬網(wǎng)殼結(jié)構(gòu)中的構(gòu)件，此法建模復(fù)雜，文中應(yīng)用UPFs創(chuàng)建了能夠考慮節(jié)點剛度影響的梁元，并與彈簧梁元模型計算進(jìn)行了對比，結(jié)果吻合較好，說明了該方法的可行性，因此，也能為其它多構(gòu)件結(jié)構(gòu)，考慮節(jié)點影響提供參考.

[1]師訪.ANSYS二次開發(fā)及應(yīng)用實例詳解[M].北京:中國水利水電出版社,2011.

[2]陳正清,楊孟剛.梁桿索結(jié)構(gòu)幾何非線性有限元[M].北京：人民交通出版社，2013.

[3]ORAN C. Tangent stiffness in space frames[J]. Journal of the Structural Engineering,1973,99(ST6):987-1001.

[4]沈祖炎,羅永峰.網(wǎng)殼結(jié)構(gòu)分析中節(jié)點大位移迭加跟蹤技術(shù)的修正[J].空間結(jié)構(gòu)，1994(1):11-16.

[5]沈祖炎,羅永峰.節(jié)點大位移條件下的梁-柱單元坐標(biāo)轉(zhuǎn)換矩陣[J].上海力學(xué),1994,15(4):13-18.

[6]羅永峰.網(wǎng)殼結(jié)構(gòu)彈塑性穩(wěn)定及承載全過程研究[D].上海：同濟(jì)大學(xué)，1991.

[7]王新敏.ANSYS工程結(jié)構(gòu)數(shù)值分析[M].北京：人民交通出版社，2007.

[8]劉海峰,羅堯治,許賢.焊接球節(jié)點剛度對網(wǎng)殼結(jié)構(gòu)有限元分析精度的影響[J].工程力學(xué)，2013,30(1):350-358.

[9]王星，董石麟，完海鷹.焊接球節(jié)點剛度的有限元分析[J].浙江大學(xué)學(xué)報,2003,34(1)：77-82.

Study on the Development of Beam Element Applying the ANSYS User Programmable Features

HUANG Bin DENG Hongzhou

(DepartmentofStructuralEngineering,TongjiUniversity,Shanghai200092,China)

The User Programmable Features (UPFs) secondary development tools provided by ANSYS platform is to explore the process of establishing the beam element considering the geometric non-linearity. Sub program files uec102, uel102 and uep102 are studied which are used in the process of establishing the beam element, and the reliability of the program are checked through the classic example. A new type of beam element is developed considering the effect of joint stiffness to improve the efficiency of modeling and the whole calculation process. Besides, this new beam element is applied to a single layer spherical reticulated shell structure. Compared with the calculating method based on the common beam element and spring element, the results show good agreement, which proves the feasibility of the method based on the new type of beam element developed using UPFs.

secondary development of ANSYS; geometric nonlinearity; beam element development; joint stiffness

2016-09-07

*國家電網(wǎng)公司科技項目(521300130CSS)、國家自然科學(xué)基金項目(51578421)資助

TU3

10.3963/j.issn.2095-3844.2016.06.010

黃斌(1985—)：男，博士生，主要研究領(lǐng)域為輸電鐵塔節(jié)點極限承載力