机械必威体育网址

 找回密码
 注册会员

QQ登录

只需一步,快速开始

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

分析一段matlab有限元程序

[复制链接]
跳转到指定楼层
1#
发表于 2013-3-24 13:41:49 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
这是关于梁单元的,可能大家觉得很简单。。。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
回复

使用道具 举报

2#
 楼主| 发表于 2013-3-24 13:44:32 | 只看该作者
那个地方怎么变成笑脸了,冒号+右括号=笑脸,改了下,应该是下面,把英文括号改为中文,就好了吧0 i# U. M9 }0 u% ~
K(1,:)=0;        K(:,1)=0;        %可以省略8 ]# @0 L9 |: \; N* d% ?: Z
K(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
3 E6 g' A% {( L* G: Y) `谢谢你啊,太感谢你了
- S* [+ G* Q. S6 A6 B0 J2 V
这谢啥呢?
回复 支持 反对

使用道具 举报

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;
. F9 a3 o# {7 N; V' J9 D+ b% C/ wv=x.^j;  是干啥的
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 注册会员

本版积分规则

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

GMT+8, 2024-11-7 19:17 , Processed in 0.058897 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4 Licensed

© 2001-2017 Comsenz Inc.

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