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

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

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

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

4、0.0 ; % 節(jié)點 6 % 單元定義 % 節(jié)點1 節(jié)點2 節(jié)點3 材料號 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é)點號 自由度號 約束值 6, 2, 0.0 ; % 集中力 % 節(jié)點號 自由度號 集中力值 gNF = 1, 2, -1 ;returnfunction SolveModel% 求解有限元模型% 輸入?yún)?shù)：無% 說明：% 該函數(shù)

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

6、 ;% step2. 計算單元剛度矩陣，并集成到整體剛度矩陣中 element_number,dummy = size( gElement ) ; for ie=1:1:element_number k = StiffnessMatrix( ie ) ; AssembleStiffnessMatrix( ie, k ) ; end % step3. 把集中力直接集成到整體節(jié)點力向量中 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é)點力向量。采用乘大數(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é)點位移向量 gDelta = gK f ;returnfunction D

8、isplayResults% 顯示計算結(jié)果% 輸入?yún)?shù)：無 global gNode gElement gMaterial gBC1 gNF gK gDelta fprintf( '節(jié)點位移n' ) ; fprintf( ' 節(jié)點號 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( '單元號%6d %16.8e %16.8e %16.8en', . ie, es(1), es(2), es(3); endreturn function k = StiffnessMatrix( ie )% 計算單元剛度矩陣% 輸入?yún)?shù):%

10、 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 ), 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 )% 把單元剛度矩陣集成到整體剛度矩陣% 輸入?yún)?shù):% 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 ; 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 )% 計算單元的應(yīng)力% 輸入?yún)?shù)% ie - 節(jié)點號% 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,