一用三结点三角形平面单元计算平面结构的应力和位移

1、一：用三结点三角形平面单元计算平面结构的应力和位移。1，设计说明书计算简图，网格划分，单元及结点的编号如下图所示。由于结构对称，去四分之一结构分析。其中E=2e10pa,mu=0.167,h=1m.开始输入节点坐标进行单元定义给出材料性质定义有限元模型引入约束条件给定荷载信息计算单刚形成总刚形成结点力求解有限元模型引进位移约束求解平衡方程显示结点位移显示计算结果显示单元应力结束变量注释：Node - 节点定义gElement - 单元定义gMaterial - 材料定义，包括弹性模量，泊松比和厚度gBC1 - 约束条件gNF - 集中力gk-总刚gDelta-结点位移子程序注释： PlaneS

2、tructualModel 定义有限元模型 SolveModel 求解有限元模型 DisplayResults 显示计算结果 k = StiffnessMatrix( ie )计算单元刚度 AssembleStiffnessMatrix( ie, k )形成总刚 es = ElementStress( ie )计算单元应力function exam1% 输入参数： 无% 输出结果： 节点位移和单元应力PlaneStructualModel ; % 定义有限元模型 SolveModel ; % 求解有限元模型 DisplayResults ; % 显示计算结果return ;function P

3、laneStructualModel% 定义平面结构的有限元模型% 输入参数：无% 说明：% 该函数定义平面结构的有限元模型数据：% gNode - 节点定义% gElement - 单元定义% gMaterial - 材料定义，包括弹性模量，泊松比和厚度% gBC1 - 约束条件% gNF - 集中力 global gNode gElement gMaterial gBC1 gNF % 节点坐标 % x y gNode = 0.0, 2.0 % 节点 1 0.0, 1.0 % 节点 2 1.0, 1.0 % 节点 3 0.0, 0.0 % 节点 4 1.0, 0.0 % 节点 5 2.0,

4、0.0 ; % 节点 6 % 单元定义 % 节点1 节点2 节点3 材料号 gElement = 3, 1, 2, 1 % 单元 1 5, 2, 4, 1 % 单元 2 2, 5, 3, 1 % 单元 3 6, 3, 5, 1; % 单元 4 % 材料性质 % 弹性模量 泊松比 厚度 gMaterial = 1e0, 0, 1 ; % 材料 1% 第一类约束条件 % 节点号 自由度号 约束值 gBC1 = 1, 1, 0.0 2, 1, 0.0 4, 1, 0.0 4, 2, 0.0 5, 2, 0.0 6, 2, 0.0 ; % 集中力 % 节点号 自由度号 集中力值 gNF = 1, 2,

5、 -1 ;returnfunction SolveModel% 求解有限元模型% 输入参数：无% 说明：% 该函数求解有限元模型，过程如下% 1. 计算单元刚度矩阵，集成整体刚度矩阵% 2. 计算单元的等效节点力，集成整体节点力向量% 3. 处理约束条件，修改整体刚度矩阵和节点力向量% 4. 求解方程组，得到整体节点位移向量 global gNode gElement gMaterial gBC1 gNF gK gDelta % step1. 定义整体刚度矩阵和节点力向量 node_number,dummy = size( gNode ) ; gK = sparse( node_number

6、* 2, node_number * 2 ) ; f = sparse( node_number * 2, 1 ) ;% step2. 计算单元刚度矩阵，并集成到整体刚度矩阵中 element_number,dummy = size( gElement ) ; for ie=1:1:element_number k = StiffnessMatrix( ie ) ; AssembleStiffnessMatrix( ie, k ) ; end % step3. 把集中力直接集成到整体节点力向量中 nf_number, dummy = size( gNF ) ; for inf=1:1:nf_n

7、umber n = gNF( inf, 1 ) ; d = gNF( inf, 2 ) ; f( (n-1)*2 + d ) = gNF( inf, 3 ) ; end% step4. 处理约束条件，修改刚度矩阵和节点力向量。采用乘大数法 bc_number,dummy = size( gBC1 ) ; for ibc=1:1:bc_number n = gBC1(ibc, 1 ) ; d = gBC1(ibc, 2 ) ; m = (n-1)*2 + d ; f(m) = gBC1(ibc, 3)* gK(m,m) * 1e15 ; gK(m,m) = gK(m,m) * 1e15 ; en

