紅-計(jì)算力學(xué)框架編程matlab

1、計(jì)算力學(xué)框架編程（）班級(jí)：號(hào)：期：計(jì)力（春）許 映 紅 1 4 2 5 6 42015-06學(xué)日目錄第一章題 目1題目選擇1題目解讀1空間框架1編程對(duì)象3第二章程 序4程序體系4基本設(shè)想4主程序5函數(shù)5程序代碼6主程序代碼6函數(shù)代碼19第三章算 例28基本信息28運(yùn)行結(jié)果29運(yùn)行截圖29其他截圖37結(jié)果文件37第一章題 目1.1 題目選擇編制框架結(jié)構(gòu)在外荷載作用下結(jié)構(gòu)的應(yīng)力、位移等響應(yīng)。1)2)3)4)5)Level 1（ Level 2（ Level 3（Level 4（60%）：平面框架+靜力荷載；70%）：平面框架+動(dòng)力荷載；80%）：空間框架+靜力荷載；90%）：空間框架+動(dòng)力荷載；L

2、evel 5（100%）：空間框架+靜力荷載+動(dòng)力荷載。選擇“Level 3（ 80%）：空間框架+靜力荷載”進(jìn)行編制。1.2 題目解讀空間框架采用空間三維有限元，一般用考慮了橫向剪切變形的 Timoshenko 梁理論進(jìn)行考慮。這種梁?jiǎn)卧扔芯植孔鴺?biāo)又有總體坐標(biāo)。1.2.1 空間框架空間桿單元和平面桿單元的區(qū)別在于：桿件除了可能承受軸力和彎矩的作用外，還可能承受扭矩的作用；而且彎矩可能同時(shí)在兩個(gè)坐標(biāo)面內(nèi)存在?？臻g桿單元的系數(shù)有彈性模量 E 、剪切彈性模量G 、橫截面積 A 、慣性矩 I x 和 I y 、極慣性矩 J 及長度 L 。每個(gè)空間桿單元有 2 個(gè)節(jié)點(diǎn)，并且相對(duì)于總體坐標(biāo)系的 X、Y

3、、Z 軸正向逆時(shí)針的傾斜角分別為xX 、yY 、zZ 。如圖 1.1 所示。圖 1.1 空間桿單元1因此，單元?jiǎng)偠染仃嚳筛鶕?jù)下面的矩陣給出： = 式中，旋轉(zhuǎn)矩陣T 和彈性模量矩陣k定義如下：（1-1） = 1112（1-2）21221 01 00020002001 02 00002 0001 00020002 000 20230000 = （1-2-1）11 002310002001 02 000020002 00 = 0 0 000 0 = 2（1-2-2）122103 0000 3 1 001 02 0 0001000 20000 = （1-2-3）22 230 0000223式中， =

4、， = ，1262126=， =， =，=，=121231232322=。3 = 000 = 式中， C0 0000 0000（1-3）（1-3-1），coxY ，coxZ ，等等。xYxZ很顯然，空間桿單元有 12 個(gè)度每個(gè)節(jié)點(diǎn)有 6 個(gè)度（3 個(gè)線位移和 3 個(gè)角位移）。因此，對(duì)一個(gè)有 n 個(gè)節(jié)點(diǎn)的空間框架結(jié)構(gòu)而言，其整體剛度矩陣 K 是 6n6n 的（因?yàn)槊總€(gè)節(jié)點(diǎn)有 6 個(gè)度）。一旦得到整體剛度矩陣 K ，就可以列出以下方程組：2（1-4） = 式中，U 是總體節(jié)點(diǎn)位移矢量，F(xiàn) 是總體節(jié)點(diǎn)載荷矢量。在這一步中，邊界條件被手動(dòng)賦值給矢量U 和 F 。然后用分解和消去法求解矩陣方程（1-4）

