|
这是关于梁单元的,可能大家觉得很简单。。。9 _7 H* {; ^! K/ Q$ A B7 y2 t9 Q
今天翻电脑里的东西时发现的,是我大三时有限元时的作业,记得当时花了很多时间研究,学了不少东西,一个简单的作业可以涉及各方面的知识。* |! {" J5 f; g. }
毕业半年了,虽然记不清弹性矩阵、刚度矩阵等,但是只要一看书(好像现在的心境不下来)应该会,学校交的学习方法和理解问题的能力,有些东西不必记,除非每天从事,那该叫熟能生巧。: l2 o! r% k+ d, g, J
记得当时编了两个程序,梁和三维实体的,我分享一下梁的程序吧,起码有亮点和创意,可以自己选择划分个数,在matlab方面花了不少功夫。
! t- d- P3 n! p9 H8 A$ O! R非常感谢当时的有限元授课老师“韩清凯”,虽然老师已经到了大工。上课的教材为韩清凯 的《弹性力学及有限元法基础教程》,书上有个梁单元的例子,在82页。+ U& Z5 ]. V/ k
--------------------------------------------------------------------------------
1 N7 L. P; P$ X8 C5 _' e% _5 L梁程序如下,名字就不写了,用昵称吧“太平岛”,这是常用的网络昵称
6 B# Y0 w- A$ f2 A' V5 M9 s; N" S( \3 \- |0 Y
& s) ?7 K4 x# r! {4 a%------------------------ BEAM2 ----------------$ ^& w' g; f$ o( g; J4 P+ |1 E. ]: W
disp('========================================');2 h. f9 \. u! q o! m6 ?
disp(' PROGRAM BEAM2 ');
; P' o7 }9 p" d8 c, G+ h! xdisp(' Beam Bending Analysis ');
5 R+ v" t& P7 Y0 ^3 `2 k- u& odisp(' The programmer:太平岛 ');
3 V* t% R" j7 b% p5 U- V% c: e: [disp('========================================');# R5 R3 r4 `9 I) \
disp(' '); %输出空行+ [7 Z2 P+ a6 h2 c
warning off all %关闭所有警告提示(求积分时,警告信息比较多)2 G* o7 U( B9 L2 K! k f6 |
Ne=input('Please input the number of elements:'); %Ne为划分的单元数
0 g& U6 }( B+ TNJ=Ne+1; %NJ为节点数
) L0 \4 ^5 E, E) T8 g+ N( m7 G2 |* }: f# yx1=0;0 v, L$ L: H- n$ l: c- Q
x2=sym('L');
1 p, K$ K1 d$ b. vx=sym('x'); # d. ?; _; h8 `1 D7 z/ _; L# {
j=0:3;
$ w% [1 D8 r/ V# @. h" J: G/ @2 jv=x.^j;2 ~$ r2 l) H" J5 c
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];
, Z( E8 O7 E/ BN=v*inv(A); %形函数
2 { N' q4 Y2 b. U2 ^# [. T%-----------------------------------------------2 l% J, y5 G" ^2 T0 s! [
% 求单元刚度矩阵
. V/ M3 Y8 k' R% {/ K2 KE=4.0e11;
$ H8 ^5 E/ M* q; A1 x/ Z( R' T! mI=5.2e-7; %I=bh^3/12=5.2e-7;
T3 k$ O1 x! C) j8 f6 fEI=4.0e11*5.2e-7;
* b1 @! q: R' g2 e% [" }& t- rB=diff(N,2);; s+ H5 O. c( M0 J6 j( D' S, H
k=B'*B;# O) F+ K# j5 k% I1 x" ^2 h
Kee=EI*int(k,x,0,'L');* w% c$ M+ H7 K, j" R) j$ f
Ke=sym(zeros(4,4,Ne)); %用三维矩阵存放单元刚度矩阵,每页存一个,并初始化0+ H4 P6 v% V6 {1 f. W+ C, N- O' h
Ke(:,:,1)=subs(Kee,'L',(10/Ne));
7 G, N; S Z- v8 s1 GT=eye(4); % 因为梁于x轴的夹角为0,所以坐标变换矩阵为T=eye(4)
& @. U3 H; Y/ Z; IKe(:,:,1)=T*Ke(:,:,1)*T';2 ^4 C, n" j) Y2 u
for ii=2:Ne
. T% O9 |/ ?+ x Ke(:,:,ii)=Ke(:,:,1);3 H2 u$ U+ R; d$ L" \. Z% u
end
) r" J( [# m' B% i( SKe=double(Ke);
* o U5 f2 c* R% q5 x%------------------------------------------------
" U9 s: O' ^/ x! w5 Q% 由矩阵装换法求整体刚度矩阵
" y$ [5 B- K' y7 ^) N3 q2 ~% 自由度Ndof=2*NJ
2 `( S4 q" |: w) l/ |Ndof=2*NJ;6 @# ~7 K7 z4 O8 F0 |1 y$ ?
K=zeros(Ndof);
* j0 V4 ~% Z) Y% X+ K7 A: x( `2 tfor ii=1:Ne* Q* w$ u. }' h
G=[zeros(4,2*(ii-1)),eye(4,Ndof-2*(ii-1))];
& F" B3 O0 g# P& M9 c8 ~! j KK=G'*Ke(:,:,Ne)*G;, a" w$ A* V" H9 I
K=K+KK;; J% w% [/ B2 j% _/ A* E( G- V
end 3 T6 P4 `' W9 j3 Y
%------------------------------------------------+ i; p& U+ t$ L+ a1 [/ e
% 约束条件,对整体刚度矩阵进行修正,以便计算
) |9 E5 `6 r! f2 n) ]3 s$ J2 FF=zeros(Ndof,1);( W2 G$ S1 ^$ W. [! P* w! L. u
F(Ndof-1)=-100;
% h4 v6 L) n7 R# M% X% u=F*inv(K),u=[v1,xta1,v2,xta2,v3,xta3,v4,xta4]',v1=0,xta1=0
; i6 L: P8 I h7 ?K(1,=0; K(:,1)=0; %可以省略0 J1 z: r, |. J
K(2,=0; K(:,2)=0; %可以省略- i1 G2 g* ^2 w2 b
KX=K(3:Ndof,3:Ndof);; A9 W" u1 j+ f$ |1 D& z; X7 b
FX=F(3:Ndof,1);
# B2 G9 @1 c3 q. v; J%------------------------------------------------
- Y0 o5 j" U8 U8 P%求整体节点位移列阵. ^4 y& ` }' i) S$ k; t' }/ ?
uu=inv(KX)*FX;2 P6 c, ?' ]+ e8 W7 _0 T) G! F$ J Q
u=[0;0;uu];
: t0 M1 D8 K( r2 W" Uii=1:2:2*NJ;( J, q/ `' o7 I( Y
v=u(ii); %各节点挠度7 i5 ~, n& r+ L" `# C4 U3 {
xta=u(ii+1); %各节点转角
& u: w7 R7 d& A- u%------------------------------------------------* g s7 @8 T7 }) x
% 后处理,计算单元应变应力3 a/ Z5 v7 U0 u+ D9 t
B=subs(B,x,'L');* M$ q5 C; v* s. P7 e, I
B=subs(B,'L',10/Ne);
7 H+ D3 T5 L2 A: S8 |, DStrain=zeros(1,Ne); %单元应变,并初始化( Q9 o& G; c/ F* B6 E( c
Stress=zeros(1,Ne); %单元应力,并初始化
* v- ^" K: e( j6 |/ x4 d, ?for ii=1:Ne
, s) {" g# A+ |! j6 b) V Strain(1,ii)=B*u(2*ii-1:2*ii+2,1);
# f2 l; V( X* x" _7 K Stress(1,ii)=E*Strain(1,ii);
; d6 j7 ]/ M. K( f* F" ]end
/ F A6 t% {, A2 M& P) o! ]: V- F B( P%--------------------------------------------------4 E- M3 z' Z. Y4 a* y
% 以下程序为屏幕输出处理
2 x; c' ^% m+ b. N2 Q( U! p; |disp(' ');
6 D' o) Z. w$ f+ C" g/ D4 L) cdisp(' Input:1-print Node displacement ');7 M0 j! S% d$ s7 @
disp(' Input:2-print strain '); w) j: C: j( m2 ?
disp(' Input:3-print stress ');
8 I7 ~4 ?$ G/ g+ S& I# tdisp(' Input:q-quit or use ''ctrl+c'' Shortcut to end the program');- Y6 l7 J/ z* q3 z; G& F4 R
* V8 l; _/ `8 z* L2 F5 X. Iwhile 1' ]4 v5 Z. H4 Y2 W( `' T& o: v
ii=input('Please input1 or 2 or 3:','s');. J4 p' p3 s2 M) z. P. E
switch(ii)
+ Q6 Y5 w1 B* F7 O4 g case '1'' \# b0 y, O8 C$ L: S
disp('Node displacement');$ j- Y" P% q; D
disp(u');
3 Y6 H& t9 w ?: r$ n case '2'7 u' ^0 [* F. f$ j$ m* ?1 L
disp('element strain');+ Q; V. g0 B4 s) c; `5 j
disp(Strain);
1 E( F; o4 f& [& v: w# v5 A6 Y1 w2 T% h case '3'
: T% P3 o: x9 n5 @6 _" D% d disp('element stress');
8 H! D( x! P8 `/ P V8 n disp(Stress);0 u# I+ y" S; I
case 'q' e3 Q& e) A) `) \3 u9 z4 _
disp('End of program');
( a, d( v; l/ J( K5 r* _ break;
* F! r9 i5 D- @. R9 i otherwise" T# ~( E, x9 c/ b$ F
disp('error!Please input again');
6 H- B, ~+ C3 S# n# t continue;
9 P/ N6 N# R2 \8 r9 D" \% F0 x end
! s8 ^# O: J* @7 |, e$ iend+ P$ h: f; }( j% P& h
( @% C4 i H. p: K. b* \( I E/ ]1 C' _- @. K" }3 s3 c
|
|