有限单元法matlab编程实例.doc
有限单元法matlab编程实例主程序E=210e6；A=2e-2;I=5e-5;L1=3；L2=4；L3=3；k1=PlaneFrameElementStiffness(E,A，I,L1,90）；k2=PlaneFrameElementStiffness（E,A,I，L2，0）;k3=PlaneFrameElementStiffness（E，A,I,L3,270);K=zeros（12，12)；K=PlaneFrameAssemble(K，k1,1,2)；K=PlaneFrameAssemble(K，k2,2，3)；K=PlaneFrameAssemble（K，k3,3,4）k=K（4:9,4：9）;f=-20;0；0;0；0;12;u=kfU=0;0;0;u;0;0；0F=K*Uu1=U（1）;U(2);U(3）;U（4）；U(5）;U(6）;u2=U(4)；U(5）；U（6)；U(7）；U（8）；U(9）；u3=U(7）；U（8）;U（9）;U（10）；U(11）；U（12);f1=PlaneFrameElementForces（E,A,I，L1,90,u1)f2=PlaneFrameElementForces（E,A，I,L2,0,u2)f3=PlaneFrameElementForces(E，A，I,L3,270,u3）需调用的函数和子程序function y=PlaneFrameAssemble（K，k，i,j)PlaneFrameAssemble This function assembles the element stiffness%matrix k of the plane frame element with nodes i and j into the global%stiffness matrix K .This function returns the global stiffness matrix K afterthe element stiffness matrix k is assembled.K（3*i2，3*i-2)=K（3i2，3i2）+k(1，1）；K(3i-2，3i-1)=K(3i2,3*i-1）+k（1,2）；K(3*i2，3i)=K(3*i-2,3*i)+k（1,3）;K（3*i-2，3*j2）=K(3i2，3*j-2)+k(1,4）；K（3i-2,3*j1)=K（3i-2，3j-1)+k（1,5）；K(3i2，3*j）=K(3*i-2,3j）+k（1,6）;K(3i1，3*i2)=K(3i-1,3i2)+k（2,1)；K（3i-1,3i1）=K(3*i-1，3*i-1）+k(2，2)；K（3i-1，3i）=K(3i-1，3i)+k（2，3)；K（3i1,3j2)=K(3*i-1，3*j-2）+k（2,4);K(3*i-1，3j-1)=K(3i-1，3*j-1）+k(2，5)；K（3*i1，3*j)=K（3i1，3*j）+k（2,6）；K(3*i，3i-2）=K(3*i,3*i-2)+k(3，1)；K（3i,3*i1)=K（3i，3*i-1)+k(3，2);K(3*i，3i)=K(3i,3i）+k(3，3）；K（3i，3*j2)=K（3*i，3*j-2)+k（3,4）;K(3i,3j1)=K（3*i，3j1)+k(3,5)；K（3*i,3*j)=K(3i，3*j)+k(3,6）;K（3*j-2,3*i2）=K（3j2,3*i2)+k（4，1);K(3j-2,3i1）=K（3j2，3i-1）+k（4，2）;K(3*j2，3*i)=K（3*j-2,3*i）+k(4,3)；K(3*j2，3j-2）=K(3j-2，3j-2)+k（4，4）；K（3j2,3*j1)=K（3*j-2,3j1）+k（4，5);K(3j2,3j）=K（3j-2,3j）+k(4，6）;K(3*j-1，3*i-2）=K（3*j1,3i-2)+k（5,1)；K（3j-1,3*i-1)=K(3*j1，3*i-1）+k(5,2);K（3j-1，3*i）=K（3j-1,3i)+k（5,3)；K（3j-1,3j-2)=K(3*j-1，3*j2)+k（5，4)；K(3j-1，3j-1）=K（3*j-1,3j1)+k(5,5)；K(3*j1,3j)=K（3j-1,3*j）+k(5,6）;K（3*j，3*i-2)=K(3*j,3*i-2）+k(6，1)；K(3j，3i1)=K(3*j,3i-1)+k（6，2)；K(3j,3i）=K(3*j，3*i)+k（6，3);K（3*j,3*j-2)=K(3*j，3*j2)+k(6,4);K（3j，3j-1)=K(3j，3j-1)+k(6,5）;K（3j，3j)=K（3*j，3*j）+k（6，6)；y=K; function y=PlaneFrameElementForces(E，A，I,L，theta，u）PlaneFrameElementforce This function returns the element force given the modulus of elasticity% E，the cross sectional area A,the moment of inetia the length L，the angle theta ，and the element nodal% displacement vector u.x=thetapi/180;C=cos（x);S=sin(x）；w1=E*A/L;w2=12EI/(L3）;w3=6E*I/(L2)；w4=4*EI/L;w5=2*E*I/L；kprime=w1 0 0 -w1 0 0；0 w2 w3 0 w2 w3;0 w3 w4 0 -w3 w5;w1 0 0 w1 0 0；0 w2 -w3 0 w2 w3；0 w3 w5 0 w3 w4；T=C S 0 0 0 0；-S C 0 0 0 0；0 0 1 0 0 0；0 0 0 C S 0;0 0 0 S C 0;0 0 0 0 0 1；y=kprimeT*u;function y=PlaneFrameElementLength(x1，y1,x2,y2）%PlaneFrameElementLength This function returns the length of the% plane frame element whose first node% has coordinates(x1，y1）and second nodes has% coordinates(x2，y2)y=sqrt(x2x1）*（x2-x1）+(y2y1)(y2y1）);function y=PlaneFrameElementStiffness（E，A,I，L，theta)%PlaneFrameElementStiffness This function returns the stiffness matrix of the plane frame element with modulus of% elasticity E，cross sectional area A ，length L，moment of inertia and angle theta. x=thetapi/180；C=cos(x）；S=sin(x）；w1=AC*C+12*IS*S/(LL)；w2=A*SS+12IC*C/（LL)；w3=（A-12*I/(L*L)）*C*S；w4=6I*S/L；w5=6IC/L;y=E/Lw1 w3 -w4 w1 -w3 w4;w3 w2 w5 -w3 -w2 w5;-w4 w5 4*I w4 w5 2*I;-w1 w3 w4 w1 w3 w4;-w3 -w2 w5 w3 w2 -w5;-w4 w5 2*I w4 w5 4*I；
