算例计算如图1所示平面刚架各结点的位移和各梁的内力及支座反力.docx

1、算例计算如图1所示平面刚架各结点的位移和各梁的内力及支座反力算例计算如图-1所示平面刚架各结点的位移和各梁的内力及支座反力。已知： 图-1受分布力和集中力的平面刚架分析：首先建立有限元模型，即定义结点坐标，定义单元的结点号和材料特性，定义约束条件，给定结点力等。把杆件的连接点和集中力的作用点取为结点，并按17编号，其中5号结点就是集中力作用点，如下图所示。为了确定结点的坐标，我们要建立一个整体坐标系，其原点为结点1，水平向右为x轴正方向，竖直向上为y轴正方向。这样就能根据框架的尺寸确定7个结点的坐标。然后将这7个结点两两组合成6个单元，并按编号，示于图-2中。该刚架有3个结点被固定，每个结点有

2、3个自由度，因此共有9个自由度被约束。另外还有1个集中力，两个单元有分布荷载。这些信息即构成了有限元模型，可编制一个函数程序PlaneFrameModel来指定。图-2 单元划分图有限元模型生成函数function PlaneFrameModel% 定义平面杆系的有限元模型% 输入参数：% 无% 返回值：% 无% 说明：% 该函数定义平面杆系的有限元模型数据：% gNode - 节点定义% gElement - 单元定义% gMaterial - 材料定义，包括弹性模量，梁的截面积和梁的抗弯惯性矩% gBC1 - 约束条件% gNF - 集中力% gDF - 分布力 global gNode

3、gElement gMaterial gBC1 gNF gDF % 节点坐标 % x y gNode = 0.0, 0.0 % 节点 1 0.0, 4.0 % 节点 2 3.0, 0.0 % 节点 3 3.0, 4.0 % 节点 4 4.5, 4.0 % 节点 5 6.0, 0.0 % 节点 6 6.0, 4.0 ; % 节点 7 % 单元定义 % 节点1 节点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 % 材料性质 %

4、弹性模量 抗弯惯性矩 截面积 gMaterial = 2.1e11, 2.0e-4, 1.0e-2 ; % 材料 1 % 第一类约束条件 % 节点号 自由度号 约束值 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 ; % 集中力 % 节点号 自由度号 集中力值 gNF = 5, 2, -80e3 ; % 分布载荷（线性分布） % 单元号 节点1载荷值 节点2载荷值 自由度号 gDF = 1 -30e3 0 2 2 -15e3 -15e3 2 ;re

5、turn有了有限元模型数据，下一步的工作就是求解。其具体流程如下图所示图-3有限元程序流程图计算分析主程序MATLAB源程序代码function plane-beam% 本程序为采用平面梁单元计算平面刚架的变形和内力% 输入参数： 无% 输出结果： 节点位移和单元节点力 PlaneFrameModel ; % 定义有限元模型 SolveModel ; % 求解有限元模型 DisplayResults ; % 显示计算结果return ;function PlaneFrameModel% 定义平面杆系的有限元模型% 输入参数：% 无% 返回值：% 无% 说明：% 该函数定义平面杆系的有限元模型数

6、据：% gNode - 节点定义% gElement - 单元定义% gMaterial - 材料定义，包括弹性模量，梁的截面积和梁的抗弯惯性矩% gBC1 - 约束条件% gNF - 集中力% gDF - 分布力 global gNode gElement gMaterial gBC1 gNF gDF % 节点坐标 % x y gNode = 0.0, 0.0 % 节点 1 0.0, 4.0 % 节点 2 3.0, 0.0 % 节点 3 3.0, 4.0 % 节点 4 4.5, 4.0 % 节点 5 6.0, 0.0 % 节点 6 6.0, 4.0 ; % 节点 7 % 单元定义 % 节点1

7、 节点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.0e-4, 1.0e-2 ; % 材料 1 % 第一类约束条件 % 节点号 自由度号 约束值 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.