5、。一旦求出未知的位移和支反力，就可以用下式求出每個(gè)單元的節(jié)點(diǎn)力矢量 f ： = （1-5）式中， f 是 61 的節(jié)點(diǎn)力矢量、u 是 121 的單元節(jié)點(diǎn)位移矢量。單元位移矢量u的第 1 個(gè)、第 2 個(gè)和第 3 個(gè)元素是第一個(gè)節(jié)點(diǎn)的位移分量，第 4 個(gè)、第 5 個(gè)和第 6 個(gè)元素是第一個(gè)節(jié)點(diǎn)的轉(zhuǎn)角分量。單元位移矢量u的第7 個(gè)、第 8 個(gè)和第 9 個(gè)元素是第二個(gè)節(jié)點(diǎn)的位移分量，第 10 個(gè)、第 11 個(gè)和第 12個(gè)元素是第二個(gè)節(jié)點(diǎn)的轉(zhuǎn)角分量。1.2.2 編程對(duì)象本程序采用進(jìn)行編程，主要針對(duì)的滿足以下要求的空間框架進(jìn)行編制：1.2.3.4.5.水平為 X 層，豎向?yàn)?Y 層，橫向?yàn)?Z 層；水平桿

6、長度為 L（m），豎桿長度為 H（m），橫桿長度為 W（m）；桿件截面面積 A（m2），截面慣性矩 Iy（m4）、Iz（m4）,扭轉(zhuǎn)慣性矩J（m4）；桿件材料彈性模量 E（MPa），剪切模量 E（MPa）；頂層節(jié)點(diǎn)受 Y 向荷載 P（kN），底層節(jié)點(diǎn)固結(jié)。其中部分函數(shù)可適用于任意形式空間框架計(jì)算。3第二章程 序2.1 程序體系2.1.1 基本設(shè)想本程序需要達(dá)到的目標(biāo)有：1.模型基本參數(shù)輸入、節(jié)點(diǎn)、單元、荷載施加、約束施加以及模型的三維顯示；對(duì)框架的節(jié)點(diǎn)及單元如圖 2.13 所示：1)節(jié)點(diǎn)以平行于 XY 面為，再以平行于 X 軸為， X 向增大，Y 向增大，Z 向增大。圖 2.1 節(jié)點(diǎn)示意圖2)

7、單元先編平行于 Y 軸的豎桿，以平行于 XY 面為向增大，Y 向增大，Z 向增大。，再編較低層， X圖 2.2 豎桿單元示意圖3)單元后編垂直于 Y 軸的水平桿和橫桿，以平行于 XY 面為，再編較4低層， X 向增大，Y 向增大，Z 向增大。圖 2.3 水平桿、橫桿單元示意圖2. 剛度矩陣計(jì)算、位移及力計(jì)算以及計(jì)算結(jié)果的保存與輸出。關(guān)鍵在于空間框架的剛度矩陣計(jì)算。2.1.2 主程序主程序共分為 7 大模塊，分別為建立模型；集成總剛；求解位移；求解內(nèi)力；保存結(jié)果；輸出結(jié)果；繪制內(nèi)力圖。見“源代碼”文件夾中的“主程序”。2.1.3 函數(shù)主程序中主要用到 3 個(gè)關(guān)鍵計(jì)算用函數(shù)以及 6 個(gè)繪圖用函數(shù)，

8、分別為：計(jì)算單元?jiǎng)偠染仃嚕豢傮w剛度矩陣；求解單元節(jié)點(diǎn)力；繪制單元軸力圖；繪制單元剪力（Y 軸）圖；5繪制單元剪力（Z 軸）圖；繪制單元扭矩圖；繪制單元彎矩（Y 軸）圖；繪制單元彎矩（Z 軸）圖。見“源代碼”文件夾中的“函數(shù)”。2.2 程序代碼2.1.1 主程序代碼clear建立模型%*%1.建立模型%*%模型基本參數(shù)fprf(請(qǐng)輸入模型參數(shù)n);X=input(水平向?qū)訑?shù)：);Y=input(豎向?qū)訑?shù)：);Z=input(橫向?qū)訑?shù)：);H=input(豎桿長度（m）：);L=input(水平桿長度（m）：);W=input(橫桿長度（m）：);A=input(桿件截面面積（m2）：);Iy=i

9、nput(在 xz 面內(nèi)截面慣性矩（m4）：);Iz=input(在 xy 面內(nèi)截面慣性矩（m4）：);J=input(扭轉(zhuǎn)慣性矩（m4）：);E=input(材料彈性模量（MPa）：);G=input(材料剪切模量（MPa）：);P=input(節(jié)點(diǎn)力 N_y（kN）：);%定義節(jié)點(diǎn)坐標(biāo)6Node=zeros(X+1)*(Y+1)*(Z+1),3);x=0;y=0;z=0;t1=1;for s=1:(Z+1)z=W*(s-1);for j=1:(Y+1)y=H*(j-1);for i=1:(X+1)x=L*(i-1);Node(t1,:)=xy z;t1=t1+1;endendend%節(jié)點(diǎn)的

