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
评论0