8、0 ; % 集中力 % 节点号 自由度号 集中力值 gNF = 5, 2, -80e3 ; % 分布载荷（线性分布） % 单元号 节点1载荷值 节点2载荷值 自由度号 gDF = 1 -30e3 0 2 2 -15e3 -15e3 2 ;returnfunction SolveModel% 求解有限元模型% 输入参数：% 无% 返回值：% 无% 说明：% 该函数求解有限元模型，过程如下% 1. 计算单元刚度矩阵，集成整体刚度矩阵% 2. 计算单元的等效节点力，集成整体节点力向量% 3. 处理约束条件，修改整体刚度矩阵和节点力向量% 4. 求解方程组，得到整体节点位移向量 global gNod

9、e gElement gMaterial gBC1 gNF gDF gK gDelta % step1. 定义整体刚度矩阵和节点力向量 node_number,dummy = size( gNode ) ; gK = sparse( node_number * 3, node_number * 3 ) ; f = sparse( node_number * 3, 1 ) ; % step2. 计算单元刚度矩阵，并集成到整体刚度矩阵中 element_number,dummy = size( gElement ) ; for ie=1:1:element_number k = Stiffness

10、Matrix( ie, 1 ) ; AssembleStiffnessMatrix( ie, k ) ; end % step3. 把集中力直接集成到整体节点力向量中 nf_number, dummy = size( gNF ) ; for inf=1:1:nf_number n = gNF( inf, 1 ) ; d = gNF( inf, 2 ) ; f( (n-1)*3 + d ) = gNF( inf, 3 ) ; end % step4. 计算分布力的等效节点力，并集成到整体节点力向量中 df_number, dummy = size( gDF ) ; for idf = 1:1:d

11、f_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( (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 ) ; e

12、nd % step5. 处理约束条件，修改刚度矩阵和节点力向量。采用乘大数法 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) = gBC1(ibc, 3)* gK(m,m) * 1e15 ; gK(m,m) = gK(m,m) * 1e15 ; end % step 6. 求解方程组，得到节点位移向量 gDelta = gK f ;returnfunction DisplayResults% 显示计算结果%

13、输入参数：% 无% 返回值：% 无 global gNode gElement gMaterial gBC1 gNF gDF gK gDelta fprintf( 节点位移n ) ; fprintf( 节点号 x方向位移 y方向位移 转角n ) ; node_number,dummy = size( gNode ) ; 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节点力n )

14、; fprintf( 轴力 剪力 弯矩n ) ; element_number, dummy = size( gElement ) ; for ie = 1:element_number enf = ElementNodeForce( ie ) ; fprintf( 单元号%6d 节点号%6d %16.8e %16.8e %16.8en, . ie, gElement(ie,1), enf(1), enf(2), enf(3) ) ; fprintf( 节点号%6d %16.8e %16.8e %16.8en, . gElement(ie,2), enf(4), enf(5), enf(6)

15、) ; endreturnfunction k = StiffnessMatrix( ie, icoord )% 计算单元刚度矩阵% 输入参数:% ie - 单元号% icoord - 坐标系参数，可以是下面两个之一% 1 - 整体坐标系% 2 - 局部坐标系% 返回值:% k - 根据icoord的值，相应坐标系下的刚度矩阵 global gNode gElement gMaterial k = zeros( 6, 6 ) ; E = gMaterial( gElement(ie, 3), 1 ) ; I = gMaterial( gElement(ie, 3), 2 ) ; A = gMa