10、總數(shù)量，用 t1n 表示t1n=t1-1;fprf(n“節(jié)點(diǎn)建立”，完成！n);%定義單元Element=zeros(3*X*Z+2*(X+Z)+1)*Y),2);%定義平行于 y 軸的豎桿單元t1=1;t2=1;for s=1:(Z+1)for j=1:(Y+1)for i=1:(X+1)if j1&i1&i=(X+1)Element(t2,:)=t1(t1+(X+1)*(Y+1);t2=t2+1;endt1=t1+1;endendendfor j=1:(Y+1)for i=1:(X+1)if j1&iziLambda=0 0 1;0 1 0;-1 0 0;elseLambda=0 0 -1

11、;0 1 0;1 0 0;endelseCxX=(xj-xi)/L;CxY=(yj-yi)/L;20CxZ=(zj-zi)/L;D=sqrt(CxX*CxX+CxY*CxY);CyX=-CxY/D;CyY=CxX/D;CyZ=0;CzX=-CxX*CxZ/D;CzY=-CxY*CxZ/D;CzZ=D;Lambda=CxX CxY CxZ;.CyX CyY CyZ;.CzX CzY CzZ;endT=Lambda zeros(3) zeros(3) zeros(3);.zeros(3)Lambda zeros(3) zeros(3);.zeros(3)zeros(3) Lambda zeros(

12、3);.zeros(3)zeros(3) zeros(3) Lambda;%計(jì)算整體坐標(biāo)系下的單元?jiǎng)偠染仃噛=T*kprime*T;end總體剛度矩陣%集成總體剛度矩陣function y=Assemble(K,k,i,j)%將第 16 行數(shù)據(jù)加入總剛for s=1:6%將第 s 行數(shù)據(jù)加入總剛t=6-s;K(6*i-t,6*i-5)=K(6*i-t,6*i-5)+k(s,1);K(6*i-t,6*i-4)=K(6*i-t,6*i-4)+k(s,2);K(6*i-t,6*i-3)=K(6*i-t,6*i-3)+k(s,3);K(6*i-t,6*i-2)=K(6*i-t,6*i-2)+k(s,4

13、);21K(6*i-t,6*i-1)=K(6*i-t,6*i-1)+k(s,5);K(6*i-t,6*i)=K(6*i-t,6*i)+k(s,6);K(6*i-t,6*j-5)=K(6*i-t,6*j-5)+k(s,7);K(6*i-t,6*j-4)=K(6*i-t,6*j-4)+k(s,8);K(6*i-t,6*j-3)=K(6*i-t,6*j-3)+k(s,9);K(6*i-t,6*j-2)=K(6*i-t,6*j-2)+k(s,10);K(6*i-t,6*j-1)=K(6*i-t,6*j-1)+k(s,11);K(6*i-t,6*j)=K(6*i-t,6*j)+k(s,12);end%將

14、第 712 行數(shù)據(jù)加入總剛fors=7:12%將第 s 行數(shù)據(jù)加入總剛t=12-s;K(6*j-t,6*i-5)=K(6*j-t,6*i-5)+k(s,1);K(6*j-t,6*i-4)=K(6*j-t,6*i-4)+k(s,2);K(6*j-t,6*i-3)=K(6*j-t,6*i-3)+k(s,3);K(6*j-t,6*i-2)=K(6*j-t,6*i-2)+k(s,4);K(6*j-t,6*i-1)=K(6*j-t,6*i-1)+k(s,5);K(6*j-t,6*i)=K(6*j-t,6*i)+k(s,6);K(6*j-t,6*j-5)=K(6*j-t,6*j-5)+k(s,7);K(6

15、*j-t,6*j-4)=K(6*j-t,6*j-4)+k(s,8);K(6*j-t,6*j-3)=K(6*j-t,6*j-3)+k(s,9);K(6*j-t,6*j-2)=K(6*j-t,6*j-2)+k(s,10);K(6*j-t,6*j-1)=K(6*j-t,6*j-1)+k(s,11);K(6*j-t,6*j)=K(6*j-t,6*j)+k(s,12);endy=K;end求解單元節(jié)點(diǎn)力%計(jì)算單元節(jié)點(diǎn)力矢量22function y=ElementForce(E,G,A,Iy,Iz,J,xi,yi,zi,xj,yj,zj,u)%計(jì)算桿件長度，本程序桿件長度均為 L 或 H 或 WL=(xj

16、-xi)2+(yj-yi)2+(zj-zi)2)(1/2);%計(jì)算中的非零項(xiàng)a1=E*A/L;a2=G*J/L;b1=12*E*Iz/(L3);b2=6*E*Iz/(L2);b3=2*E*Iz/L;c1=12*E*Iy/(L3);c2=6*E*Iy/(L2);c3=2*E*Iy/L;%計(jì)算局部坐標(biāo)系下的單元?jiǎng)偠染仃噆prime=a1 0 0 0 0 0 -a1 0 0 0 0 0;.0b1 0 00 b2 0 -b1 0 0 0 b2;.00c1 0-c2 000-c1 0 -c2 0;.000 a20 0 000-a2 0 0;.00-c2 0 2*c3000 c2 0 c3 0;.0b2

