計算如圖-1所示平面剛架各結點的位移和各梁的內力及支

1、算例計算如圖-1所示平面剛架各結點的位移和各梁的內力及支座反力。已知：圖-1受分布力和集中力的平面剛架分析：首先建立有限元模型，即定義結點坐標，定義單元的結點號和材料特性，定義約束條件，給定結點力等。把桿件的連接點和集中力的作用點取為結點，并按17編號，其中5號結點就是集中力作用點，如下圖所示。為了確定結點的坐標，我們要建立一個整體坐標系，其原點為結點1，水平向右為x軸正方向，豎直向上為y軸正方向。這樣就能根據(jù)框架的尺寸確定7個結點的坐標。然后將這7個結點兩兩組合成6個單元，并按編號，示于圖-2中。該剛架有3個結點被固定，每個結點有3個自由度，因此共有9個自由度被約束。另外還有1個集中力，兩個

2、單元有分布荷載。這些信息即構成了有限元模型，可編制一個函數(shù)程序PlaneFrameModel來指定。圖-2 單元劃分圖有限元模型生成函數(shù)function PlaneFrameModel% 定義平面桿系的有限元模型% 輸入?yún)?shù)：% 無% 返回值：% 無% 說明：% 該函數(shù)定義平面桿系的有限元模型數(shù)據(jù)：% gNode - 節(jié)點定義% gElement - 單元定義% gMaterial - 材料定義，包括彈性模量，梁的截面積和梁的抗彎慣性矩% gBC1 - 約束條件% gNF - 集中力% gDF - 分布力 global gNode gElement gMaterial gBC1 gNF gDF

3、 % 節(jié)點坐標 % x y gNode = 0.0, 0.0 % 節(jié)點 1 0.0, 4.0 % 節(jié)點 2 3.0, 0.0 % 節(jié)點 3 3.0, 4.0 % 節(jié)點 4 4.5, 4.0 % 節(jié)點 5 6.0, 0.0 % 節(jié)點 6 6.0, 4.0 ; % 節(jié)點 7 % 單元定義 % 節(jié)點1 節(jié)點2 材料號 gElement = 1, 2, 1 % 單元 1 2, 4, 1 % 單元 2 3, 4, 1 % 單元 3 4, 5, 1 % 單元 4 5, 7, 1 % 單元 5 6, 7, 1 ; % 單元 6 % 材料性質 % 彈性模量 抗彎慣性矩 截面積 gMaterial = 2.1e

4、11, 2.0e-4, 1.0e-2 ; % 材料 1 % 第一類約束條件 % 節(jié)點號 自由度號 約束值 gBC1 = 1, 1, 0.0 1, 2, 0.0 1, 3, 0.0 3, 1, 0.0 3, 2, 0.0 3, 3, 0.0 6, 1, 0.0 6, 2, 0.0 6, 3, 0.0 ; % 集中力 % 節(jié)點號 自由度號 集中力值 gNF = 5, 2, -80e3 ; % 分布載荷（線性分布） % 單元號 節(jié)點1載荷值 節(jié)點2載荷值 自由度號 gDF = 1 -30e3 0 2 2 -15e3 -15e3 2 ;return有了有限元模型數(shù)據(jù)，下一步的工作就是求解。其具體流程如

5、下圖所示建立有限元模型（PlaneFrameModel）求解有限元模型（Solve Model）顯示計算結果(Display Results）結束計算單元剛度矩陣（Stiffness Matrix）集成整體剛度矩陣（Assemble Stiffness Matrix）計算等效結點力（Equivalent Node Force）處理約束條件求解方程組開始圖-3有限元程序流程圖計算分析主程序MATLAB源程序代碼function plane-beam% 本程序為采用平面梁單元計算平面剛架的變形和內力% 輸入?yún)?shù)： 無% 輸出結果： 節(jié)點位移和單元節(jié)點力 PlaneFrameModel ; % 定義

6、有限元模型 SolveModel ; % 求解有限元模型 DisplayResults ; % 顯示計算結果return ;function PlaneFrameModel% 定義平面桿系的有限元模型% 輸入?yún)?shù)：% 無% 返回值：% 無% 說明：% 該函數(shù)定義平面桿系的有限元模型數(shù)據(jù)：% gNode - 節(jié)點定義% gElement - 單元定義% gMaterial - 材料定義，包括彈性模量，梁的截面積和梁的抗彎慣性矩% gBC1 - 約束條件% gNF - 集中力% gDF - 分布力 global gNode gElement gMaterial gBC1 gNF gDF % 節(jié)點坐

