机械必威体育网址

 找回密码
 注册会员

QQ登录

只需一步,快速开始

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

分析一段matlab有限元程序

[复制链接]
跳转到指定楼层
1#
发表于 2013-3-24 13:41:49 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
这是关于梁单元的,可能大家觉得很简单。。。+ B  g8 [- h# L
今天翻电脑里的东西时发现的,是我大三时有限元时的作业,记得当时花了很多时间研究,学了不少东西,一个简单的作业可以涉及各方面的知识。6 X0 D/ R) s( ]; s1 X$ ~9 Q
毕业半年了,虽然记不清弹性矩阵、刚度矩阵等,但是只要一看书(好像现在的心境不下来)应该会,学校交的学习方法和理解问题的能力,有些东西不必记,除非每天从事,那该叫熟能生巧。
! [' \, D* f7 F% Y记得当时编了两个程序,梁和三维实体的,我分享一下梁的程序吧,起码有亮点和创意,可以自己选择划分个数,在matlab方面花了不少功夫。7 C/ @: W, \: [: e6 Y& A" P* _& D
非常感谢当时的有限元授课老师“韩清凯”,虽然老师已经到了大工。上课的教材为韩清凯 的《弹性力学及有限元法基础教程》,书上有个梁单元的例子,在82页。
1 Q/ Q, f& s: U5 ^& E- A% v--------------------------------------------------------------------------------
" @* w; w2 A3 d% V; S梁程序如下,名字就不写了,用昵称吧“太平岛”,这是常用的网络昵称( ^2 G: q* [/ z
$ B( E* ]5 [7 g8 I7 Q
8 ~$ w$ ~# @  K$ j* K
%------------------------ BEAM2  ----------------( }  A) Q4 |0 @+ M& {& a+ V
disp('========================================');5 t( B7 c0 A( N) e
disp('            PROGRAM BEAM2               ');
( j' x+ _" `2 K9 Ddisp('        Beam Bending Analysis           ');
: T( W% D& i. l( C- tdisp('        The programmer:太平岛           ');
1 m* N& u4 t$ k0 f. |disp('========================================');# q0 ^2 R: g: [- m
disp(' ');                        %输出空行
! I2 x1 e+ I' p% p& `warning off all                        %关闭所有警告提示(求积分时,警告信息比较多)5 `! w; {6 i5 w9 C+ q
Ne=input('Please input the number of elements:');        %Ne为划分的单元数, `( {) Y( ]: n* v5 b
NJ=Ne+1;                        %NJ为节点数: s# j  B$ x2 y/ i1 ]
x1=0;& t  o, l0 C: h3 p1 H! X* [8 A
x2=sym('L');" d* Y0 R: H) P6 r) W8 ^& Y
x=sym('x');                        
3 |! a0 K- X* y9 {7 g9 ?j=0:3;
% f' s# a: G/ P3 u3 X3 P/ h- Jv=x.^j;3 ]# b! o, s& J5 L+ q* B
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];) a$ q$ t; P& _
N=v*inv(A);                        %形函数
& q& y2 H* i! v/ ]6 |& ^1 v%-----------------------------------------------
' r, l# H9 O" Q+ B: t0 ^% 求单元刚度矩阵) s( j. C+ N& r
E=4.0e11;$ D% |" R  M' n* \1 m9 y
I=5.2e-7;                        %I=bh^3/12=5.2e-7;. t. V  h9 x$ F+ G9 O
EI=4.0e11*5.2e-7;
$ F2 @" l. |& u; B: _- t: sB=diff(N,2);
: E8 v( J( @5 Bk=B'*B;. I. t) R5 e+ b- V) [
Kee=EI*int(k,x,0,'L');& G* Y  d8 Z8 _9 j) c% n5 W
Ke=sym(zeros(4,4,Ne));        %用三维矩阵存放单元刚度矩阵,每页存一个,并初始化0
1 o8 D' n. |& N/ z: t' B+ Y7 KKe(:,:,1)=subs(Kee,'L',(10/Ne));8 q" }' Y' X- u- {) o
T=eye(4);                % 因为梁于x轴的夹角为0,所以坐标变换矩阵为T=eye(4)
6 k+ T$ Z1 N2 j4 H5 e* p% A( b. wKe(:,:,1)=T*Ke(:,:,1)*T';5 k2 v% j* i# t6 r& O+ f9 W
for ii=2:Ne6 ^# t6 Q: M1 Y! J. r3 L- U- x
    Ke(:,:,ii)=Ke(:,:,1);4 ]: q3 Q/ O' t9 A- o- `
end " M- e4 v$ y- W8 ]
Ke=double(Ke);2 k& x0 h# ?8 K9 ?9 j  a
%------------------------------------------------) f5 o. Y+ S5 b! u" B) |! L+ v
% 由矩阵装换法求整体刚度矩阵- r6 e9 g8 D8 x- G0 h% s, O& E
% 自由度Ndof=2*NJ' T4 u8 E- `" b- M& q4 T: N0 k) I
Ndof=2*NJ;
/ t- O) ]  ^5 x  X+ S/ X1 j) LK=zeros(Ndof);) Q6 J1 n6 N# W
for ii=1:Ne
: B6 l$ Q8 d0 t: n  T' P+ Y    G=[zeros(4,2*(ii-1)),eye(4,Ndof-2*(ii-1))];' T: a) n4 ?2 ]0 z) Z# y4 _$ z
    KK=G'*Ke(:,:,Ne)*G;
