|
這是關(guān)于梁單元的,可能大家覺得很簡單。。。
4 ?2 }4 O4 \+ X, H今天翻電腦里的東西時發(fā)現(xiàn)的,是我大三時有限元時的作業(yè),記得當(dāng)時花了很多時間研究,學(xué)了不少東西,一個簡單的作業(yè)可以涉及各方面的知識。3 h8 ?( i/ t) L) P8 H9 m
畢業(yè)半年了,雖然記不清彈性矩陣、剛度矩陣等,但是只要一看書(好像現(xiàn)在的心境不下來)應(yīng)該會,學(xué)校交的學(xué)習(xí)方法和理解問題的能力,有些東西不必記,除非每天從事,那該叫熟能生巧。6 K% ?' O5 Y/ T U/ [
記得當(dāng)時編了兩個程序,梁和三維實體的,我分享一下梁的程序吧,起碼有亮點和創(chuàng)意,可以自己選擇劃分個數(shù),在matlab方面花了不少功夫。
# D& M" l+ O6 m非常感謝當(dāng)時的有限元授課老師“韓清凱”,雖然老師已經(jīng)到了大工。上課的教材為韓清凱 的《彈性力學(xué)及有限元法基礎(chǔ)教程》,書上有個梁單元的例子,在82頁。# o/ r6 Z% n* R: d2 h- A
--------------------------------------------------------------------------------, k$ E5 \0 d9 u0 I
梁程序如下,名字就不寫了,用昵稱吧“太平島”,這是常用的網(wǎng)絡(luò)昵稱
7 u: X: w; v' R1 d, ]: q- _! ~- L" R" H0 _% F
# F) D" H2 T7 T ^ `, a
%------------------------ BEAM2 ----------------
8 f0 P! e- p. D( g( t Adisp('========================================');
! n& `- b/ `; L C* zdisp(' PROGRAM BEAM2 ');9 z6 |( O, S1 K+ ~ A& s+ B
disp(' Beam Bending Analysis ');
! A6 s; X9 ^+ udisp(' The programmer:太平島 ');
4 X1 ^* N1 Z" P* t# F' l/ O! ndisp('========================================');
) U3 k" R- v% ~' Fdisp(' '); %輸出空行7 H! a4 v- o& G6 ~% P
warning off all %關(guān)閉所有警告提示(求積分時,警告信息比較多)
( N' h1 S6 P, KNe=input('Please input the number of elements:'); %Ne為劃分的單元數(shù)
/ v6 `1 d' p& \NJ=Ne+1; %NJ為節(jié)點數(shù)9 S: e9 `. N0 o
x1=0;
+ Y* T, d& A% w& Y+ y" C1 Vx2=sym('L');
! S: R. O# x5 q% H2 F' |0 Cx=sym('x'); 8 d' [$ b; z( J4 e. H8 R
j=0:3;9 h1 P! |% Z$ S& d6 W- k# h5 Z. c+ `" Q
v=x.^j;1 a, y0 {/ s! o1 T# I4 E4 ?7 H8 \
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];
( K. S2 D; ]5 z$ d8 @- A( y7 DN=v*inv(A); %形函數(shù)
& }" G( l! P3 N1 t%-----------------------------------------------$ F( H8 h5 `3 S5 E
% 求單元剛度矩陣
* s* Q& D& u8 J/ h) lE=4.0e11;3 _% ?% A& q$ G. I9 b/ [, u8 c
I=5.2e-7; %I=bh^3/12=5.2e-7;
( c3 _1 n% I; ^0 |EI=4.0e11*5.2e-7;
! K3 A" q" G# J, \, k6 ^7 |1 wB=diff(N,2);
) i" ?2 X9 V( O0 t1 Z9 ?k=B'*B;
* E. t3 M0 X' t" [Kee=EI*int(k,x,0,'L');( @3 e4 G! t1 o+ X, c
Ke=sym(zeros(4,4,Ne)); %用三維矩陣存放單元剛度矩陣,每頁存一個,并初始化0
{3 ^3 n; z' @/ y! ]& B- iKe(:,:,1)=subs(Kee,'L',(10/Ne));
) @ i( h4 S8 @" E4 U; x, [$ C2 sT=eye(4); % 因為梁于x軸的夾角為0,所以坐標(biāo)變換矩陣為T=eye(4)
% J* _9 H J# Y* ^% `Ke(:,:,1)=T*Ke(:,:,1)*T';
: \& P! O s' f6 ?2 _4 afor ii=2:Ne) _2 [; V) C3 E. {; D8 A; A2 X/ W
Ke(:,:,ii)=Ke(:,:,1);8 U) v, Z# P5 T: j
end
/ `# B8 e& ?+ a& J) o( ?( SKe=double(Ke);* W2 V( }. b' r: M: b, h
%------------------------------------------------6 ^/ \+ Z8 `7 B `
% 由矩陣裝換法求整體剛度矩陣# g* \9 Y7 [, Q$ B2 T! W) T
% 自由度Ndof=2*NJ
+ c: P; G: n1 i( \# [ ^Ndof=2*NJ;4 X5 D/ U; r/ S" j# h7 E
K=zeros(Ndof);$ T- `2 {; X) @/ R4 Z
for ii=1:Ne
* ]9 H/ A o6 L- O G=[zeros(4,2*(ii-1)),eye(4,Ndof-2*(ii-1))];; v; u) v I1 h, w/ ]# `
KK=G'*Ke(:,:,Ne)*G;
# N$ f9 [9 {. Z4 n$ l K=K+KK;' T! ]9 N; C0 a/ ]8 C/ O; @
end
9 k! |. J) O- w6 k%------------------------------------------------
( s! s) ~' V' p0 d+ O8 E% 約束條件,對整體剛度矩陣進行修正,以便計算
, c$ U E) Y& iF=zeros(Ndof,1);! i( N5 U( {1 k3 j
F(Ndof-1)=-100;
~% z4 v, O+ D% t* z% u=F*inv(K),u=[v1,xta1,v2,xta2,v3,xta3,v4,xta4]',v1=0,xta1=0; w1 U+ h1 u# T0 t( p% t! p
K(1, =0; K(:,1)=0; %可以省略
( e# [ Z6 V L- w) \/ \K(2, =0; K(:,2)=0; %可以省略
9 l9 M/ e; z, V4 o; |0 c0 ]# ZKX=K(3:Ndof,3:Ndof);7 ^+ G8 X1 R2 h% d( P+ P
FX=F(3:Ndof,1);
' {3 V4 o; u |) }%------------------------------------------------# M; j; U* ^' q, L$ ~$ G
%求整體節(jié)點位移列陣
! ^, l, Z5 @& l5 e* g3 zuu=inv(KX)*FX;
5 ~& [3 w' B+ H7 d. yu=[0;0;uu];% K7 O E2 Q( y+ M7 c- b9 \' K
ii=1:2:2*NJ;+ y2 b$ O* w- B0 a' e! r
v=u(ii); %各節(jié)點撓度+ C, k5 b$ d; |" E
xta=u(ii+1); %各節(jié)點轉(zhuǎn)角 K9 Z+ ?) R6 p2 B+ J' u
%------------------------------------------------
6 F) b* }* e- t1 V4 t0 w2 {9 i% 后處理,計算單元應(yīng)變應(yīng)力. q3 P! K9 y% q
B=subs(B,x,'L');
3 Y- ]. I' `! |3 M( o6 U) hB=subs(B,'L',10/Ne);0 v$ t7 |: h0 j3 Q
Strain=zeros(1,Ne); %單元應(yīng)變,并初始化# z- `' E$ ^* r
Stress=zeros(1,Ne); %單元應(yīng)力,并初始化
9 J3 q @& @- p1 B" V A7 L; f* Ifor ii=1:Ne; m. F9 Y* Y5 r" }5 d0 P3 p; }
Strain(1,ii)=B*u(2*ii-1:2*ii+2,1);
$ N5 G0 q7 _! U( Y% J' b Stress(1,ii)=E*Strain(1,ii);
; m( Y6 @1 }+ f1 q# \% cend5 ^; z9 _- c, I, _- L/ ?& [
%--------------------------------------------------
! B" Y% F h5 K" V* \( g- U; O% 以下程序為屏幕輸出處理2 a& M+ T5 d# v$ v( R5 _& D3 u1 A. ?
disp(' ');' h8 w' l- C0 S! j2 e
disp(' Input:1-print Node displacement ');
- w& h3 `% f8 f* K. udisp(' Input:2-print strain ');
5 ~' U& P* K% \% E9 Q- ~disp(' Input:3-print stress ');
5 k/ E/ G# r2 l+ C8 l8 Q" edisp(' Input:q-quit or use ''ctrl+c'' Shortcut to end the program');# F9 v. e* F2 ?: k
9 z, K0 z2 n+ N# Uwhile 1
7 b6 E- D4 \7 C# V0 o* \4 u; g ii=input('Please input1 or 2 or 3:','s');1 v9 k! ], S9 N! d8 v
switch(ii)
/ l- v; f; f7 ]7 l case '1'. |/ I& ]" X7 |: T4 y
disp('Node displacement'); x! i. k& U* }5 G# Y* q
disp(u');0 q, M- v. L# c- Z4 B
case '2'& ]6 {" R' G& s$ ^4 a9 h
disp('element strain');
) I. \5 f( C5 ]& U( d disp(Strain);% T* s; T0 [- C; j( T
case '3'
1 y0 i! O. w/ ~5 e6 Z$ p. u: ] disp('element stress');' l y% P# Q6 U! O3 s+ H7 @( k7 R
disp(Stress);* y$ f" O. K2 S6 H5 u1 j, S2 B
case 'q' P" S$ _( z: w8 b% \3 p: L
disp('End of program');! @3 l- U% ]0 c! V, Y5 @
break;
$ m+ E. i& w$ S2 ~: x otherwise( y) \% O5 i7 [, D4 W+ h! E
disp('error!Please input again');9 M2 ^2 h& A# p" H6 e' v0 l
continue;& h( l/ K3 q7 U( m
end$ v3 Z1 r8 D8 ]
end) O7 g5 R& b4 e. m
( B4 Y2 \0 X8 b2 {' W' O
; J2 l; L4 t! a- _, e |
|