有限元法15年作业

1、有限元法作业班级：学号：姓名：作业提交至：femhomework, 标题：班级+学号+姓名+有限元作业（不含+）1、已知弹性模量E=29.5x104N/mm2，A=100mm2，求结点位移、单元应力及支座反力。解：1） 计算程序1. 主程序prep;kw=fun_kw(XY,ELB,EA);P=fun_load(XY,NPQ);kw,P=fun_restrain1(kw,P,SU);U=fun_mguass(kw,P);stress=fun_interf(XY,ELB,U,EA);posp;2. 前处理程序clear;clc;XY=1,0,0 2,400e-3,0 3,400e-3,300e-

2、3 4,0,300e-3;EA=295e9,1e-4;ELB=1,1,2,1 2,2,3,1 3,1,3,1 4,3,4,1;NPQ=3,20000 6,-25000;SU=1,0 2,0 4,0 7,0 8,0;3. 后处理程序N,m=size(XY);U=reshape(U,2,N);U=U;U1=1:N; U=U1,UMB,m=size(ELB);stress1=1:MB;stress=stress1,stress 2）计算结果3） 功能函数1单元刚度矩阵函数function ke,t = fun_ke_t( dxy,E,A )%UNTITLED2 Summary of this fun

3、ction goes here% Detailed explanation goes herebm=dxy(2,2)-dxy(1,2);cm=dxy(2,1)-dxy(1,1);d=sqrt(bm2+cm2);s=bm/d;c=cm/d;ke=(E*A/d)*1,0,-1,0 0,0,0,0 -1,0,1,0 0,0,0,0;t=c,s,0,0 -s,c,0,0 0,0,c,s 0,0,-s,c; End2总体刚度矩阵函数function kw = fun_kw(XY,ELB,EA )%UNTITLED3 Summary of this function goes here% Detailed

