机械必威体育网址
标题:
分析一段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' ?& k
disp('========================================');
; t! F% f, Q9 \" ^* [ h
disp(' PROGRAM BEAM2 ');
) f- r1 @' i3 W
disp(' Beam Bending Analysis ');
3 o4 f' t( e8 }& C6 J# r& Y
disp(' 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" n
Ne=input('Please input the number of elements:'); %Ne为划分的单元数
3 w. Q& \; K; W! u. S
NJ=Ne+1; %NJ为节点数
4 v* p1 x0 d; a
x1=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 ?" S
v=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 w
E=4.0e11;
/ v, t/ r6 m- o, }4 q8 m- v
I=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 u
Kee=EI*int(k,x,0,'L');
6 L& B. v* u8 M. M3 I+ Z4 C
Ke=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( T
for ii=2:Ne
! f( T5 \5 n$ W
Ke(:,:,ii)=Ke(:,:,1);
- V! m0 S* ~0 K: \) U8 L
end
8 b+ c" _( g3 S; b: ?( U
Ke=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) m
Ndof=2*NJ;
6 u9 l6 m6 g$ _. m
K=zeros(Ndof);
3 G, H5 w3 ?, u
for 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( O
F(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 x
K(1,
=0; K(:,1)=0; %可以省略
+ T O7 k$ a" c& H1 s, K8 f
K(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% `( E
uu=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' [) n
xta=u(ii+1); %各节点转角
. g5 Y8 Q+ }3 i$ b; B3 f( t
%------------------------------------------------
/ [" k, `# I7 n6 m
% 后处理,计算单元应变应力
' {* J3 C& n& B: p& R4 q! O
B=subs(B,x,'L');
" N$ ?; n0 @! y1 Q
B=subs(B,'L',10/Ne);
7 M; o3 ?1 s5 Y# z5 _* K
Strain=zeros(1,Ne); %单元应变,并初始化
0 ^- y0 n4 J5 V6 C
Stress=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: o
end
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* g
disp(' Input:1-print Node displacement ');
2 K' a& G+ y1 y& `0 E
disp(' Input:2-print strain ');
# S; Q1 [( y/ A0 P
disp(' 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 l
K(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+ J
v=x.^j; 是干啥的
欢迎光临 机械必威体育网址 (//www.szfco.com/)
Powered by Discuz! X3.4