189 8069 5689
你提出的问题我之前刚好做过,使用有限元方法来进行桁架结构分析.
目前创新互联建站已为上1000家的企业提供了网站建设、域名、网站空间、绵阳服务器托管、企业网站设计、丰台网站维护等服务，公司将坚持客户导向、应用为本的策略，正道将秉承"和谐、参与、激情"的文化，与客户和合作伙伴齐心协力一起成长，共同发展。
Matlab编程实现平面杆单元分析
首先,明确Matlab程序要实现的5个重要模块分别为：单元刚度矩阵的求解、单元组装、节点位移的求解、单元应力的求解、节点力的求解.下面给出这5个模块的实现.
1．\x09单元刚度矩阵求解
定义函数Bar2D2Node_Stiffness,该函数计算单元的刚度矩阵,输入弹性模量E,横截面积A,两个节点坐标输出单元刚度矩阵k(4X4).具体代码如下：
function k=Bar2D2Node_Stiffness(E,A,x1,y1,x2,y2)
L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
x=acos((x2-x1)/L);
C=cos(x);
S=sin(x);
k=E*A/L*[C*C C*S -C*C -C*S; C*S S*S -C*S -S*S;
-C*C -C*S C*C C*S; -C*S -S*S C*S S*S];
2．\x09单元组装
定义函数Bar2D2Node_Assembly,该函数进行单元刚度矩阵的组装,输入单元刚度矩阵k,单元的节点编号i、j.输出整体刚度矩阵KK,具体代码如下：
function z = Bar2D2Node_Assembly(KK,k,i,j)
DOF(1)=2*i-1;
DOF(2)=2*i;
DOF(3)=2*j-1;
DOF(4)=2*j;
for n1=1:4
for n2=1:4
KK(DOF(n1),DOF(n2))= KK(DOF(n1),DOF(n2))+k(n1,n2);
end
end
z=KK;
3．\x09节点位移的求解
定义函数Bar2D2Node_Disp(KK,num,p),该函数输入KK为总体刚度矩阵；num为活动自由度编号数组；p为活动自由度方向上的节点力；输出节点位移列阵.具体代码如下：
function u = Bar2D2Node_Disp(KK,num,p)
k=KK(num,num)
u=k\p
4．\x09单元应力的求解
定义函数函数Bar2D2Node_Stress(E,x1,y1,x2,y2,u),该函数计算单元的应力输入弹性模量E,第一个节点坐标（x1,y1）,第二个节点坐标（x2,y2）单元节点位移矢量u,返回单元应力标量stress .具体代码如下：
L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
x=acos((x2-x1)/L);
C=cos(x);
S=sin(x);
stress=E/L*[-C -S C S]*u;
5．\x09计算节点力
定义函数Bar2D2Node_Forces(KK,q),该函数用于计算节点力,KK为刚度矩阵,q为节点位移阵列
function P= Bar2D2Node_Forces(KK,q)
q=zeros(8,1);
q(num)=u;
P=KK*q;
至此,基于Matlab的杆单元有限元分析的程序设计已经完成,遇到实际问题时可以直接调用这些函数就可以解决问.
经典算例
如图所示的结构,各个杆的弹性模量和横截面积都为 , .试基于MATLAB平台求解该结构的节点位移、单元应力以及支反力.
四杆桁架结构
对该问题进行有限元分析的过程如下
（1） 结构的离散化与编号
\x09对该结构进行自然离散,节点编号和单元编号如上图所示
（2）计算各单元的刚度矩阵(基于国际标准单位)
\x09\x09输入弹性模量E、横截面积A,各点坐标.然后分别针对单元1,2,3和4,调用4次Bar2D2Node_Stiffness,就可以得到单元的刚度矩阵.
对应的主程序中代码：
E=2.95e11;A=0.0001;x1=0;y1=0;x2=0.4;y2=0;x3=0.4;y3=0.3;x4=0;y4=0.3;
k1=Bar2D2Node_Stiffness (E,A,x1,y1,x2,y2)
k2=Bar2D2Node_Stiffness (E,A,x2,y2,x3,y3)
k3=Bar2D2Node_Stiffness (E,A,x1,y1,x3,y3)
k4=Bar2D2Node_Stiffness (E,A,x4,y4,x3,y3)
（3） 建立整体刚度方程
由于该结构共有4个节点,因此,设置结构总的刚度矩阵为KK(8×8),先对KK清零,然后四次调用函数Bar2D2Node _Assembly进行刚度矩阵的组装.相关主程序代码为：
KK=zeros(8,8);
KK=Bar2D2Node_Assembly (KK,k1,1,2);
KK=Bar2D2Node_Assembly (KK,k2,2,3);
KK=Bar2D2Node_Assembly (KK,k3,1,3);
KK=Bar2D2Node_Assembly (KK,k4,4,3)
（4）边界条件的处理及刚度方程的求解
由图可以看出,被约束的自由度有：节点1的x,y方向自由度,节点2的y方向自由度,4节点的x、y方向两个自由度.则活动自由度编号为3,5,6.活动自由度对应的节点载荷F3=20000N,F5=0N,F6=25000N,采用高斯消去法进行求解,对应的代码为：
num=[3,5,6];%可活动的自由度编号
p=[20000;0;-25000];
u=Bar2D2Node_Disp(KK,num,p)
（5）支反力的计算
在得到整个结构的节点位移后,由原整体刚度方程就可以计算出对应的支反力.这部分对应的主程序的代码如下：
q=zeros(8,1);
q(num)=u;%节点位移阵列
P=Bar2D2Node_Forces(KK,q)
（6）单元应力的计算
先从整体位移列阵q中提取出单元的位移列阵,然后,调用计算单元应力的函数Bar2D2Node_ElementStress,就可以得到各个单元的应力分量.
u1=[q(1);q(2);q(3);q(4)]
stress1=Bar2D2Node_Stress(E,x1,y1,x2,y2,u1)
u2=[q(3);q(4);q(5);q(6)]
stress2=Bar2D2Node_Stress(E,x2,y2,x3,y3,u2)
u3=[q(1);q(2);q(5);q(6)]
stress3=Bar2D2Node_Stress(E,x1,y1,x3,y3,u3)
u4=[q(7);q(8);q(5);q(6)]
stress4=Bar2D2Node_Stress(E,x4,y4,x3,y3,u4)
（7）计算结果的整理
通过主程序的运行得计算结果.
主程序
%计算各单元的刚度矩阵(以国际标准单位)
E=2.95e11;
A=0.0001;
x1=0;
y1=0;
x2=0.4;
y2=0;
x3=0.4;
y3=0.3;
x4=0;
y4=0.3;
k1=Bar2D2Node_Stiffness (E,A,x1,y1,x2,y2)
k2=Bar2D2Node_Stiffness (E,A,x2,y2,x3,y3)
k3=Bar2D2Node_Stiffness (E,A,x1,y1,x3,y3)
k4=Bar2D2Node_Stiffness (E,A,x4,y4,x3,y3)
%建立整体刚度方程
%由于该结构共有4个节点,因此,结构总的刚度矩阵为KK(8×8),先对K清零,然后四次调用函数Bar2D2Node _Assembly进行刚度矩阵的组装.
KK=zeros(8,8);
KK=Bar2D2Node_Assembly (KK,k1,1,2);
KK=Bar2D2Node_Assembly (KK,k2,2,3);
KK=Bar2D2Node_Assembly (KK,k3,1,3);
KK=Bar2D2Node_Assembly (KK,k4,4,3)
%边界条件的处理及刚度方程求解
num=[3,5,6];%可活动的自由度编号
p=[20000;0;-25000];
u=Bar2D2Node_Disp(KK,num,p)
%支反力的计算
q=zeros(8,1);
q(num)=u;%节点位移阵列
P=Bar2D2Node_Forces(KK,q)
%各单元的应力计算
u1=[q(1);q(2);q(3);q(4)];
stress1=Bar2D2Node_Stress(E,x1,y1,x2,y2,u1)
u2=[q(3);q(4);q(5);q(6)];
stress2=Bar2D2Node_Stress(E,x2,y2,x3,y3,u2)
u3=[q(1);q(2);q(5);q(6)];
stress3=Bar2D2Node_Stress(E,x1,y1,x3,y3,u3)
u4=[q(7);q(8);q(5);q(6)];
stress4=Bar2D2Node_Stress(E,x4,y4,x3,y3,u4)
说的可能有些罗嗦,注意其中有5个function,和最后一个主程序,计算的时候直接运行主程序就可以了.希望能帮助到你.
给你说了也白搭，你大概知道就行了，这个需要代码量。。。简单地说就是，没有多态，那么等号左边是啥右边就得是啥，这就叫耦合，有了多态，左边是父类（或者接口），右边是子类（或实现类），我只管调用接口里面的方法就是了，至于你实现类怎么去实现，那是你的事，你要修改一下实现，你只管去把实现类换掉，我这边一行代码都不用变，这就解耦了
特征向量都求出来了，用哪一种归一化还不跟玩儿一样吗？okok已经回答你了，这里再贴一下，主要因为已经消失几天了，出来透口气：div class="blockcode"复制内容到剪贴板h5代码:/h5code id="code0"function ziyouzhendong1(k,m)br /% k=600*[1 -1 0;-1 3 -2;0 -2 5];m=diag([1 1.5 2]);br /clcbr /[C,B]=eig(inv(m)*k);br /n=size(m,2);br /for i=1:nbr / Ct=C(:,i);br / zhenxing(:,i)=Ct/Ct(1);br /endbr /w=sqrt(diag(B));br /B1=diag(B)';br /Eigen=inv(m)*k;br /for i=1:nbr / t1(:,i)=B1(i)*zhenxing(:,i);br / t2(:,i)=Eigen*zhenxing(:,i);br /endbr /B,zhenxing,t1,t2,w/code/div其实完全可以不用循环，这是上学初学MATLAB一个月左右时写的程序，目的是帮助同班美女从繁重的结构动力学作业中解脱出来并有充分的时间买化妆品and keep us amused...，当时还没人学MATLAB这玩意儿，全拿计算器算动力学矩阵的问题，显然即使最多只有4×4的K矩阵也足够让人受折磨了...br /
br /
转载
在BDF文件中SOL之下加入以下语句，可以输出刚度矩阵
compile semg
alter ’kjjz.*stiffness’ $
matprn kjjz// $