4、 explanation goes hereN,m=size(XY);kw=zeros(2*N,2*N);MB,m=size(ELB);for k=1:MB ii=ELB(k,2); jj=ELB(k,3); LTB=ELB(k,4); dxy=XY(ii,2),XY(ii,3) XY(jj,2),XY(jj,3); E=EA(LTB,1); A=EA(LTB,2); ke,t=fun_ke_t(dxy,E,A); cn=2*ii-1,2*ii,2*jj-1,2*jj; ke=(t)*ke*t; for i=1:4 for j=1:4 kw(cn(i),cn(j)=kw(cn(i),cn(j)

5、+ke(i,j); end endend end3载荷列阵函数 function P = fun_load( XY,NPQ )%UNTITLED4 Summary of this function goes here% Detailed explanation goes hereN,m=size(XY);P=zeros(2*N,1);PQ,m=size(NPQ);for k=1:PQ i=NPQ(k,1); P(i)=NPQ(k,2);endend4约束处理函数function kw,P = fun_restrain1( kw,P,SU)%UNTITLED6 Summary of this f

6、unction goes here% Detailed explanation goes hereLR,m=size(SU);for k=1:LR i=SU(k,1); kw(i,:)=0; kw(i,i)=1; P(i)=SU(k,2);end5结点位移函数function x = fun_mguass( a,b )%UNTITLED7 Summary of this function goes here% Detailed explanation goes heren,m=size(a);for (k=1:n-1) for(i=k+1:n) if abs(a(i,k)abs(a(k,k)

7、w=a(k,:); a(k,:)=a(i,:); a(i,:)=w; w=b(k); b(k)=b(i); b(i)=w; end endendx=ab; end6. 单元内力函数function stress = fun_interf( XY,ELB,U,EA )%UNTITLED8 Summary of this function goes here% Detailed explanation goes hereMB,m=size(ELB);stress=zeros(MB,5);for k=1:MB ii=ELB(k,2); jj=ELB(k,3); LTB=ELB(k,4); dxy=X

8、Y(ii,2),XY(ii,3) XY(jj,2),XY(jj,3); E=EA(LTB,1); A=EA(LTB,2); ke,t=fun_ke_t(dxy,E,A); d=U(2*ii-1) U(2*ii) U(2*jj-1) U(2*jj); ke1=(t)*ke*t; F1=ke1*d; F=t*F1; stress(k,1:4)=F; stress(k,5)=stress(k,3)/A;endstress=stress(1:MB,3),stress(1:MB,5); end 2、计算图示桁架的结点位移和支反力。解：1） 计算程序1.主程序prep;kw=fun_kw(XY,ELB,E

9、A);P=fun_load(XY,NPQ);kw,P=fun_restrain1(kw,P,SU);U=fun_mguass(kw,P);stress=fun_interf(XY,ELB,U,EA);posp;2. 前处理程序clear;clc;XY=1,0,0 2,3,0 3,6,3 4,3,3 5,0,3;EA=190e9,8e-4;ELB=1,1,2,1 2,2,5,1 3,4,5,1 4,2,4,1 5,3,4,1 6,2,3,1;NPQ=6,-5000 8,-5000;SU=1,0 2,0 9,0 10,0;3. 后处理程序N,m=size(XY);U=reshape(U,2,N);

10、U=U;U1=1:N;U=U1,UMB,m=size(ELB);stress1=1:MB;stress=stress1,stress 2）计算结果3）功能函数与前一题相同。3、求如图所示刚架的结点位移和支座反力。解：1） 计算程序1.主程序GJINPUT;KZ=kztb(XY,ELB,EAIZ);P=ydx(XY,ELB,EAIZ,ELPQ);U=weiyi(KZ,P,SU);PD=gdn1(U,XY,ELB,EAIZ,ELPQ);GJOUTPUT;2. 前处理程序clear;clc;XY=1,0,0.9 2,1,0.9 3,1,0;ELB=1,1,2,2 2,2,3,1;A1=7.65e-4

11、;IZ1=204e-8;A2=7.65e-4;IZ2=204e-8;EAIZ=30e9,A1,IZ1 30e9,A2,IZ2;ELPQ=1,1,8,1SU=1,0 2,0 7,0 8,0 9,0;3. 后处理程序clc;N,M=size(XY);DU=zeros(N,4);for i=1:N DU(i,1)=i; DU(i,2)=U(3*i-2); DU(i,3)=U(3*i-1); DU(i,4)=U(3*i);endMB,M=size(ELB);for i=1:MB P1(i,1)=i;endDUDP=P1,PD2）计算结果3） 功能函数1.单元刚度矩阵函数function KE,T =

12、ketb( dxy,E,A,IZ )%UNTITLED3 Summary of this function goes here% Detailed explanation goes hereDY=dxy(2,2)-dxy(1,2);DX=dxy(2,1)-dxy(1,1);L=sqrt(DX2+DY2);S=DY/L;C=DX/L;a1=IZ/L;a2=a1/L;a3=E/L;KE=a3*A,0,0,-A,0,0 0,12*a2,6*a1,0,-12*a2,6*a1 0,6*a1,4*IZ,0,-6*a1,2*IZ -A,0,0,A,0,0 0,-12*a2,-6*a1,0,12*a2,-6*

13、a1 0,6*a1,2*IZ,0,-6*a1,4*IZ; t=C,S,0;-S,C,0;0,0,1; t1=zeros(3,3); T=t,t1;t1,t;End2. 总体刚度矩阵函数function KZ = kztb( XY,ELB,EAIZ )%UNTITLED4 Summary of this function goes here% Detailed explanation goes hereN,m=size(XY);KZ=zeros(3*N,3*N);MB,m=size(ELB);for k=1:MB ii=ELB(k,2); jj=ELB(k,3); LTB=ELB(k,4); d

14、xy=XY(ii,2),XY(ii,3) XY(jj,2),XY(jj,3); E=EAIZ(LTB,1); A=EAIZ(LTB,2); IZ=EAIZ(LTB,3); KE,T=ketb(dxy,E,A,IZ); CN=3*ii-2,3*ii-1,3*ii,3*jj-2,3*jj-1,3*jj; KE=(T)*KE*T; for i=1:6 for j=1:6 KZ(CN(i),CN(j)=KZ(CN(i),CN(j)+KE(i,j); end end end 3. 载荷简化函数function PO,ii,jj = dxjd( ELB,XY,ELPQ1 )%UNTITLED5 Summa

15、ry of this function goes here% Detailed explanation goes herePO=zeros(6,1);k=ELPQ1(1);ii=ELB(k,2);jj=ELB(k,3);dxy=XY(ii,2),XY(ii,3) XY(jj,2),XY(jj,3);DY=dxy(2,2)-dxy(1,2);DX=dxy(2,1)-dxy(1,1);L=sqrt(DX2+DY2);C=ELPQ1(2);Q=ELPQ1(3);type=ELPQ1(4);switch type case 1 PO(2)=PO(2)+0.5*Q*C*(2-2*C2/L2+C3/L3)

16、; PO(3)= PO(3)+Q*C2*(6-8*C/L+3*C2/L2)/12; PO(5)= PO(5)+Q*C- PO(2); PO(6)= PO(6)-Q*C3*(4-3*C/L)/12/L; case 2 D=L-C; PO(2)= PO(2)+Q*(L+2*C)*D2/L3; PO(3)= PO(3)+Q*C*D2/L2; PO(5)= PO(5)+Q- PO(2); PO(6)= PO(6)+Q*D*C2/L2; case 3 D=L-C; PO(1)= PO(1)+Q*D/L; PO(4)= PO(4)+Q*C/L; case 4 PO(2)= PO(2)+7*Q*L/20;

17、PO(3)= PO(3)+Q*L2/20; PO(5)= PO(5)+3*Q*L/20; PO(6)= PO(6)+Q*L2/30;end end 4. 载荷列阵函数function P = ydx( XY,ELB,EAIZ,ELPQ)%UNTITLED6 Summary of this function goes here% Detailed explanation goes hereN,m=size(XY)P=zeros(3*N,1)LPQ,m=size(ELPQ)for j=1:LPQ ELPQ1=ELPQ(j,:) k=ELPQ(j,1) PO,ii,jj=dxjd(ELB,XY,EL

18、PQ1) LTB=ELB(k,4); dxy=XY(ii,2),XY(ii,3) XY(jj,2),XY(jj,3); E=EAIZ(LTB,1); A=EAIZ(LTB,2); IZ=EAIZ(LTB,3); KE,T=ketb(dxy,E,A,IZ); PO=(T)*PO; CN=3*ii-2,3*ii-1,3*ii,3*jj-2,3*jj-1,3*jj for i=1:6 P(CN(i)=P(CN(i)+PO(i); endendend 5. 结点位移函数function U = weiyi( KZ,P,SU )%UNTITLED7 Summary of this function goes here% Detailed explanation goes hereLR,m=size(SU);for k=1:LR i=SU(k,1)