机械必威体育网址

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

作者: yinzengguang    时间: 2013-3-24 13:41
标题: 分析一段matlab有限元程序
这是关于梁单元的,可能大家觉得很简单。。。
" J/ \; x4 b4 T6 r' L今天翻电脑里的东西时发现的,是我大三时有限元时的作业,记得当时花了很多时间研究,学了不少东西,一个简单的作业可以涉及各方面的知识。
, v$ p  D- B9 m* x0 j  C9 W: X9 R. A毕业半年了,虽然记不清弹性矩阵、刚度矩阵等,但是只要一看书(好像现在的心境不下来)应该会,学校交的学习方法和理解问题的能力,有些东西不必记,除非每天从事,那该叫熟能生巧。" l  G8 S7 O% K1 e, ?
记得当时编了两个程序,梁和三维实体的,我分享一下梁的程序吧,起码有亮点和创意,可以自己选择划分个数,在matlab方面花了不少功夫。
9 S  m' x; h% u非常感谢当时的有限元授课老师“韩清凯”,虽然老师已经到了大工。上课的教材为韩清凯 的《弹性力学及有限元法基础教程》,书上有个梁单元的例子,在82页。) r* y: h- [# a7 T
--------------------------------------------------------------------------------" n1 L  d9 m9 a# E# _! k( y4 a2 [# O, V
梁程序如下,名字就不写了,用昵称吧“太平岛”,这是常用的网络昵称
, i2 I. W9 ^& a, j
; K9 y0 \) t4 M: p' T6 E/ R. }( ~' s2 n2 b/ v
%------------------------ BEAM2  ----------------
6 m: G2 T' S2 s% q: M. T6 b' ?& kdisp('========================================');
; t! F% f, Q9 \" ^* [  hdisp('            PROGRAM BEAM2               ');) f- r1 @' i3 W
disp('        Beam Bending Analysis           ');
3 o4 f' t( e8 }& C6 J# r& Ydisp('        The programmer:太平岛           ');" }3 ]& @, C4 u  y* q8 ]) O1 c
disp('========================================');4 t4 c, [* T& H! {$ f4 U
disp(' ');                        %输出空行/ _- c9 E! L: A+ z( h: x! o) r( I5 E
warning off all                        %关闭所有警告提示(求积分时,警告信息比较多)
) m9 M/ F% d, R" nNe=input('Please input the number of elements:');        %Ne为划分的单元数
3 w. Q& \; K; W! u. SNJ=Ne+1;                        %NJ为节点数
4 v* p1 x0 d; ax1=0;. k& ?. j+ j+ |, u
x2=sym('L');# M5 Q- a9 T2 B% ?& `
x=sym('x');                        6 A4 ]8 s9 t' y* b7 K0 ~. ~" M# `
j=0:3;
7 a4 |0 {6 E  T% U3 ?" Sv=x.^j;  w- A4 L  I1 R! \$ o$ ?1 G5 h
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];/ W$ ^9 H& \) }9 Q5 T4 _
N=v*inv(A);                        %形函数$ }# N! b5 G' C' @
%-----------------------------------------------
- M; z/ L6 z$ j( ^8 \; I6 a  |8 `% 求单元刚度矩阵
0 k% i9 a$ ]' E1 \9 [/ u. Y+ \3 wE=4.0e11;
/ v, t/ r6 m- o, }4 q8 m- vI=5.2e-7;                        %I=bh^3/12=5.2e-7;: ^1 ~( v! z) z
EI=4.0e11*5.2e-7;4 r& ~2 {* ?# f
B=diff(N,2);2 Y" _7 o- U/ W4 d  ?0 D" A
k=B'*B;
7 ?/ }$ I: {5 e' I3 S  u8 uKee=EI*int(k,x,0,'L');
6 L& B. v* u8 M. M3 I+ Z4 CKe=sym(zeros(4,4,Ne));        %用三维矩阵存放单元刚度矩阵,每页存一个,并初始化0* I8 b# b" {6 d5 V  G
Ke(:,:,1)=subs(Kee,'L',(10/Ne));% t7 \2 _2 y/ N
T=eye(4);                % 因为梁于x轴的夹角为0,所以坐标变换矩阵为T=eye(4); J% g9 W" V) d( v% ]8 ~" }
Ke(:,:,1)=T*Ke(:,:,1)*T';
4 l! {( o0 l+ W6 o( Tfor ii=2:Ne! f( T5 \5 n$ W
    Ke(:,:,ii)=Ke(:,:,1);
- V! m0 S* ~0 K: \) U8 Lend
8 b+ c" _( g3 S; b: ?( UKe=double(Ke);
) X4 S3 ^1 o4 l  m! x' C! F% p- F%------------------------------------------------4 p) }4 x* H$ P+ d6 b  H$ _: G3 e
% 由矩阵装换法求整体刚度矩阵
* t& R; C4 H2 ]5 I  v% 自由度Ndof=2*NJ
# a* r  t# r) W% E8 K) mNdof=2*NJ;
6 u9 l6 m6 g$ _. mK=zeros(Ndof);
3 G, H5 w3 ?, ufor ii=1:Ne
6 k( _: K- S. a# {    G=[zeros(4,2*(ii-1)),eye(4,Ndof-2*(ii-1))];/ J) s% S- K. m4 L3 _! _8 I* z
    KK=G'*Ke(:,:,Ne)*G;: o+ ]/ i/ c( o" N2 Q& S: _& q
    K=K+KK;# g' f$ o. A! S) |
end  
5 Y4 l8 D8 t% s6 v%------------------------------------------------
" o/ G: z! N( Q: r% 约束条件,对整体刚度矩阵进行修正,以便计算- j4 `7 Q' |3 m* N
F=zeros(Ndof,1);
8 l1 t) q+ w( OF(Ndof-1)=-100;
$ u% x4 `/ }9 j/ J% u=F*inv(K),u=[v1,xta1,v2,xta2,v3,xta3,v4,xta4]',v1=0,xta1=0
( z/ W7 @& e7 xK(1,=0;        K(:,1)=0;        %可以省略
+ T  O7 k$ a" c& H1 s, K8 fK(2,=0;        K(:,2)=0;        %可以省略& N1 r7 r9 C* ^+ p+ u
KX=K(3:Ndof,3:Ndof);( X1 P. j! V& v- F; [
FX=F(3:Ndof,1);7 M" G2 O" D& E- V5 G* c+ J# m
%------------------------------------------------
: C6 M: `4 E$ R+ ^4 P+ T. M0 _. j%求整体节点位移列阵
8 `( f7 V: t- M% `( Euu=inv(KX)*FX;5 k2 p! \9 a1 D; \* n
u=[0;0;uu];$ u. A" O5 t, a! U3 @% p; y) n: r
ii=1:2:2*NJ;+ a& B+ u3 A1 I: \
v=u(ii);                        %各节点挠度
8 A: a8 e3 W4 o( X4 p; L' [) nxta=u(ii+1);                        %各节点转角. g5 Y8 Q+ }3 i$ b; B3 f( t
%------------------------------------------------/ [" k, `# I7 n6 m
% 后处理,计算单元应变应力
' {* J3 C& n& B: p& R4 q! OB=subs(B,x,'L');" N$ ?; n0 @! y1 Q
B=subs(B,'L',10/Ne);
7 M; o3 ?1 s5 Y# z5 _* KStrain=zeros(1,Ne);                %单元应变,并初始化
0 ^- y0 n4 J5 V6 CStress=zeros(1,Ne);                %单元应力,并初始化' y) P6 Q. E; A7 p# C7 L
for ii=1:Ne
) _8 M$ n1 B- g, ^3 G    Strain(1,ii)=B*u(2*ii-1:2*ii+2,1);) k0 f$ ?! K- e& u# t0 i  V
    Stress(1,ii)=E*Strain(1,ii);
: x0 {7 W- ]* z8 Y' d8 z, j: oend
3 t% \3 |: X, }) X9 J* u5 U, m  `%--------------------------------------------------, c5 i0 H8 }* G& k& p8 D4 ?$ Y0 W- i
% 以下程序为屏幕输出处理* e/ n  e( _; d% t
disp(' ');
* G$ e( Y6 r3 V) I/ n* gdisp(' Input:1-print Node displacement ');2 K' a& G+ y1 y& `0 E
disp(' Input:2-print strain ');
# S; Q1 [( y/ A0 Pdisp(' Input:3-print stress ');
3 p) ]5 H% I% ~0 }disp(' Input:q-quit or use ''ctrl+c'' Shortcut to end the program');
5 M7 s2 f! N# l- U) o( `6 q; R6 o- G4 j5 ~$ V% b( v6 z
while 1- G4 K5 b+ f, g( [
    ii=input('Please input1 or 2 or 3:','s');9 A6 m# j8 b; x3 k5 R+ w
        switch(ii)
& H' h: n. Q" U1 z: \            case '1'5 W' C: \: C' F  P5 M- C+ h- ~
                disp('Node displacement');+ {' O) i. W2 u7 @' a
                disp(u');
0 Y0 U. q  Q" w% ?- X- y            case '2'% `7 z+ v6 ?  _; F
                disp('element strain');7 h7 m$ K# ^: S" \0 g% h0 N3 k2 T
                disp(Strain);2 h" \* c$ t1 j. ~' J5 r- n
            case '3'1 c! {' {! h1 y6 m$ I. b# p
                disp('element stress');) c1 Z, }& \: s7 f( r
                disp(Stress);
! n& F) N3 t& _$ M/ }8 f* n            case 'q'5 ~7 B" F) c2 Y
                disp('End of program');- I# }! h" h  v' d, I8 F5 J0 F5 |
                break;
* `) x9 q/ {& |' V2 [" N            otherwise) a0 i: |! F1 I$ w  D
                disp('error!Please input again');
3 Q/ e( s  B3 J, _# K( A                continue;  p3 s! a8 |! p5 U$ V1 U
        end# B3 J9 M2 |/ z& b* S
end
$ K! m: V4 R- o* L! Q/ b

% d8 n  u' S' E( b" X1 E3 P8 p  \$ T

作者: yinzengguang    时间: 2013-3-24 13:44
那个地方怎么变成笑脸了,冒号+右括号=笑脸,改了下,应该是下面,把英文括号改为中文,就好了吧
1 K' C& s) }1 G& j0 lK(1,:)=0;        K(:,1)=0;        %可以省略
- t5 G& f7 k0 }. k& \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
谢谢你啊,太感谢你了

  i( w$ {$ r, B) u& j这谢啥呢?
作者: 怀念昔日熊猫    时间: 2013-3-24 21:21
我以前编过一个5体展开的多体姿态动力学计算程序。。。可惜早就忘记了,sigh
作者: d調旳華孋    时间: 2013-11-7 20:39
这程序也运行不了啊
作者: d調旳華孋    时间: 2013-11-7 20:44
j=0:3;
* A3 u: a! A5 O+ Jv=x.^j;  是干啥的




欢迎光临 机械必威体育网址 (//www.szfco.com/) Powered by Discuz! X3.4