|
这是关于梁单元的,可能大家觉得很简单。。。' C/ ^0 T/ K/ P" ?
今天翻电脑里的东西时发现的,是我大三时有限元时的作业,记得当时花了很多时间研究,学了不少东西,一个简单的作业可以涉及各方面的知识。5 ^: k' { w3 n( w0 v6 G9 f
毕业半年了,虽然记不清弹性矩阵、刚度矩阵等,但是只要一看书(好像现在的心境不下来)应该会,学校交的学习方法和理解问题的能力,有些东西不必记,除非每天从事,那该叫熟能生巧。2 ^ g0 B: U. q+ f
记得当时编了两个程序,梁和三维实体的,我分享一下梁的程序吧,起码有亮点和创意,可以自己选择划分个数,在matlab方面花了不少功夫。! b) O+ f6 w/ O/ {+ w. Y
非常感谢当时的有限元授课老师“韩清凯”,虽然老师已经到了大工。上课的教材为韩清凯 的《弹性力学及有限元法基础教程》,书上有个梁单元的例子,在82页。
) _: B$ p+ b1 F) y/ Z& X--------------------------------------------------------------------------------- x/ T# s2 p k" c' L3 y
梁程序如下,名字就不写了,用昵称吧“太平岛”,这是常用的网络昵称
@0 m) H2 P+ B M' s9 U# n5 c$ Y$ G- j4 i& ?( Q5 N
) s# q2 M i1 ]
%------------------------ BEAM2 ----------------
" x0 n6 K' W$ n$ b' d4 @. F# |disp('========================================');' u. T7 L% I2 U& `2 f$ @
disp(' PROGRAM BEAM2 '); ^: V3 _# _! [" T$ l
disp(' Beam Bending Analysis ');
; }( Q5 O+ J4 c7 g: ldisp(' The programmer:太平岛 ');
6 G$ p( [2 W1 u, y1 S. Ldisp('========================================');
) @* n1 n/ V! f' Y) ldisp(' '); %输出空行6 m- P. v* s+ P9 [' w5 F2 q; }
warning off all %关闭所有警告提示(求积分时,警告信息比较多)
2 c) c% [. [% m& N, hNe=input('Please input the number of elements:'); %Ne为划分的单元数7 a* G( D8 L e( g) y
NJ=Ne+1; %NJ为节点数' A8 f" v' J) s; t* G7 T/ ?( h
x1=0;
" A$ Q. ~- `- b7 ux2=sym('L');
1 \7 q' P! X, t& y6 A" rx=sym('x'); 3 v% }* I# i' L- Z: b' P) g& ~
j=0:3;/ ~- z8 g# o6 `7 j1 H
v=x.^j; w6 h. L! d' b( w. g( g2 L& A
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];
7 h# Z# Y) W uN=v*inv(A); %形函数5 ?0 ~4 G6 |$ G$ k# k f
%-----------------------------------------------
& _, i- f. p0 \% 求单元刚度矩阵$ h6 ~' J4 I+ E- l" O) e
E=4.0e11;
8 p8 {: R+ B6 s4 T& b: UI=5.2e-7; %I=bh^3/12=5.2e-7;4 H/ q' z! f( t1 x8 X7 g3 K" P
EI=4.0e11*5.2e-7;1 Z7 C6 C/ `( E. E" m
B=diff(N,2);1 N) F9 |1 D2 S$ t: h& b
k=B'*B;
/ {9 K4 }4 k: {. M5 ZKee=EI*int(k,x,0,'L');( C( [, k# G. O
Ke=sym(zeros(4,4,Ne)); %用三维矩阵存放单元刚度矩阵,每页存一个,并初始化00 y8 t+ G! H" j
Ke(:,:,1)=subs(Kee,'L',(10/Ne));
# ^4 R; j5 n! jT=eye(4); % 因为梁于x轴的夹角为0,所以坐标变换矩阵为T=eye(4)5 x9 N; |8 g$ L8 t. j# F @
Ke(:,:,1)=T*Ke(:,:,1)*T';! Y) G! f4 j. `8 R6 M
for ii=2:Ne
6 @/ w# t" P6 e( W M2 ~/ J1 c I Ke(:,:,ii)=Ke(:,:,1);
( V3 [4 E# f; D1 p& rend $ j) ^4 t/ T2 w0 W7 E
Ke=double(Ke);8 L+ c' {! c* G& s2 Z
%------------------------------------------------! l- U b9 g2 j; B& w+ Z
% 由矩阵装换法求整体刚度矩阵
6 [! X) [ b3 I% 自由度Ndof=2*NJ
( U8 h: B9 q9 c' \7 F4 VNdof=2*NJ;$ z3 q, V: ~) q @) c
K=zeros(Ndof);
+ k# K; n) z b% [5 pfor ii=1:Ne
4 ~: J9 [3 V- Q6 T G=[zeros(4,2*(ii-1)),eye(4,Ndof-2*(ii-1))];# f; F1 c, b5 c
KK=G'*Ke(:,:,Ne)*G;
6 U& Y/ T+ H9 L4 Y) h `6 o% _ K=K+KK;
4 S* p9 r. a6 p! f2 }end $ ]3 ?2 O1 l3 t. G* p: H1 P
%------------------------------------------------
( u+ E0 }7 ?/ A4 N, N% 约束条件,对整体刚度矩阵进行修正,以便计算
& C: b& d. }7 |2 ~. u) n, oF=zeros(Ndof,1);
* V' V% q6 z6 w& b& Z. rF(Ndof-1)=-100;$ n3 t9 G8 d# L( P9 F* f W/ U
% u=F*inv(K),u=[v1,xta1,v2,xta2,v3,xta3,v4,xta4]',v1=0,xta1=05 m9 E. c% G; C) a/ }; m
K(1,=0; K(:,1)=0; %可以省略* s7 H5 n b5 y$ i
K(2,=0; K(:,2)=0; %可以省略
x7 e f2 ^2 Q; @( A1 C9 w- mKX=K(3:Ndof,3:Ndof);
% E* L& `( F7 ]- u2 r4 bFX=F(3:Ndof,1);3 n0 H: U. D+ E, r( L
%------------------------------------------------) ^+ v/ l; `" G/ D8 B7 d4 c
%求整体节点位移列阵
" F. }2 S8 t! I. muu=inv(KX)*FX;
4 m2 S- b% Z- o. eu=[0;0;uu];+ m" l5 x& j/ F- t2 d3 t4 O% b
ii=1:2:2*NJ;; n) T5 X4 r+ {8 s9 t1 V/ B
v=u(ii); %各节点挠度" W& g& X: u( E4 R9 Z
xta=u(ii+1); %各节点转角; k- R$ ]& v- ?6 m4 ^# _7 U; H
%------------------------------------------------
* b6 ^9 H1 r) H: I9 S% 后处理,计算单元应变应力; i2 u3 u4 C7 M
B=subs(B,x,'L');
7 w4 Y0 g6 A6 Z- O. M+ XB=subs(B,'L',10/Ne);
, ^; e9 B* F/ a* IStrain=zeros(1,Ne); %单元应变,并初始化+ p. b9 q( r' q6 x$ p
Stress=zeros(1,Ne); %单元应力,并初始化
. k( ]6 T+ C0 L! s/ \) Lfor ii=1:Ne
6 i# h+ |$ s3 Y) A4 ?% l6 z Strain(1,ii)=B*u(2*ii-1:2*ii+2,1);
7 q1 z; F9 [0 W# J; ^7 v8 t( I Stress(1,ii)=E*Strain(1,ii);
- y( D @1 M. t8 Z$ b' ^end
- D$ R* G& X+ Y% J* c%--------------------------------------------------
9 }" d2 [1 g* Y' f0 d% 以下程序为屏幕输出处理; C% U% ]% e" F" s0 S0 F. \
disp(' ');
- b, b9 b5 z8 K9 _; i2 W _# ydisp(' Input:1-print Node displacement ');
& K$ h/ p8 H/ L+ M% v, w& q0 ddisp(' Input:2-print strain ');, w Y, P& i1 \9 I& p
disp(' Input:3-print stress ');) y7 y( B7 J0 Z t- s/ o
disp(' Input:q-quit or use ''ctrl+c'' Shortcut to end the program');
$ c( q) G3 \2 z2 Y- E: x. c3 F5 y
8 p$ d) x, u) |% Owhile 12 T% X* N% f# U) K
ii=input('Please input1 or 2 or 3:','s');( D) h4 t; _' Y4 {
switch(ii)6 `+ `# ~, [. @/ j8 |2 m
case '1'
9 j" n8 l$ J- a disp('Node displacement');6 s+ a# L7 P- h# o7 B
disp(u');
9 ?$ j$ @( K M1 x& J" A case '2'
9 [! K2 S( J% q4 F& N& C1 Y disp('element strain');
- V$ S2 N/ I( [/ T0 t1 L$ p% Q" B0 U disp(Strain);0 c$ l0 H) M8 Y$ j3 G8 _
case '3'( R3 [( C+ K& e, N6 C
disp('element stress');) J5 N- l9 I& X* ~8 n
disp(Stress);5 o+ t! M4 }! B3 B# Q/ m; X: k
case 'q'
# M6 @9 c8 s& v9 e0 A% ~ disp('End of program');
0 k6 u8 ~; i9 U2 e2 `; y* E break;, m5 K3 ]9 F3 x( d: u
otherwise2 Y* m0 U, \1 X5 g) s% ?
disp('error!Please input again');+ ]% |, H Q2 P9 x* o; J7 r$ a
continue;' v+ K! u2 w/ |" L; z3 a2 ?
end
+ E5 c ~0 ?$ \2 Z7 _: G. J0 }% ~end C2 I: i4 O/ e. s5 K W. y) K
" m0 `* `) Y; X, S: X4 O; e' ?" s+ b( }! Q' a0 t& r
|
|