17、0 0 0 2*b3 0 -b2 0 0 0 b3;.-a1 0 000 0 a1 0 0 0 0 0;.0-b1 000 -b2 0 b1 0 0 0 -b2;.00-c10c2 0 0 0 c1 0 c2 0;.000 -a20 0 0 0 0 a2 0 0;.00-c2 0c3 0 0 0 c2 0 2*c3 0;.0b2 0 0 0 b3 0 -b2 0 0 0 2*b3;%計(jì)算由局部坐標(biāo)系 x、y、z 軸轉(zhuǎn)換到整體坐標(biāo)系 X、Y、Z 軸的坐標(biāo)轉(zhuǎn)換矩陣 Tif xi=xj&yi=yjif zjziLambda=0 0 1;0 1 0;-1 0 0;else23Lambda=0 0 -1

18、;0 1 0;100;endelseCxX=(xj-xi)/L;CxY=(yj-yi)/L;CxZ=(zj-zi)/L;D=sqrt(CxX*CxX+CxY*CxY);CyX=-CxY/D;CyY=CxX/D;CyZ=0;CzX=-CxX*CxZ/D;CzY=-CxY*CxZ/D;CzZ=D;Lambda=CxX CxY CxZ;.CyX CyY CyZ;.CzX CzY CzZ;endT=Lambda zeros(3) zeros(3) zeros(3);.zeros(3)Lambda zeros(3) zeros(3);.zeros(3)zeros(3) Lambda zeros(3);.

19、zeros(3)zeros(3) zeros(3) Lambda;%計(jì)算整體坐標(biāo)系下的單元?jiǎng)偠染仃噛=kprime*T*u;end繪制單元軸力圖%繪制節(jié)點(diǎn)力矢量沿桿長變化的軸向力曲線圖function y=ElementAxialDiagram(f,)x=0;L;z=-f(1);f(7);hold on;24title(軸力圖-Element ,num2str(ie);xlabel(L);ylabel(Nx);plot(x,z);y1=0;0;plot(x,y1,k);end繪制單元剪力（Y 軸）圖%繪制節(jié)點(diǎn)力矢量沿桿長變化的剪力曲線圖（沿 Z 軸）function y=ElementShea

20、rYDiagram(f,)x=0;L;z=f(2);-f(8);hold on;title(剪力圖(Y)-Element ,num2str(ie);xlabel(L);ylabel(Ny);plot(x,z);y1=0;0;plot(x,y1,k);end繪制單元剪力（Z 軸）圖%繪制節(jié)點(diǎn)力矢量沿桿長變化的剪力曲線圖（沿 Z 軸）function y=ElementShearZDiagram(f,)x=0;L;z=f(3);-f(9);hold on;title(剪力圖(Z)-Element ,num2str(ie);xlabel(L);ylabel(Nz);plot(x,z);25y1=0;0;plot(x,y1,k);end繪制單元扭矩圖%繪制節(jié)點(diǎn)力矢量沿桿長變化的扭矩曲線圖function y=ElementTorDiagram(f,)x=0;L;z=f(4);-f(10);hold on;title(扭矩圖-Element ,num2str(ie);xlabel(L);ylabel(Mx);plot(x,z);y1=0;0;plot(x,y1,k);end繪制單元彎矩（Y 軸）圖%繪制節(jié)點(diǎn)力矢量沿桿長變化的彎矩曲線圖（沿 Y 軸）functi