机械必威体育网址

找回密码
注册会员

QQ登录

只需一步,快速开始

搜索
查看: 3554 | 回复: 7
打印 上一主题 下一主题

分析一段matlab有限元程序

[复制链接]
跳转到指定楼层
1#
发表于 2013-3-24 13:41:49 | 只看该作者 回帖奖励 | 倒序浏览 | 阅读模式
这是关于梁单元的,可能大家觉得很简单。。。/ S! q5 k9 p8 M" ^9 Z* S8 |" G+ g J
今天翻电脑里的东西时发现的,是我大三时有限元时的作业,记得当时花了很多时间研究,学了不少东西,一个简单的作业可以涉及各方面的知识。8 z1 b9 B5 p- C4 Z# s0 x F. O3 x
毕业半年了,虽然记不清弹性矩阵、刚度矩阵等,但是只要一看书(好像现在的心境不下来)应该会,学校交的学习方法和理解问题的能力,有些东西不必记,除非每天从事,那该叫熟能生巧。! [. n4 R6 l, S: Y5 L5 _6 N' J' ]
记得当时编了两个程序,梁和三维实体的,我分享一下梁的程序吧,起码有亮点和创意,可以自己选择划分个数,在matlab方面花了不少功夫。) ?& `0 d) P; ^1 O( J6 z9 Z* M: @
非常感谢当时的有限元授课老师“韩清凯”,虽然老师已经到了大工。上课的教材为韩清凯 的《弹性力学及有限元法基础教程》,书上有个梁单元的例子,在82页。
' b/ ^2 ?- {) P% D! {--------------------------------------------------------------------------------
1 N- ]6 n- U& r梁程序如下,名字就不写了,用昵称吧“太平岛”,这是常用的网络昵称5 Z5 f" ~- @# |8 \6 X- q: @

6 m: m( z0 u5 T( q- t& y* w
' @8 A. \6 W. P, X%------------------------ BEAM2 ----------------
% V3 r6 J6 \, K$ z+ S& }disp('========================================');
2 u! R3 L9 Q' f5 P: y+ ?disp(' PROGRAM BEAM2 ');7 ^! ? \' `4 R) e5 L! ^ k4 A5 p/ c. P
disp(' Beam Bending Analysis ');' w; i/ {( O4 g1 W8 G
disp(' The programmer:太平岛 ');/ J8 s) ]2 @$ P6 e! B
disp('========================================');
" F4 \8 J6 A5 \4 a W: y [: pdisp(' '); %输出空行
. k/ L$ Z" V* m' q- \8 O7 Wwarning off all %关闭所有警告提示(求积分时,警告信息比较多)
4 C7 G/ y: f; V# \/ Y5 p* V( LNe=input('Please input the number of elements:'); %Ne为划分的单元数7 G+ b6 q! I, m& U
NJ=Ne+1; %NJ为节点数% `+ A0 z$ u6 g9 A* v1 v: O
x1=0;
! F. b* z0 ?/ C8 yx2=sym('L');- a& ?8 B- E; E1 ?. J8 Y" N
x=sym('x');
! _ [) u% S* Z+ H, a6 T. Sj=0:3;2 [0 }/ B! g# s- y* u' w
v=x.^j;
0 V& Y! T; V( 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];
6 }8 c; d& p6 N! @5 \7 tN=v*inv(A); %形函数
) C5 {8 M( q. `, t1 ?2 _& a%-----------------------------------------------
- C3 f5 N; _* v% 求单元刚度矩阵
3 z3 e8 @. I, s8 z/ AE=4.0e11;
( H7 j8 D# I) h$ pI=5.2e-7; %I=bh^3/12=5.2e-7;
* v9 [% @ Y3 Z! p. OEI=4.0e11*5.2e-7;
X: J' v1 x' P" N0 {. O7 K6 A! fB=diff(N,2);6 j1 { o v5 ?# Q2 [, X
k=B'*B;8 z9 U" c& C+ X- @- @. M1 d+ b/ V) @7 e
Kee=EI*int(k,x,0,'L');
0 Y* R( A) C+ b- _2 j7 bKe=sym(zeros(4,4,Ne)); %用三维矩阵存放单元刚度矩阵,每页存一个,并初始化0) L f7 { L5 C; I# _( A
Ke(:,:,1)=subs(Kee,'L',(10/Ne));
# `' C; Y+ t* j% l1 J- X( Q- dT=eye(4); % 因为梁于x轴的夹角为0,所以坐标变换矩阵为T=eye(4)
5 k, A" \7 Q7 j5 A- XKe(:,:,1)=T*Ke(:,:,1)*T';
F4 I1 K. |7 d5 Y" O8 v6 Y i+ mfor ii=2:Ne
; l2 A4 u/ I' `9 I, j" H, zKe(:,:,ii)=Ke(:,:,1);" O; |0 D) w$ k" ~% A
end1 y, K6 | ^0 _3 R
Ke=double(Ke);# @% |5 w% f. x7 t- f
%------------------------------------------------
8 d) t3 G' x8 k; |( W7 W% 由矩阵装换法求整体刚度矩阵' F5 g- P* C7 e7 r( d' p
% 自由度Ndof=2*NJ
4 z" P/ W1 C$ m+ @$ pNdof=2*NJ;
3 r7 e* e5 O% u2 }: RK=zeros(Ndof);0 ?/ I5 u% k# s8 Z
for ii=1:Ne
9 M7 E Y1 N: ?3 m, UG=[zeros(4,2*(ii-1)),eye(4,Ndof-2*(ii-1))];
4 ?- w8 Q! U1 n AKK=G'*Ke(:,:,Ne)*G;9 S& @" e! A/ R6 T5 ~- J2 j
K=K+KK;; Z, j6 Q" s/ P2 E1 E
end: l) C9 _! N8 g6 _! d; o
%------------------------------------------------
+ L" ]. g! }* g8 h @% 约束条件,对整体刚度矩阵进行修正,以便计算/ ~8 M( x3 {: k) l9 `" C
F=zeros(Ndof,1);
3 D, f. @( `7 E7 OF(Ndof-1)=-100;
" k6 h' ]! r6 I' |' m! k s x% u=F*inv(K),u=[v1,xta1,v2,xta2,v3,xta3,v4,xta4]',v1=0,xta1=0
+ Y0 ~) U. o6 m; T0 l- RK(1,=0; K(:,1)=0; %可以省略. s: p/ J3 K+ C* v
K(2,=0; K(:,2)=0; %可以省略
" V" f. [4 P" Q: C# e, X9 gKX=K(3:Ndof,3:Ndof);" Z; {- X0 C) k8 J$ v( @
FX=F(3:Ndof,1);# ?. m! P/ ~6 k/ F$ u* q: _
%------------------------------------------------" x# H+ W$ H0 ?8 g7 \
%求整体节点位移列阵
K, }: ]2 p( E! duu=inv(KX)*FX;
a/ l0 ^, d8 a1 }) R5 ?$ Xu=[0;0;uu];/ g& X! {* h* l
ii=1:2:2*NJ;
% V0 p% K6 F5 B* qv=u(ii); %各节点挠度- Z5 r) w+ E5 g- N
xta=u(ii+1); %各节点转角& c( g; y! e7 M
%------------------------------------------------
, ?6 b& ?9 ?$ `% 后处理,计算单元应变应力$ _# j3 z6 a1 w/ i
B=subs(B,x,'L');
& w" O0 M% l1 Q# xB=subs(B,'L',10/Ne);
/ }, @+ p7 @4 v& OStrain=zeros(1,Ne); %单元应变,并初始化$ i4 w+ X, u+ B+ c7 e
Stress=zeros(1,Ne); %单元应力,并初始化
' q8 z8 z7 [1 f8 Wfor ii=1:Ne
" x. [3 O8 n. F& TStrain(1,ii)=B*u(2*ii-1:2*ii+2,1);# T! G4 n/ \+ ^" ?+ z* g# I
Stress(1,ii)=E*Strain(1,ii);+ `& k1 z7 |" U$ L
end
3 e, z8 n$ D E2 R1 o" j%--------------------------------------------------
* U3 G' v! ?1 z/ I/ x% 以下程序为屏幕输出处理2 R/ l6 C6 M# r! q/ P# n8 U B
disp(' ');! c% w, b$ m8 m. d/ D; I
disp(' Input:1-print Node displacement ');
3 _7 O/ T$ A0 bdisp(' Input:2-print strain ');
! D+ L( ~' p. G! V2 F" Vdisp(' Input:3-print stress ');4 n6 J( O7 _+ r5 G
disp(' Input:q-quit or use ''ctrl+c'' Shortcut to end the program');1 m. U9 ?; R* b7 ]- f

. ?+ w0 W: t" g7 G. l+ e3 e5 L9 Swhile 1' v# ]+ y, ^9 Y7 R2 a6 ]& S
ii=input('Please input1 or 2 or 3:','s');
) ]5 T" h7 |: Bswitch(ii)# g& ^! c3 b- P: f9 |. a
case '1'
% U; h# d/ l* O: k n4 \disp('Node displacement');: f. e# v% x0 w6 p
disp(u');
0 J3 p+ Y; [) P2 qcase '2'
. ~" _" _1 G5 a" Xdisp('element strain');
2 s* g7 A$ z0 g9 K1 |' `disp(Strain);4 d+ }, {' ~3 J% H
case '3') i( e4 [1 J2 [% N' Y9 p: l
disp('element stress');# C) r# L8 W2 l( i. Z U
disp(Stress);3 L0 M: Y p6 J6 }
case 'q'
, X& N0 }: ?6 t+ j# vdisp('End of program');
3 F0 }7 n9 n3 abreak;
& S+ L% ]; R# y; o kotherwise: k, m9 w$ o3 r1 |5 N
disp('error!Please input again');) L0 S, [5 \# w6 G8 i9 C
continue;
+ V/ H2 I8 [2 `& G7 \$ b8 `end8 W2 M5 f5 G4 c# X4 m5 G% f; P
end
4 I4 i2 N, {; A1 V" ^8 j

# K, m$ b, T4 Z/ r( `
9 d/ |2 Z0 z4 ~! z$ i' {- G6 Z3 Q9 S
2#
楼主 | 发表于 2013-3-24 13:44:32 | 只看该作者
那个地方怎么变成笑脸了,冒号+右括号=笑脸,改了下,应该是下面,把英文括号改为中文,就好了吧/ d" A6 n0 E, n \
K(1,:)=0; K(:,1)=0; %可以省略
t/ I+ {7 ?# D$ ^" eK(2,:)=0; K(:,2)=0; %可以省略
3#
发表于 2013-3-24 14:55:52 | 只看该作者
没看懂
4#
发表于 2013-3-24 16:53:09 | 只看该作者
谢谢你啊,太感谢你了
5#
楼主 | 发表于 2013-3-24 18:35:59 | 只看该作者
jiaweicz 发表于 2013-3-24 16:53 ' v% p2 |$ `) `9 W; k, }3 Z
谢谢你啊,太感谢你了
/ x9 N& T) l+ P$ V% A+ [
这谢啥呢?
6#
发表于 2013-3-24 21:21:19 | 只看该作者
我以前编过一个5体展开的多体姿态动力学计算程序。。。可惜早就忘记了,sigh

点评

是啊,很多东西不用就忘记啦 发表于 2013-3-26 12:34
7#
发表于 2013-11-7 20:39:02 | 只看该作者
这程序也运行不了啊
8#
发表于 2013-11-7 20:44:34 | 只看该作者
j=0:3;
. L5 O9 D7 \; {. m& `) q* O# E1 x7 pv=x.^j; 是干啥的
您需要登录后才可以回帖 登录| 注册会员

本版积分规则

小黑屋|手机版|Archiver|机械必威体育网址(京ICP备10217105号-1,京ICP证050210号,浙公网安备33038202004372号)

GMT+8, 2024-7-29 02:18, Processed in 0.058868 second(s), 23 queries , Gzip On.

Powered byDiscuz!X3.4Licensed

? 2001-2017Comsenz Inc.

快速回复 返回顶部 返回列表