机械必威体育网址

 找回密码
 注册会员

QQ登录

只需一步,快速开始

搜索
查看: 3725|回复: 7
打印 上一主题 下一主题

分析一段matlab有限元程序

[复制链接]
跳转到指定楼层
1#
发表于 2013-3-24 13:41:49 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
这是关于梁单元的,可能大家觉得很简单。。。. X& Z2 a6 C) D- F( w( u' G
今天翻电脑里的东西时发现的,是我大三时有限元时的作业,记得当时花了很多时间研究,学了不少东西,一个简单的作业可以涉及各方面的知识。
9 l0 f  o+ d( h5 f4 D毕业半年了,虽然记不清弹性矩阵、刚度矩阵等,但是只要一看书(好像现在的心境不下来)应该会,学校交的学习方法和理解问题的能力,有些东西不必记,除非每天从事,那该叫熟能生巧。
! n/ m4 Y& }- c记得当时编了两个程序,梁和三维实体的,我分享一下梁的程序吧,起码有亮点和创意,可以自己选择划分个数,在matlab方面花了不少功夫。9 C; T7 u& ^. b' Z
非常感谢当时的有限元授课老师“韩清凯”,虽然老师已经到了大工。上课的教材为韩清凯 的《弹性力学及有限元法基础教程》,书上有个梁单元的例子,在82页。
, D  d+ G& ^+ g2 H* R8 k- M--------------------------------------------------------------------------------9 v$ \! _+ o% @' A9 z; w
梁程序如下,名字就不写了,用昵称吧“太平岛”,这是常用的网络昵称3 G3 x% d, p$ p
7 a  z& N+ K9 {5 {3 J3 V/ S1 w
" X# }4 @7 K' F
%------------------------ BEAM2  ----------------
% p; N$ n6 F7 _7 sdisp('========================================');( N  N5 k( _6 h
disp('            PROGRAM BEAM2               ');
8 N) r4 F8 _! `( l% b. I( O, ddisp('        Beam Bending Analysis           ');
% m1 e; C' i% d/ rdisp('        The programmer:太平岛           ');0 y, h$ T) ?( L9 x
disp('========================================');
9 [6 r* F/ {( D9 Z( A# S) qdisp(' ');                        %输出空行3 r2 T2 ^" A3 x" i1 W  A8 {
warning off all                        %关闭所有警告提示(求积分时,警告信息比较多)
- m" C1 T1 p- t9 l4 aNe=input('Please input the number of elements:');        %Ne为划分的单元数
; Y% Z! q' {* @  n1 {NJ=Ne+1;                        %NJ为节点数
8 |# k4 p6 e; }8 i, a+ r$ nx1=0;
: d  X( L9 e9 B! F3 r2 jx2=sym('L');
- @' D, r% Z# @x=sym('x');                        
! R7 m% I! o" Q% Kj=0:3;/ X5 x9 N1 x) S  _& y) P" Q
v=x.^j;# v- U9 x# Q* e2 E  [
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];5 J% e' q" q. m1 l5 L
N=v*inv(A);                        %形函数6 v! H( c3 Q  Z; V0 x% c
%-----------------------------------------------3 b4 _; Y3 E/ C% W2 y" y; t( @
% 求单元刚度矩阵9 r! r6 \) @$ k8 I% K- J2 \; R& i. F
E=4.0e11;
' ]7 r0 a; R" G" TI=5.2e-7;                        %I=bh^3/12=5.2e-7;
; H6 t7 Z* L6 p5 m5 [% B+ R8 b" KEI=4.0e11*5.2e-7;8 B! r: l' H' v, c- P2 A, ?& L
B=diff(N,2);
/ w8 d% J! x* A: t6 Ak=B'*B;7 G- J& i. B4 s1 I
Kee=EI*int(k,x,0,'L');3 @& h% {. e! ~, V' `
Ke=sym(zeros(4,4,Ne));        %用三维矩阵存放单元刚度矩阵,每页存一个,并初始化0
2 v: r  @5 o1 q+ UKe(:,:,1)=subs(Kee,'L',(10/Ne));
" ?% |. c' m  V: `; e- c6 S, Z- Y$ YT=eye(4);                % 因为梁于x轴的夹角为0,所以坐标变换矩阵为T=eye(4)8 t: O3 B& L8 l& ~2 g
Ke(:,:,1)=T*Ke(:,:,1)*T';7 f# c- |  q# Y
for ii=2:Ne
% c& U$ t) }! ]" B7 ^! d: c    Ke(:,:,ii)=Ke(:,:,1);
; `5 i3 _5 v, `+ `end
7 n# G! o+ \0 T; H& `2 h8 `Ke=double(Ke);% T* f8 Y) D/ N) t- y4 }
%------------------------------------------------' R; D% z. D: R0 |- r9 Z& Z
% 由矩阵装换法求整体刚度矩阵
' ~+ W$ u- q! n2 ?' Y! \% 自由度Ndof=2*NJ
$ W7 B. i: G$ Z1 L/ uNdof=2*NJ;& d$ C, ^! R: N) V( z6 s
K=zeros(Ndof);& `  T! I" B4 T/ K9 |; I8 i
for ii=1:Ne
1 ~. {9 o6 u% _% B$ i# ~. R0 a7 q    G=[zeros(4,2*(ii-1)),eye(4,Ndof-2*(ii-1))];
) g- i6 X) ~$ }1 z7 w2 @& L( v    KK=G'*Ke(:,:,Ne)*G;4 v& ~6 l% \, e6 M4 ^
    K=K+KK;
0 U6 B; @" ~$ u6 S" m( ~) H  j6 kend  
: V+ j- U% m) z0 p%------------------------------------------------
( L/ l5 N) ?3 S& Z6 H8 |, T6 @% 约束条件,对整体刚度矩阵进行修正,以便计算' _0 F: _* E0 u! t# r, V& T# t8 _
F=zeros(Ndof,1);
+ y5 `/ H- }% O+ B7 ?0 aF(Ndof-1)=-100;
9 s% }+ v) V/ [: U5 v# E% u=F*inv(K),u=[v1,xta1,v2,xta2,v3,xta3,v4,xta4]',v1=0,xta1=0& N4 Y1 h9 B+ c
K(1,=0;        K(:,1)=0;        %可以省略) e* ^3 Q3 B& T3 N9 I" p- w
K(2,=0;        K(:,2)=0;        %可以省略
$ A8 A% l; D  M8 x& bKX=K(3:Ndof,3:Ndof);
8 q! y3 x2 y+ s' z2 KFX=F(3:Ndof,1);
  p2 [, U5 Q& @%------------------------------------------------
9 Y: L2 u' X: A1 [7 h; m# ^%求整体节点位移列阵  X$ w5 ~) k8 C
uu=inv(KX)*FX;
" K5 j/ E0 l6 qu=[0;0;uu];
. Z. S5 d) ~6 b0 j* B: Uii=1:2:2*NJ;/ ~2 v! T2 c' [  E
v=u(ii);                        %各节点挠度
$ w! R9 }8 p6 w) N2 rxta=u(ii+1);                        %各节点转角
! q. C; d3 h  k/ c" A9 {%------------------------------------------------6 x+ n7 D& k- K2 b! h, J1 S
% 后处理,计算单元应变应力
0 o: T( d9 z8 P+ ?B=subs(B,x,'L');
7 }* N7 h% p' `% X, ]B=subs(B,'L',10/Ne);9 ]* U. N" I, o8 S$ i0 D
Strain=zeros(1,Ne);                %单元应变,并初始化+ n/ W" h7 ?, g: F0 c( a5 y: K7 d
Stress=zeros(1,Ne);                %单元应力,并初始化
9 u5 A# U9 G# R. Q; Rfor ii=1:Ne
2 f2 f5 H, T. C    Strain(1,ii)=B*u(2*ii-1:2*ii+2,1);
1 p6 {3 a  Z, A1 v8 r6 h    Stress(1,ii)=E*Strain(1,ii);
3 L6 N) C: j* v( N! ^- _, I! N% aend
  |7 j/ u# c' L- L%--------------------------------------------------
  ~$ T- r% K2 l, ?' Q% 以下程序为屏幕输出处理, [; E4 J9 h$ d6 D) S5 O
disp(' ');4 c7 w1 e$ B( U3 t7 b" I! S* w1 _
disp(' Input:1-print Node displacement ');
0 P: O5 r" k1 n* [9 Q6 bdisp(' Input:2-print strain ');
) B) g% D% F4 ?9 T# s( vdisp(' Input:3-print stress ');; F# j+ h& M9 S2 U
disp(' Input:q-quit or use ''ctrl+c'' Shortcut to end the program');- Z1 A% p) S! `/ A* i3 ^

" P: M1 _, \( \while 1/ ^; H6 O+ n. ]% p) B, p) f! d! ?
    ii=input('Please input1 or 2 or 3:','s');
) ^$ }/ |0 k' ~: B        switch(ii)
" J( B1 n: o) Q" X            case '1'% B6 m  b8 Y1 m/ L, a1 z5 n- a7 G
                disp('Node displacement');) t/ T+ u( x1 V: I/ v' H  t8 @
                disp(u');/ ^: o: f1 t6 x
            case '2'$ E3 B; s5 ^; O8 @7 B' r: _: C% c
                disp('element strain');
$ ]- g" F; l- \. S                disp(Strain);4 d' V  U. D- n8 K8 d- L
            case '3'
) U' M1 j0 B5 x6 o$ c                disp('element stress');$ v) V2 w' l! @. ?
                disp(Stress);
6 _6 w% t' m+ S! I$ H5 c+ A            case 'q'# a  R2 ]1 m& G! ^, X/ h
                disp('End of program');" P0 n& V6 F9 Q3 ]* w2 z
                break;
' ~/ k/ b) n1 }$ r! }0 R            otherwise6 P9 X/ {' b* ^" O6 O
                disp('error!Please input again');
& O8 l. ~: Y, s                continue;- x* w/ U" K' U5 a9 {  Z
        end
- ^# J3 F+ Y. M+ m" Tend2 B7 ]8 d+ |$ r. z; C! ]) u4 |' M. G% B
9 `1 I0 Q7 U; H0 ?7 _
! B+ F! J) h4 H0 D* N) A8 N
回复

