搜索
cad2007下载
软件下载
solidworks下载
proe下载
机械标准
机械招聘
机械图纸
机械百科
机械交易网
网站建设
机械设计手册
proe视频教程
cad2013视频教程
solidworks2008视频教程
CAD2004视频教程

分析一段matlab有限元程序

[复制链接]
查看: 467|回复: 7

该用户从未签到

发表于 2013-11-15 18:34:59 | 显示全部楼层 |阅读模式
这是关于梁单元的,可能大家觉得很简单。。。
今天翻电脑里的东西时发现的,是我大三时有限元时的作业,记得当时花了很多时间研究,学了不少东西,一个简单的作业可以涉及各方面的知识。
毕业半年了,虽然记不清弹性矩阵、刚度矩阵等,但是只要一看书(好像现在的心境不下来)应该会,学校交的学习方法和理解问题的能力,有些东西不必记,除非每天从事,那该叫熟能生巧。
记得当时编了两个程序,梁和三维实体的,我分享一下梁的程序吧,起码有亮点和创意,可以自己选择划分个数,在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
回复

使用道具 举报

该用户从未签到

 楼主| 发表于 2013-11-15 18:35:46 | 显示全部楼层
那个地方怎么变成笑脸了,冒号+右括号=笑脸,改了下,应该是下面,把英文括号改为中文,就好了吧
K(1,:)=0;  K(:,1)=0;  %可以省略
K(2,:)=0;  K(:,2)=0;  %可以省略
回复 支持 反对

使用道具 举报

该用户从未签到

发表于 2013-11-15 18:36:19 | 显示全部楼层
没看懂
回复 支持 反对

使用道具 举报

该用户从未签到

发表于 2013-11-15 18:36:53 | 显示全部楼层
谢谢你啊,太感谢你了
回复 支持 反对

使用道具 举报

该用户从未签到

 楼主| 发表于 2013-11-15 18:37:38 | 显示全部楼层

jiaweicz 发表于 2013-3-24 16:53

谢谢你啊,太感谢你了

这谢啥呢?
回复 支持 反对

使用道具 举报

该用户从未签到

发表于 2013-11-15 18:38:30 | 显示全部楼层
我以前编过一个5体展开的多体姿态动力学计算程序。。。可惜早就忘记了,sigh
回复 支持 反对

使用道具 举报

该用户从未签到

发表于 2013-11-15 18:39:25 | 显示全部楼层
这程序也运行不了啊
回复 支持 反对

使用道具 举报

该用户从未签到

发表于 2013-11-15 18:40:19 | 显示全部楼层
j=0:3;
v=x.^j;是干啥的
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 注册会员

本版积分规则

Copyright © 2012-2035 厦门鑫时器科技有限公司 版权所有
闽ICP备2023009579号-1 技术支持:机械网站建设  Powered by Discuz! X3.4
快速回复 返回顶部 返回列表