|
这是关于梁单元的,可能大家觉得很简单。。。+ B g8 [- h# L
今天翻电脑里的东西时发现的,是我大三时有限元时的作业,记得当时花了很多时间研究,学了不少东西,一个简单的作业可以涉及各方面的知识。6 X0 D/ R) s( ]; s1 X$ ~9 Q
毕业半年了,虽然记不清弹性矩阵、刚度矩阵等,但是只要一看书(好像现在的心境不下来)应该会,学校交的学习方法和理解问题的能力,有些东西不必记,除非每天从事,那该叫熟能生巧。
! [' \, D* f7 F% Y记得当时编了两个程序,梁和三维实体的,我分享一下梁的程序吧,起码有亮点和创意,可以自己选择划分个数,在matlab方面花了不少功夫。7 C/ @: W, \: [: e6 Y& A" P* _& D
非常感谢当时的有限元授课老师“韩清凯”,虽然老师已经到了大工。上课的教材为韩清凯 的《弹性力学及有限元法基础教程》,书上有个梁单元的例子,在82页。
1 Q/ Q, f& s: U5 ^& E- A% v--------------------------------------------------------------------------------
" @* w; w2 A3 d% V; S梁程序如下,名字就不写了,用昵称吧“太平岛”,这是常用的网络昵称( ^2 G: q* [/ z
$ B( E* ]5 [7 g8 I7 Q
8 ~$ w$ ~# @ K$ j* K
%------------------------ BEAM2 ----------------( } A) Q4 |0 @+ M& {& a+ V
disp('========================================');5 t( B7 c0 A( N) e
disp(' PROGRAM BEAM2 ');
( j' x+ _" `2 K9 Ddisp(' Beam Bending Analysis ');
: T( W% D& i. l( C- tdisp(' The programmer:太平岛 ');
1 m* N& u4 t$ k0 f. |disp('========================================');# q0 ^2 R: g: [- m
disp(' '); %输出空行
! I2 x1 e+ I' p% p& `warning off all %关闭所有警告提示(求积分时,警告信息比较多)5 `! w; {6 i5 w9 C+ q
Ne=input('Please input the number of elements:'); %Ne为划分的单元数, `( {) Y( ]: n* v5 b
NJ=Ne+1; %NJ为节点数: s# j B$ x2 y/ i1 ]
x1=0;& t o, l0 C: h3 p1 H! X* [8 A
x2=sym('L');" d* Y0 R: H) P6 r) W8 ^& Y
x=sym('x');
3 |! a0 K- X* y9 {7 g9 ?j=0:3;
% f' s# a: G/ P3 u3 X3 P/ h- Jv=x.^j;3 ]# b! o, s& J5 L+ q* B
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];) a$ q$ t; P& _
N=v*inv(A); %形函数
& q& y2 H* i! v/ ]6 |& ^1 v%-----------------------------------------------
' r, l# H9 O" Q+ B: t0 ^% 求单元刚度矩阵) s( j. C+ N& r
E=4.0e11;$ D% |" R M' n* \1 m9 y
I=5.2e-7; %I=bh^3/12=5.2e-7;. t. V h9 x$ F+ G9 O
EI=4.0e11*5.2e-7;
$ F2 @" l. |& u; B: _- t: sB=diff(N,2);
: E8 v( J( @5 Bk=B'*B;. I. t) R5 e+ b- V) [
Kee=EI*int(k,x,0,'L');& G* Y d8 Z8 _9 j) c% n5 W
Ke=sym(zeros(4,4,Ne)); %用三维矩阵存放单元刚度矩阵,每页存一个,并初始化0
1 o8 D' n. |& N/ z: t' B+ Y7 KKe(:,:,1)=subs(Kee,'L',(10/Ne));8 q" }' Y' X- u- {) o
T=eye(4); % 因为梁于x轴的夹角为0,所以坐标变换矩阵为T=eye(4)
6 k+ T$ Z1 N2 j4 H5 e* p% A( b. wKe(:,:,1)=T*Ke(:,:,1)*T';5 k2 v% j* i# t6 r& O+ f9 W
for ii=2:Ne6 ^# t6 Q: M1 Y! J. r3 L- U- x
Ke(:,:,ii)=Ke(:,:,1);4 ]: q3 Q/ O' t9 A- o- `
end " M- e4 v$ y- W8 ]
Ke=double(Ke);2 k& x0 h# ?8 K9 ?9 j a
%------------------------------------------------) f5 o. Y+ S5 b! u" B) |! L+ v
% 由矩阵装换法求整体刚度矩阵- r6 e9 g8 D8 x- G0 h% s, O& E
% 自由度Ndof=2*NJ' T4 u8 E- `" b- M& q4 T: N0 k) I
Ndof=2*NJ;
/ t- O) ] ^5 x X+ S/ X1 j) LK=zeros(Ndof);) Q6 J1 n6 N# W
for ii=1:Ne
: B6 l$ Q8 d0 t: n T' P+ Y G=[zeros(4,2*(ii-1)),eye(4,Ndof-2*(ii-1))];' T: a) n4 ?2 ]0 z) Z# y4 _$ z
KK=G'*Ke(:,:,Ne)*G;
4 K6 A; ~/ I7 H K=K+KK;
# I- X) z0 r6 V6 M1 f$ J0 L& Uend
8 U* k0 e t+ a! y%------------------------------------------------
1 [$ l6 X; o2 t8 S3 A% 约束条件,对整体刚度矩阵进行修正,以便计算
) n: C) s! ]. _F=zeros(Ndof,1);2 C& L0 w1 E8 _0 M8 p* n
F(Ndof-1)=-100;1 U' R$ j- l8 I4 R5 Z/ K* R2 [
% u=F*inv(K),u=[v1,xta1,v2,xta2,v3,xta3,v4,xta4]',v1=0,xta1=0
( ~+ M+ c! x1 R0 ~/ k0 K1 FK(1, =0; K(:,1)=0; %可以省略% F% ~9 c2 U( v7 c! y
K(2, =0; K(:,2)=0; %可以省略) Y8 |6 f! ^* t6 U, y4 P* d7 g$ p4 ~
KX=K(3:Ndof,3:Ndof);# y" Y% f8 r; a- d7 p$ ~" V
FX=F(3:Ndof,1);
- I0 b3 b7 d; n$ W8 _4 x%------------------------------------------------$ r# [' u5 t) L g* i
%求整体节点位移列阵
3 p0 f1 _+ t/ f9 c1 p; Juu=inv(KX)*FX;' D, F* \# I; H1 F
u=[0;0;uu];' A, G/ U$ k8 B( z4 Y
ii=1:2:2*NJ;
}! r+ s, j/ X9 z* [0 v) Jv=u(ii); %各节点挠度& X7 w* l8 `, S. f2 |6 T( S
xta=u(ii+1); %各节点转角
* x6 _5 M: Z# W9 t0 K! ~. J: [3 S( U%------------------------------------------------8 U) E* c( ~! a) f9 d* V( `
% 后处理,计算单元应变应力* t. _3 Q v3 l! c
B=subs(B,x,'L');
6 Y0 c( b3 W/ X) | ^B=subs(B,'L',10/Ne);; `" u6 ^3 F4 X2 U+ O* a9 ~ G
Strain=zeros(1,Ne); %单元应变,并初始化
9 u) V1 H% v z$ vStress=zeros(1,Ne); %单元应力,并初始化7 \! }+ {8 v" A: P( j. M5 L% {/ j8 d
for ii=1:Ne
- V- M$ V' M+ R/ i2 m: Y: S/ G Strain(1,ii)=B*u(2*ii-1:2*ii+2,1);: s' _' p: T1 A6 D( e5 U) q- a
Stress(1,ii)=E*Strain(1,ii);. }' t$ p6 ^7 S* A
end; t& @4 I+ ]% Z. I7 ]' Z
%--------------------------------------------------
; O2 E3 y7 }; K/ c5 Y% 以下程序为屏幕输出处理
2 [" }0 V/ M' q3 ]; T8 Wdisp(' ');1 T% N7 ^- Q! x8 e
disp(' Input:1-print Node displacement ');- K" ?6 N2 P/ u9 J3 p
disp(' Input:2-print strain ');
. i3 n6 {2 d" j. |7 I% K% Bdisp(' Input:3-print stress ');
6 r$ H% b9 l4 M; V* F* a9 Q Wdisp(' Input:q-quit or use ''ctrl+c'' Shortcut to end the program');. x" v+ z$ N; U
3 W" q3 O9 ], a: G( e. bwhile 1
9 U# V8 r6 U1 M7 h ii=input('Please input1 or 2 or 3:','s');. x, r5 c: d3 S9 ?
switch(ii)
7 f3 I, J# I1 Q7 J& P$ e' x case '1'
: S+ @+ i m, ~: a! D/ V! H$ d disp('Node displacement');
! |6 i$ T6 C+ C$ C" d disp(u');
7 u7 q# y& ?* |& ~# U2 J case '2'" |9 w/ H& I+ }. j
disp('element strain');( @' a- |1 k$ _$ e7 h5 }0 M# N
disp(Strain);
y) A: ^5 i4 l' O7 c+ A8 R7 r case '3'
1 A2 w2 t4 e# V6 _ disp('element stress');
3 T" A5 [+ m, r8 z6 r disp(Stress);
W% y7 S3 b7 v/ _7 {5 @ case 'q'0 [* d: E$ K5 M& n; Q
disp('End of program');
! F- d3 F, z$ y) X4 p( V5 o! E0 k break;
; ?1 |3 Y) p' x2 P! n otherwise
( ]4 j, |3 s/ i) P: U( e- @; Z disp('error!Please input again');
; f. N. P/ }9 h( w" }& K0 d2 c) _# T% a continue;
0 r; O7 m }9 r' \/ o end" ~2 S4 U1 l8 j1 d- M' l
end) h$ e& `3 y8 z& {; b
. F2 S' [/ Y0 i& H: `2 z2 {& M
5 `3 I% @2 `6 a3 @5 @. m |
|