|
这是关于梁单元的,可能大家觉得很简单。。。. 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
|
|