一用三結(jié)點(diǎn)三角形平面單元計(jì)算平面結(jié)構(gòu)的應(yīng)力和位移

1、一：用三結(jié)點(diǎn)三角形平面單元計(jì)算平面結(jié)構(gòu)的應(yīng)力和位移。1，設(shè)計(jì)說(shuō)明書(shū)計(jì)算簡(jiǎn)圖，網(wǎng)格劃分，單元及結(jié)點(diǎn)的編號(hào)如下圖所示。由于結(jié)構(gòu)對(duì)稱，去四分之一結(jié)構(gòu)分析。其中E=2e10pa,mu=0.167,h=1m.開(kāi)始輸入節(jié)點(diǎn)坐標(biāo)進(jìn)行單元定義給出材料性質(zhì)定義有限元模型引入約束條件給定荷載信息計(jì)算單剛形成總剛形成結(jié)點(diǎn)力求解有限元模型引進(jìn)位移約束求解平衡方程顯示結(jié)點(diǎn)位移顯示計(jì)算結(jié)果顯示單元應(yīng)力結(jié)束變量注釋：Node - 節(jié)點(diǎn)定義gElement - 單元定義gMaterial - 材料定義，包括彈性模量，泊松比和厚度gBC1 - 約束條件gNF - 集中力gk-總剛gDelta-結(jié)點(diǎn)位移子程序注釋： PlaneS

2、tructualModel 定義有限元模型 SolveModel 求解有限元模型 DisplayResults 顯示計(jì)算結(jié)果 k = StiffnessMatrix( ie )計(jì)算單元?jiǎng)偠?AssembleStiffnessMatrix( ie, k )形成總剛 es = ElementStress( ie )計(jì)算單元應(yīng)力function exam1% 輸入?yún)?shù)： 無(wú)% 輸出結(jié)果： 節(jié)點(diǎn)位移和單元應(yīng)力PlaneStructualModel ; % 定義有限元模型 SolveModel ; % 求解有限元模型 DisplayResults ; % 顯示計(jì)算結(jié)果return ;function P

3、laneStructualModel% 定義平面結(jié)構(gòu)的有限元模型% 輸入?yún)?shù)：無(wú)% 說(shuō)明：% 該函數(shù)定義平面結(jié)構(gòu)的有限元模型數(shù)據(jù)：% gNode - 節(jié)點(diǎn)定義% gElement - 單元定義% gMaterial - 材料定義，包括彈性模量，泊松比和厚度% gBC1 - 約束條件% gNF - 集中力 global gNode gElement gMaterial gBC1 gNF % 節(jié)點(diǎn)坐標(biāo) % x y gNode = 0.0, 2.0 % 節(jié)點(diǎn) 1 0.0, 1.0 % 節(jié)點(diǎn) 2 1.0, 1.0 % 節(jié)點(diǎn) 3 0.0, 0.0 % 節(jié)點(diǎn) 4 1.0, 0.0 % 節(jié)點(diǎn) 5 2.0,

4、0.0 ; % 節(jié)點(diǎn) 6 % 單元定義 % 節(jié)點(diǎn)1 節(jié)點(diǎn)2 節(jié)點(diǎn)3 材料號(hào) gElement = 3, 1, 2, 1 % 單元 1 5, 2, 4, 1 % 單元 2 2, 5, 3, 1 % 單元 3 6, 3, 5, 1; % 單元 4 % 材料性質(zhì) % 彈性模量 泊松比 厚度 gMaterial = 1e0, 0, 1 ; % 材料 1% 第一類約束條件 % 節(jié)點(diǎn)號(hào) 自由度號(hào) 約束值 6, 2, 0.0 ; % 集中力 % 節(jié)點(diǎn)號(hào) 自由度號(hào) 集中力值 gNF = 1, 2, -1 ;returnfunction SolveModel% 求解有限元模型% 輸入?yún)?shù)：無(wú)% 說(shuō)明：% 該函數(shù)

