材料力学求解器-刚架与桁架杆系的计算机求解（附matlab代码）

1 刚架的计算机求解
- 1.1位移法与刚度矩阵
- 1.2 matlab程序
2 桁架的计算机求解

材料力学是一门非常成熟的学科，里面有大量的已经成熟的计算方法、计算模型。然而在掌握了基本原理之后，求解问题的主要时间都花在了计算上面。这时，则完全可以利用计算机采用数值方法进行快速求解。

本文的参考文献很简单，就是单辉祖编写的材料力学第3版下册（或者说是第II册），第18章的内容。

还是惯例声明：本人没有相关的工程应用经验，只是纯粹对相关算法感兴趣才写此博客。所以如果有错误，欢迎在评论区指正，不胜感激。本文主要关注于算法的实现，对于实际应用等问题本人没有任何经验，所以也不再涉及。

1 刚架的计算机求解

1.1位移法与刚度矩阵

计算机求解刚架的原理主要是位移法和刚度矩阵。
基本思路为：

1 定义划分单元与节点
2 列出每个节点的载荷F
3 求解每个单元的刚度矩阵Km
4 把所有单元的刚度矩阵相加，得到总的刚度矩阵K
5 利用位移法与胡克定律F=K*D，求解出位移D
6 根据位移，反推出其余受力

以书上的例18-4，举例说明

1 定义划分单元与节点

上面的刚架，可以很轻易的划分出单元与节点，包含4个节点与3和单元。此外，刚度矩阵的方法要求把受力点算作节点。

2 列出每个节点的载荷F

对于刚架来说，每个节点受力包含Fx、Fy、Me三个力。所以F是一个3*4=12个元素的列向量。
外加的力载荷只有2号节点处的三个力。所以，列向量F可以写作：
[000Fx−Fy−Me000000]′\left[ \begin{matrix} 0&0&0&Fx&-Fy&-Me&0&0&0&0&0&0\\ \end{matrix} \right] ' [000Fx−Fy−Me000000]′
这里’代表转置。
外加力的方向定义为与x、y轴方向相同为正。力矩定义为逆时针为正。

3 求解每个单元的刚度矩阵Km
首先是局部坐标系下，列出刚度矩阵
之后，以③号单元为例，3号单元的两个节点分别为i=2和j=4。
单独某个单元的刚度矩阵Km如下：
[EAL00−EAL00012EIL36EIL20−12EIL36EIL206EIL24EIL0−6EIL22EIL−EAL00EAL000−12EIL3−6EIL2012EIL3−6EIL206EIL22EIL0−6EIL24EIL]\left[ \begin{matrix} \frac{EA}{L}& 0& 0& -\frac{EA}{L}& 0& 0\\ 0& \frac{12EI}{L^3}& \frac{6EI}{L^2}& 0& -\frac{12EI}{L^3}& \frac{6EI}{L^2}\\ 0& \frac{6EI}{L^2}& \frac{4EI}{L}& 0& -\frac{6EI}{L^2}& \frac{2EI}{L}\\ -\frac{EA}{L}& 0& 0& \frac{EA}{L}& 0& 0\\ 0& -\frac{12EI}{L^3}& -\frac{6EI}{L^2}& 0& \frac{12EI}{L^3}& -\frac{6EI}{L^2}\\ 0& \frac{6EI}{L^2}& \frac{2EI}{L}& 0& -\frac{6EI}{L^2}& \frac{4EI}{L}\\ \end{matrix} \right] ⎣⎢⎢⎢⎢⎢⎢⎡LEA00−LEA000L312EIL26EI0−L312EIL26EI0L26EIL4EI0−L26EIL2EI−LEA00LEA000−L312EI−L26EI0L312EI−L26EI0L26EIL2EI0−L26EIL4EI⎦⎥⎥⎥⎥⎥⎥⎤
其中E是弹性模量，A是横截面积，L是单元的长度，I是惯性矩。
之后旋转单元，将单元转换为全局坐标系下
Km=T′∗Km∗TKm=T'*Km*T Km=T′∗Km∗T
其中旋转矩阵T为：
[cos⁡αsin⁡α0000−sin⁡αcos⁡α0000001000000cos⁡αsin⁡α0000−sin⁡αcos⁡α0000001]\left[ \begin{matrix} \cos \alpha& \sin \alpha& 0& 0& 0& 0\\ -\sin \alpha& \cos \alpha& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& \cos \alpha& \sin \alpha& 0\\ 0& 0& 0& -\sin \alpha& \cos \alpha& 0\\ 0& 0& 0& 0& 0& 1\\ \end{matrix} \right] ⎣⎢⎢⎢⎢⎢⎢⎡cosα−sinα0000sinαcosα0000001000000cosα−sinα0000sinαcosα0000001⎦⎥⎥⎥⎥⎥⎥⎤
之后还需要把这个6*6的局部刚度矩阵，拓展到12*12的全局刚度矩阵。方法就是：将局部刚度矩阵的1-3行移动到12*12全局矩阵中的3*i-2到3*i行（也就是4-6行），将4-6行移动到12*12全局矩阵中的3*j-2到3*j行（也就是10-12行）。列也同理。

