75 lines
2.7 KiB
Matlab
75 lines
2.7 KiB
Matlab
clc,clear,close all
|
||
tic
|
||
nodeCoord=[0 0;0 3;0 6 ;0 9; 5 0;5 3;5 6;5 9];
|
||
EleNode=[1 2 ;2 3; 3 4 ;5 6; 6 7; 7 8 ;2 6 ;3 7;4 8];
|
||
x=nodeCoord(:,1)'; % 节点x轴方向坐标
|
||
y=nodeCoord(:,2)'; % 节点y轴方向坐标
|
||
A=0.08;E=3e10;I=0.0128/12; % 定义截面面积和弹性模量
|
||
%单元信息：编号，节点1编号，节点2编号，截面面积，弹性模量,惯性矩
|
||
ele=[1 1 2 A E I;2 2 3 A E I; 3 3 4 A E I;4 5 6 A E I; 5 6 7 A E I;6 7 8 A E I;7 2 6 A E I;8 3 7 A E I;9 4 8 A E I];;
|
||
%载荷信息：节点编号，x方向力，y方向力，z方向力
|
||
Load=[4 2e5 0 0];
|
||
%约束：节点编号，x方向约束，y方向约束，z方向约束
|
||
Constr=[1 0 0 0;5 0 0 0];
|
||
Dofs=3*size(x,2); %总自由度数
|
||
EleCount=size(ele,1); %单元总数
|
||
K=zeros(Dofs,Dofs); %总体刚度矩阵
|
||
F=zeros(Dofs,1); %总体载荷列阵
|
||
U=zeros(Dofs,1); %总体位移列阵
|
||
BarLength=BarsLength(x,y,ele);
|
||
figure('Name','Undeformed Truss')
|
||
RenderFrame(ele,Load,Constr,x,y,U,1,'-k') %绘制框架
|
||
figure('Name','Deformed Truss')
|
||
RenderFrame(ele,Load,Constr,x,y,U,0.5,'-k') %绘制未变形后框架
|
||
hold on
|
||
%遍历所有单元，将各单元刚度阵分块组装到总体刚度阵
|
||
for iEle =1:EleCount
|
||
%该单元的两个节点的编号
|
||
n1=ele(iEle,2);n2=ele(iEle,3);
|
||
%计算坐标变换矩阵
|
||
R=CoordTransform([x(n1) x(n2)],[y(n1) y(n2)],BarLength(iEle));
|
||
%计算单元刚度矩阵 Ke=R'*ke*R;局部坐标系下的单元刚度阵转换为全局坐标下的单元刚度阵
|
||
ke= FrameElementKe(ele(iEle,4),ele(iEle,5),ele(iEle,6),R,BarLength(iEle));
|
||
%将各单元刚度分块组装到总刚相应位置
|
||
eleDof=[n1*3-2:n1*3,n2*3-2:n2*3];
|
||
K(eleDof,eleDof)=K(eleDof,eleDof)+ke;
|
||
|
||
end
|
||
%形成载荷列阵--1、2、3、4自由度赋载荷值
|
||
for LoadNum=1:size(Load,1)
|
||
for i=1:3
|
||
F(3*Load(LoadNum)+i-3,1)=Load(LoadNum,i+1);
|
||
end
|
||
end
|
||
%施加约束--乘大数法
|
||
for iConstr=1:size(Constr,1)
|
||
for j=2:4
|
||
if ~isnan(Constr(iConstr,j))
|
||
K(3*Constr(iConstr,1 )+j-4,3*Constr(iConstr,1 )+j-4)=...
|
||
1e12*K(3*Constr(iConstr,1 )+j-4,3*Constr(iConstr,1 )+j-4);
|
||
F(3*Constr(iConstr,1)+j-4)=Constr(iConstr,j)*...
|
||
K(3*Constr(iConstr,1 )+j-4,3*Constr(iConstr,1)+j-4);
|
||
end
|
||
end
|
||
end
|
||
U=K\F; %全局坐标系下位移,此处原来为逆但是矩阵奇异，故改成伪逆
|
||
%%%%%%%%%%%%
|
||
for iEle =1:EleCount
|
||
%计算杆局部坐标下的位移
|
||
n1=ele(iEle,2);n2=ele(iEle,3);
|
||
R=CoordTransform([x(n1) x(n2)],[y(n1) y(n2)],BarLength(iEle));
|
||
localU = R*[U(3*n1-2:3*n1,1);U(3*n2-2:3*n2,1)];
|
||
[Ke_local] = FrameElementKeLocal(ele(iEle,4),ele(iEle,5),ele(iEle,6),R,BarLength(iEle));
|
||
EleForce(:, iEle)=Ke_local*localU; %力
|
||
end
|
||
|
||
%保存位移，应力与轴力到文本文件
|
||
|
||
MomentForce=EleForce([3,6],:);
|
||
MomentForce(2,:)=-MomentForce(2,:);%体现出正负弯矩的效果
|
||
RenderFrame(ele,Load,Constr,x,y,U,1,'-.b') %绘制变形后的框架
|
||
figure
|
||
ForcePlot(ele,Load,Constr,x,y,U,3,'-b',MomentForce)%绘制弯矩图
|
||
toc
|
||
|