有限元的matlab編程

有限元編程示例有限元的matlab編程題目描述:如下圖所示的平面桁架，桿件長度、彈性模量、截面積以及所受節(jié)點力P的大小可以自行定義。求節(jié)點位移及桿件軸力。例一：桁架有限元的matlab編程解題思路：建立模型集成總剛求解位移求解桿件軸力輸出結(jié)果有限元的matlab編程建立模型：定義節(jié)點坐標(biāo)Node=zeros(10,2);x=-1*L;%L為橫桿長度fori=1:2:10x=x+L;Node(i,:)=[x0];endx=-1*L;fori=2:2:10x=x+L;Node(i,:)=[xH];%H為豎桿長度end模型相關(guān)參數(shù)輸入H=input('豎桿長度(m):');

L=input('水平桿長度(m):');E=input('桿件彈性模量(Gpa):');A=input('桿件截面積(m^2):');a=input('節(jié)點力P(kN)：');節(jié)點編號方式有限元的matlab編程定義單元，即儲存單元兩端的節(jié)點號Element=zeros(21,2);fori=1:2:7Element(5/2*i-3/2,:)=[i,i+1];Element(5/2*i-1/2,:)=[i,i+2];Element(5/2*i+1/2,:)=[i,i+3];endfori=2:2:8Element(5*i/2-1,:)=[i,i+1];Element(5*i/2,:)=[i,i+2];endElement(21,:)=[9,10];加下劃線的為單元編號有限元的matlab編程集成總剛：xi=Node(Element(ie,1),1);%ie為單元號，以下相同yi=Node(Element(ie,1),2);xj=Node(Element(ie,2),1);yj=Node(Element(ie,2),2);獲取單元兩端節(jié)點坐標(biāo)L=((xj-xi)^2+(yj-yi)^2)^(1/2);計算桿件長度形成等效荷載列陣f=[0;0;0;a;0;0;0;a;0;0;0;a;0;0;0;a;0;0;0;a];%每個節(jié)點兩個自由度，a為之前輸入的節(jié)點力有限元的matlab編程計算從局部坐標(biāo)到整體坐標(biāo)的坐標(biāo)轉(zhuǎn)換矩陣TfunctionT=TransformMatrix(ie)%ie為單元號c=(xj-xi)/L;s=(yj-yi)/L;T=[c-s00sc0000c-s00sc];計算單元剛度矩陣kk=[E*A/L0-E*A/L00000-E*A/L0E*A/L00000];T=TransformMatrix(ie);k=T*k*transpose(T);%transpose(T)為T的轉(zhuǎn)置矩陣2有限元的matlab編程集成整體剛度矩陣Kforie=1:1:21%按單元順序進行循環(huán)k=PlaneTrussElementStiffness(ie);%計算第ie個單元的單剛m=Element(ie,1);%ie單元的首節(jié)點號n=Element(ie,2);%ie單元的末節(jié)點號K(2*m-1,2*n-1)=k(1,3);K(2*m-1,2*n)=k(1,4);K(2*m,2*n-1)=k(2,3);K(2*m,2*n)=k(2,4);K=zeros(20,20);%用來存儲整體剛度矩陣集成總剛的非對角線元素(這里的元素指2*2的小矩陣）在下面的集成中，將總剛看成10*10的矩陣，每個元素為2*2的小矩陣m=Element(ie,2);%ie單元的末節(jié)點號n=Element(ie,1);%ie單元的首節(jié)點號K(2*m-1,2*n-1)=k(3,1);K(2*m-1,2*n)=k(3,2);K(2*m,2*n-1)=k(4,1);K(2*m,2*n)=k(4,2);end有限元的matlab編程集成總剛的對角線元素（這里的元素指2*2的小矩陣）fori=1:1:10%按節(jié)點的順序循環(huán)forj=1:1:21%對于每個節(jié)點，再按單元的順序循環(huán)k=PlaneTrussElementStiffness(j);ifElement(j,1)==I

%如果i節(jié)點為j單元的首節(jié)點K(2*i-1,2*i-1)=K(2*i-1,2*i-1)+k(1,1);K(2*i-1,2*i)=K(2*i-1,2*i)+k(1,2);K(2*i,2*i-1)=K(2*i,2*i-1)+k(2,1);K(2*i,2*i)=K(2*i,2*i)+k(2,2);endifElement(j,2)==i

%如果i節(jié)點為j單元的末節(jié)點K(2*i-1,2*i-1)=K(2*i-1,2*i-1)+k(3,3);K(2*i-1,2*i)=K(2*i-1,2*i)+k(3,4);K(2*i,2*i-1)=K(2*i,2*i-1)+k(4,3);K(2*i,2*i)=K(2*i,2*i)+k(4,4);endendend有限元的matlab編程求解位移：u=zeros(20);根據(jù)約束情況修改總剛，采用對角元素置1法fori=1:1:20

K(1,i)=0;K(2,i)=0;K(18,i)=0;K(i,1)=0;K(i,2)=0;K(i,18)=0;end

%自由度1、2、18被約束了,所在的行和列的其他元素都改為0K=K*1e15;%乘以一個大數(shù)，減小計算誤差f=f*1e15;u=K\f;求解K(1,1)=1;%對角線元素置1K(2,2)=1;K(18,18)=1;有限元的matlab編程求解軸力：獲取單元兩端的節(jié)點號i=Element(ie,1);%ie為單元號j=Element(ie,2);獲取單元兩端的節(jié)點位移uElement=zeros(4,1);uElement(1:2)=u((i-1)*2+1:(i-1)*2+2);uElement(3:4)=u((j-1)*2+1:(j-1)*2+2);k=PlaneTrussElementStiffness(ie);nodef=k*uElement;%整體坐標(biāo)下的節(jié)點力T=TransformMatrix(ie);%計算坐標(biāo)轉(zhuǎn)換矩陣nodef=transpose(T)*nodef;%從整體坐標(biāo)轉(zhuǎn)換到局部坐標(biāo)計算單元的節(jié)點力有限元的matlab編程輸出求解結(jié)果：輸出位移fprintf('節(jié)點位移\n');fori=1:1:10disp(['節(jié)點號',num2str(i),'x方向位移：',num2str(u(2*i-1,1)),'y方向位移：',num2str(u(2*i,1))]);end輸出節(jié)點力fprintf('\n\n節(jié)點力\n');forie=1:1:21nodef=NodeForce(ie);disp(['單元號:',num2str(ie),'節(jié)點號：',num2str(Element(ie,1)),'節(jié)點號：',num2str(Element(ie,2)),'軸力:',num2str(nodef(1))]);end有限元的matlab編程例二：網(wǎng)架有限元的matlab編程思路分析網(wǎng)架是由多根桿件按照一定的網(wǎng)格形式通過節(jié)點連結(jié)而成的空間結(jié)構(gòu)。構(gòu)成網(wǎng)架的基本單元有三角錐，三棱體，正方體，截頭四角錐等。鑒于網(wǎng)架的形式較多，本程序提供一種通用的網(wǎng)架輸入方法，但錄入較為繁瑣，同時提供一種正放四角錐網(wǎng)架的簡易輸入方法作為典型?？紤]幾何非線性。本程序采用荷載分級加載的方式考慮網(wǎng)架的幾何非線性。將總荷載分成1000份分步施加，求解各荷載步下的節(jié)點位移，修改網(wǎng)架相應(yīng)節(jié)點坐標(biāo)以及剛度矩陣，依次迭代求出網(wǎng)架的總位移。本程序的網(wǎng)架位移求解函數(shù)附在主程序后面，主程序運行時調(diào)用該函數(shù)。幾點說明有限元的matlab編程用戶自定義輸入幾何建模正放四角錐網(wǎng)架簡易輸入定義荷載定義邊界條件網(wǎng)架分析位移應(yīng)變應(yīng)力位移求解函數(shù)單剛矩陣荷載矩陣約束矩陣總剛矩陣求解位移&&&分級加載，通過修改節(jié)點坐標(biāo)，迭代求解幾何非線性有限元的matlab編程e=input('選擇網(wǎng)架類型，0代表自由定義網(wǎng)架，1代表四角錐網(wǎng)架')%網(wǎng)架類型的選擇網(wǎng)架類型的選擇用戶自定義網(wǎng)架（網(wǎng)架信息的錄入，包括節(jié)點、單元、截面、彈性模量等）ife==0%選擇自定義網(wǎng)架Node=input(‘定義節(jié)點編號及對應(yīng)坐標(biāo)，按[1x1y1z1;2x2y2z2;...]輸入’);%形成節(jié)點儲存矩陣Men=input(‘定義單元與節(jié)點的關(guān)系，按[1node1node2;2node3node4;...]輸入,node1<node2,依次類推’);%形成單元儲存矩陣Msum=length(Men);%查找網(wǎng)架錄入的單元數(shù)Cont1=input('定義單元實常數(shù)，若所有桿件截面面積和彈性模量不變，則輸入0，否則輸入1');定義單元屬性的輸入方式有限元的matlab編程ifCont1==0AE1=input('請輸入統(tǒng)一的截面面積與彈性模量，按[AE]輸入');AE=zeros(Msum,3);AE(:,1)=1:Msum;AE(:,2)=AE1(1,1);AE(:,3)=AE1(1,2);elseAE=input('請輸入相應(yīng)單元的截面面積與彈性模量，按[1,A1E1;2,A2E2;...]輸入');endP=input(‘定義節(jié)點荷載，按[node1P1;node2P2;...]輸入’);%網(wǎng)架荷載輸入BC=input(‘定義邊界約束，按[node1ConxConyConz;node2ConxConyConz);...]輸入,Con代表x、y、z方向約束，取0為約束，取1無約束’);%網(wǎng)架邊界條件end單元屬性相同單元屬性不同荷載及邊界條件有限元的matlab編程正放四角錐網(wǎng)架定義ife==1hu=input('輸入網(wǎng)架上層節(jié)點行數(shù)');%定義網(wǎng)架上層節(jié)點的行數(shù)lu=input('輸入網(wǎng)架上層節(jié)點列數(shù)');%定義網(wǎng)架上層節(jié)點的列數(shù)dis_xu=input('輸入網(wǎng)架上層節(jié)點列間距');%定義網(wǎng)架上層的行間距dis_yu=input('輸入網(wǎng)架上層節(jié)點行間距');%定義網(wǎng)架上層的列間距hd=hu-1;%網(wǎng)架下層節(jié)點的行數(shù)ld=lu-1;%網(wǎng)架下層節(jié)點的列數(shù)dis_xd=dis_xu;%網(wǎng)架下層的行間距dis_yd=dis_yu;%網(wǎng)架下層的行間距dis_z=input('輸入網(wǎng)架上下層間距');%網(wǎng)架上下層間距定義網(wǎng)架上層節(jié)點定義網(wǎng)架下層節(jié)點定義網(wǎng)架高度有限元的matlab編程fori=1:huforj=1:luNode((i-1)*lu+j,2)=(j-1)*dis_xu;Node((i-1)*lu+j,3)=(i-1)*dis_yu;Node((i-1)*lu+j,4)=dis_z;endendfori=1:hdforj=1:ldNode((i-1)*ld+j+hu*lu,2)=(j-1+0.5)*dis_xd;Node((i-1)*ld+j+hu*lu,3)=(i-1+0.5)*dis_yd;Node((i-1)*ld+j+hu*lu,4)=0;endend網(wǎng)架上層節(jié)點編號與對應(yīng)坐標(biāo)網(wǎng)架下層節(jié)點編號與對應(yīng)坐標(biāo)有限元的matlab編程Nsum=length(Node);%查詢網(wǎng)架的節(jié)點數(shù)fori=1:Nsum%將節(jié)點編號錄入節(jié)點矩陣Node(i,1)=i;endfori=1:huforj=1:lu-1Men((i-1)*(lu-1)+j,2)=(i-1)*lu+j;Men((i-1)*(lu-1)+j,3)=(i-1)*lu+j+1;endendfori=1:luforj=1:hu-1Men((i-1)*(hu-1)+j+(lu-1)*hu,2)=(j-1)*lu+i;Men((i-1)*(hu-1)+j+(lu-1)*hu,3)=j*lu+i;endend節(jié)點編號的錄入網(wǎng)架上層橫向單元的拓?fù)渚W(wǎng)架上層縱向單元的拓?fù)溆邢拊膍atlab編程fori=1:hdforj=1:ld-1Men((i-1)*(ld-1)+(lu-1)*hu+(hu-1)*lu+j,2)=(i-1)*ld+j+hu*lu;Men((i-1)*(ld-1)+(lu-1)*hu+(hu-1)*lu+j,3)=(i-1)*ld+j+hu*lu+1;endendfori=1:ldforj=1:hd-1Men((i-1)*(hd-1)+(ld-1)*hd+(lu-1)*hu+(hu-1)*lu+j,2)=(j-1)*ld+i+hu*lu;Men((i-1)*(hd-1)+(ld-1)*hd+(lu-1)*hu+(hu-1)*lu+j,3)=j*ld+i+hu*lu;endend網(wǎng)架下層縱向單元的拓?fù)渚W(wǎng)架下層橫向單元的拓?fù)溆邢拊膍atlab編程網(wǎng)架腹桿單元的拓?fù)鋐ori=1:hdforj=1:ldMen(((i-1)*ld+j-1)*4+(hu-1)*lu+(lu-1)*hu+(hd-1)*ld+(ld-1)*hd+1,2)=(i-1)*lu+j;Men(((i-1)*ld+j-1)*4+(hu-1)*lu+(lu-1)*hu+(hd-1)*ld+(ld-1)*hd+1,3)=(i-1)*ld+hu*lu+j;Men(((i-1)*ld+j-1)*4+(hu-1)*lu+(lu-1)*hu+(hd-1)*ld+(ld-1)*hd+2,2)=(i-1)*lu+j+1;Men(((i-1)*ld+j-1)*4+(hu-1)*lu+(lu-1)*hu+(hd-1)*ld+(ld-1)*hd+2,3)=(i-1)*ld+j+hu*lu;Men(((i-1)*ld+j-1)*4+(hu-1)*lu+(lu-1)*hu+(hd-1)*ld+(ld-1)*hd+3,2)=i*lu+j;Men(((i-1)*ld+j-1)*4+(hu-1)*lu+(lu-1)*hu+(hd-1)*ld+(ld-1)*hd+3,3)=(i-1)*ld+j+hu*lu;Men(((i-1)*ld+j-1)*4+(hu-1)*lu+(lu-1)*hu+(hd-1)*ld+(ld-1)*hd+4,2)=i*lu+j+1;Men(((i-1)*ld+j-1)*4+(hu-1)*lu+(lu-1)*hu+(hd-1)*ld+(ld-1)*hd+4,3)=(i-1)*ld+j+hu*lu;endend腹桿N腹桿N+1腹桿N+2腹桿N+3有限元的matlab編程單元編號錄入單元儲存矩陣Msum=length(Men);%查詢網(wǎng)架單元數(shù)