7、標 % x y gNode = 0.0, 0.0 % 節(jié)點 1 0.0, 4.0 % 節(jié)點 2 3.0, 0.0 % 節(jié)點 3 3.0, 4.0 % 節(jié)點 4 4.5, 4.0 % 節(jié)點 5 6.0, 0.0 % 節(jié)點 6 6.0, 4.0 ; % 節(jié)點 7 % 單元定義 % 節(jié)點1 節(jié)點2 材料號 gElement = 1, 2, 1 % 單元 1 2, 4, 1 % 單元 2 3, 4, 1 % 單元 3 4, 5, 1 % 單元 4 5, 7, 1 % 單元 5 6, 7, 1 ; % 單元 6 % 材料性質 % 彈性模量 抗彎慣性矩 截面積 gMaterial = 2.1e11, 2.

8、0e-4, 1.0e-2 ; % 材料 1 % 第一類約束條件 % 節(jié)點號 自由度號 約束值 gBC1 = 1, 1, 0.0 1, 2, 0.0 1, 3, 0.0 3, 1, 0.0 3, 2, 0.0 3, 3, 0.0 6, 1, 0.0 6, 2, 0.0 6, 3, 0.0 ; % 集中力 % 節(jié)點號 自由度號 集中力值 gNF = 5, 2, -80e3 ; % 分布載荷（線性分布） % 單元號 節(jié)點1載荷值 節(jié)點2載荷值 自由度號 gDF = 1 -30e3 0 2 2 -15e3 -15e3 2 ;returnfunction SolveModel% 求解有限元模型% 輸入?yún)?/p>

9、數(shù)：% 無% 返回值：% 無% 說明：% 該函數(shù)求解有限元模型，過程如下% 1. 計算單元剛度矩陣，集成整體剛度矩陣% 2. 計算單元的等效節(jié)點力，集成整體節(jié)點力向量% 3. 處理約束條件，修改整體剛度矩陣和節(jié)點力向量% 4. 求解方程組，得到整體節(jié)點位移向量 global gNode gElement gMaterial gBC1 gNF gDF gK gDelta % step1. 定義整體剛度矩陣和節(jié)點力向量 node_number,dummy = size( gNode ) ; gK = sparse( node_number * 3, node_number * 3 ) ; f =

10、sparse( node_number * 3, 1 ) ; % step2. 計算單元剛度矩陣，并集成到整體剛度矩陣中 element_number,dummy = size( gElement ) ; for ie=1:1:element_number k = StiffnessMatrix( ie, 1 ) ; AssembleStiffnessMatrix( ie, k ) ; end % step3. 把集中力直接集成到整體節(jié)點力向量中 nf_number, dummy = size( gNF ) ; for inf=1:1:nf_number n = gNF( inf, 1 ) ;

11、 d = gNF( inf, 2 ) ; f( (n-1)*3 + d ) = gNF( inf, 3 ) ; end % step4. 計算分布力的等效節(jié)點力，并集成到整體節(jié)點力向量中 df_number, dummy = size( gDF ) ; for idf = 1:1:df_number enf = EquivalentNodeForce( gDF(idf,1), gDF(idf, 2), gDF( idf, 3), gDF( idf, 4 ) ) ; i = gElement( gDF(idf,1), 1 ) ; j = gElement( gDF(idf,1), 2 ) ; f

12、( (i-1)*3+1 : (i-1)*3+3 ) = f( (i-1)*3+1 : (i-1)*3+3 ) + enf( 1:3 ) ; f( (j-1)*3+1 : (j-1)*3+3 ) = f( (j-1)*3+1 : (j-1)*3+3 ) + enf( 4:6 ) ; end % step5. 處理約束條件，修改剛度矩陣和節(jié)點力向量。采用乘大數(shù)法 bc_number,dummy = size( gBC1 ) ; for ibc=1:1:bc_number n = gBC1(ibc, 1 ) ; d = gBC1(ibc, 2 ) ; m = (n-1)*3 + d ; f(m) =

13、 gBC1(ibc, 3)* gK(m,m) * 1e15 ; gK(m,m) = gK(m,m) * 1e15 ; end % step 6. 求解方程組，得到節(jié)點位移向量 gDelta = gK f ;returnfunction DisplayResults% 顯示計算結果% 輸入?yún)?shù)：% 無% 返回值：% 無 global gNode gElement gMaterial gBC1 gNF gDF gK gDelta fprintf( 節(jié)點位移n ) ; fprintf( 節(jié)點號 x方向位移 y方向位移 轉角n ) ; node_number,dummy = size( gNode )

