基于MTLAB的3節(jié)點三角形單元求解平面彈性問題

1、基于MTLAB的3節(jié)點三角形單元求解平面彈性問題一 引言 本程序相比其他參考資料的程序，在操作方面得到大大簡化，只需要一次輸入相關力學參量、單元總數(shù)、節(jié)點總數(shù)、節(jié)點坐標、單元節(jié)點編號，利用此程序就能夠很快得到所有單元的單元剛度矩陣，以及系統(tǒng)的整體剛度矩陣。為了便于后面計算單元的應力，在計算單元剛度矩陣時，MATLAB程序輸出了彈性系數(shù)矩陣D和幾何矩陣B。 此外，本程序在引入位移邊界條件求解剛度方程時，采用了化零置1法，避免了手工化簡整體剛度矩陣的麻煩。當然，這在一定程度上加重了計算機的負擔，增加了計算機的計算時間，但是這個延長的時間相對于手工化簡剛度矩陣所需的時間來說，是微不足道的。程序需要提

2、前輸入的參數(shù)有：（1） 節(jié)點坐標node_gcoord（2） 單元節(jié)點編號element_node_number（3） 節(jié)點總數(shù)node_number（4） 單元總數(shù)element_number（5） 彈性模量E（6） 泊松比NU（7） 厚度t（8） 平面問題性質指示參數(shù)ID，其中1為平面應力問題，2為平面應變問題該程序主要包括以下部分：1. 單元剛度矩陣程序：用于計算各單元的單元剛度矩陣；2. 整體剛度矩陣程序：用于組裝整體剛度矩陣；3. 位移及支反力求解程序：用于求解未知位移和未知支反力4. 單元應力計算程序：用于計算各單元的單元應力二 問題描述為了驗證程序的正確性，我們選取了參考文獻1中

3、的一個簡單實例進行驗證。問題如下：如圖所示為一矩形薄平板，在右端部受集中力F=100kN的作用，材料常數(shù)為：彈性模量E=1e7Pa，泊松比NU=1/3，板的厚度t=0.1m。試按平面應力問題計算各節(jié)點位移及支座反力。三 問題求解（1） 結構的離散化與編號1） 節(jié)點描述 2) 單元描述 （2） 計算單元剛度矩陣KE(:,:,1) = 1.0e+006 * 0.2813 0 0 0.1875 -0.2813 -0.1875 0 0.0938 0.1875 0 -0.1875 -0.0938 0 0.1875 0.3750 0 -0.3750 -0.1875 0.1875 0 0 1.1250 -0

4、.1875 -1.1250 -0.2813 -0.1875 -0.3750 -0.1875 0.6563 0.3750 -0.1875 -0.0938 -0.1875 -1.1250 0.3750 1.2188KE(:,:,2) = 1.0e+006 * 0.2813 0 0 0.1875 -0.2813 -0.1875 0 0.0938 0.1875 0 -0.1875 -0.0938 0 0.1875 0.3750 0 -0.3750 -0.1875 0.1875 0 0 1.1250 -0.1875 -1.1250 -0.2813 -0.1875 -0.3750 -0.1875 0.65

5、63 0.3750 -0.1875 -0.0938 -0.1875 -1.1250 0.3750 1.2188（3） 組裝整體剛度矩陣KK = 1.0e+006 * 0.6563 0.3750 -0.3750 -0.1875 -0.2813 -0.1875 0 0 0.3750 1.2188 -0.1875 -1.1250 -0.1875 -0.0938 0 0 -0.3750 -0.1875 0.6563 0 0 0.3750 -0.2813 -0.1875 -0.1875 -1.1250 0 1.2188 0.3750 0 -0.1875 -0.0938 -0.2813 -0.1875 0

6、 0.3750 0.6563 0 -0.3750 -0.1875 -0.1875 -0.0938 0.3750 0 0 1.2188 -0.1875 -1.1250 0 0 -0.2813 -0.1875 -0.3750 -0.1875 0.6563 0.3750 0 0 -0.1875 -0.0938 -0.1875 -1.1250 0.3750 1.2188（4） 求解位移及支反力delta = 0.1877 -0.8992 -0.1497 -0.8422 0 0 0 0F = 1.0e+004 * 0 -0.5000 0 -0.5000 -2.0000 -0.0702 2.0000 1.

7、0702（5） 計算單元應力stress(:,:,1) = 1.0e+005 * -0.8419 -0.2806 -1.5791stress(:,:,2) = 1.0e+004 * 8.4187 -2.8953 -4.2094四 結果分析可以看出，計算得到的單元1的應力分量為x =-84190Pa,y=-28060Pa,xy=-157910Pa,單元2的應力分量為x =84187Pa,y=-28953Pa,xy=-42094Pa.此結果與參考文獻1上給出的結果完全吻合。五 程序優(yōu)點 很多參考資料中的程序在求解單元剛度矩陣時，往往需要人工進行重復操作，例如一個系統(tǒng)具有n個單元，那么就需要人為的運

8、行n次程序，并且每次都需要手動更改所需參數(shù)來求解單元剛度矩陣。在組裝整體剛度矩陣時，也需要人為的運行n次程序，并且每次都需要手動更改所需參數(shù)來組裝整體剛度矩陣。這樣一來，使得操作極為不便，浪費大量時間。本程序相比其他參考資料的程序，在操作方面得到大大簡化，只需要一次輸入相關力學參量、單元總數(shù)、節(jié)點總數(shù)、節(jié)點坐標、單元節(jié)點編號，利用此程序就能夠很快得到所有單元的單元剛度矩陣，以及系統(tǒng)的整體剛度矩陣。為了便于后面計算單元的應力，在計算單元剛度矩陣時，MATLAB程序輸出了彈性系數(shù)矩陣D和幾何矩陣B。 此外，本程序在引入位移邊界條件求解剛度方程時，采用了化零置1法，避免了手工化簡整體剛度矩陣的麻煩。

