机械必威体育网址

标题: 分析一段matlab有限元程序 [打印本页]

作者: yinzengguang 时间: 2013-3-24 13:41
标题: 分析一段matlab有限元程序
这是关于梁单元的，可能大家觉得很简单。。。
今天翻电脑里的东西时发现的，是我大三时有限元时的作业，记得当时花了很多时间研究，学了不少东西，一个简单的作业可以涉及各方面的知识。
毕业半年了，虽然记不清弹性矩阵、刚度矩阵等，但是只要一看书（好像现在的心境不下来）应该会，学校交的学习方法和理解问题的能力，有些东西不必记，除非每天从事，那该叫熟能生巧。
记得当时编了两个程序，梁和三维实体的，我分享一下梁的程序吧，起码有亮点和创意，可以自己选择划分个数，在matlab方面花了不少功夫。
非常感谢当时的有限元授课老师“韩清凯”，虽然老师已经到了大工。上课的教材为韩清凯的《弹性力学及有限元法基础教程》，书上有个梁单元的例子，在82页。
--------------------------------------------------------------------------------
梁程序如下，名字就不写了，用昵称吧“太平岛”，这是常用的网络昵称

%------------------------ BEAM2  ----------------
disp('========================================');
disp('          PROGRAM BEAM2             ');
disp('       Beam Bending Analysis          ');
disp('       The programmer:太平岛          ');
disp('========================================');
disp(' ');                      %输出空行
warning off all                      %关闭所有警告提示(求积分时，警告信息比较多)
Ne=input('Please input the number of elements:');       %Ne为划分的单元数
NJ=Ne+1;                      %NJ为节点数
x1=0;
x2=sym('L');
x=sym('x');
j=0:3;
v=x.^j;
A=[1,x1,x1^2,x1^3;0,1,2*x1,3*x1^2;1,x2,x2^2,x2^3;0,1,2*x2,3*x2^2];
N=v*inv(A);                      %形函数
%-----------------------------------------------
% 求单元刚度矩阵
E=4.0e11;
I=5.2e-7;                      %I=bh^3/12=5.2e-7；
EI=4.0e11*5.2e-7;
B=diff(N,2);
k=B'*B;
Kee=EI*int(k,x,0,'L');
Ke=sym(zeros(4,4,Ne));       %用三维矩阵存放单元刚度矩阵,每页存一个，并初始化0
Ke(:,:,1)=subs(Kee,'L',(10/Ne));
T=eye(4);             % 因为梁于x轴的夹角为0，所以坐标变换矩阵为T=eye(4)
Ke(:,:,1)=T*Ke(:,:,1)*T';
for ii=2:Ne
Ke(:,:,ii)=Ke(:,:,1);
end
Ke=double(Ke);
%------------------------------------------------
% 由矩阵装换法求整体刚度矩阵
% 自由度Ndof=2*NJ
Ndof=2*NJ;
K=zeros(Ndof);
for ii=1:Ne
G=[zeros(4,2*(ii-1)),eye(4,Ndof-2*(ii-1))];
KK=G'*Ke(:,:,Ne)*G;
K=K+KK;
end
%------------------------------------------------
% 约束条件,对整体刚度矩阵进行修正，以便计算
F=zeros(Ndof,1);
F(Ndof-1)=-100;
% u=F*inv(K),u=[v1,xta1,v2,xta2,v3,xta3,v4,xta4]',v1=0,xta1=0
K(1,

=0; K(:,1)=0; %可以省略
K(2,

=0;       K(:,2)=0;       %可以省略
KX=K(3:Ndof,3:Ndof);
FX=F(3:Ndof,1);
%------------------------------------------------
%求整体节点位移列阵
uu=inv(KX)*FX;
u=[0;0;uu];
ii=1:2:2*NJ;
v=u(ii);                      %各节点挠度
xta=u(ii+1);                      %各节点转角
%------------------------------------------------
% 后处理，计算单元应变应力
B=subs(B,x,'L');
B=subs(B,'L',10/Ne);
Strain=zeros(1,Ne);             %单元应变，并初始化
Stress=zeros(1,Ne);             %单元应力，并初始化
for ii=1:Ne
Strain(1,ii)=B*u(2*ii-1:2*ii+2,1);
Stress(1,ii)=E*Strain(1,ii);
end
%--------------------------------------------------
% 以下程序为屏幕输出处理
disp(' ');
disp(' Input:1-print Node displacement ');
disp(' Input:2-print strain ');
disp(' Input:3-print stress ');
disp(' Input:q-quit or use ''ctrl+c'' Shortcut to end the program');

while 1
ii=input('Please input1 or 2 or 3:','s');
      switch(ii)
         case '1'
            disp('Node displacement');
            disp(u');
         case '2'
            disp('element strain');
            disp(Strain);
         case '3'
            disp('element stress');
            disp(Stress);
         case 'q'
            disp('End of program');
            break;
         otherwise
            disp('error!Please input again');
            continue;
      end
end

作者: yinzengguang 时间: 2013-3-24 13:44
那个地方怎么变成笑脸了，冒号+右括号=笑脸，改了下，应该是下面，把英文括号改为中文，就好了吧
K(1,:）=0; K(:,1)=0; %可以省略
K(2,:）=0; K(:,2)=0; %可以省略

作者: 孤若流星 时间: 2013-3-24 14:55
没看懂

作者: jiaweicz 时间: 2013-3-24 16:53
谢谢你啊，太感谢你了

作者: yinzengguang 时间: 2013-3-24 18:35

jiaweicz 发表于 2013-3-24 16:53 / j3 m( v2 P) P2 S
谢谢你啊，太感谢你了

这谢啥呢？

作者: 怀念昔日熊猫 时间: 2013-3-24 21:21
我以前编过一个5体展开的多体姿态动力学计算程序。。。可惜早就忘记了，sigh

作者: d調旳華孋 时间: 2013-11-7 20:39
这程序也运行不了啊

作者: d調旳華孋 时间: 2013-11-7 20:44
j=0:3;
v=x.^j; 是干啥的

欢迎光临机械必威体育网址 (//www.szfco.com/)