|
這是關于梁單元的,可能大家覺得很簡單。。。
; [+ a. X& e2 C5 k' v( l今天翻電腦里的東西時發現的,是我大三時有限元時的作業,記得當時花了很多時間研究,學了不少東西,一個簡單的作業可以涉及各方面的知識。
0 o- o" Y% F( L$ S K畢業半年了,雖然記不清彈性矩陣、剛度矩陣等,但是只要一看書(好像現在的心境不下來)應該會,學校交的學習方法和理解問題的能力,有些東西不必記,除非每天從事,那該叫熟能生巧。+ _. |0 Z# ], Z, j" S. }
記得當時編了兩個程序,梁和三維實體的,我分享一下梁的程序吧,起碼有亮點和創意,可以自己選擇劃分個數,在matlab方面花了不少功夫。. }+ G/ t' z6 {7 H
非常感謝當時的有限元授課老師“韓清凱”,雖然老師已經到了大工。上課的教材為韓清凱 的《彈性力學及有限元法基礎教程》,書上有個梁單元的例子,在82頁。
. S9 E, {. b% |--------------------------------------------------------------------------------
# V- L5 Q/ `! y* {8 c5 n! N% F* Q梁程序如下,名字就不寫了,用昵稱吧“太平島”,這是常用的網絡昵稱
, }8 Q+ w+ Y: A& A5 h0 ~/ ~- L2 [( r: R& ]( F" }' S; q% h. P! [
7 t5 ]. d; q3 O0 D! s& b%------------------------ BEAM2 ----------------
9 f3 w& n, T( a9 \- H6 ndisp('========================================');: I) e. [5 g% J5 J5 h
disp(' PROGRAM BEAM2 ');
+ p. f( ^3 S8 @9 vdisp(' Beam Bending Analysis '); v8 X7 x- i1 P5 m
disp(' The programmer:太平島 ');
; L/ l1 \5 M" J3 C, b( H8 cdisp('========================================');6 q& K+ U& |3 n- ]
disp(' '); %輸出空行3 N/ J- z/ z; t$ d* L2 |) S2 G
warning off all %關閉所有警告提示(求積分時,警告信息比較多)
# G1 D' [/ N0 L* z1 tNe=input('Please input the number of elements:'); %Ne為劃分的單元數6 j! t* A# u% b- s! ]* c/ T7 Z0 W
NJ=Ne+1; %NJ為節點數
2 K0 u; _. J( X, ?5 Hx1=0;
5 ]5 o$ ^- ^1 \x2=sym('L');8 J/ k8 O/ z+ b5 |
x=sym('x'); % H/ z1 L9 H9 H) o
j=0:3;' L- @' `" e3 ?- F8 L, g, D3 }
v=x.^j;
$ c; ]& }* d8 d1 WA=[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];
% [2 T, G6 Y- L5 KN=v*inv(A); %形函數8 o( a0 T e! H$ l% E
%-----------------------------------------------
) i ~2 E; E }0 L$ @6 v" R, g% 求單元剛度矩陣
6 D6 J8 b3 ?$ u3 `1 oE=4.0e11;6 S# \5 P; L- F$ H0 c
I=5.2e-7; %I=bh^3/12=5.2e-7;
. m L* ^, H3 e" s4 T) `/ \EI=4.0e11*5.2e-7;1 Y. z& @6 v' o r, ?; n' X
B=diff(N,2);; A% R4 `6 j# T" g3 u6 }, g
k=B'*B;+ a+ x" B! I5 |+ Z D. \) c! T0 ^6 h
Kee=EI*int(k,x,0,'L');- h; t7 u7 \ n k
Ke=sym(zeros(4,4,Ne)); %用三維矩陣存放單元剛度矩陣,每頁存一個,并初始化0
6 ^5 e& V" B% d! V7 g |Ke(:,:,1)=subs(Kee,'L',(10/Ne));$ s4 C7 g9 w5 X7 {* o
T=eye(4); % 因為梁于x軸的夾角為0,所以坐標變換矩陣為T=eye(4)
1 x7 I9 I3 ~! ?" M8 uKe(:,:,1)=T*Ke(:,:,1)*T';
" U5 _1 ]" n* ^! P8 K1 L2 c5 ]for ii=2:Ne
* _0 n/ U& Q! M: h k Ke(:,:,ii)=Ke(:,:,1);
! J8 v8 B. p' {( M: ~end $ F1 x; M6 o' g" \
Ke=double(Ke);) |: X7 B+ @3 y, a
%------------------------------------------------
7 p7 N$ a8 i% x$ F7 b% 由矩陣裝換法求整體剛度矩陣
. n! ^* F2 k/ Q6 N" f% 自由度Ndof=2*NJ* ]( y4 ~# J3 g
Ndof=2*NJ;6 q6 ~! V$ W; ?& a2 I4 W2 M8 l8 S* w
K=zeros(Ndof); _% H5 w5 }/ w7 f: P! K
for ii=1:Ne
& v, r- Q& e9 u2 c' N9 \ G=[zeros(4,2*(ii-1)),eye(4,Ndof-2*(ii-1))];7 D S6 D, R: G+ I N' m$ ^
KK=G'*Ke(:,:,Ne)*G;
- l9 M7 i" x0 x K=K+KK;$ _) J" H j. W
end
" r2 z7 @9 `& s0 z% e%------------------------------------------------
7 n1 X- D, `/ `0 F; N$ }+ `$ a8 w* c1 L% 約束條件,對整體剛度矩陣進行修正,以便計算
9 f' Q N& H! V: ^+ L6 z: V7 |5 k9 ~F=zeros(Ndof,1);* W4 g, k% h8 J2 J
F(Ndof-1)=-100;5 C/ c+ n# t( t! P7 o' _
% u=F*inv(K),u=[v1,xta1,v2,xta2,v3,xta3,v4,xta4]',v1=0,xta1=0- R+ A$ h. J/ I* ?
K(1, =0; K(:,1)=0; %可以省略
7 r6 I9 e8 q* v7 e* oK(2, =0; K(:,2)=0; %可以省略2 o$ L- ]: t1 {$ n5 D
KX=K(3:Ndof,3:Ndof);3 [8 t/ w' e* d2 ?; A
FX=F(3:Ndof,1);
3 Q9 A" ~0 @ Q$ e3 _- h%------------------------------------------------
' A! l! R+ A7 w* N%求整體節點位移列陣2 Q7 w4 y5 i9 J$ K) P( @$ A
uu=inv(KX)*FX;
# \- }' T' W7 r ~u=[0;0;uu];; |" f2 u: {/ T1 ~1 d
ii=1:2:2*NJ;
: U& R5 g9 O* ^" K* Vv=u(ii); %各節點撓度
" I+ m% |" D5 cxta=u(ii+1); %各節點轉角( x4 @* L$ w. e" _ ]' {7 {
%------------------------------------------------
8 t K, b s3 J( \; V% ~2 w' ?% 后處理,計算單元應變應力
3 ~% w2 q4 D0 A$ B1 sB=subs(B,x,'L');6 l9 y+ x4 t3 @1 n! G* i
B=subs(B,'L',10/Ne);
: }( |/ a# c; a- {. b" Q9 Y; bStrain=zeros(1,Ne); %單元應變,并初始化, X$ K6 T2 a5 ~
Stress=zeros(1,Ne); %單元應力,并初始化. \! E. t; h# \+ x+ v3 A: @
for ii=1:Ne7 p r" U1 N, G& q( X: w
Strain(1,ii)=B*u(2*ii-1:2*ii+2,1);
# u# I) R4 O6 }; F7 { Stress(1,ii)=E*Strain(1,ii);, Y, ?7 C9 P1 k
end( j: g5 ?6 H$ z( R
%--------------------------------------------------/ c3 X" s3 \- W9 _$ e4 d
% 以下程序為屏幕輸出處理
/ V5 Z0 \& X: o0 `! Ldisp(' ');8 f( V7 a7 r% K- T
disp(' Input:1-print Node displacement ');
) O- r- e3 {8 ], D# H0 D- f6 pdisp(' Input:2-print strain ');
$ y' t3 ~5 z8 _/ l4 f- vdisp(' Input:3-print stress ');
- {; ]5 ]/ e( M* O6 g* gdisp(' Input:q-quit or use ''ctrl+c'' Shortcut to end the program');
V1 r% ^! }) o/ v/ Z5 V) s, P
G( C8 T) |% wwhile 1
4 N1 h U* M' b3 i) F8 [+ _6 i. I ii=input('Please input1 or 2 or 3:','s');- h2 H/ L6 N7 T# \& l3 r- E* |$ Z
switch(ii)8 @# L$ N; G1 ^! k! J* V9 D
case '1'
7 D/ g% T! x7 E. E6 v disp('Node displacement');
. F3 B/ s3 k$ o8 o) P disp(u');
7 O& u% V* H* r) p' ~3 d* i case '2'
2 \2 w5 {5 F" [+ e# V disp('element strain');3 w4 a, K2 A; s* h
disp(Strain);
Y$ X2 B, [! s# |/ j' L: Y/ C case '3'+ M, @. K* U$ s. B# e# D
disp('element stress');
5 W) c9 b$ k, i; C" ^% w0 }. T% b5 g disp(Stress);
& l* ?; X1 G% J" a% c% [ case 'q'
" I& R( _2 h5 N# {: q5 g disp('End of program');% W( A$ L( M0 x+ T
break;
" ^6 R4 q* p/ K9 V6 I otherwise# E+ S3 d! h' \% k" B: B, M5 W
disp('error!Please input again');$ K# z D7 S4 `( X! s6 @
continue;
6 K2 x8 ]. Z5 j0 W7 t: t, ~ end3 P% N$ }7 K1 c/ k' O* K* z4 h$ \
end$ v5 j. {' b$ D9 C5 C# j$ _& d
+ V; I7 ?! I a* h- a
% z) t( k/ ~- q( K
|
|