8、d% step 5. 求解方程组，得到节点位移向量 gDelta = gK f ;returnfunction DisplayResults% 显示计算结果% 输入参数：无 global gNode gElement gMaterial gBC1 gNF gK gDelta fprintf( '节点位移n' ) ; fprintf( ' 节点号 x方向位移 y方向位移n' ) ; node_number,dummy = size( gNode ) ; for i=1:node_number fprintf( '%6d %16.8e %16.8en'

9、;,. i, gDelta(i-1)*2+1), gDelta(i-1)*2+2) ; end fprintf( 'nn单元应力n' ) ; fprintf( ' X-STR Y-STR XY-STRn' ) ; element_number, dummy = size( gElement ) ; for ie = 1:element_number es = ElementStress( ie ) ; fprintf( '单元号%6d %16.8e %16.8e %16.8en', . ie, es(1), es(2), es(3); endre

10、turn function k = StiffnessMatrix( ie )% 计算单元刚度矩阵% 输入参数:% ie - 单元号 global gNode gElement gMaterial k = zeros( 6, 6 ) ; E = gMaterial( gElement(ie, 4), 1 ) ; mu = gMaterial( gElement(ie, 4), 2 ) ; h = gMaterial( gElement(ie, 4), 3 ) ; xi = gNode( gElement( ie, 1 ), 1 ) ; yi = gNode( gElement( ie, 1 )

11、, 2 ) ; xj = gNode( gElement( ie, 2 ), 1 ) ; yj = gNode( gElement( ie, 2 ), 2 ) ; xm = gNode( gElement( ie, 3 ), 1 ) ; ym = gNode( gElement( ie, 3 ), 2 ) ; ai=xj*ym-xm*yj; aj=xm*yi-xi*ym; am=xi*yj-xj*yi; bi=yj-ym; bj=ym-yi; bm=yi-yj; ci=-(xj-xm); cj=-(xm-xi); cm=-(xi-xj); A=(ai+aj+am)/2; B=bi 0 bj 0

12、 bm 0 0 ci 0 cj 0 cm ci bi cj bj cm bm; B=B/2/A; D=1 mu 0 mu 1 0 0 0 (1-mu)/2; D=D*E/(1-mu2); k=transpose(B)*D*B*h*abs(A);returnfunction AssembleStiffnessMatrix( ie, k )% 把单元刚度矩阵集成到整体刚度矩阵% 输入参数:% ie - 单元号% k - 单元刚度矩阵 global gElement gK for i=1:1:3 for j=1:1:3 for p=1:1:2 for q =1:1:2 m = (i-1)*2+p ;

13、 n = (j-1)*2+q ; M = (gElement(ie,i)-1)*2+p ; N = (gElement(ie,j)-1)*2+q ; gK(M,N) = gK(M,N) + k(m,n) ; end end end endreturnfunction es = ElementStress( ie )% 计算单元的应力% 输入参数% ie - 节点号% es - 单元应力 global gElement gNode gDelta gMaterial es=zeros(1,6); de=zeros(6,1); for j=1:1:3 de(2*j-1)=gDelta(2*gElem

14、ent(ie,j)-1); de(2*j)=gDelta(2*gElement(ie,j); end E = gMaterial( gElement(ie, 4), 1 ) ; mu = gMaterial( gElement(ie, 4), 2 ) ; h = gMaterial( gElement(ie, 4), 3 ) ; xi = gNode( gElement( ie, 1 ), 1 ) ; yi = gNode( gElement( ie, 1 ), 2 ) ; xj = gNode( gElement( ie, 2 ), 1 ) ; yj = gNode( gElement( i

15、e, 2 ), 2 ) ; xm = gNode( gElement( ie, 3 ), 1 ) ; ym = gNode( gElement( ie, 3 ), 2 ) ; ai=xj*ym-xm*yj; aj=xm*yi-xi*ym; am=xi*yj-xj*yi; bi=yj-ym; bj=ym-yi; bm=yi-yj; ci=-(xj-xm); cj=-(xm-xi); cm=-(xi-xj); A=(ai+aj+am)/2; B=bi 0 bj 0 bm 0 0 ci 0 cj 0 cm ci bi cj bj cm bm; B=B/2/A; D=1 mu 0 mu 1 0 0 0

16、 (1-mu)/2; D=D*E/(1-mu2); S=D*B; es(1:3)=S*de; es(6)=0.5*sqrt(es(1)-es(2)2+4*es(3)2); es(4)=0.5*(es(1)+es(2)+es(6); es(5)=0.5*(es(1)+es(2)-es(6);return3，数据文件：输入数据：gNode = 0.0, 2.0 0.0, 1.0 1.0, 1.0 0.0, 0.0 1.0, 0.0 2.0, 0.0 ; gElement = 3, 1, 2, 1 5, 2, 4, 1 2, 5, 3, 1 6, 3, 5, 1; gMaterial = 1e0,

17、0, 1 ;gBC1 = 1, 1, 0.0 2, 1, 0.0 4, 1, 0.0 4, 2, 0.0 5, 2, 0.0 6, 2, 0.0 ;gNF = 1, 2, -1 ;输出数据：节点位移 节点号 x方向位移 y方向位移 1 -8.79120879e-016 -3.25274725e+000 2 8.79120879e-017 -1.25274725e+000 3 -8.79120879e-002 -3.73626374e-001 4 1.17216117e-016 -8.35164835e-016 5 1.75824176e-001 -2.93040293e-016 6 1.758