目录
1、此文是基于matlab的受轴弯剪的梁单元的结构求解,代码如下:
1、此文是基于matlab的受轴弯剪的梁单元的结构求解,代码如下:
% 基于matlab的受轴弯剪的梁单元的结构求解;
% 输入量为:节点数目,坐标; 单元节点编号,*度编号; 已知力,位移; 未知位移;
clc
clear all
num_nod=3; %节点数目;
ele_nod=[1 2;2 3]; %单元节点编号;
nod_coor=[0 0;12 0; 24.73 12.73] ; %节点坐标,向右为x轴,向下为y轴;
ele_dof=[1 2 3 4 5 6;4 5 6 7 8 9]; %单元*度编号,对应于子块搬家对应的家;
force=zeros(3*num_nod,1);
displacement=zeros(3*num_nod,1);
stiffness=zeros(3*num_nod);
ele_num=size(ele_nod,1);
force(5) = 1000000000.0; %已知节点力;
displacement(1)=0.0; %约束条件;
displacement(2)=0.0;
displacement(3)=0.0;
displacement(7)=0.0;
displacement(8)=0.0;
displacement(9)=0.0;
unknown_f_a=[4;5;6]; %未知节点位移;
% 截面基本参数
b=0.5;
h=0.5;
E=2.1e11;
%初使化所要的矩阵,F K 位移矩阵;
A=b*h;
I=b*h^3/12;
for e = 1 : ele_num
L(e) = sqrt( (nod_coor(ele_nod(e,1),1)-nod_coor(ele_nod(e,2),1))^2+ ...
(nod_coor(ele_nod(e,1),2)-nod_coor(ele_nod(e,2),2))^2 );
C=(nod_coor(ele_nod(e,2),1)-nod_coor(ele_nod(e,1),1))/L(e);
S=(nod_coor(ele_nod(e,2),2)-nod_coor(ele_nod(e,1),2))/L(e);
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]; %坐标转换矩阵
k11=A*E/L(e);
k22=12*E*I/L(e)^3;
k23=6*E*I/L(e)^2;
k33=4*E*I/L(e);
k36=2*E*I/L(e);
k_local = [k11 0 0 -k11 0 0;0 k22 k23 0 -k22 k23;0 k23 k33 0 -k23 k36;-k11 0 0 k11 0 0;0 -k22 -k23 0 k22 -k23; 0 k23 k36 0 -k23 k33];
k=t'*k_local*t;
ele_dof_a = ele_dof(e,:);
for i = 1 : 6
for j = 1 : 6
stiffness(ele_dof_a(1,i),ele_dof_a(1,j))=...
stiffness(ele_dof_a(1,i),ele_dof_a(1,j))+k(i,j);
end
end
end
%因为此时的K是奇异的,那么要消除对应的值;
for i = 1:size(unknown_f_a)
dis_new(i,1)=displacement(unknown_f_a(i,1),1);
force_new(i,1)=force(unknown_f_a(i,1),1);
end
for i = 1:size(unknown_f_a)
for j =1:size(unknown_f_a)
stiff_new(i,j) = stiffness(unknown_f_a(i,1),unknown_f_a(j,1));
end
end
%消除K的奇异性后,求解所对应*度处的位移;
dis_new = inv(stiff_new)*force_new;
for i = 1 : size(unknown_f_a)
fprintf('对应的第%d个*度的未知位移大小为 %d;\n',unknown_f_a(i),dis_new(i))
end
for i = 1 :size(unknown_f_a)
displacement(unknown_f_a(i,1),1) = dis_new(i,1);
end
for e = 1 :ele_num
x = [nod_coor(ele_nod(e,1),1) nod_coor(ele_nod(e,2),1)];
y = [nod_coor(ele_nod(e,1),2) nod_coor(ele_nod(e,2),2)];
plot(x,-y,'b')
hold on
end
for i=1:num_nod
nod_coor_def(i,1)=nod_coor(i,1)+displacement(3*i-2,1);
nod_coor_def(i,2)=nod_coor(i,2)+displacement(3*i-1,1);
end
for e=1:ele_num
x=[nod_coor_def(ele_nod(e,1),1) nod_coor_def(ele_nod(e,2),1)];
y=[nod_coor_def(ele_nod(e,1),2) nod_coor_def(ele_nod(e,2),2)];
plot(x,-y,'r')
hold on
end