fori=1:Msum%將單元編號錄入單元儲存矩陣Men(i,1)=i;end定義截面屬性E=2.1e11;%默認(rèn)材料為steelA1=input('請輸入網(wǎng)架上層單元的截面面積');%默認(rèn)網(wǎng)架上層單元截面尺寸相同A2=input('請輸入網(wǎng)架下層單元的截面面積');%默認(rèn)網(wǎng)架下層單元截面尺寸相同A3=input('請輸入網(wǎng)架腹桿單元的截面面積');%默認(rèn)網(wǎng)架腹桿單元截面尺寸相同AE=zeros(Msum,3);%定義單元屬性矩陣m1=(hu-1)*lu+(lu-1)*hu;%上層單元截止編號m2=m1+(hd-1)*ld+(ld-1)*hd;%下層單元截止編號有限元的matlab編程AE(1:m1,2)=A1;%將上層單元尺寸錄入AE矩陣AE((m1+1):m2,2)=A2;%將下層單元尺寸錄入AE矩陣AE((m2+1):Msum,2)=A3;%將腹桿單元尺寸錄入AE矩陣AE(:,1)=1:Msum;%將單元編號錄入AE矩陣AE(:,3)=E;%將材料彈性模量錄入AE矩陣定義荷載cont2=input('定義節(jié)點荷載，若網(wǎng)架上層節(jié)點力與下層節(jié)點力均布，則輸入0，否則輸入1');ifcont2==0P1=input('請輸入網(wǎng)架上層節(jié)點荷載');P2=input('請輸入網(wǎng)架下層節(jié)點荷載');m3=hu*lu;P(1:Nsum,1)=1:Nsum;P(1:m3,2)=P1;P((m3+1):Nsum,2)=P2;elseP=input('定義節(jié)點荷載，按[node1P1;node2P2;...]輸入');end有限元的matlab編程定義邊界條件cont3=input('定義邊界約束，若網(wǎng)架上層周邊節(jié)點全約束，則輸入0，若下層周邊節(jié)點全約束，輸入1，否則輸入2');ifcont3==0n1=2*(hu+lu-2);BC=zeros(n1,4);BC(1:lu-2,1)=2:lu-1;BC((lu-1):(2*lu-4),1)=lu*(hu-1)+2:lu*hu-1;BC((2*lu-3):(2*lu-4+hu),1)=1:lu:lu*(hu-1)+1;BC((2*lu-3+hu):n1,1)=lu:lu:hu*lu;有限元的matlab編程elseifcont3==1n1=2*(hd+ld-2);BC=zeros(n1,4);BC(1:ld-2,1)=2:ld-1;BC((ld-1):(2*ld-4),1)=ld*(hd-1)+2:ld*hd-1;BC((2*ld-3):(2*ld-4+hd),1)=1:ld:ld*(hd-1)+1;BC((2*ld-3+hd):n1,1)=ld:ld:hd*ld;fori=1:n1BC(i,1)=BC(i,1)+hu*lu;endelseBC=input('定義邊界約束，按[node1ConxConyConz;node2ConxConyConz);...]輸入,Con代表x、y、z方向約束，取0為約束，取1無約束');endend有限元的matlab編程Nsum=length(Node);Msum=length(Men);Psum=length(P);BCsum=length(BC);%提取各矩陣的行數(shù)考慮幾何非線性分析網(wǎng)架fori=1:Psum%將力分為1000份P(i,2)=P(i,2)/1000;endU=zeros(3*Nsum,1);%總位移矩陣fori=1:1000[u,L1,Kz]=grid(Node,Men,AE,P,BC,Nsum,Msum,Psum,BCsum);forj=1:Nsum%根據(jù)節(jié)點位移修改網(wǎng)架的節(jié)點坐標(biāo)Node(j,2)=Node(j,2)+u(3*j-2,1);Node(j,3)=Node(j,3)+u(3*j-1,1);Node(j,4)=Node(j,4)+u(3*j,1);endU=U+u;%每次迭代位移的疊加end迭代法修正剛度矩陣和網(wǎng)架位移有限元的matlab編程求解網(wǎng)架桿件的應(yīng)力應(yīng)變L0=zeros(Msum,1);

