|
這是關(guān)于梁單元的,可能大家覺得很簡單。。。) P! K2 m" F+ [: T* K" B
今天翻電腦里的東西時(shí)發(fā)現(xiàn)的,是我大三時(shí)有限元時(shí)的作業(yè),記得當(dāng)時(shí)花了很多時(shí)間研究,學(xué)了不少東西,一個(gè)簡單的作業(yè)可以涉及各方面的知識(shí)。$ V4 i3 B) [+ K# y3 i' o! j* s
畢業(yè)半年了,雖然記不清彈性矩陣、剛度矩陣等,但是只要一看書(好像現(xiàn)在的心境不下來)應(yīng)該會(huì),學(xué)校交的學(xué)習(xí)方法和理解問題的能力,有些東西不必記,除非每天從事,那該叫熟能生巧。6 u# K- B2 H7 H0 g
記得當(dāng)時(shí)編了兩個(gè)程序,梁和三維實(shí)體的,我分享一下梁的程序吧,起碼有亮點(diǎn)和創(chuàng)意,可以自己選擇劃分個(gè)數(shù),在matlab方面花了不少功夫。
' m: c; L. J7 V3 K1 P/ _) v1 \非常感謝當(dāng)時(shí)的有限元授課老師“韓清凱”,雖然老師已經(jīng)到了大工。上課的教材為韓清凱 的《彈性力學(xué)及有限元法基礎(chǔ)教程》,書上有個(gè)梁單元的例子,在82頁。
% w2 R' z/ F2 P3 f# t--------------------------------------------------------------------------------
( Q0 H- j3 r' Q梁程序如下,名字就不寫了,用昵稱吧“太平島”,這是常用的網(wǎng)絡(luò)昵稱
4 C, ?* V. B" a. a% p
) I/ F5 I0 d7 l( Q9 [- g; S# t4 U: v# T ?$ k3 n
%------------------------ BEAM2 ----------------
" E) J/ b% S l* I9 v, Xdisp('========================================');
- @ n) [5 J/ I+ f/ S) Ndisp(' PROGRAM BEAM2 ');: s1 I" J5 ~8 F: [; }4 ] i* E
disp(' Beam Bending Analysis ');
! s: v. h" P6 \6 L& S7 M5 G9 ldisp(' The programmer:太平島 ');! u( E( h( Y5 D3 X* I2 q! `% B8 s- {: c
disp('========================================');
2 _' B" ]- [8 f- F$ P& Xdisp(' '); %輸出空行
. i6 R4 ]. p& T& ^: F4 x6 zwarning off all %關(guān)閉所有警告提示(求積分時(shí),警告信息比較多)
& t2 D+ V5 O+ i; ^- u7 ONe=input('Please input the number of elements:'); %Ne為劃分的單元數(shù)
7 e0 s& v. |4 H' h. ]" n0 U) q5 ENJ=Ne+1; %NJ為節(jié)點(diǎn)數(shù)
9 Q; N% r% ~+ n" N: v, D& w) G# ex1=0;
: i5 ~2 p3 H9 V# G9 }x2=sym('L');' E" u2 t; R) A) G
x=sym('x');
) l4 u* q% Z7 w7 Z2 {# |j=0:3;- t" i0 W3 @& h7 A( J
v=x.^j;
, j9 \! g8 W/ m+ f3 i2 ]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];
$ J3 Y. }. k& JN=v*inv(A); %形函數(shù); S7 f# ~- d" q- ~
%-----------------------------------------------
! o* G, f$ B8 p) Q% 求單元?jiǎng)偠染仃?font class="jammer">5 w4 i+ t, V+ K+ s
E=4.0e11;7 |" k; m y; h; G
I=5.2e-7; %I=bh^3/12=5.2e-7;
: i. z: `8 S7 ^4 N+ t9 ~" z+ `EI=4.0e11*5.2e-7;( V+ |/ ^& j r+ |9 U6 t1 _$ K
B=diff(N,2);
% N1 t( i8 v H: w- N7 f. ~k=B'*B;
% T0 L; W* a5 F9 i- J* b1 G. xKee=EI*int(k,x,0,'L');
6 m& u- d7 r3 y3 y7 {( n% I, TKe=sym(zeros(4,4,Ne)); %用三維矩陣存放單元?jiǎng)偠染仃?每頁存一個(gè),并初始化0
u! y, L9 D, h7 HKe(:,:,1)=subs(Kee,'L',(10/Ne));
" B4 v3 u$ d: OT=eye(4); % 因?yàn)榱河趚軸的夾角為0,所以坐標(biāo)變換矩陣為T=eye(4)' |; ~- \1 p0 v* q$ R0 ?( U
Ke(:,:,1)=T*Ke(:,:,1)*T';/ R* t" g& l0 Z# |( z4 }
for ii=2:Ne! }! @0 Z, i; Y3 @5 _. u
Ke(:,:,ii)=Ke(:,:,1);7 I' q, A: W9 ?
end
) ?# ~" x. |; H S% [( XKe=double(Ke);& S2 F) `8 P- {" W
%------------------------------------------------/ q2 @) Z: T. P2 T5 r6 t8 P
% 由矩陣裝換法求整體剛度矩陣
8 s4 x' t5 A6 Q. |, j% 自由度Ndof=2*NJ2 L) o ~& V/ g) T% d% _1 u% N
Ndof=2*NJ;
3 S+ ]- e& E$ P+ l2 x4 \% i$ oK=zeros(Ndof);
) j6 v: B) F! N+ `& i) Pfor ii=1:Ne
( P4 y# s( f4 Z% g7 H9 }1 O6 e, |2 C G=[zeros(4,2*(ii-1)),eye(4,Ndof-2*(ii-1))];* y/ H6 @/ C7 f1 T1 [/ e6 d" a4 V
KK=G'*Ke(:,:,Ne)*G;4 {1 k# i2 v) M& j1 I; N1 H4 \
K=K+KK;
1 W5 B1 E) i+ A# F& Vend 6 [0 N2 A0 u6 k. ?& y$ h% A) Z `
%------------------------------------------------
. ?% f9 a- t; Q% 約束條件,對整體剛度矩陣進(jìn)行修正,以便計(jì)算
+ e( N- m# I, A0 z7 L2 q( x; k3 \F=zeros(Ndof,1);. \# P; [$ r" q8 D V0 E
F(Ndof-1)=-100;
1 [- W9 Z2 | ~1 a3 Z8 ~2 C% u=F*inv(K),u=[v1,xta1,v2,xta2,v3,xta3,v4,xta4]',v1=0,xta1=0) Y# [2 @. V( ?$ b- K% J8 Z
K(1, =0; K(:,1)=0; %可以省略
i8 v8 \1 ? {1 g, D+ DK(2, =0; K(:,2)=0; %可以省略
9 `( e5 T4 Z; q+ m( mKX=K(3:Ndof,3:Ndof);
7 H0 n) O, w" x6 NFX=F(3:Ndof,1);
! ~0 u ]1 \" g' q9 j* ~%------------------------------------------------
* Y; p1 h" V$ U+ c7 O, C" B%求整體節(jié)點(diǎn)位移列陣5 `; J2 X. c: M _! _
uu=inv(KX)*FX;/ Z( N u4 [& R+ }; v* M
u=[0;0;uu];
& X8 V+ r) P$ k9 f. w4 Lii=1:2:2*NJ;* y+ p: Z+ v! o2 j
v=u(ii); %各節(jié)點(diǎn)撓度
" c. I. M' L @0 f( Cxta=u(ii+1); %各節(jié)點(diǎn)轉(zhuǎn)角" d a9 A* T& a( E
%------------------------------------------------ ^$ e& l: N" l4 s4 p
% 后處理,計(jì)算單元應(yīng)變應(yīng)力
1 U& Y/ C( T* C: H+ kB=subs(B,x,'L');" L. d! X; F: g" |
B=subs(B,'L',10/Ne);2 d, G. k; w" R( B* Q- B
Strain=zeros(1,Ne); %單元應(yīng)變,并初始化- `- V! R9 I0 e9 Q
Stress=zeros(1,Ne); %單元應(yīng)力,并初始化
/ x4 s: A4 i% Rfor ii=1:Ne0 q c; V, V/ c
Strain(1,ii)=B*u(2*ii-1:2*ii+2,1);
& l0 Z2 t% N l# I# h$ _9 L/ d Stress(1,ii)=E*Strain(1,ii);: B1 M: d" B5 L" A5 h
end
1 c0 ^2 {: y( B8 l; D%--------------------------------------------------
0 U. x! v$ v( ^" a1 r$ e* U2 u% 以下程序?yàn)槠聊惠敵鎏幚?br />
" @( q' p- \& r& K/ j5 Odisp(' ');9 Z2 y6 ?; c& R' i
disp(' Input:1-print Node displacement ');% G/ A6 v7 z# y4 Q4 J+ Y" ?- X. `
disp(' Input:2-print strain ');& H/ ]8 U! U0 l9 s' R
disp(' Input:3-print stress ');* n @- h+ s1 `/ w# B" C
disp(' Input:q-quit or use ''ctrl+c'' Shortcut to end the program');
& ^: P! M! r: Z, p% i3 l8 ?5 S
) R u3 e4 R$ \- |4 C' `; y7 u/ Swhile 1/ x# L. B9 g* Z8 {* Y7 Z. m# K' n
ii=input('Please input1 or 2 or 3:','s');
4 n6 d, } f [4 [8 I2 u1 G1 {, `$ x switch(ii)0 {# m; f$ @% ? E% [/ s' \, P
case '1'0 P! n x, U% T' ]4 H. E+ \/ g+ J
disp('Node displacement');% z# i6 S1 H# R9 }3 X7 K# a
disp(u');! e2 P' L- K5 w& j; J
case '2'; V. |1 R( A! `7 i* z+ q( J; h$ D- s
disp('element strain');" K5 V3 k6 z, w$ I8 G v
disp(Strain);) C& g' `$ f% x
case '3'- |# w8 |' A9 n1 T5 ?
disp('element stress');- b) [4 a% V6 e ~* }4 {/ `
disp(Stress);# ]* N6 U; |' o! z& [; @
case 'q'8 ?" W/ F. r- W Q% W* r
disp('End of program');
8 `8 g# b) C5 J, M/ F% O break;
. F9 X% x& w8 @ n( s" C otherwise! R. x: E* B& D0 a) s9 C
disp('error!Please input again');( F* E% D' E2 f3 r, x8 n( G
continue;
- X0 n; O M3 j Q end- F4 O4 K6 f: a+ }0 M+ M
end. O8 P% l7 Q0 x( v5 E% C( {7 r8 I# h, f
( `; I' i7 g: l3 {9 k
. S; z% q6 R) ]; F7 H4 q& A |
|