14、 ; for i=1:node_number fprintf( %6d %16.8e %16.8e %16.8en,. i, gDelta(i-1)*3+1), gDelta(i-1)*3+2), gDelta(i-1)*3+3) ) ; end fprintf( nn節(jié)點力n ) ; fprintf( 軸力 剪力 彎矩n ) ; element_number, dummy = size( gElement ) ; for ie = 1:element_number enf = ElementNodeForce( ie ) ; fprintf( 單元號%6d 節(jié)點號%6d %16.8e %16

15、.8e %16.8en, . ie, gElement(ie,1), enf(1), enf(2), enf(3) ) ; fprintf( 節(jié)點號%6d %16.8e %16.8e %16.8en, . gElement(ie,2), enf(4), enf(5), enf(6) ) ; endreturnfunction k = StiffnessMatrix( ie, icoord )% 計算單元剛度矩陣% 輸入?yún)?shù):% ie - 單元號% icoord - 坐標系參數(shù)，可以是下面兩個之一% 1 - 整體坐標系% 2 - 局部坐標系% 返回值:% k - 根據(jù)icoord的值，相應坐標系

16、下的剛度矩陣 global gNode gElement gMaterial k = zeros( 6, 6 ) ; E = gMaterial( gElement(ie, 3), 1 ) ; I = gMaterial( gElement(ie, 3), 2 ) ; A = gMaterial( gElement(ie, 3), 3 ) ; xi = gNode( gElement( ie, 1 ), 1 ) ; yi = gNode( gElement( ie, 1 ), 2 ) ; xj = gNode( gElement( ie, 2 ), 1 ) ; yj = gNode( gEle

17、ment( ie, 2 ), 2 ) ; L = ( (xj-xi)2 + (yj-yi)2 )(1/2) ; k = E*A/L 0 0 -E*A/L 0 0 0 12*E*I/L3 6*E*I/L2 0 -12*E*I/L3 6*E*I/L2 0 6*E*I/L2 4*E*I/L 0 -6*E*I/L2 2*E*I/L -E*A/L 0 0 E*A/L 0 0 0 -12*E*I/L3 -6*E*I/L2 0 12*E*I/L3 -6*E*I/L2 0 6*E*I/L2 2*E*I/L 0 -6*E*I/L2 4*E*I/L ; if icoord = 1 T = TransformMat

18、rix( ie ) ; k = T*k*transpose(T) ; endreturnfunction AssembleStiffnessMatrix( ie, k )% 把單元剛度矩陣集成到整體剛度矩陣% 輸入?yún)?shù):% ie - 單元號% k - 單元剛度矩陣% 返回值:% 無 global gElement gK for i=1:1:2 for j=1:1:2 for p=1:1:3 for q =1:1:3 m = (i-1)*3+p ; n = (j-1)*3+q ; M = (gElement(ie,i)-1)*3+p ; N = (gElement(ie,j)-1)*3+q ;

19、gK(M,N) = gK(M,N) + k(m,n) ; end end end endreturnfunction enf = EquivalentNodeForce( ie, p1, p2, idof )% 計算線性分布荷載的等效節(jié)點力% 輸入?yún)?shù):% ie - 單元號% p1 - 第一個節(jié)點上的分布力集度值% p2 - 第二個節(jié)點上的分布力集度值% idof - 分布力的種類，它可以是下面幾種% 1 - 分布軸向力% 2 - 分布橫向力% 3 - 分布彎矩% 返回值:% enf - 整體坐標系下等效節(jié)點力向量 global gElement gNode enf = zeros( 6, 1

20、 ) ; % 定義 6x1 的等效節(jié)點力向量 xi = gNode( gElement( ie, 1 ), 1 ) ; yi = gNode( gElement( ie, 1 ), 2 ) ; xj = gNode( gElement( ie, 2 ), 1 ) ; yj = gNode( gElement( ie, 2 ), 2 ) ; L = sqrt( (xj-xi)2 + (yj-yi)2 ) ; switch idof case 1 % 分布軸向力 enf( 1 ) = (2*p1+p2)*L/6 ; enf( 4 ) = (p1+2*p2)*L/6 ; case 2 % 分布橫向力

