應用ANSYS用戶可編程特性實現梁單元開發研究*

黃 斌 鄧洪洲

(同濟大學建筑工程系 上海 200092)

應用ANSYS用戶可編程特性實現梁單元開發研究*

黃 斌 鄧洪洲

(同濟大學建筑工程系 上海 200092)

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

ANSYS二次開發；幾何非線性；梁單元開發；節點剛度

0 引 言

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

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

UPFs是對ANSYS程序提供的可供用戶修改子程序文件(uex100～uex105)進行調整或修改，用于建立新型材料的本構關系，定義復合材料的失效準則，或根據用戶要求開發具有特定力學性能的單元等[1].

1 梁單元的開發過程探索

1.1 單元形函數

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

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

采用量綱一的量坐標

1.2 幾何非線性分析

1.2.1 空間大轉動處理

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

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

(1)

式中：ax，ay,az為旋轉軸的方向余弦.

1.2.2 單元坐標轉換矩陣

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

圖1 單元坐標系示意圖

梁單元坐標系與整體坐標系的轉換矩陣為

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

(2)

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

單元矩陣的y軸與z軸在整體坐標系中相對位置的確定比較復雜.設節點坐標系的任意時刻的坐標矩陣為

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

初始構形下節點坐標系分別平行于3根整體坐標系的坐標軸，則初始節點坐標系矩陣為

假設在t時刻αt已經求得t+Δt時刻的轉角增量為Δθ(Δθx，Δθy，Δθz)，則根據式(1)可以得到t+Δt時刻的節點定向矩陣αt+Δt.

αt+Δt=TΔθαt

(3)

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

Sit+Δt=αit+Δtr0

Sjt+Δt=αjt+Δtr0

(4)

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

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

(5)

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

由此可根據兩端截面的定向矩陣得到兩單元坐標系的定向矩陣

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

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

(6)

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

(7)

1.3 UPFs子程序文件說明

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

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

2 經典算例檢驗

為檢驗程序的可靠性及有效性，與經典算例進行了對比分析.

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

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

圖2 William框架結構簡圖(尺寸單位：mm)

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

圖3 點A處P-δA曲線

圖4 P-H曲線

圖5～6分別對應h=8.128 mm時點A處荷載p-豎向位移δA曲線及荷載p-水平支座反力H曲線.單元劃分數量均為4單元.

圖5 點A處P-δA曲線

圖6 P-H曲線

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

六角星形穹頂結構平衡路徑分析也是幾何非線性分析的經典算例，結構布置見圖7.

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

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

圖8～9為結構點1處的荷載位移曲線.圖8單元劃分數量均為1，圖9單元劃分數量為4.圖10為4單元時，荷載與點2處的x向水平位移ux關系曲線，圖11為4單元時點2處的z向位移uz關系曲線.

圖8 點1處P-δ1曲線

圖9 點1處P-δ1曲線

圖10 點2處P-ux曲線

圖11 點2處P-uz曲線

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

2.3 六邊形12構件空間框架幾何非線性分析

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

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

計算結果見圖13～14，圖13為單元劃分數量均為1單元，圖14為單元劃分數量均為4單元.兩者的荷載位移曲線均較好的吻合.

圖13 點A處P-δA曲線

圖14 點A處P-δA曲線

3 新型梁元的網殼分析應用

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

當不考慮節點尺寸的影響，僅考慮節點半剛性時，計算簡圖見圖15.

圖15 考慮節點半剛性單元計算簡圖

彈簧的力矩-轉角關系可以表示為

Ms=Ksθs

(8)

Ms=[MyaMzaMyaiMzaiMyb

MzbMybiMzbi]T

(9)

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

(10)

式中：θyai，θzai，θybi，θzbi為增加的內部節點自由度；Myai，Mzai，Mybi，Mzbi為與其對應的力矩.

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

構件采用普通梁元模擬每個單元有12個自由度，考慮彈簧之后，增加了4個內部自由度，此時，構件的自由度變為16個.將原構件的12×12剛度矩陣與彈簧的8×8剛度矩陣均擴展為16×16階矩陣，并且進行合并，則有：

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)

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

(16)

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

(17)

式(17)中

(18)

(19)

(20)

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

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

根據向量運算法則，a點的位移向量和荷載向量與A點的位移向量和荷載向量有如下關系

(21)

(22)

(23)

(24)

(25)

Kab與KAB之間的關系為

(26)

式(26)中，轉換矩陣的表達式為

(27)

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

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

圖17 單層球面網殼結構布置圖

圖18 單層球面網殼1/8結構

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

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

表1 工況1時的構件應力對比 ×108 N/m2

表2 工況2下的構件應力對比 ×108 N/m2

4 結 論

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

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

3) 在網殼結構中常采用彈簧與梁元組合模擬網殼結構中的構件，此法建模復雜，文中應用UPFs創建了能夠考慮節點剛度影響的梁元，并與彈簧梁元模型計算進行了對比，結果吻合較好，說明了該方法的可行性，因此，也能為其它多構件結構，考慮節點影響提供參考.

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

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

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

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

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

[6]羅永峰.網殼結構彈塑性穩定及承載全過程研究[D].上海：同濟大學，1991.

[7]王新敏.ANSYS工程結構數值分析[M].北京：人民交通出版社，2007.

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

[9]王星，董石麟，完海鷹.焊接球節點剛度的有限元分析[J].浙江大學學報,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

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

TU3

10.3963/j.issn.2095-3844.2016.06.010

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