|
这是关于梁单元的,可能大家觉得很简单。。。/ S! q5 k9 p8 M" ^9 Z* S8 |" G+ g J 今天翻电脑里的东西时发现的,是我大三时有限元时的作业,记得当时花了很多时间研究,学了不少东西,一个简单的作业可以涉及各方面的知识。8 z1 b9 B5 p- C4 Z# s0 x F. O3 x 毕业半年了,虽然记不清弹性矩阵、刚度矩阵等,但是只要一看书(好像现在的心境不下来)应该会,学校交的学习方法和理解问题的能力,有些东西不必记,除非每天从事,那该叫熟能生巧。! [. n4 R6 l, S: Y5 L5 _6 N' J' ] 记得当时编了两个程序,梁和三维实体的,我分享一下梁的程序吧,起码有亮点和创意,可以自己选择划分个数,在matlab方面花了不少功夫。) ?& `0 d) P; ^1 O( J6 z9 Z* M: @ 非常感谢当时的有限元授课老师“韩清凯”,虽然老师已经到了大工。上课的教材为韩清凯 的《弹性力学及有限元法基础教程》,书上有个梁单元的例子,在82页。 ' b/ ^2 ?- {) P% D! {-------------------------------------------------------------------------------- 1 N- ]6 n- U& r梁程序如下,名字就不写了,用昵称吧“太平岛”,这是常用的网络昵称5 Z5 f" ~- @# |8 \6 X- q: @
6 m: m( z0 u5 T( q- t& y* w ' @8 A. \6 W. P, X%------------------------ BEAM2 ---------------- % V3 r6 J6 \, K$ z+ S& }disp('========================================'); 2 u! R3 L9 Q' f5 P: y+ ?disp(' PROGRAM BEAM2 ');7 ^! ? \' `4 R) e5 L! ^ k4 A5 p/ c. P disp(' Beam Bending Analysis ');' w; i/ {( O4 g1 W8 G disp(' The programmer:太平岛 ');/ J8 s) ]2 @$ P6 e! B disp('========================================'); " F4 \8 J6 A5 \4 a W: y [: pdisp(' '); %输出空行 . k/ L$ Z" V* m' q- \8 O7 Wwarning off all %关闭所有警告提示(求积分时,警告信息比较多) 4 C7 G/ y: f; V# \/ Y5 p* V( LNe=input('Please input the number of elements:'); %Ne为划分的单元数7 G+ b6 q! I, m& U NJ=Ne+1; %NJ为节点数% `+ A0 z$ u6 g9 A* v1 v: O x1=0; ! F. b* z0 ?/ C8 yx2=sym('L');- a& ?8 B- E; E1 ?. J8 Y" N x=sym('x'); ! _ [) u% S* Z+ H, a6 T. Sj=0:3;2 [0 }/ B! g# s- y* u' w v=x.^j; 0 V& Y! T; V( A+ {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]; 6 }8 c; d& p6 N! @5 \7 tN=v*inv(A); %形函数 ) C5 {8 M( q. `, t1 ?2 _& a%----------------------------------------------- - C3 f5 N; _* v% 求单元刚度矩阵 3 z3 e8 @. I, s8 z/ AE=4.0e11; ( H7 j8 D# I) h$ pI=5.2e-7; %I=bh^3/12=5.2e-7; * v9 [% @ Y3 Z! p. OEI=4.0e11*5.2e-7; X: J' v1 x' P" N0 {. O7 K6 A! fB=diff(N,2);6 j1 { o v5 ?# Q2 [, X k=B'*B;8 z9 U" c& C+ X- @- @. M1 d+ b/ V) @7 e Kee=EI*int(k,x,0,'L'); 0 Y* R( A) C+ b- _2 j7 bKe=sym(zeros(4,4,Ne)); %用三维矩阵存放单元刚度矩阵,每页存一个,并初始化0) L f7 { L5 C; I# _( A Ke(:,:,1)=subs(Kee,'L',(10/Ne)); # `' C; Y+ t* j% l1 J- X( Q- dT=eye(4); % 因为梁于x轴的夹角为0,所以坐标变换矩阵为T=eye(4) 5 k, A" \7 Q7 j5 A- XKe(:,:,1)=T*Ke(:,:,1)*T'; F4 I1 K. |7 d5 Y" O8 v6 Y i+ mfor ii=2:Ne ; l2 A4 u/ I' `9 I, j" H, zKe(:,:,ii)=Ke(:,:,1);" O; |0 D) w$ k" ~% A end1 y, K6 | ^0 _3 R Ke=double(Ke);# @% |5 w% f. x7 t- f %------------------------------------------------ 8 d) t3 G' x8 k; |( W7 W% 由矩阵装换法求整体刚度矩阵' F5 g- P* C7 e7 r( d' p % 自由度Ndof=2*NJ 4 z" P/ W1 C$ m+ @$ pNdof=2*NJ; 3 r7 e* e5 O% u2 }: RK=zeros(Ndof);0 ?/ I5 u% k# s8 Z for ii=1:Ne 9 M7 E Y1 N: ?3 m, UG=[zeros(4,2*(ii-1)),eye(4,Ndof-2*(ii-1))]; 4 ?- w8 Q! U1 n AKK=G'*Ke(:,:,Ne)*G;9 S& @" e! A/ R6 T5 ~- J2 j K=K+KK;; Z, j6 Q" s/ P2 E1 E end: l) C9 _! N8 g6 _! d; o %------------------------------------------------ + L" ]. g! }* g8 h @% 约束条件,对整体刚度矩阵进行修正,以便计算/ ~8 M( x3 {: k) l9 `" C F=zeros(Ndof,1); 3 D, f. @( `7 E7 OF(Ndof-1)=-100; " k6 h' ]! r6 I' |' m! k s x% u=F*inv(K),u=[v1,xta1,v2,xta2,v3,xta3,v4,xta4]',v1=0,xta1=0 + Y0 ~) U. o6 m; T0 l- RK(1,=0; K(:,1)=0; %可以省略. s: p/ J3 K+ C* v K(2,=0; K(:,2)=0; %可以省略 " V" f. [4 P" Q: C# e, X9 gKX=K(3:Ndof,3:Ndof);" Z; {- X0 C) k8 J$ v( @ FX=F(3:Ndof,1);# ?. m! P/ ~6 k/ F$ u* q: _ %------------------------------------------------" x# H+ W$ H0 ?8 g7 \ %求整体节点位移列阵 K, }: ]2 p( E! duu=inv(KX)*FX; a/ l0 ^, d8 a1 }) R5 ?$ Xu=[0;0;uu];/ g& X! {* h* l ii=1:2:2*NJ; % V0 p% K6 F5 B* qv=u(ii); %各节点挠度- Z5 r) w+ E5 g- N xta=u(ii+1); %各节点转角& c( g; y! e7 M %------------------------------------------------ , ?6 b& ?9 ?$ `% 后处理,计算单元应变应力$ _# j3 z6 a1 w/ i B=subs(B,x,'L'); & w" O0 M% l1 Q# xB=subs(B,'L',10/Ne); / }, @+ p7 @4 v& OStrain=zeros(1,Ne); %单元应变,并初始化$ i4 w+ X, u+ B+ c7 e Stress=zeros(1,Ne); %单元应力,并初始化 ' q8 z8 z7 [1 f8 Wfor ii=1:Ne " x. [3 O8 n. F& TStrain(1,ii)=B*u(2*ii-1:2*ii+2,1);# T! G4 n/ \+ ^" ?+ z* g# I Stress(1,ii)=E*Strain(1,ii);+ `& k1 z7 |" U$ L end 3 e, z8 n$ D E2 R1 o" j%-------------------------------------------------- * U3 G' v! ?1 z/ I/ x% 以下程序为屏幕输出处理2 R/ l6 C6 M# r! q/ P# n8 U B disp(' ');! c% w, b$ m8 m. d/ D; I disp(' Input:1-print Node displacement '); 3 _7 O/ T$ A0 bdisp(' Input:2-print strain '); ! D+ L( ~' p. G! V2 F" Vdisp(' Input:3-print stress ');4 n6 J( O7 _+ r5 G disp(' Input:q-quit or use ''ctrl+c'' Shortcut to end the program');1 m. U9 ?; R* b7 ]- f
. ?+ w0 W: t" g7 G. l+ e3 e5 L9 Swhile 1' v# ]+ y, ^9 Y7 R2 a6 ]& S ii=input('Please input1 or 2 or 3:','s'); ) ]5 T" h7 |: Bswitch(ii)# g& ^! c3 b- P: f9 |. a case '1' % U; h# d/ l* O: k n4 \disp('Node displacement');: f. e# v% x0 w6 p disp(u'); 0 J3 p+ Y; [) P2 qcase '2' . ~" _" _1 G5 a" Xdisp('element strain'); 2 s* g7 A$ z0 g9 K1 |' `disp(Strain);4 d+ }, {' ~3 J% H case '3') i( e4 [1 J2 [% N' Y9 p: l disp('element stress');# C) r# L8 W2 l( i. Z U disp(Stress);3 L0 M: Y p6 J6 } case 'q' , X& N0 }: ?6 t+ j# vdisp('End of program'); 3 F0 }7 n9 n3 abreak; & S+ L% ]; R# y; o kotherwise: k, m9 w$ o3 r1 |5 N disp('error!Please input again');) L0 S, [5 \# w6 G8 i9 C continue; + V/ H2 I8 [2 `& G7 \$ b8 `end8 W2 M5 f5 G4 c# X4 m5 G% f; P end 4 I4 i2 N, {; A1 V" ^8 j # K, m$ b, T4 Z/ r( ` 9 d/ |2 Z0 z4 ~! z$ i' {- G6 Z3 Q9 S |
|