%所有根桿的最初長度fori=1:Msum%單元兩端的節(jié)點編號p=Men(i,2);q=Men(i,3);X1=Node(p,2);%單元端節(jié)點的坐標(biāo)Y1=Node(p,3);Z1=Node(p,4);X2=Node(q,2);Y2=Node(q,3);Z2=Node(q,4);L0(i,1)=sqrt((X2-X1)^2+(Y2-Y1)^2+(Z2-Z1)^2);%網(wǎng)架桿件的初始長度end有限元的matlab編程Lt=L1-L0;%所有桿件長度的增量strain=zeros(Msum,1);%定義應(yīng)變矩陣stress=zeros(Msum,1);%定義應(yīng)力矩陣fori=1:MsumE=AE(i,3);strain(i,1)=Lt(i,1)/L0(i,1);%第i根桿件應(yīng)變stress(i,1)=E*strain(i,1);%第i根桿件應(yīng)力enddisp(U);%輸出網(wǎng)架節(jié)點位移disp(stress);%輸出網(wǎng)架桿件應(yīng)力有限元的matlab編程網(wǎng)架節(jié)點位移求解函數(shù)function[u,L1,Kz]=grid(Node,Men,AE,P,BC,Nsum,Msum,Psum,BCsum)Kz=zeros(3*Nsum,3*Nsum);%定義剛度矩陣L=zeros(Msum,1);fori=1:Msum

%單元兩端的節(jié)點編號p=Men(i,2);q=Men(i,3);A=AE(i,2);E=AE(i,3);%單元兩端節(jié)點坐標(biāo)X1=Node(p,2);Y1=Node(p,3);Z1=Node(p,4);X2=Node(q,2);Y2=Node(q,3);Z2=Node(q,4);

%單元長度L(i,1)=sqrt((X2-X1)^2+(Y2-Y1)^2+(Z2-Z1)^2);%單元的方向余弦l=(X2-X1)/L(i,1);m=(Y2-Y1)/L(i,1);n=(Z2-Z1)/L(i,1);有限元的matlab編程%定義轉(zhuǎn)化矩陣t=sqrt(l^2+n^2);ift==0r=[0m0;-m00;001];elser=[lmn;-l*m/tt-m*n/t;-n/t0l/t];endO=zeros(3,3);T=[rO;Or];%整體坐標(biāo)下的單剛矩陣ke=E*A/L(i,1)*[100-100;0000