|
這是關于梁單元的,可能大家覺得很簡單。。。2 Z$ c8 \% {) a# l; d
今天翻電腦里的東西時發現的,是我大三時有限元時的作業,記得當時花了很多時間研究,學了不少東西,一個簡單的作業可以涉及各方面的知識。
. o, ~9 g, G6 t p9 B4 Q畢業半年了,雖然記不清彈性矩陣、剛度矩陣等,但是只要一看書(好像現在的心境不下來)應該會,學校交的學習方法和理解問題的能力,有些東西不必記,除非每天從事,那該叫熟能生巧。
4 |! Z7 d! e7 p8 R記得當時編了兩個程序,梁和三維實體的,我分享一下梁的程序吧,起碼有亮點和創意,可以自己選擇劃分個數,在matlab方面花了不少功夫。
' x" m: i4 m$ v6 R' x8 ?非常感謝當時的有限元授課老師“韓清凱”,雖然老師已經到了大工。上課的教材為韓清凱 的《彈性力學及有限元法基礎教程》,書上有個梁單元的例子,在82頁。7 u) W2 ~$ ?) I) B( D
--------------------------------------------------------------------------------( E8 w% ?* b4 X/ K/ ^
梁程序如下,名字就不寫了,用昵稱吧“太平島”,這是常用的網絡昵稱( \2 N* Y) x/ E$ L
: A7 L; f9 ~# u* [& w( ]" @5 P6 R4 H% O6 @8 \' { r& E
%------------------------ BEAM2 ----------------6 x$ C# d. m8 t
disp('========================================');+ U0 x. h# f9 u/ H
disp(' PROGRAM BEAM2 ');
1 ~6 ~' _- A2 l* N+ _' a+ Zdisp(' Beam Bending Analysis ');2 l2 b$ M p3 n' Z0 n2 ]# Y
disp(' The programmer:太平島 ');, o3 m4 t$ s2 i8 y d7 n
disp('========================================');
, v, g1 l6 `8 g7 V% u* Z1 N' fdisp(' '); %輸出空行
' A: V5 M( W& b! Zwarning off all %關閉所有警告提示(求積分時,警告信息比較多)
, b) s! f$ U5 s) p) e$ vNe=input('Please input the number of elements:'); %Ne為劃分的單元數
2 Z. f$ n3 D- X6 p. RNJ=Ne+1; %NJ為節點數
3 o8 y/ D) e- z/ n: _& \x1=0;$ f/ t% c4 p( f2 J* t! }, N
x2=sym('L');: k. K: K) C5 z9 D. U X
x=sym('x');
/ H" P! |/ U5 A& \1 n& e, }j=0:3;
9 ]/ t4 y9 C* j/ v" S, p! m7 Nv=x.^j;
) Q) g- }; U; n2 M7 s3 c7 b: lA=[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];. @0 Z. t! f6 l9 Q) m" _) W* r
N=v*inv(A); %形函數; E) S8 D) T: X' s& ~
%-----------------------------------------------
) M+ a4 T5 d1 O7 Z+ u% 求單元剛度矩陣
' ~; {' g' b/ ]1 S2 L0 q! d, KE=4.0e11;& Z0 S2 C4 M% Z: {7 d# U
I=5.2e-7; %I=bh^3/12=5.2e-7;
" ~. U7 c) ^' V8 S$ qEI=4.0e11*5.2e-7;
- N# C2 u# S* H y. U. H) X* o4 X7 }; rB=diff(N,2);0 K3 |" ~$ U! o
k=B'*B;+ F6 o* P% j9 Y3 |
Kee=EI*int(k,x,0,'L');
) y3 r9 v7 _ E- ~3 e8 C& }0 kKe=sym(zeros(4,4,Ne)); %用三維矩陣存放單元剛度矩陣,每頁存一個,并初始化0
4 f4 H& I' _$ l& |2 u; o& m2 fKe(:,:,1)=subs(Kee,'L',(10/Ne));
B) w7 T: q6 ^+ d1 cT=eye(4); % 因為梁于x軸的夾角為0,所以坐標變換矩陣為T=eye(4)+ H: M1 e! x6 G8 b/ ]
Ke(:,:,1)=T*Ke(:,:,1)*T';
/ l4 \7 x! `" D& B1 x, U/ B6 D7 tfor ii=2:Ne+ {2 M5 n8 k% t/ P1 r
Ke(:,:,ii)=Ke(:,:,1);: I0 r0 F4 B) u9 a' E- H
end 3 |/ K, t5 D, U
Ke=double(Ke);& `$ Y. P0 c' }$ ^6 X' g9 \ c
%------------------------------------------------; e* X& E3 L* q. {
% 由矩陣裝換法求整體剛度矩陣
/ M% n0 g/ R8 \7 z% 自由度Ndof=2*NJ
& }3 ?' f) j9 F& |Ndof=2*NJ;
. D' V. X) l! K- |5 s: `* M8 hK=zeros(Ndof);$ L" h- V1 C: ?2 T) \4 P- q
for ii=1:Ne u8 F3 z& ~" h5 d( U& W' e4 _
G=[zeros(4,2*(ii-1)),eye(4,Ndof-2*(ii-1))];
s3 e2 U2 Y8 p; S5 G KK=G'*Ke(:,:,Ne)*G;5 V8 w$ i. R/ l, H& P' o1 {
K=K+KK;) U6 C4 }7 ~5 d. y: l
end
5 F2 ^7 D. D' p: V- [%------------------------------------------------
; H" V, F) k: U4 Y& P7 N2 W3 `& y$ g% 約束條件,對整體剛度矩陣進行修正,以便計算
7 U( C# D8 K _& S4 H4 hF=zeros(Ndof,1);
7 c$ T; n2 ]( Z/ L( L& r6 h1 dF(Ndof-1)=-100;, {. }7 V0 |5 ?6 ] \7 M: ?3 V, i
% u=F*inv(K),u=[v1,xta1,v2,xta2,v3,xta3,v4,xta4]',v1=0,xta1=0
6 K! f2 u" y* W6 p/ p" r' EK(1, =0; K(:,1)=0; %可以省略
9 i- k& T$ W$ g6 z) e5 uK(2, =0; K(:,2)=0; %可以省略
# g2 I/ {! C# d, C: l# P# a0 zKX=K(3:Ndof,3:Ndof);% w) L# Z9 r& k/ {5 C1 X& F: T
FX=F(3:Ndof,1);
, ^- A+ f3 t7 `* U4 d%------------------------------------------------
% L8 C3 F- r! ~( y%求整體節點位移列陣
5 a$ K9 O- ]' I' i: L( Auu=inv(KX)*FX;+ j" `3 F9 F" k a( ?$ ~
u=[0;0;uu];
0 a; W) `, u3 y6 [6 @ii=1:2:2*NJ;
9 p0 s! B' n0 }4 X. U( ^4 w7 G) mv=u(ii); %各節點撓度
! d+ P, O$ U( `0 s4 V# `xta=u(ii+1); %各節點轉角
& N$ j$ Z) b2 ]2 d! g7 D9 D- L' m. L%------------------------------------------------
5 y2 ~, ]1 K/ M0 N- K% 后處理,計算單元應變應力 ]) _! B: g6 A: V- a" W& x/ w3 K/ ?% Q
B=subs(B,x,'L');4 A' b( v2 b/ q, N2 Y
B=subs(B,'L',10/Ne);
# B! D4 ^7 I) [Strain=zeros(1,Ne); %單元應變,并初始化6 P0 ]- w p: K4 F8 q, v# n' z$ S
Stress=zeros(1,Ne); %單元應力,并初始化
4 W4 v+ h- z! ?- R1 Dfor ii=1:Ne; F \# m* J! J! w
Strain(1,ii)=B*u(2*ii-1:2*ii+2,1);
& R- k: b, l& {5 S Stress(1,ii)=E*Strain(1,ii);
+ k! Q4 x) x: n0 E" Hend
6 Q" _& E! @3 k" H7 g) g8 Q%--------------------------------------------------
, U" c% E& ^' U- g0 o+ {% 以下程序為屏幕輸出處理
( J; J- g5 G1 xdisp(' ');: K) N' \1 e4 W: i
disp(' Input:1-print Node displacement ');
. G7 o3 h M6 u% }' n% B8 odisp(' Input:2-print strain ');
$ k; a. d3 Q1 n5 rdisp(' Input:3-print stress ');& W8 A/ l5 P& A- \4 A
disp(' Input:q-quit or use ''ctrl+c'' Shortcut to end the program');) i8 c; s) w9 D9 c
: s C4 V! ^ i) l; K0 e6 lwhile 1+ E6 S) W% U4 O* q/ n
ii=input('Please input1 or 2 or 3:','s');* p. D7 v% E, A2 {+ [3 E- D
switch(ii)* r1 c; _1 f/ s/ b4 H# `4 [
case '1'4 n2 K, R p }4 o
disp('Node displacement');3 c- [& t, ^! Q7 W: ]# u) u
disp(u');* Y+ z4 @( t3 p
case '2'
8 x3 ] H% i: N5 ` disp('element strain');9 K& e" F( P l g
disp(Strain);
% Z8 V2 E- J& b& C case '3'
( N; L: P8 G: U! _3 t% t disp('element stress');% a( D! @, o' d: Z
disp(Stress);
L [" ]) Q; ]% i1 d r& q case 'q'
' R" v5 ?. S5 o" R disp('End of program');- M1 z. j5 R! Z2 k' ?
break;; f D9 G/ v2 V7 m3 U. s
otherwise/ F* ?" H0 F2 F& j+ T
disp('error!Please input again');
' \8 {1 D! i v8 p# K$ S continue;6 M0 F& h+ `( }
end
' d+ u4 `5 F; U1 h7 ]1 eend
! J- x9 |& @7 f. \6 }
6 q7 _7 J" b$ d
! q8 D/ j+ F& E" j |
|