youxianyuan.rar_位移MATLAB程序

clear NJ=8 %单元个数 Ne=6 %单元总数 XY=[0 0 0.08 0 0 0.04 0.08 0.04 0 0.08 0.08 0.08 0 0.12 0.08 0.12] Code=[1 2 3 3 2 4 3 4 5 5 4 6 5 6 7 7 6 8]; E=2.06e11 %材料参数 E弹性模量 Nu=0.3 %泊松比 t=0.001 %%%计算单元刚度矩阵 D=E/(1-Nu*Nu)*[1 Nu 0;Nu 1 0;0 0 (1-Nu)/2]; Kz=zeros(2*NJ,2*NJ); for e=1:Ne I=Code(e,1); J=Code(e,2); M=Code(e,3); x1=XY(I,1); x2=XY(J,1); x3=XY(M,1); y1=XY(I,2); y2=XY(J,2); y3=XY(M,2); A=0.5*det([1 x1 y1;1 x2 y2;1 x3 y3]); b1=y2-y3;b2=y3-y1;b3=y1-y2; c1=-(x2-x3);c2=x1-x3;c3=x2-x1; B=[b1 0 b2 0 b3 0 0 c1 0 c2 0 c3 c1 b1 c2 b2 c3 b3]/(2*A); Ke=t*A*B'*D*B; %%%%%% Kz(2*I-1:2*I,2*I-1:2*I)=Kz(2*I-1:2*I,2*I-1:2*I)+Ke(1:2,1:2); Kz(2*I-1:2*I,2*J-1:2*J)=Kz(2*I-1:2*I,2*J-1:2*J)+Ke(1:2,3:4); Kz(2*I-1:2*I,2*M-1:2*M)=Kz(2*I-1:2*I,2*M-1:2*M)+Ke(1:2,5:6); %%%%%%%%%%%% Kz(2*J-1:2*J,2*I-1:2*I)=Kz(2*J-1:2*J,2*I-1:2*I)+Ke(3:4,1:2); Kz(2*J-1:2*J,2*J-1:2*J)=Kz(2*J-1:2*J,2*J-1:2*J)+Ke(3:4,3:4); Kz(2*J-1:2*J,2*M-1:2*M)=Kz(2*J-1:2*J,2*M-1:2*M)+Ke(3:4,5:6); %%%%%%%%%%% Kz(2*M-1:2*M,2*I-1:2*I)=Kz(2*M-1:2*M,2*I-1:2*I)+Ke(5:6,1:2); Kz(2*M-1:2*M,2*J-1:2*J)=Kz(2*M-1:2*M,2*J-1:2*J)+Ke(5:6,3:4); Kz(2*M-1:2*M,2*M-1:2*M)=Kz(2*M-1:2*M,2*M-1:2*M)+Ke(5:6,5:6); end Kz %%%Kz整体刚度矩阵 %%%%%%引入载荷F1，f1，f2 %F=Kz*U % F1=zeros(2*NJ,1); % F1(15)=100;F1(16)=50; % %引入约束条件u1=v1=0;u2=v2=0 % Kz(1,:)=0;Kz(:,1)=0;Kz(1,1)=1; % Kz(2,:)=0;Kz(:,2)=0;Kz(2,2)=1; % Kz(3,:)=0;Kz(:,3)=0;Kz(3,3)=1; % Kz(4,:)=0;Kz(:,4)=0;Kz(4,4)=1; % Kz%%% 新的刚度矩阵 % F(1)=0;F(2)=0;F(3)=0;F(4)=0; % F %%新的载荷矩阵 % U=inv(Kz)*F