浮 濤
(鄭州市交通規劃勘察設計研究院,河南 鄭州 450000)
有限元法是當今工程分析中應用最為廣泛的數值分析方法。在以往的結構有限元分析編程實踐中,通常采用BASIC、FORTRAN、C、C++語言等進行編程,相比較而言,MATLAB 被視為有效快捷的工程分析工具,其優點是可將編程計算、數據分析及數據圖形化成功地集成在工作環境中,被廣泛應用于多領域。使用者可編譯M 文件或者運用軟件自身提供的強大工具箱來進行動態仿真,可求解出一系列工程問題,且求出的數值結果具有可視化、圖形功能完善等特點。
在求解普通結構的工程分析中,使用較多的二維有限元是平面梁單元、平面桁架單元和平面剛架單元,三種單元都具有兩種坐標系,分別為總體坐標系和局部坐標系。
平面梁單元與平面剛架單元的參數有E、I、A、L(E代表彈性模量,I代表截面慣性矩,A代表截面面積,L代表單元長度),每個平面梁單元與平面剛架單元均有2個節點,并且相對于總體坐標系的X軸正向逆時針的傾斜角為θ。其中平面梁單元有4 個自由度(每個節點有2 個自由度即1 個位移自由度和1 個轉角自由度);平面剛架單元有6個自由度(每個節點有3個自由度即2個位移自由度和一個轉角自由度)[1-3]。
平面桁架單元的參數有E、A、L(E代表彈性模量,A代表截面面積,L代表單元長度),每個平面桁架單元有2個節點,并且相對于總體坐標系的X軸正向逆時針的傾斜角為θ。平面桁架單元有4 個自由度(每個節點有2個自由度即2個位移自由度)。
上述平面單元均約定位移向上為正,轉角逆時針方向為正。
本文采用MATLAB 進行有限元分析的方法,并對某平面桁架結構、平面剛架結構進行實例分析研究。
結構的離散化,寫出單元剛度矩陣,集成整體剛度矩陣,引入邊界條件,求解方程組,結果處理。
1.2.1 平面桁架單元的M函數文件
(1)PTEStiffness(E,A,L,θ)—計算平面桁架單元(E,A,L,θ)的單剛矩陣。
(2)PTAssemble(K,k,i,j)—將連接節點i、j的平面桁架單元的單剛矩陣k集成到整剛矩陣K。
(3)PTEForce(E,A,L,θ,u)—計算單元節點力矢量。(4)PTEStress(E,L,θ,u)—計算單元應力。
1.2.2 平面剛架單元的M函數文件
(5)PFEStiffness(E,A,I,L,θ)—計算平面剛架單元(E,A,I,L,theta)的單剛矩陣。
(6)PFAssemble(K,k,i,j)—將連接節點i、j的平面剛架單元的單剛矩陣k集成到整剛矩陣K。
(7)PFEForces(E,A,I,L,θ,u)—計算平面剛架單元的節點力矢量。
(8)PFEADiagram(f,L)—繪制節點力矢量為f和長度為L的單元軸力圖。
(9)PFESDiagram(f,L)—繪制節點力矢量為f和長度為L的單元剪力圖。
(10)PFEMDiagram(f,L)—繪制節點力矢量為f和長度為L的單元彎矩圖。
在MATLAB軟件集成環境中,若進行平面梁單元結構分析,使用(1)~(4)這4個M函數文件,根據實際結構給出函數所需要的參數直接調用即可;若進行平面桁架單元結構分析,使用(5)~(10)這6個M函數文件,根據實際結構給出函數所需要的參數直接調用即可[4]。
MATLAB 軟件的使用方法,本文不再詳細敘述,用MATLAB 軟件寫出矩陣方程和在MATLAB 集成環境中編寫并使用M函數文件請參閱文獻[5]。
如圖1 的桁架,已知各桿E及A均相同,且E=210GPa,A=1E-2m2,L1=400mm,L2=300mm,L3=500mm,L4=400mm,求解各桿件內力。
圖1 桁架結構圖
根據M函數文件PTEStiffness()求出單元的剛度矩陣k1、k2、k3、k4,該結構有四個節點,所以整體剛度矩陣為8×8。由于該結構有四個桁架單元,故需要四次調用函數文件PTAssemble(),求得整體剛度矩陣。
被約束的自由度有:U1x=U1y=U2y=U4x=U4y=0。則活動自由度編號為3,5,6。活動自由度對應的節點載荷F3=20000N,F5=0,F6=-25000N,采用高斯消去法進行求解,最后求得節點2的水平位移3.81e-3(mm),節點3的水平位移7.94e-4(mm),節點3的豎向位移3.125e-3(mm)。建立結構的節點位移矢量U,然后計算結構的節點力矢量F。由此可求得節點1的水平支反力-15.833kN(方向向左),垂直支反力3.125kN(方向向上),節點2的水平支反力20kN(方向向右),垂直支反力21.875kN(方向向上),節點4的水平支反力4.167kN(方向向左),垂直支反力0,滿足力的平衡條件。
建 立 單 元 位 移 矢 量u1、u2、u3和u4,然 后 調 用PTEForces()求解各桿軸力,調用PTEStress()求解各桿應力,見表1所示。
表1 各單元軸力及應力[4]
如圖2的剛架,已知各桿E、I及A均相同,且A=1E-2m2,I=3E-5m4,L1=6m,L2=8m,L3=6m,E=210GPa,求節點2和3的位移和轉角,并采用MATLAB繪出單元2的剪力圖和彎矩圖。
圖2 剛架結構圖
根據M函數文件PFEStiffness()求出單元的剛度矩陣k1、k2、k3,該結構有四個節點,所以整體剛度矩陣為12×12。由于該結構有三個剛架單元,故需要三次調用函數文件PFAssemble(),求得整體剛度矩陣。
被約束的自由度有:U1x=U1y=R1z=0;U4x=U4y=R4z=0。則活動自由度編號為4,5,6,7,8,9。活動自由度對應的節點載荷F4=-30000N,F5=0,F6=0,F7=0,F8=0,F9=-20000kN·m,采用高斯消去法進行求解,最后求得節點2 的水平位移-0.0611m(向左),豎向位移0,轉角位移0.0078rad(逆時針);節點3 的水平位移-0.0610m(向左),豎向位移0,轉角位移0.0043rad(逆時針)。
建立結構的節點位移矢量U,然后計算結構的節點力矢量F。由此可求得節點1 的水平支反力13.187kN(方向向右),垂直支反力7.158kN(方向向上),彎矩為47.753kN·m(順時針)節點4 的水平支反力16.813kN(方向向右),垂直支反力7.158kN(方向向下),彎矩為54.983kN·m(順時針),滿足力的平衡條件。
表2 建立單元位移矢量u1、u2和u3,然后調用PFE Forces()函數求出單元矢量f1、f2和f3。
表2 單元軸力、剪力、彎矩值[4](單位:kN)
最后調用PFEADiagram()、PFESDiagram()以及PFEMDiagram()可分別繪出單元的軸力圖、剪力圖和彎矩圖。單元2 的軸力圖、剪力圖和彎矩圖見圖3、圖4、圖5。
圖3 單元2軸力圖(單位:N)
圖4 單元2剪力圖(單位:N)
圖5 單元2彎矩圖(單位:N)
本文以平面桁架結構和平面剛架結構為例,利用MATLAB 進行有限元結構分析,經實踐證明,該方法簡便、可行、實用。本文雖僅從平面的角度為例驗證其可行性,但其原理和方法可推廣到其他單元,如板殼、實體單元的剛度矩陣,甚至三維空間和非線性剛度矩陣之中。同時,該方法可為相關工程人員繼續深入研究提供一定參考依據。