9、當然，這在一定程度上加重了計算機的負擔，增加了計算機的計算時間，但是這個延長的時間相對于手工化簡剛度矩陣所需的時間來說，是微不足道的。附錄1. 計算單元節(jié)點坐標%triangle2d3nodeeng%function element_node_gcoord=triangle2d3nodeeng(node_gcoord,element_node_number,element_number)%-%element_node_gcoord=triangle2d3nodeeng(node_gcoord,element_node_number)%該函數(shù)用于生成單元的節(jié)點坐標 %element_node_g

10、coord返回單元的節(jié)點坐標%node_gcoord節(jié)點坐標%element_node_number單元節(jié)點編號%-for loopi=1:element_number element_node_gcoord(loopi,1)=node_gcoord(element_node_number(loopi,1),1); element_node_gcoord(loopi,2)=node_gcoord(element_node_number(loopi,1),2); element_node_gcoord(loopi,3)=node_gcoord(element_node_number(loopi,

11、2),1); element_node_gcoord(loopi,4)=node_gcoord(element_node_number(loopi,2),2); element_node_gcoord(loopi,5)=node_gcoord(element_node_number(loopi,3),1); element_node_gcoord(loopi,6)=node_gcoord(element_node_number(loopi,3),2); end%end%2. 計算單元剛度矩陣%triangle2d3nodestiffness%function KE,B,D,A=triangle

12、2d3nodestiffness(E,NU,t,element_number,element_node_gcoord,ID)%-%KE,B,D=triangle2d3nodestiffness(E,NU,t,element_number,element_node_gcoord,ID)%該函數(shù)用于計算單元的剛度矩陣%KE返回單元剛度矩陣%B返回各單元的幾何矩陣%D返回平面彈性系數(shù)矩陣%輸入彈性模量E，泊松比NU，厚度t%輸入單元節(jié)點坐標文件element_node_gcooord%輸入單元總數(shù)element_number%輸入平面問題性質指示參數(shù)ID，其中1為平面應力問題，2為平面應變問題%-f

13、or loopi=1:element_numberxi=element_node_gcoord(loopi,1);yi=element_node_gcoord(loopi,2);xj=element_node_gcoord(loopi,3);yj=element_node_gcoord(loopi,4);xm=element_node_gcoord(loopi,5);ym=element_node_gcoord(loopi,6);A(loopi)=1/2*det(1 xi yi;1 xj yj;1 xm ym); %計算各三角形單元的面積b1=yj-ym;b2=ym-yi;b3=yi-yj;c

14、1=xm-xj;c2=xi-xm;c3=xj-xi;%計算各單元的幾何矩陣BB(:,:,loopi)=1/(2*A(loopi)*b1 0 b2 0 b3 0;0 c1 0 c2 0 c3;c1 b1 c2 b2 c3 b3;if ID=1 E=E; NU=NU;elseif ID=2 E=E/(1-NU2); NU=NU/(1-NU);end%形成平面問題的彈性系數(shù)矩陣D D=E/(1-NU2)*1 NU 0;NU 1 0;0 0 (1-NU)/2;%形成單元的剛度矩陣KE(:,:,loopi)=B(:,:,loopi)'*D*B(:,:,loopi)*A(loopi)*t;end%

15、end%3. 形成整體剛度矩陣%triangle2d3nodeassembly%function KK=triangle2d3nodeassembly(KE,node_number,element_number,element_node_number)%-%KK=triangle2d3nodeassembly(KE,node_number,element_number,element_node_number)%該函數(shù)用于組裝整體剛度矩陣%KK返回系統(tǒng)的整體剛度矩陣%KE單元剛度矩陣%node_numbre節(jié)點總數(shù)%element_number單元總數(shù)%element_node_number單元

16、的節(jié)點編號%-KK=zeros(2*node_number,2*node_number);for loopi=1:element_numberi=element_node_number(loopi,1);j=element_node_number(loopi,2);m=element_node_number(loopi,3);dof(1)=2*i-1;dof(2)=2*i;dof(3)=2*j-1;dof(4)=2*j;dof(5)=2*m-1;dof(6)=2*m;k=KE(:,:,loopi);for n1=1:6 for n2=1:6 KK(dof(n1),dof(n2)=KK(dof(

17、n1),dof(n2)+k(n1,n2); endendend%end%4. 求解未知位移和未知支反力%triangle2d3noded_f%function delta,F=triangle2d3noded_f(KK,deltas,FS,node_number)%-%該函數(shù)用于引入位移邊界條件化簡整體剛度矩陣和節(jié)點載荷列陣，并且求出未知位移和支反力%delta返回求解出的節(jié)點位移列陣%F返回求解出的節(jié)點力列陣%KK原整體剛度矩陣%deltas節(jié)點位移列陣，里面含有未知節(jié)點位移%FS節(jié)點力列陣，里面含有未知節(jié)點力%node_number節(jié)點總數(shù)%-KK1=KK;for loopi=1:2*no

18、de_number if (deltas(loopi,1)=0) KK1(:,loopi)=0; KK1(loopi,:)=0; KK1(loopi,loopi)=1; FS(loopi,1)=0; endend KKS=KK1;delta=double(KKSFS);F=KK*delta;%end%5. 計算單元應變和單元應力%triangle2d3nodestrainstress%function strain,stress=triangle2d3nodestrainstress(D,B,delta,element_node_number,element_number)%-%strain,stress=triangle2d3nodestrainstress(D,B,delta,element_node_number,element_number)%該函數(shù)用于計算各