5、求解有限元模型，過(guò)程如下% 1. 計(jì)算單元?jiǎng)偠染仃?，集成整體剛度矩陣% 2. 計(jì)算單元的等效節(jié)點(diǎn)力，集成整體節(jié)點(diǎn)力向量% 3. 處理約束條件，修改整體剛度矩陣和節(jié)點(diǎn)力向量% 4. 求解方程組，得到整體節(jié)點(diǎn)位移向量 global gNode gElement gMaterial gBC1 gNF gK gDelta % step1. 定義整體剛度矩陣和節(jié)點(diǎn)力向量 node_number,dummy = size( gNode ) ; gK = sparse( node_number * 2, node_number * 2 ) ; f = sparse( node_number * 2, 1 )

6、 ;% step2. 計(jì)算單元?jiǎng)偠染仃嚕⒓傻秸w剛度矩陣中 element_number,dummy = size( gElement ) ; for ie=1:1:element_number k = StiffnessMatrix( ie ) ; AssembleStiffnessMatrix( ie, k ) ; end % step3. 把集中力直接集成到整體節(jié)點(diǎn)力向量中 nf_number, dummy = size( gNF ) ; for inf=1:1:nf_number n = gNF( inf, 1 ) ; d = gNF( inf, 2 ) ; f( (n-1)*2

7、+ d ) = gNF( inf, 3 ) ; end% step4. 處理約束條件，修改剛度矩陣和節(jié)點(diǎn)力向量。采用乘大數(shù)法 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 ; end% step 5. 求解方程組，得到節(jié)點(diǎn)位移向量 gDelta = gK f ;returnfunction D

8、isplayResults% 顯示計(jì)算結(jié)果% 輸入?yún)?shù)：無(wú) global gNode gElement gMaterial gBC1 gNF gK gDelta fprintf( '節(jié)點(diǎn)位移n' ) ; fprintf( ' 節(jié)點(diǎn)號(hào) x方向位移 y方向位移n' ) ; node_number,dummy = size( gNode ) ; for i=1:node_number fprintf( '%6d %16.8e %16.8en',. i, gDelta(i-1)*2+1), gDelta(i-1)*2+2) ; end fprintf(

9、'nn單元應(yīng)力n' ) ; fprintf( ' X-STR Y-STR XY-STRn' ) ; element_number, dummy = size( gElement ) ; for ie = 1:element_number es = ElementStress( ie ) ; fprintf( '單元號(hào)%6d %16.8e %16.8e %16.8en', . ie, es(1), es(2), es(3); endreturn function k = StiffnessMatrix( ie )% 計(jì)算單元?jiǎng)偠染仃? 輸入?yún)?shù):%

10、 ie - 單元號(hào) 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 ), 2 ) ; xj = gNode( gElement( ie, 2 ), 1 ) ; yj = gNode(

11、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 bm 0 0 ci 0 cj 0 cm ci bi cj bj cm bm; B=B/2/A; D=1 mu 0

12、 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 )% 把單元?jiǎng)偠染仃嚰傻秸w剛度矩陣% 輸入?yún)?shù):% ie - 單元號(hào)% k - 單元?jiǎng)偠染仃?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 ; n = (j-1)*2+q ; M = (gElement(ie,i)-1)*2+p ; N = (gEleme

13、nt(ie,j)-1)*2+q ; gK(M,N) = gK(M,N) + k(m,n) ; end end end endreturnfunction es = ElementStress( ie )% 計(jì)算單元的應(yīng)力% 輸入?yún)?shù)% ie - 節(jié)點(diǎn)號(hào)% es - 單元應(yīng)力 global gElement gNode gDelta gMaterial es=zeros(1,6); de=zeros(6,1); for j=1:1:3 de(2*j-1)=gDelta(2*gElement(ie,j)-1); de(2*j)=gDelta(2*gElement(ie,j); end E = gM

14、aterial( 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( ie, 2 ), 2 ) ; xm = gNode( gElement( ie, 3 ), 1 ) ; ym = g

15、Node( 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 (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，數(shù)據(jù)文件：輸入數(shù)據(jù)：gNode = 0.0, 2.0 0.0, 1.0 1.0, 1.0 0.0, 0.0 1.0, 0.0 2.0,