16、terial( gElement(ie, 3), 3 ) ; 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 = ( (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

17、 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 = TransformMatrix( ie ) ; k = T*k*transpose(T) ; endreturnfunction AssembleStiffnessMatrix( ie, k )% 把单元刚度矩阵集成到整体刚度矩阵% 输入参数:% ie - 单元号% k - 单元刚度矩阵% 返回值:% 无

18、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 ; gK(M,N) = gK(M,N) + k(m,n) ; end end end endreturnfunction enf = EquivalentNodeForce( ie, p1, p2, idof )% 计算线性分布荷载的等效节点力% 输入参数:% ie - 单元号% p1

19、- 第一个节点上的分布力集度值% p2 - 第二个节点上的分布力集度值% idof - 分布力的种类，它可以是下面几种% 1 - 分布轴向力% 2 - 分布横向力% 3 - 分布弯矩% 返回值:% enf - 整体坐标系下等效节点力向量 global gElement gNode enf = zeros( 6, 1 ) ; % 定义 6x1 的等效节点力向量 xi = gNode( gElement( ie, 1 ), 1 ) ; yi = gNode( gElement( ie, 1 ), 2 ) ; xj = gNode( gElement( ie, 2 ), 1 ) ; yj = gNo

20、de( 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 % 分布横向力 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(

21、 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 = TransformMatrix( ie ) ; % 计算单元的转换矩阵 enf = T * enf ; % 把等效节点力转换到整体坐标下returnfunction T = TransformMatrix( ie )% 计算单元的坐标转换矩阵( 局部坐标 - 整体坐标 )% 输入参数% ie

22、 - 节点号% 返回值% T - 从局部坐标到整体坐标的坐标转换矩阵 global gElement gNode 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 ) ; c = (xj-xi)/L ; s = (yj-yi)/L ; T= c -s 0 0 0 0 s c 0 0 0 0 0

23、 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 )% 计算单元的节点力% 输入参数% ie - 节点号% 返回值% enf - 单元局部坐标系下的节点力 global gElement 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 )

24、= 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 ) enf = enf - EquivalentNodeForce( gDF(idf,1), . gDF(idf, 2), gDF( idf, 3), gDF( idf, 4 ) ) ; break ; end end T = TransformMatrix( ie ) ; enf =

25、 transpose( T ) * enf ;return程序计算机如果如下节点位移 节点号 x方向位移 y方向位移 转角 1 7.19450697e-019 -3.10908802e-020 -2.78411083e-019 2 7.88387267e-004 -3.10908802e-005 -3.44682851e-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

26、-004 2.15592934e-005 6 1.29964769e-018 -7.70240078e-020 -4.19430144e-019 7 7.56145765e-004 -7.70240078e-005 2.71750964e-004节点力 轴力 剪力 弯矩单元号 1 节点号 1 1.63227121e+004 4.76656742e+004 3.56932655e+004 节点号 2 -1.63227121e+004 1.23343258e+004 -5.03056852e+003单元号 2 节点号 2 1.23343258e+004 1.63227121e+004 5.0305

27、6852e+003 节点号 4 -1.23343258e+004 2.86772879e+004 -2.35624322e+004单元号 3 节点号 3 6.82396838e+004 2.09960018e+003 6.84599262e+003 节点号 4 -6.82396838e+004 -2.09960018e+003 1.55240811e+003单元号 4 节点号 4 1.02347256e+004 3.95623959e+004 2.20100241e+004 节点号 5 -1.02347256e+004 -3.95623959e+004 3.73335698e+004单元号 5

28、 节点号 5 1.02347256e+004 -4.04376041e+004 -3.73335698e+004 节点号 7 -1.02347256e+004 4.04376041e+004 -2.33228363e+004单元号 6 节点号 6 4.04376041e+004 1.02347256e+004 1.76160660e+004 节点号 7 -4.04376041e+004 -1.02347256e+004 2.33228363e+004ANSYS计算结果对比* POST1 NODAL DEGREE OF FREEDOM LISTING * THE FOLLOWING DEGREE

29、 OF FREEDOM RESULTS ARE IN GLOBAL COORDINATES NODE UX UY UZ ROTX ROTY ROTZ 1 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 2 0.78839E-03 -0.31091E-04 0.0000 0.0000 0.0000 -0.34468E-04 3 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 4 0.77077E-03 -0.12998E-03 0.0000 0.0000 0.0000 -0.25208E-03 5 0.76346E-03 -0.56779E-03 0.0000 0.0000 0.0000 0.21559E-04 6 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 7 0.75615E-03 -0.77024E-04 0.0000 0.0000 0.0000 0.27175E-03* POST1 NODAL TOTAL FORCE SUMMATION * THE FOLLOWING X