21、 enf( 2 ) = (7*p1+3*p2)*L/20 ; enf( 3 ) = (3*p1+2*p2)*L2/60 ; enf( 5 ) = (3*p1+7*p2)*L/20 ; enf( 6 ) = -(2*p1+3*p2)*L2/60 ; case 3 % 分布彎矩 enf( 2 ) = -(p1+p2)/2 ; enf( 3 ) = (p1-p2)*L/12 ; enf( 5 ) = (p1+p2)/2 ; enf( 6 ) = -(p1-p2)*L/12 ; otherwise disp( sprintf( 分布力的種類錯誤，單元號:%d,ie ) ) ; end T = Tran

22、sformMatrix( ie ) ; % 計算單元的轉換矩陣 enf = T * enf ; % 把等效節(jié)點力轉換到整體坐標下returnfunction T = TransformMatrix( ie )% 計算單元的坐標轉換矩陣( 局部坐標 - 整體坐標 )% 輸入?yún)?shù)% ie - 節(jié)點號% 返回值% T - 從局部坐標到整體坐標的坐標轉換矩陣 global gElement gNode xi = gNode( gElement( ie, 1 ), 1 ) ; yi = gNode( gElement( ie, 1 ), 2 ) ; xj = gNode( gElement( ie, 2

23、 ), 1 ) ; yj = gNode( gElement( ie, 2 ), 2 ) ; L = sqrt( (xj-xi)2 + (yj-yi)2 ) ; c = (xj-xi)/L ; s = (yj-yi)/L ; T= c -s 0 0 0 0 s c 0 0 0 0 0 0 1 0 0 0 0 0 0 c -s 0 0 0 0 s c 0 0 0 0 0 0 1 ;returnfunction enf = ElementNodeForce( ie )% 計算單元的節(jié)點力% 輸入?yún)?shù)% ie - 節(jié)點號% 返回值% enf - 單元局部坐標系下的節(jié)點力 global gElemen

24、t gNode gDelta gDF i = gElement( ie, 1 ) ; j = gElement( ie, 2 ) ; de = zeros( 6, 1 ) ; de( 1:3 ) = gDelta( (i-1)*3+1:(i-1)*3+3 ) ; de( 4:6 ) = gDelta( (j-1)*3+1:(j-1)*3+3 ) ; k = StiffnessMatrix( ie, 1 ) ; enf = k * de ; df_number, dummy = size( gDF ) ; for idf = 1:1:df_number if ie = gDF( idf, 1 )

25、 enf = enf - EquivalentNodeForce( gDF(idf,1), . gDF(idf, 2), gDF( idf, 3), gDF( idf, 4 ) ) ; break ; end end T = TransformMatrix( ie ) ; enf = transpose( T ) * enf ;return程序計算機如果如下節(jié)點位移 節(jié)點號 x方向位移 y方向位移 轉角 1 7.19450697e-019 -3.10908802e-020 -2.78411083e-019 2 7.88387267e-004 -3.10908802e-005 -3.446828

26、51e-005 3 2.66615896e-019 -1.29980350e-019 -1.62999824e-019 4 7.70766801e-004 -1.29980350e-004 -2.52075453e-004 5 7.63456283e-004 -5.67794228e-004 2.15592934e-005 6 1.29964769e-018 -7.70240078e-020 -4.19430144e-019 7 7.56145765e-004 -7.70240078e-005 2.71750964e-004節(jié)點力 軸力 剪力 彎矩單元號 1 節(jié)點號 1 1.63227121e

27、+004 4.76656742e+004 3.56932655e+004 節(jié)點號 2 -1.63227121e+004 1.23343258e+004 -5.03056852e+003單元號 2 節(jié)點號 2 1.23343258e+004 1.63227121e+004 5.03056852e+003 節(jié)點號 4 -1.23343258e+004 2.86772879e+004 -2.35624322e+004單元號 3 節(jié)點號 3 6.82396838e+004 2.09960018e+003 6.84599262e+003 節(jié)點號 4 -6.82396838e+004 -2.0996001

28、8e+003 1.55240811e+003單元號 4 節(jié)點號 4 1.02347256e+004 3.95623959e+004 2.20100241e+004 節(jié)點號 5 -1.02347256e+004 -3.95623959e+004 3.73335698e+004單元號 5 節(jié)點號 5 1.02347256e+004 -4.04376041e+004 -3.73335698e+004 節(jié)點號 7 -1.02347256e+004 4.04376041e+004 -2.33228363e+004單元號 6 節(jié)點號 6 4.04376041e+004 1.02347256e+004 1.76160660e+004 節(jié)點號 7 -4.04376041e+004 -1.02347256e+004 2.33228363e+004ANSYS計算結果對比* POST1 NODAL DEGREE OF FREEDOM LISTING * THE FOLLOWING DEGREE OF FREE