4 把所有单元的刚度矩阵相加，得到总的刚度矩阵K
这一步就是将所有单元得到的刚度矩阵相加就行。

5 利用位移法与胡克定律F=K*D，求解出位移D
这里用到matlab里的
D=K\F;
就可以直接求解出D。

此外，如果有固定约束（比如固定端约束了x，y，theta三个自由度），需要用到乘大数法。
比如位移D列向量中的第i个分量被约束住了，规定要移动Di，则将刚度矩阵K的i行i列乘以一个很大的数（比如a=10E12），记为a*Kii。并且将载荷向量F中的第i行，更改为a*Kii*Di。

6 根据位移，反推出其余受力
这个可以根据梁弯曲受力的挠度公式，待定系数反求出受力。也可以根据之前得到的局部坐标系下的刚度矩阵，求解受力。

1.2 matlab程序

clear
clc
close all%材料力学求解器
%平面刚架%% 1前处理，划分定义节点与单元，输入初始条件
Node=4;%节点数量
Element=3;%单元数量
%节点位置，每一行分别为x、y
N_xy=[0,1;1,1;2,1;1,0];
%单元属性【节点i,节点j,弹性模量E,横截面积A,惯性矩I】
El_ijEAI=[1 ,2 ,210E9 ,2E3*10^(-3*2) ,3E6*10^(-3*4) ;2 ,3 ,210E9 ,2E3*10^(-3*2) ,3E6*10^(-3*4) ;2 ,4 ,210E9 ,2E3*10^(-3*2) ,3E6*10^(-3*4) ;];
%定义节点处施加的力的大小，顺序为Fx1,Fy1,Me1, Fx2,Fy2,Me2, Fx3,...
%这里例子中，F=[0;0;0;20E3;-20E3;-40E3;0;0;0;0;0;0];
F=zeros(3*Node,1);
F(4)=20E3;
F(5)=-20E3;
F(6)=-40E3;
%定义位移约束，顺序为X1,Y1,T1, X2,Y2,T2, X3,...
%后面为约束具体的值
%这里例子中，1,2,3,7,8,9,10,11,12都是0位移约束
Move_Limit=[1,2,3,7,8,9,10,11,12;0,0,0,0,0,0,0,0,0];%% 2采用刚度矩阵方法计算D=SolveStiffnessMethod(Node,Element,N_xy,El_ijEAI,F,Move_Limit);%% 3后处理%3.1求解变形后的图形
Scale_X=1000;%位移变形放大系数(宁小勿大）
Scale_T=0.5/max(abs( D(3:3:3*Node) ));%角度变形放大系数(宁小勿大）
%Scale_T=100;%角度变形放大系数(宁小勿大）
figure(1)
hold on
%绘制变形前的图像
for m=1:Element%原先节点的坐标位置xyI=N_xy(El_ijEAI(m,1),:);xyJ=N_xy(El_ijEAI(m,2),:);plot([xyI(1),xyJ(1)],[xyI(2),xyJ(2)],'b')
end
%绘制变形后的图像
for m=1:ElementNode_I=El_ijEAI(m,1);Node_J=El_ijEAI(m,2);%原先节点的坐标位置xyI=N_xy(Node_I,:);xyJ=N_xy(Node_J,:);%变形后的节点坐标位置xyIJ_New=[xyI'+Scale_X*D(3*Node_I-2:3*Node_I-1),xyJ'+Scale_X*D(3*Node_J-2:3*Node_J-1)];%原先的节点对应的夹角Theta_IJ=angle((xyJ(1)-xyI(1))+(xyJ(2)-xyI(2))*1i);%变形后的节点对应的角度DirI=[cos( Theta_IJ+Scale_T*D(3*Node_I) ) ; sin( Theta_IJ+Scale_T*D(3*Node_I) )];DirJ=[cos( Theta_IJ+Scale_T*D(3*Node_J) ) ; sin( Theta_IJ+Scale_T*D(3*Node_J) )];%生成样条曲线cs3 = spline([1,2],[DirI,xyIJ_New,DirJ]);yy=ppval(cs3,linspace(1,2));%绘制plot(yy(1,:),yy(2,:),'r')
end
hold off
axis equal%3.2求解变形后的弯矩图
Scale_M=1E-5;%弯矩变形放大系数(默认1）
figure()
hold on
%绘制原先的图像
for m=1:Element%原先节点的坐标位置xyI=N_xy(El_ijEAI(m,1),:);xyJ=N_xy(El_ijEAI(m,2),:);plot([xyI(1),xyJ(1)],[xyI(2),xyJ(2)],'b')
end
%绘制力矩图
for m=1:Element%旋转杆件，回到局部坐标系Node_I=El_ijEAI(m,1);Node_J=El_ijEAI(m,2);E=El_ijEAI(m,3);A=El_ijEAI(m,4);I=El_ijEAI(m,5);xyI=N_xy(Node_I,:);xyJ=N_xy(Node_J,:);L=sqrt(sum((xyI-xyJ).^2));Theta_IJ=angle((xyJ(1)-xyI(1))+(xyJ(2)-xyI(2))*1i);%计算水平梁的变形T=RotateMatrixT(xyI,xyJ);D_IJ=[D(3*Node_I-2:3*Node_I);D(3*Node_J-2:3*Node_J)];DT=T'*D_IJ;%根据弯矩公式反推弯矩Me_abcd=[0,0,1,0;0.5*L^2,L,1,0;0,0,0,1;L^3/6,L^2/2,L,1]/E/I;abcd=Me_abcd\[DT(3);DT(6);DT(4);DT(5)];%绘图x1=linspace(0,L);y1=abcd(1)*x1+abcd(2);y1=Scale_M*y1;%旋转xyxy2=[cos(Theta_IJ),-sin(Theta_IJ);sin(Theta_IJ),cos(Theta_IJ)]*[x1;y1];%将0点平移到xyI上xy2=xy2+xyI';%plot(xy2(1,:),xy2(2,:),'r')patch([xyI(1),xy2(1,:),xyJ(1)],[xyI(2),xy2(2,:),xyJ(2)],'k','FaceAlpha',0.5)
end
hold off
axis equal
%3.3 输出每个单元的受力，每个单元计算6个力，分别为I端点的Fx,Fy,Me和J端点的Fx,Fy,Me
FxFyMeIJ=zeros(Element,6);
for m=1:Element
Node_I=El_ijEAI(m,1);
Node_J=El_ijEAI(m,2);
E=El_ijEAI(m,3);
A=El_ijEAI(m,4);
I=El_ijEAI(m,5);
xyI=N_xy(Node_I,:);
xyJ=N_xy(Node_J,:);
L=sqrt(sum((xyI-xyJ).^2));
Theta_IJ=angle((xyJ(1)-xyI(1))+(xyJ(2)-xyI(2))*1i);
%旋转杆件，回到局部坐标系
T=RotateMatrixT(xyI,xyJ);
D_IJ=[D(3*Node_I-2:3*Node_I);D(3*Node_J-2:3*Node_J)];
DT=T*D_IJ;
%计算受力
FxI=-E*A/L*DT(4)+E*A/L*DT(1);
FxJ=E*A/L*DT(4)-E*A/L*DT(1);
FyI=-12*E*I/L^3*DT(5)+6*E*I/L^2*DT(6)+12*E*I/L^3*DT(2)+6*E*I/L^2*DT(3);
FyJ=12*E*I/L^3*DT(5)-6*E*I/L^2*DT(6)-12*E*I/L^3*DT(2)-6*E*I/L^2*DT(3);
MeI=-6*E*I/L^2*DT(5)+2*E*I/L*DT(6)+6*E*I/L^2*DT(2)+4*E*I/L*DT(3);
MeJ=-6*E*I/L^2*DT(5)+4*E*I/L*DT(6)+6*E*I/L^2*DT(2)+2*E*I/L*DT(3);
FxFyMeIJ(m,:)=[FxI,FyI,MeI,FxJ,FyJ,MeJ];
end
FxFyMeIJ( abs(FxFyMeIJ) < (max(abs(FxFyMeIJ))/1E10) )=0;
D( abs(D) < (max(abs(D))/1E10) )=0;
%输出
Str_Output1=char(zeros(Node+1,60)+32);%32\12288
Str_Output1(1,1:3)='节 点';
Str_Output1(1,14:17)='水平位移';
Str_Output1(1,34:37)='竖直位移';
Str_Output1(1,54:57)='转  角';
for m=2:Node+1Str_Output1(m,2:3)=num2str(m-1,'%02.f');Str_Output1(m,4:5)=char(12288);Str_Output1(m,12:23)=num2str(D(3*m-5),'%+10.5e');Str_Output1(m,24:28)=char(12288);Str_Output1(m,30:41)=num2str(D(3*m-4),'%+10.5e');Str_Output1(m,42:46)=char(12288);Str_Output1(m,48:59)=num2str(D(3*m-3),'%+10.5e');
end
Str_Output1Str_Output2=char(zeros(2*Element+1,60)+32);
Str_Output2(1,1:3)='单 元';
Str_Output2(1,14:16)='横向力';
Str_Output2(1,34:36)='切向力';
Str_Output2(1,54:56)='力 矩';
for m=1:ElementStr_Output2(2*m,2:3)=num2str(m-1,'%02.f');Str_Output2(2*m,4:5)=char(12288);Str_Output2(2*m,12:23)=num2str(FxFyMeIJ(m,1),'%+10.5e');Str_Output2(2*m,24:28)=char(12288);Str_Output2(2*m,30:41)=num2str(FxFyMeIJ(m,2),'%+10.5e');Str_Output2(2*m,42:46)=char(12288);Str_Output2(2*m,48:59)=num2str(FxFyMeIJ(m,3),'%+10.5e');Str_Output2(2*m+1,4:5)=char(12288);Str_Output2(2*m+1,12:23)=num2str(FxFyMeIJ(m,4),'%+10.5e');Str_Output2(2*m+1,24:28)=char(12288);Str_Output2(2*m+1,30:41)=num2str(FxFyMeIJ(m,5),'%+10.5e');Str_Output2(2*m+1,42:46)=char(12288);Str_Output2(2*m+1,48:59)=num2str(FxFyMeIJ(m,6),'%+10.5e');
end
Str_Output2%% 一些自定义函数
function T=RotateMatrixT(xyI,xyJ)
xI=xyI(1);
yI=xyI(2);
xJ=xyJ(1);
yJ=xyJ(2);
L=sqrt((xJ-xI)^2+(yJ-yI)^2);
T=eye(6);
T(1,1)=(xJ-xI)/L;T(1,2)=(yJ-yI)/L;T(2,1)=-(yJ-yI)/L;T(2,2)=(xJ-xI)/L;
T(4,4)=(xJ-xI)/L;T(4,5)=(yJ-yI)/L;T(5,4)=-(yJ-yI)/L;T(5,5)=(xJ-xI)/L;
endfunction D=SolveStiffnessMethod(Node,Element,N_xy,El_ijEAI,F,Move_Limit)
K_El=zeros(3*Node,3*Node,Element);
for m=1:Element%每个单元的属性Node_I=El_ijEAI(m,1);Node_J=El_ijEAI(m,2);E=El_ijEAI(m,3);A=El_ijEAI(m,4);I=El_ijEAI(m,5);%两端节点的坐标xyI=N_xy(Node_I,:);xyJ=N_xy(Node_J,:);L=sqrt(sum((xyI-xyJ).^2));%每个单元对应的刚度矩阵%计算系数KFxu=E*A/L;KFyv=12*E*I/L^3;KFyt=6*E*I/L^2;KMev=6*E*I/L^2;KMet=4*E*I/L;%刚度矩阵%Node_I=1;Node_J=2;K_El_m=zeros(6);K_El_m(1,1)=KFxu;K_El_m(2,2)=KFyv;K_El_m(2,3)=KFyt;K_El_m(3,2)=KMev;K_El_m(3,3)=KMet;K_El_m(1,4)=-KFxu;K_El_m(2,5)=-KFyv;K_El_m(2,6)=KFyt;K_El_m(3,5)=-KMev;K_El_m(3,6)=KMet;K_El_m(4,1)=-KFxu;K_El_m(5,2)=-KFyv;K_El_m(5,3)=-KFyt;K_El_m(6,2)=KMev;K_El_m(6,3)=KMet;K_El_m(4,4)=KFxu;K_El_m(5,5)=KFyv;K_El_m(5,6)=-KFyt;K_El_m(6,5)=-KMev;K_El_m(6,6)=KMet;%旋转杆件T=RotateMatrixT(xyI,xyJ);K_El_m_T=T'*K_El_m*T;%扩展到全局方程中，II、IJ、JI、JJ节点方程与K_El_m_T中的相同，其余全是0K_El_m_S=zeros(3*Node,3*Node);K_El_m_S(Node_I*3-2:Node_I*3-0,Node_I*3-2:Node_I*3-0)=K_El_m_T(1:3,1:3);K_El_m_S(Node_I*3-2:Node_I*3-0,Node_J*3-2:Node_J*3-0)=K_El_m_T(1:3,4:6);K_El_m_S(Node_J*3-2:Node_J*3-0,Node_I*3-2:Node_I*3-0)=K_El_m_T(4:6,1:3);K_El_m_S(Node_J*3-2:Node_J*3-0,Node_J*3-2:Node_J*3-0)=K_El_m_T(4:6,4:6);%保存到K_El里K_El(:,:,m)=K_El_m_S;
end%整体的包含各个元素的刚度矩阵
K=sum(K_El,3);a=1E12;%放大系数
%利用乘大数法，模拟位移约束的影响
for m=1:size(Move_Limit,2)n=Move_Limit(1,m);K(n,n)=a*K(n,n);F(n)=K(n,n)*Move_Limit(2,m);
end
%求解位移结果
D=K\F;
end

输出得到的变形图如下。这里变形为了展示效果有所夸张，实际变形基本看不出来。

输出的力矩分布图如下：

输出的计算结果如下：

'节 点          水平位移                竖直位移                转  角   '' 01　　      +0.00000e+00　　　　　 +0.00000e+00　　　　　 +0.00000e+00 '' 02　　      +4.72998e-05　　　　　 -4.59643e-05　　　　　 -5.31466e-03 '' 03　　      +0.00000e+00　　　　　 +0.00000e+00　　　　　 +0.00000e+00 '' 04　　      +0.00000e+00　　　　　 +0.00000e+00　　　　　 +0.00000e+00 ''单 元          横向力                 切向力                 力 矩    '' 00　　      -1.98659e+04　　　　　 -1.97419e+04　　　　　 -6.52272e+03 ''   　　      +1.98659e+04　　　　　 +1.97419e+04　　　　　 -1.32192e+04 '' 01　　      +1.98659e+04　　　　　 -2.04369e+04　　　　　 -1.35667e+04 ''   　　      -1.98659e+04　　　　　 +2.04369e+04　　　　　 -6.87021e+03 '' 02　　      +1.93050e+04　　　　　 -1.97318e+04　　　　　 -1.32141e+04 ''   　　      -1.93050e+04　　　　　 +1.97318e+04　　　　　 -6.51767e+03 '

2 桁架的计算机求解

求解思路与之前的刚架是一致的。因为这里只涉及到二力杆，所以刚度方程简化了不少。

求解顺序还是如下：

1 定义划分单元与节点
2 列出每个节点的载荷F
3 求解每个单元的刚度矩阵Km
4 把所有单元的刚度矩阵相加，得到总的刚度矩阵K
5 利用位移法与胡克定律F=K*D，求解出位移D
6 根据位移，反推出其余受力

1 定义划分单元与节点

这个只需要把每一个杆定义为一个单元即可。

2 列出每个节点的载荷F

对于二力杆来说，每个节点受力包含Fx、Fy两个力，不存在力矩。相比较刚架少了一个力，总的力只有2m个。

外加力的方向定义为与x、y轴方向相同为正。

3 求解每个单元的刚度矩阵Km
这里局部坐标系的刚度矩阵变成了4*4的矩阵。分别对应Fxi,Fyi,Fxj,Fyj。
[EAL0−EAL00000−EAL0EAL00000]\left[ \begin{matrix} \frac{EA}{L}& 0& -\frac{EA}{L}& 0 \\ 0& 0& 0& 0 \\ -\frac{EA}{L}& 0& \frac{EA}{L}& 0\\ 0& 0& 0& 0\\ \end{matrix} \right] ⎣⎢⎢⎡LEA0−LEA00000−LEA0LEA00000⎦⎥⎥⎤
其中E是弹性模量，A是横截面积，L是单元的长度。
之后旋转单元，将单元转换为全局坐标系下
Km=T′∗Km∗TKm=T'*Km*T Km=T′∗Km∗T
其中旋转矩阵T为：
[cos⁡αsin⁡α00−sin⁡αcos⁡α0000cos⁡αsin⁡α00−sin⁡αcos⁡α]\left[ \begin{matrix} \cos \alpha& \sin \alpha& 0& 0\\ -\sin \alpha& \cos \alpha& 0& 0\\ 0& 0& \cos \alpha& \sin \alpha\\ 0& 0& -\sin \alpha& \cos \alpha\\ \end{matrix} \right] ⎣⎢⎢⎡cosα−sinα00sinαcosα0000cosα−sinα00sinαcosα⎦⎥⎥⎤
之后还需要把这个4*4的局部刚度矩阵，拓展到全局刚度矩阵。方法同之前的刚架。

4 把所有单元的刚度矩阵相加，得到总的刚度矩阵K
这一步就是将所有单元得到的刚度矩阵相加就行。

5 利用位移法与胡克定律F=K*D，求解出位移D
这里用到matlab里的
D=K\F;
就可以直接求解出D。

此外，如果有固定约束（比如固定端约束了x，y两个自由度），需要用到乘大数法。
比如位移D列向量中的第i个分量被约束住了，规定要移动Di，则将刚度矩阵K的i行i列乘以一个很大的数（比如a=10E12），记为a*Kii。并且将载荷向量F中的第i行，更改为a*Kii*Di。

6 根据位移，反推出其余受力
对于二力杆来说，方程很简单，只需要带入到局部坐标系中求解出压力就行。

还是以书上例题为例：

matlab代码如下：

clear
clc
close all%材料力学求解器
%平面二力杆系%% 1前处理，划分定义节点与单元，输入初始条件
Node=4;%节点数量
Element=6;%单元数量
%节点位置，每一行分别为x、y
N_xy=[0,0;0,2;1,2;1,0];
%单元属性【节点i,节点j,弹性模量E,横截面积A】
El_ijEA=[1 ,2 ,200E9 ,100*10^(-3*2) ;2 ,3 ,200E9 ,100*10^(-3*2) ;3 ,4 ,200E9 ,100*10^(-3*2) ;1 ,4 ,200E9 ,100*10^(-3*2) ;1 ,3 ,200E9 ,100*10^(-3*2) ;2 ,4 ,200E9 ,100*10^(-3*2) ;];
%定义节点处施加的力的大小，顺序为Fx1,Fy1, Fx2,Fy2, Fx3,...
F=zeros(2*Node,1);
F(4)=-3E3;
F(5)= 1E3;
F(6)=-2E3;
%定义位移约束，顺序为X1,Y1, X2,Y2, X3,...
%后面为约束具体的值
%这里例子中，1,2,8都是0位移约束
Move_Limit=[1,2,8;0,0,0];
%检测
if size(N_xy,1)~=Node || size(El_ijEA,1)~=Elementerror('输入有误')
end
%% 2采用刚度矩阵方法计算
D=SolveStiffnessMethod(Node,Element,N_xy,El_ijEA,F,Move_Limit);%% 3后处理
%3.1计算轴力
F_IJ=zeros(Element,1);
for m=1:Element
Node_I=El_ijEA(m,1);
Node_J=El_ijEA(m,2);
E=El_ijEA(m,3);
A=El_ijEA(m,4);
xyI=N_xy(Node_I,:);
xyJ=N_xy(Node_J,:);
L=sqrt(sum((xyI-xyJ).^2));
Theta_IJ=angle((xyJ(1)-xyI(1))+(xyJ(2)-xyI(2))*1i);
%旋转杆件，回到局部坐标系
T=RotateMatrixT(xyI,xyJ);
D_IJ=[D(2*Node_I-1:2*Node_I);D(2*Node_J-1:2*Node_J)];
DT=T*D_IJ;
%计算受力
FxI=-E*A/L*DT(1)+E*A/L*DT(3);F_IJ(m,:)=FxI;
end%3.2求解变形后的图形
Scale_X=200;%位移变形放大系数(宁小勿大）
%绘制颜色条
N_color=6;
mcp=[[linspace(1,1/(N_color-1),N_color-1)',zeros(N_color-1,1),zeros(N_color-1,1)];...[zeros(N_color,1),zeros(N_color,1),linspace(0,1,N_color)']];
mcp=flipud(mcp);
Max_F=max(abs(F_IJ));
ColorIndex=round(F_IJ/Max_F*(N_color-1))+N_color;
%颜色索引
figure(1)
hold on
%绘制变形前的图像
for m=1:Element%原先节点的坐标位置xyI=N_xy(El_ijEA(m,1),:);xyJ=N_xy(El_ijEA(m,2),:);plot([xyI(1),xyJ(1)],[xyI(2),xyJ(2)],'color',[0.5,0.5,0.5])
end
%绘制变形后的图像
for m=1:ElementNode_I=El_ijEA(m,1);Node_J=El_ijEA(m,2);%原先节点的坐标位置xyI=N_xy(Node_I,:);xyJ=N_xy(Node_J,:);%变形后的节点坐标位置xyIJ_New=[xyI'+Scale_X*D(2*Node_I-1:2*Node_I),xyJ'+Scale_X*D(2*Node_J-1:2*Node_J)];%绘制plot(xyIJ_New(1,:),xyIJ_New(2,:),'color',mcp(ColorIndex(m),:),'linewidth',1.5)end
hold off
axis equal
caxis([-Max_F,Max_F])
colormap(mcp)
colorbar()%% 一些自定义函数
function T=RotateMatrixT(xyI,xyJ)
xI=xyI(1);
yI=xyI(2);
xJ=xyJ(1);
yJ=xyJ(2);
L=sqrt((xJ-xI)^2+(yJ-yI)^2);
T=eye(4);
T(1,1)=(xJ-xI)/L;T(1,2)=(yJ-yI)/L;T(2,1)=-(yJ-yI)/L;T(2,2)=(xJ-xI)/L;
T(3,3)=(xJ-xI)/L;T(3,4)=(yJ-yI)/L;T(4,3)=-(yJ-yI)/L;T(4,4)=(xJ-xI)/L;
endfunction D=SolveStiffnessMethod(Node,Element,N_xy,El_ijEA,F,Move_Limit)
K_El=zeros(2*Node,2*Node,Element);
for m=1:Element%每个单元的属性Node_I=El_ijEA(m,1);Node_J=El_ijEA(m,2);E=El_ijEA(m,3);A=El_ijEA(m,4);%两端节点的坐标xyI=N_xy(Node_I,:);xyJ=N_xy(Node_J,:);L=sqrt(sum((xyI-xyJ).^2));%每个单元对应的刚度矩阵%计算系数KFxu=E*A/L;%刚度矩阵%Node_I=1;Node_J=2;K_El_m=zeros(4);K_El_m(1,1)=KFxu;K_El_m(1,3)=-KFxu;K_El_m(3,1)=-KFxu;K_El_m(3,3)=KFxu;%旋转杆件T=RotateMatrixT(xyI,xyJ);K_El_m_T=T'*K_El_m*T;%扩展到全局方程中，II、IJ、JI、JJ节点方程与K_El_m_T中的相同，其余全是0K_El_m_S=zeros(2*Node,2*Node);K_El_m_S(Node_I*2-1:Node_I*2-0,Node_I*2-1:Node_I*2-0)=K_El_m_T(1:2,1:2);K_El_m_S(Node_I*2-1:Node_I*2-0,Node_J*2-1:Node_J*2-0)=K_El_m_T(1:2,3:4);K_El_m_S(Node_J*2-1:Node_J*2-0,Node_I*2-1:Node_I*2-0)=K_El_m_T(3:4,1:2);K_El_m_S(Node_J*2-1:Node_J*2-0,Node_J*2-1:Node_J*2-0)=K_El_m_T(3:4,3:4);%保存到K_El里K_El(:,:,m)=K_El_m_S;
end%整体的包含各个元素的刚度矩阵
K=sum(K_El,3);a=1E11;%放大系数
%利用乘大数法，模拟位移约束的影响
for m=1:size(Move_Limit,2)n=Move_Limit(1,m);K(n,n)=a*K(n,n);F(n)=K(n,n)*Move_Limit(2,m);
end
%求解位移结果
D=K\F;
end

输出的图像结果：

变形同样被放大了，实际变形几乎看不出来。颜色显示了每根杆的受力情况。

材料力学求解器-刚架与桁架杆系的计算机求解（附matlab代码）相关推荐

【路径规划】基于遗传算法求解多车多类型车辆的车辆路径优化问题附matlab代码
1 内容介绍多车辆多路线的交通路线优化涉及到排序问题,是一个N-P难题,高效精确的算法存在的可能性不大.提出了基于遗传算法的求解方法,给出了实例来证明如何利用遗传算法解决多车辆多路线的优化问题.结果 ...
【智能优化求解】基于粒子群算法实现综合能源系统优化附matlab代码
1 简介为了解决现有冷热电联供型综合能源系统大多只单一考虑系统机组投资成本或系统环境污染,影响系统整体优化运行的问题,以系统经济性和环保性为目标,对冷热电联供系统进行研究分析.构建含燃气轮机.燃气锅 ...
【路径规划】基于遗传算法求解带时间窗多电动车充电路径规划问题附matlab代码
1 简介电动车在物流领域中取代燃油车是一个广泛的发展趋势.但电动车的电池利用率低,充电时间长,相关充电配套设施建设不完善,存在"续驶里程焦虑"等现象成为了电动车推广和应用的重要制 ...
CST入门——求解器简介与时域、频域和积分求解器设置
目录 1. 高频电磁仿真求解器简介 1.1. 时域求解器 Time Domain Solver(主) 1.2. 频域求解器 Frequency Domain Solver(主) 1.3. 本征模求解器 ...
惩罚函数外点matlab,禁忌搜索算法求解带时间窗的车辆路径问题(惩罚函数版附MATLAB代码)...
本周应小伙伴要求继续学习TS求VRPTW,不过这次通过使用惩罚约束的形式允许解违反时间窗约束和容量约束,不过要给违反约束的解加以惩罚. 这次我们的目标函数就不单单只有车辆总行驶距离了,还要包括当前解中 ...
智能优化算法-蚁狮优化器Ant Lion Optimizer(附Matlab代码)
引言蚁狮优化器(Ant Lion Optimizer,ALO)模拟了自然界蚁狮的捕猎机制.实施了蚁群随机行走.设置陷阱.设陷阱.捕捉猎物.重建陷阱等5个主要捕猎步骤.于2015年发表在Seyedal ...
【智能优化算法-灰狼算法】基于贪婪非分级灰狼优化器求解单目标优化问题附matlab代码
1 内容介绍灰狼优化(GWO)算法是一种新兴的算法,它基于灰狼的社会等级以及它们的狩猎和合作策略. 该算法于 2014 年推出,已被大量研究人员和设计人员使用,原始论文的引用次数超过了许多其他算法. ...
【智能优化算法-蒲公英优化器】基于蒲公英优化器求解单目标优化问题附matlab代码
1 内容介绍群智能优化算法作为当前优化算法中的一个主要研究热点,经过近年的发展,已经发展为较为新颖的演化计算技术,受到越来越多不同领域研究工作者的关注.群智能优化算法比传统优化方法求解各种复杂优化问 ...
【智能优化算法】基于沙猫群优化算法求解单目标优化问题附matlab代码
1 内容介绍这项研究提出了一种新的元启发式算法,称为沙猫群优化 (SCSO),它模仿试图在自然界中生存的沙猫行为.这些猫能够探测到低于 2 kHz 的低频,并且具有难以置信的挖掘猎物的能力.受这两个 ...

材料力学求解器-刚架与桁架杆系的计算机求解（附matlab代码）

材料力学求解器-刚架与桁架杆系的计算机求解（附matlab代码）

1 刚架的计算机求解

1.1位移法与刚度矩阵

1.2 matlab程序

2 桁架的计算机求解

材料力学求解器-刚架与桁架杆系的计算机求解（附matlab代码）相关推荐

最新文章

热门文章