《有限元分析》課程作業(yè)

《有限元分析》課程作業(yè)《有限元分析》課程作業(yè)《有限元分析》課程作業(yè)《有限元分析》課程作業(yè)編制僅供參考審核批準(zhǔn)生效日期地址：電話：傳真:郵編:《有限元分析》課程作業(yè)任課教師：徐亞蘭學(xué)生姓名：陳新杰學(xué)號：班級：1304012時間：2016-01-05

一、問題描述及分析問題：如圖1所示，有一矩形平板，在右側(cè)受到P=10KN/m的分布力，材料常數(shù)為：彈性模量；泊松比；板的厚度為t=；試按平面應(yīng)力問題利用三角形與矩形單元分別計算各個節(jié)點位移及支座反力。P=10KN/m1m1P=10KN/m1m1m圖1平面矩形結(jié)構(gòu)的有限元分析分析：使用兩種方案：一、基于3節(jié)點三角形單元的有限元建模，將矩形劃分為兩個3節(jié)點三角形單元；二、基于4節(jié)點矩形單元的有限元建模，使用一個4節(jié)點矩形單元。利用MATLAB軟件計算出各要求量，再將兩種方案的計算結(jié)果進行比較、分析、得出結(jié)論。二、有限元建模及分析1、基于3節(jié)點三角形單元的有限元建模及分析（1）結(jié)構(gòu)的離散化與編號如圖2所示，將平面矩形結(jié)構(gòu)分為兩個3節(jié)點三角形單元。單元①三個節(jié)點的編號為1，2，4，單元②三個節(jié)點的編號為3，4，2，各個節(jié)點的位置坐標(biāo)為，各個節(jié)點的位移（分別沿方向和方向）為。X②①y1432X②①y1432圖2方案一：使用兩個3節(jié)點三角形單元（2）各單元的剛度矩陣及剛度方程a.單元的幾何和節(jié)點描述單元①有6個節(jié)點位移自由度（DOF）。將所有節(jié)點上的位移組成一個列陣，記作；同樣，將所有節(jié)點上的各個力也組成一個列陣，記作，則有同理，對于單元②，有b.單元的位移場描述對于單元①，設(shè)位移函數(shù)(1-1)由節(jié)點條件，在處，有(1-2)將式(1-1)代入節(jié)點條件式(1-2)中，可求出式（1-1）中待定系數(shù)，即(1-3)(1-4)(1-5)(1-6)(1-7)(1-8)在式(1-3)~式(1-8)中(1-9)(1-10)上式中的符號(1,2,3)表示下標(biāo)輪換，如同時更換。將單元①各節(jié)點的位置坐標(biāo)代入得將系數(shù)式(1-3)~式(1-8)式代入(1-1)中，重寫位移函數(shù)，并以節(jié)點位移的形式進行表示，有(1-11)令，則有形狀函數(shù)矩陣(1-12)位移函數(shù)式(1-11)寫成矩陣形式，有(1-13)對于單元②，過程同上，有形狀函數(shù)矩陣(1-14)位移函數(shù)(1-15)c.單元的應(yīng)變場描述對于單元①，應(yīng)變函數(shù)(1-16)其中幾何矩陣(1-17)對于單元②，應(yīng)變方程(1-18)其中幾何矩陣(1-19)d.單元的應(yīng)力場描述(1-20)其中，彈性系數(shù)矩陣D(1)為(1-21)(1-22)其中應(yīng)力函數(shù)矩陣S(1)為(1-23)應(yīng)力方程為(1-24)對于單元②，過程同上。彈性系數(shù)矩陣D(2)為(1-25)應(yīng)力函數(shù)矩陣S(2)為(1-26)應(yīng)力方程為(1-27)e.單元的勢能表達是單元剛度矩陣，即(1-28)其中薄板厚度。將式(1-17)、式(1-21)代入式(1-29)，得到單元①的剛度陣計算得同理，得到單元②的剛度陣為將兩個單元按節(jié)點位移所對應(yīng)的位置進行組裝，得到總體剛度矩陣為節(jié)點力系統(tǒng)的勢能(計算結(jié)果在下面呈現(xiàn))(4)邊界條件的處理及方程求解邊界條件為。因此，將針對節(jié)點2和節(jié)點3的位移求解，節(jié)點2和節(jié)點3對應(yīng)總體剛度陣KK中的第3行到第6行、第3列到第6列，則需從KK中提出，置給k，然后生成對應(yīng)的載荷列陣p，再采用高斯消去法進行求解。>>k=KK(3:6,3:6);>>p=[500;0;500;0];>>u=k\pu=[將列排成了行]再計算支反力。在得到整個結(jié)構(gòu)的節(jié)點位移后，由原整體剛度方程就可以計算出對應(yīng)的支反力；先將上面得到的位移結(jié)果與邊界條件的節(jié)點位移進行組合，得到整體的位移列陣U，再代回原整體剛度方程，計算出所有的節(jié)點力，按照位置關(guān)系找出對應(yīng)的支反力。>>U=[0;0;u;0;0]U=0000[將列排成了行]>>P=KK*UP=-50050005000-500[將列排成了行]所以，節(jié)點1的支反力為，節(jié)點2的支反力為。根據(jù)已求得的位移和支反力計算系統(tǒng)的勢能。>>A=*U'*KK*U-P'*UA=(5）結(jié)果分析上述支反力計算結(jié)果滿足靜力平衡，驗證了以上求解過程及MATLAB算法的正確性。2、基于四節(jié)點四邊形單元的有限元建模及分析（1）結(jié)構(gòu)的離散化與編號如圖3所示一個4節(jié)點矩形單元，單元的節(jié)點位移共有8個自由度（DOF）。節(jié)點編號為1,2,3,4，各自的位置坐標(biāo)為，各個節(jié)點的位移（分別沿方向和方向）為。Xy①4321Xy①4321圖3方案二：使用一個4節(jié)點矩形單元（2）局部坐標(biāo)系下單元的描述a.單元的幾何和節(jié)點描述采用無量綱坐標(biāo)其中。則單元四個節(jié)點的幾何位置為將所有節(jié)點上的位移組成一個列陣，記作q；同樣，將所有節(jié)點上的各個力也組成一個列陣，記作F，則有b.單元的位移場描述設(shè)位移函數(shù)為由節(jié)點條件，在處，有將位移試函數(shù)代入節(jié)點條件中，求出待定系數(shù)和，再代入位移函數(shù)中，整理后得其中如以無量綱坐標(biāo)來表達，可寫成將帶入上式得到形狀函數(shù)矩陣寫成矩陣形式，有c.單元的應(yīng)變場描述單元應(yīng)變?yōu)槠渲袔缀尉仃嚍閐.單元的應(yīng)力場描述應(yīng)力表達式為其中，應(yīng)力函數(shù)矩陣。e.單元的勢能表達以上已將單元的三大基本變量用基于節(jié)點的位移列陣來表示；將其代入單元勢能表達式中，有，其中為4節(jié)點矩形單元的剛度矩陣，即其中，t為薄板的厚度，，上式的各個字塊矩陣為f.單元剛度陣及剛度方程單元剛度陣在上面已經(jīng)列出。將單元的勢能對節(jié)點位移取一階極值，可得到單元的剛度方程(4)邊界條件的處理及方程求解處理方法與3節(jié)點三角形單元一致，利用上述求解程序具有的可移植性，簡化了求解過程。>>k=K(3:6,3:6);>>p=[500;0;500;0];>>u=k\pu=*[將列排成了行]再計算支反力。同樣注意按照位置關(guān)系找出對應(yīng)的支反力。>>U=[0;0;u;0;0]U=*0000[將列排成了行]>>P=K*UP=-50050005000-500[將列排成了行]所以，節(jié)點1的支反力為，節(jié)點2的支反力為。根據(jù)已求得的位移和支反力計算系統(tǒng)的勢能。>>A=*U'*K*U-P'*UA=(5）結(jié)果分析基于4節(jié)點矩形單元計算出的勢能小于基于3節(jié)點三角形單元計算出的結(jié)果，若將該系統(tǒng)分為更多的單元，計算精度也將提高。3.兩種方案的比較與分析從以上計算可以看出，用三角形單元計算時，由于形函數(shù)是完全一次式，因而其應(yīng)變場在單元內(nèi)均為常數(shù)；而對于四邊形單元，其形函數(shù)帶有二次式，計算得到的應(yīng)變場和應(yīng)力場都是坐標(biāo)的一次函數(shù)，但不是完全的一次函數(shù)，對提高精度有一定作用；根據(jù)最小勢能原理，勢能越小，則整體計算精度越高，比較兩種單元計算得到的系統(tǒng)勢能，可看出，在相同的節(jié)點自由度情況下，矩形單元的計算精度要比三角形單元高。三、基于MATLAB的編程實現(xiàn)1.基于3節(jié)點三角形單元的有限元編程實現(xiàn)（1）程序編寫說明Triangle2D3Node_Stiffness(E,NU,t,xi,yi,xj,yj,xm,ym,ID)該函數(shù)計算單元的剛度矩陣，輸入模量E，泊松比NU，厚度t，三個節(jié)點i，j，m，平面問題性質(zhì)指示參數(shù)（1為平面應(yīng)力，2為平面應(yīng)變），輸出單元剛度矩陣k（66）。Triangle2D3Node_Assemble(KK,k,i,j,m)該函數(shù)進行單元剛度矩陣的組裝，輸入單元剛度矩陣k，單元的節(jié)點編號i，j，m，輸出整體剛度矩陣KK。Triangle2D3Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,u,ID)該函數(shù)計算單元的應(yīng)力，輸入彈性模量E，泊松比NU，厚度t，三個節(jié)點i，j，m，平面問題性質(zhì)指示參數(shù)（1為平面應(yīng)力，2為平面應(yīng)變），單元的位移列陣u(61)，輸出單元的應(yīng)力，由于它為常應(yīng)力單元，則單元的應(yīng)力分量為Sx，Sy，Sxy。（2）程序清單%%%%%Triangle2D3Node%%%begin%%%%%%%functionk=Triangle2D3Node_Stiffness(E,NU,t,xi,yi,xj,yj,xm,ym,ID)%該函數(shù)計算單元的剛度矩陣%輸入彈性模量E、泊松比NU和厚度t%輸入3個節(jié)點i，j，m的坐標(biāo)xi,yi,xj,yj,xm,ym%輸入平面問題性質(zhì)指示參數(shù)ID（1位平面應(yīng)力，2為平面應(yīng)變）%輸入單元剛度矩陣k（6*6)%-----------------------------------------A=(xi*(yj-ym)+xj(ym-yi)+xm*(yi-yj))/2;betai=yj-ym;betaj=ym-yi;betam=yi-yj;gammai=xm-xj;gammaj=xi-xm;gammam=xj-xi;B=[betai0betaj0betam0;0gammai0gammaj0gammam;gammaibetaigammajbetajgammambetam]/(2*A);ifID==1D=(E/(1-NU*NU))*[1NU0;NU10;00(1-NU)/2];elseifID==2D=(E/(1+NU)/(1-2*NU))*[1-NUNU0;NU1-NU0;00(1-2*NU)/2];endk=t*A*B'*D*B;end%%%%%%%%%%%%%%%%%%%%%%%%%%%%functionz=Triangle2D3Node_Assemble(KK,k,i,j,m)%該函數(shù)進行單元剛度矩陣的組裝%輸入單元剛度矩陣k%輸入單元的節(jié)點編號i，j，m%輸入整體剛度矩陣KKrf%------------------------------------------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;forn1=1:6forn2=1:6KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);endendz=KK;%%%%%%%%%%%%%%%%%%%%%%%%%%%%functionstress=Triangle2D3Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,u,ID)%該函數(shù)計算單元的應(yīng)力%輸入彈性模量E、泊松比NU和厚度t%輸入平面問題性質(zhì)指示參數(shù)ID（1位平面應(yīng)力，2為平面應(yīng)變），單元的位移列陣u（6*1）%輸出單元的應(yīng)力stress（3*1），由于它為常應(yīng)力單元，則單元的應(yīng)力分量Sx，Sy，Sxy%--------------------------------------------------A=(xi*xj(ym-yi)+xm*(yi-yj))/2;betai=yj-ym;betaj=ym-yi;betam=yi-yj;gammai=xm-xj;gammaj=xi-xm;gammam=xj-xi;B=[batai0bataj0betam0;0gammai0gammaj0gammam;gammaibetaigammajbetajgammambetam]/(2*A);ifID==1D=(E/(1-NU*NU))*[1NU0;NU10;00(1-NU)/2];elseifID==2D=(E/(1+NU)/(1-2*NU))*[1-NUNU0;NU1-NU0;00(1-2*NU)/2];endstress=D*B*u;%%%%%%%%%%%%%%%%%%%%%%%%%%>>E=1E7;>>NU=1/3;>>t=;>>c=Triangle2D3Node_Stiffness(E,NU,t,0,0,1,0,0,1,2)>>CC=zeros(8,8);>>CC=Triangle2D3Node_Assemble(KK,k1,1,2,4)；>>CC=Triangle2D3Node_Assemble(KK,k1,3,4,2)>>k1=Triangle2D3Node_Stiffness(E,NU,t,0,0,1,0,0,1,1)>>KK=zeros(8,8);>>KK=Triangle2D3Node_Assemble(KK,k1,1,2,4)；>>KK=Triangle2D3Node_Assemble(KK,k1,3,4,2)>>k=KK(3:6,3:6);>>p=[500;0;500;0];>>u=k\p>>U=[0;0;u;0;0]>>P=KK*U>>A=*U'*KK*U-P'*U（3)計算結(jié)果①應(yīng)變CC=+05*0000000000000000②位移U=0000③節(jié)點力P=00其中，節(jié)點1的支反力為，節(jié)點2的支反力為。④勢能A=⑤單元剛度陣KK=+05*00000000000000002.基于四節(jié)點四邊形單元的有限元建模及分析（1）程序編寫說明Quad2D4Node_Stiffness(E,NU,h,xi,yi,xj,yj,xm,ym,xp,yp,ID)該函數(shù)計算單元的剛度矩陣，輸入模量E，泊松比NU，厚度h，4個節(jié)點i，j，m，p，平面問題性質(zhì)指示參數(shù)ID（1為平面應(yīng)力，2為平面應(yīng)變），輸出單元剛度矩陣k（88）。Quad2D4Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,xp,yp,u,ID)該函數(shù)計算單元的應(yīng)力，輸入彈性模量E，泊松比NU，厚度h，4個節(jié)點i，j，m，p，平面問題性質(zhì)指示參數(shù)ID（1為平面應(yīng)力，2為平面應(yīng)變），單元的位移列陣u(８1)，輸出單元的應(yīng)力，由于它為常應(yīng)力單元，則單元的應(yīng)力分量為Sx，Sy，Sxy。（2）程序清單%%%%%%%Quad2D4Node%%%begin%%%%%%%%functionk=Quad2D4Node_Stiffness(E,NU,h,xi,yi,xj,yj,xm,ym,xp,yp,ID)%該函數(shù)計算單元的剛度矩陣%輸入彈性模量E、泊松比NU和厚度h%輸入4個節(jié)點i，j，m，p的坐標(biāo)xi,yi,xj,yj,xm,ym,xp,yp%輸入平面問題性質(zhì)指示參數(shù)ID（1位平面應(yīng)力，2為平面應(yīng)變）%輸入單元剛度矩陣k（8*8)%-----------------------------------------symsst;a=(yi*(s-1)+yj*(-1-s)+ym*(1+s)+yp*(1-s))/4;b=(yi*(t-1)+yj*(1-t)+ym*(1+t)+yp*(-1-t))/4;c=(xi*(t-1)+xj*(1-t)+xm*(1+t)+xp*(-1-t))/4;d=(xi*(s-1)+xj*(-1-s)+xm*(1+s)+xp*(1-s))/4;B1=[a*(t-1)/4-b*(s-1)/40;0c*(s-1)/4-d*(t-1)/4;c*(-1+s)/4-d*(t-1)/4a*(t-1)/4-b*(s-1)/4];B2=[a*(1-t)/4-b*(-1-s)/40;0c*(-1-s)/4-d*(1-t)/4;c*(-1-s)/4-d*(1-t)/4a*(1-t)/4-b*(-1-s)/4];B3=[a*(t+1)/4-b*(s+1)/40;0c*(s+1)/4-d*(t+1)/4;c*(s+1)/4-d*(t+1)/4a*(t+1)/4-b*(s+1)/4];B4=[a*(-t-1)/4-b*(1-s)/40;0c*(1-s)/4-d*(-t-1)/4;c*(1-s)/4-d*(-t-1)/4a*(-t-1)/4-b*(1-s)/4];Bfirst=[B1B2B3B4];Jfirst=[01-tt-ss-1;t-10s+1-s-t;s-t-s-10t+1;1-ss+t-t-10];J=[xixjxmxp]*[yi;yj;ym;yp]/8;B=Bfirst/J;ifID==1D=(E/(1-NU*NU))*[1NU0;NU10;00(1-NU)/2];elseifID==2D=(E/(1+NU)/(1-2*NU))*[1-NUNU0;NU1-NU0;00(1-2*NU)/2];endBD=J*transpose(B)*D*B;r=int(int(BD,t,-1,1),s,-1,1);z=h*r;k=double(z);end%%%%%%%%%%%%%%%%%%%%%%%%%%%%functionstress=Quad2D4Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,xp,yp,u,ID)%該函數(shù)計算單元的應(yīng)力%輸入彈性模量E、泊松比NU和厚度h%輸入平面問題性質(zhì)指示參數(shù)ID（1位平面應(yīng)力，2為平面應(yīng)變）%輸入單元的位移列陣u（8*1）%輸出單元的應(yīng)力stress（3*1），由于它為常應(yīng)力單元，則單元的應(yīng)力分量Sx，Sy，Sxy%--------------------------------------------------symst;a=(yi*(s-1)+yi*(-1-s)+ym*(1+S)+yp*(1-s))/4;b=(yi*(t-1)+yj*(1-t)+ym*(1+t)+yp*(-1-t))/4c=(xi*(t-1)+xj*(1-t)+xm*(1+t)+xp*(-1-t))/4d=(xi*(s-1)+xj*(-1-s)+xm*(1+s)+xp*(1-s))/4B1=[a*(t-1)/4-b*(s-1)/40;0c*(s-1)/4-d*(t-1)/4;c*(-1+s)/4-d*(t-1)/4a*(t-1)/4-b*(s-1)/4];B2=[a*(1-t)/4-b*(-1-s)/40;0c*(-1-S)/4-d*(1-t)/4;c*(-1-s)/4-d*(1-t)/4a*(1-t)/4-b*(-1-s)/4];B3=[a*(t+1)/4-b*(s+1)/40;0c*(s+1)/4-d*(t+1)/4;c*(s+1)/4-d*(t+1)/4a*(t+1)/4-b*(s+1)/4];B4=[a*(-t-1)/4-b*(1-s)/40;0c*(1-s)/4-d*(-t-1)/4;c*(1-s)/4-d*(-t-1)/4a*(-t-1)/4-b*(1-s)/4];Bfirst=[B1B2B3B4];Jfirst=[01-tt-ss-1;t-10s+1-s-t;s-t-s-10t+1;1-ss+t-t-10];J=[xixjxmxp]*[yi;yj;ym;yp]/8;B=Bfirst/J;ifID==1D=(E/(1-NU*NU))*[1NU0;NU10;00(1-NU)/2];elseifID==2D=(E/(1+NU)/(1-2*NU))*[1-NUNU0;NU1-NU0;00(1-2*NU)/2];endstr1=D*B*u;str2=subs(str1,{s,t},{0,0});stress=double(str2);end>>Ｃ=Quad2D4Node_Stiffness(E,NU,h,0,0,1,0,1,1,0,1，２)>>K=Quad