使用道具 举报

2#
 楼主| 发表于 2013-3-24 13:44:32 | 只看该作者
那个地方怎么变成笑脸了,冒号+右括号=笑脸,改了下,应该是下面,把英文括号改为中文,就好了吧
$ g# d- |  @' F! p4 v' EK(1,:)=0;        K(:,1)=0;        %可以省略( ]& e1 i9 |' r3 Z: N; _6 B1 }
K(2,:)=0;        K(:,2)=0;        %可以省略
回复 支持 反对

使用道具 举报

3#
发表于 2013-3-24 14:55:52 | 只看该作者
没看懂
回复 支持 反对

使用道具 举报

4#
发表于 2013-3-24 16:53:09 | 只看该作者
谢谢你啊,太感谢你了
回复 支持 反对

使用道具 举报

5#
 楼主| 发表于 2013-3-24 18:35:59 | 只看该作者
jiaweicz 发表于 2013-3-24 16:53
  D; E# |% m1 N4 i+ |+ S3 ~* n谢谢你啊,太感谢你了
/ J$ c: g8 Z5 }; L$ y( w1 _: M
这谢啥呢?
回复 支持 反对

使用道具 举报

6#
发表于 2013-3-24 21:21:19 | 只看该作者
我以前编过一个5体展开的多体姿态动力学计算程序。。。可惜早就忘记了,sigh

点评

是啊,很多东西不用就忘记啦  发表于 2013-3-26 12:34
回复 支持 反对

使用道具 举报

7#
发表于 2013-11-7 20:39:02 | 只看该作者
这程序也运行不了啊
回复 支持 反对

使用道具 举报

8#
发表于 2013-11-7 20:44:34 | 只看该作者
j=0:3;
; w! X; @- G6 I3 cv=x.^j;  是干啥的
回复 支持 反对

使用道具 举报

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

本版积分规则

小黑屋|手机版|Archiver|机械必威体育网址 ( 京ICP备10217105号-1,京ICP证050210号,浙公网安备33038202004372号 )

GMT+8, 2024-11-26 02:31 , Processed in 0.056699 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4 Licensed

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表