4 K6 A; ~/ I7 H    K=K+KK;
# I- X) z0 r6 V6 M1 f$ J0 L& Uend  
8 U* k0 e  t+ a! y%------------------------------------------------
1 [$ l6 X; o2 t8 S3 A% 约束条件,对整体刚度矩阵进行修正,以便计算
) n: C) s! ]. _F=zeros(Ndof,1);2 C& L0 w1 E8 _0 M8 p* n
F(Ndof-1)=-100;1 U' R$ j- l8 I4 R5 Z/ K* R2 [
% u=F*inv(K),u=[v1,xta1,v2,xta2,v3,xta3,v4,xta4]',v1=0,xta1=0
( ~+ M+ c! x1 R0 ~/ k0 K1 FK(1,=0;        K(:,1)=0;        %可以省略% F% ~9 c2 U( v7 c! y
K(2,=0;        K(:,2)=0;        %可以省略) Y8 |6 f! ^* t6 U, y4 P* d7 g$ p4 ~
KX=K(3:Ndof,3:Ndof);# y" Y% f8 r; a- d7 p$ ~" V
FX=F(3:Ndof,1);
- I0 b3 b7 d; n$ W8 _4 x%------------------------------------------------$ r# [' u5 t) L  g* i
%求整体节点位移列阵
3 p0 f1 _+ t/ f9 c1 p; Juu=inv(KX)*FX;' D, F* \# I; H1 F
u=[0;0;uu];' A, G/ U$ k8 B( z4 Y
ii=1:2:2*NJ;
  }! r+ s, j/ X9 z* [0 v) Jv=u(ii);                        %各节点挠度& X7 w* l8 `, S. f2 |6 T( S
xta=u(ii+1);                        %各节点转角
* x6 _5 M: Z# W9 t0 K! ~. J: [3 S( U%------------------------------------------------8 U) E* c( ~! a) f9 d* V( `
% 后处理,计算单元应变应力* t. _3 Q  v3 l! c
B=subs(B,x,'L');
6 Y0 c( b3 W/ X) |  ^B=subs(B,'L',10/Ne);; `" u6 ^3 F4 X2 U+ O* a9 ~  G
Strain=zeros(1,Ne);                %单元应变,并初始化
9 u) V1 H% v  z$ vStress=zeros(1,Ne);                %单元应力,并初始化7 \! }+ {8 v" A: P( j. M5 L% {/ j8 d
for ii=1:Ne
- V- M$ V' M+ R/ i2 m: Y: S/ G    Strain(1,ii)=B*u(2*ii-1:2*ii+2,1);: s' _' p: T1 A6 D( e5 U) q- a
    Stress(1,ii)=E*Strain(1,ii);. }' t$ p6 ^7 S* A
end; t& @4 I+ ]% Z. I7 ]' Z
%--------------------------------------------------
; O2 E3 y7 }; K/ c5 Y% 以下程序为屏幕输出处理
2 [" }0 V/ M' q3 ]; T8 Wdisp(' ');1 T% N7 ^- Q! x8 e
disp(' Input:1-print Node displacement ');- K" ?6 N2 P/ u9 J3 p
disp(' Input:2-print strain ');
. i3 n6 {2 d" j. |7 I% K% Bdisp(' Input:3-print stress ');
6 r$ H% b9 l4 M; V* F* a9 Q  Wdisp(' Input:q-quit or use ''ctrl+c'' Shortcut to end the program');. x" v+ z$ N; U

3 W" q3 O9 ], a: G( e. bwhile 1
9 U# V8 r6 U1 M7 h    ii=input('Please input1 or 2 or 3:','s');. x, r5 c: d3 S9 ?
        switch(ii)
7 f3 I, J# I1 Q7 J& P$ e' x            case '1'
: S+ @+ i  m, ~: a! D/ V! H$ d                disp('Node displacement');
! |6 i$ T6 C+ C$ C" d                disp(u');
7 u7 q# y& ?* |& ~# U2 J            case '2'" |9 w/ H& I+ }. j
                disp('element strain');( @' a- |1 k$ _$ e7 h5 }0 M# N
                disp(Strain);
  y) A: ^5 i4 l' O7 c+ A8 R7 r            case '3'
1 A2 w2 t4 e# V6 _                disp('element stress');
3 T" A5 [+ m, r8 z6 r                disp(Stress);
  W% y7 S3 b7 v/ _7 {5 @            case 'q'0 [* d: E$ K5 M& n; Q
                disp('End of program');
! F- d3 F, z$ y) X4 p( V5 o! E0 k                break;
; ?1 |3 Y) p' x2 P! n            otherwise
( ]4 j, |3 s/ i) P: U( e- @; Z                disp('error!Please input again');
; f. N. P/ }9 h( w" }& K0 d2 c) _# T% a                continue;
0 r; O7 m  }9 r' \/ o        end" ~2 S4 U1 l8 j1 d- M' l
end) h$ e& `3 y8 z& {; b
. F2 S' [/ Y0 i& H: `2 z2 {& M

5 `3 I% @2 `6 a3 @5 @. m
回复

使用道具 举报

2#
 楼主| 发表于 2013-3-24 13:44:32 | 只看该作者
那个地方怎么变成笑脸了,冒号+右括号=笑脸,改了下,应该是下面,把英文括号改为中文,就好了吧
1 ~$ g5 C' \, oK(1,:)=0;        K(:,1)=0;        %可以省略
# B) N2 |! [5 W' @& Q- }; m% ~, xK(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
1 `0 Z* }7 m# X5 }, v" p谢谢你啊,太感谢你了
( y! b- p1 d( u0 O1 V2 y
这谢啥呢?
回复 支持 反对

使用道具 举报

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;, S+ L, @7 a" Z' t3 _
v=x.^j;  是干啥的
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2025-2-19 06:50 , Processed in 0.076493 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4 Licensed

© 2001-2017 Comsenz Inc.

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