有限单元法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=k\f

U=[0;0;0;u;0;0;0]

F=K*U

u1=[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 after

%the element stiffness matrix k is assembled.

K(3*i-2,3*i-2)=K(3*i-2,3*i-2)+k(1,1);

K(3*i-2,3*i-1)=K(3*i-2,3*i-1)+k(1,2);

K(3*i-2,3*i)=K(3*i-2,3*i)+k(1,3);

K(3*i-2,3*j-2)=K(3*i-2,3*j-2)+k(1,4);

K(3*i-2,3*j-1)=K(3*i-2,3*j-1)+k(1,5);

K(3*i-2,3*j)=K(3*i-2,3*j)+k(1,6);

K(3*i-1,3*i-2)=K(3*i-1,3*i-2)+k(2,1);

K(3*i-1,3*i-1)=K(3*i-1,3*i-1)+k(2,2);

K(3*i-1,3*i)=K(3*i-1,3*i)+k(2,3);

K(3*i-1,3*j-2)=K(3*i-1,3*j-2)+k(2,4);

K(3*i-1,3*j-1)=K(3*i-1,3*j-1)+k(2,5);

K(3*i-1,3*j)=K(3*i-1,3*j)+k(2,6);

K(3*i,3*i-2)=K(3*i,3*i-2)+k(3,1);

K(3*i,3*i-1)=K(3*i,3*i-1)+k(3,2);

K(3*i,3*i)=K(3*i,3*i)+k(3,3);

K(3*i,3*j-2)=K(3*i,3*j-2)+k(3,4);

K(3*i,3*j-1)=K(3*i,3*j-1)+k(3,5);

K(3*i,3*j)=K(3*i,3*j)+k(3,6);

K(3*j-2,3*i-2)=K(3*j-2,3*i-2)+k(4,1);

K(3*j-2,3*i-1)=K(3*j-2,3*i-1)+k(4,2);

K(3*j-2,3*i)=K(3*j-2,3*i)+k(4,3);

K(3*j-2,3*j-2)=K(3*j-2,3*j-2)+k(4,4);

K(3*j-2,3*j-1)=K(3*j-2,3*j-1)+k(4,5);

K(3*j-2,3*j)=K(3*j-2,3*j)+k(4,6);

K(3*j-1,3*i-2)=K(3*j-1,3*i-2)+k(5,1);

K(3*j-1,3*i-1)=K(3*j-1,3*i-1)+k(5,2);

K(3*j-1,3*i)=K(3*j-1,3*i)+k(5,3);

K(3*j-1,3*j-2)=K(3*j-1,3*j-2)+k(5,4);

K(3*j-1,3*j-1)=K(3*j-1,3*j-1)+k(5,5);

K(3*j-1,3*j)=K(3*j-1,3*j)+k(5,6);

K(3*j,3*i-2)=K(3*j,3*i-2)+k(6,1);

K(3*j,3*i-1)=K(3*j,3*i-1)+k(6,2);

K(3*j,3*i)=K(3*j,3*i)+k(6,3);

K(3*j,3*j-2)=K(3*j,3*j-2)+k(6,4);

K(3*j,3*j-1)=K(3*j,3*j-1)+k(6,5);

K(3*j,3*j)=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=theta*pi/180;

C=cos(x);

S=sin(x);

w1=E*A/L;

w2=12*E*I/(L^3);

w3=6*E*I/(L^2);

w4=4*E*I/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=kprime*T*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((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));

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=theta*pi/180;

C=cos(x);

S=sin(x);

w1=A*C*C+12*I*S*S/(L*L);

w2=A*S*S+12*I*C*C/(L*L);

w3=(A-12*I/(L*L))*C*S;