久久久国产一区二区_国产精品av电影_日韩精品中文字幕一区二区三区_精品一区二区三区免费毛片爱

 找回密碼
 注冊會員

QQ登錄

只需一步,快速開始

搜索
查看: 4220|回復: 7

分析一段matlab有限元程序

[復制鏈接]
1#
發表于 2013-3-24 13:41:49 | 只看該作者 |倒序瀏覽 |閱讀模式
這是關于梁單元的,可能大家覺得很簡單。。。2 Z$ c8 \% {) a# l; d
今天翻電腦里的東西時發現的,是我大三時有限元時的作業,記得當時花了很多時間研究,學了不少東西,一個簡單的作業可以涉及各方面的知識。
. o, ~9 g, G6 t  p9 B4 Q畢業半年了,雖然記不清彈性矩陣、剛度矩陣等,但是只要一看書(好像現在的心境不下來)應該會,學校交的學習方法和理解問題的能力,有些東西不必記,除非每天從事,那該叫熟能生巧。
4 |! Z7 d! e7 p8 R記得當時編了兩個程序,梁和三維實體的,我分享一下梁的程序吧,起碼有亮點和創意,可以自己選擇劃分個數,在matlab方面花了不少功夫。
' x" m: i4 m$ v6 R' x8 ?非常感謝當時的有限元授課老師“韓清凱”,雖然老師已經到了大工。上課的教材為韓清凱 的《彈性力學及有限元法基礎教程》,書上有個梁單元的例子,在82頁。7 u) W2 ~$ ?) I) B( D
--------------------------------------------------------------------------------( E8 w% ?* b4 X/ K/ ^
梁程序如下,名字就不寫了,用昵稱吧“太平島”,這是常用的網絡昵稱( \2 N* Y) x/ E$ L

: A7 L; f9 ~# u* [& w( ]" @5 P6 R4 H% O6 @8 \' {  r& E
%------------------------ BEAM2  ----------------6 x$ C# d. m8 t
disp('========================================');+ U0 x. h# f9 u/ H
disp('            PROGRAM BEAM2               ');
1 ~6 ~' _- A2 l* N+ _' a+ Zdisp('        Beam Bending Analysis           ');2 l2 b$ M  p3 n' Z0 n2 ]# Y
disp('        The programmer:太平島           ');, o3 m4 t$ s2 i8 y  d7 n
disp('========================================');
, v, g1 l6 `8 g7 V% u* Z1 N' fdisp(' ');                        %輸出空行
' A: V5 M( W& b! Zwarning off all                        %關閉所有警告提示(求積分時,警告信息比較多)
, b) s! f$ U5 s) p) e$ vNe=input('Please input the number of elements:');        %Ne為劃分的單元數
2 Z. f$ n3 D- X6 p. RNJ=Ne+1;                        %NJ為節點數
3 o8 y/ D) e- z/ n: _& \x1=0;$ f/ t% c4 p( f2 J* t! }, N
x2=sym('L');: k. K: K) C5 z9 D. U  X
x=sym('x');                        
/ H" P! |/ U5 A& \1 n& e, }j=0:3;
9 ]/ t4 y9 C* j/ v" S, p! m7 Nv=x.^j;
) Q) g- }; U; n2 M7 s3 c7 b: lA=[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];. @0 Z. t! f6 l9 Q) m" _) W* r
N=v*inv(A);                        %形函數; E) S8 D) T: X' s& ~
%-----------------------------------------------
) M+ a4 T5 d1 O7 Z+ u% 求單元剛度矩陣
' ~; {' g' b/ ]1 S2 L0 q! d, KE=4.0e11;& Z0 S2 C4 M% Z: {7 d# U
I=5.2e-7;                        %I=bh^3/12=5.2e-7;
" ~. U7 c) ^' V8 S$ qEI=4.0e11*5.2e-7;
- N# C2 u# S* H  y. U. H) X* o4 X7 }; rB=diff(N,2);0 K3 |" ~$ U! o
k=B'*B;+ F6 o* P% j9 Y3 |
Kee=EI*int(k,x,0,'L');
) y3 r9 v7 _  E- ~3 e8 C& }0 kKe=sym(zeros(4,4,Ne));        %用三維矩陣存放單元剛度矩陣,每頁存一個,并初始化0
4 f4 H& I' _$ l& |2 u; o& m2 fKe(:,:,1)=subs(Kee,'L',(10/Ne));
  B) w7 T: q6 ^+ d1 cT=eye(4);                % 因為梁于x軸的夾角為0,所以坐標變換矩陣為T=eye(4)+ H: M1 e! x6 G8 b/ ]
Ke(:,:,1)=T*Ke(:,:,1)*T';
/ l4 \7 x! `" D& B1 x, U/ B6 D7 tfor ii=2:Ne+ {2 M5 n8 k% t/ P1 r
    Ke(:,:,ii)=Ke(:,:,1);: I0 r0 F4 B) u9 a' E- H
end 3 |/ K, t5 D, U
Ke=double(Ke);& `$ Y. P0 c' }$ ^6 X' g9 \  c
%------------------------------------------------; e* X& E3 L* q. {
% 由矩陣裝換法求整體剛度矩陣
/ M% n0 g/ R8 \7 z% 自由度Ndof=2*NJ
& }3 ?' f) j9 F& |Ndof=2*NJ;
. D' V. X) l! K- |5 s: `* M8 hK=zeros(Ndof);$ L" h- V1 C: ?2 T) \4 P- q
for ii=1:Ne  u8 F3 z& ~" h5 d( U& W' e4 _
    G=[zeros(4,2*(ii-1)),eye(4,Ndof-2*(ii-1))];
  s3 e2 U2 Y8 p; S5 G    KK=G'*Ke(:,:,Ne)*G;5 V8 w$ i. R/ l, H& P' o1 {
    K=K+KK;) U6 C4 }7 ~5 d. y: l
end  
5 F2 ^7 D. D' p: V- [%------------------------------------------------
; H" V, F) k: U4 Y& P7 N2 W3 `& y$ g% 約束條件,對整體剛度矩陣進行修正,以便計算
7 U( C# D8 K  _& S4 H4 hF=zeros(Ndof,1);
7 c$ T; n2 ]( Z/ L( L& r6 h1 dF(Ndof-1)=-100;, {. }7 V0 |5 ?6 ]  \7 M: ?3 V, i
% u=F*inv(K),u=[v1,xta1,v2,xta2,v3,xta3,v4,xta4]',v1=0,xta1=0
6 K! f2 u" y* W6 p/ p" r' EK(1,=0;        K(:,1)=0;        %可以省略
9 i- k& T$ W$ g6 z) e5 uK(2,=0;        K(:,2)=0;        %可以省略
# g2 I/ {! C# d, C: l# P# a0 zKX=K(3:Ndof,3:Ndof);% w) L# Z9 r& k/ {5 C1 X& F: T
FX=F(3:Ndof,1);
, ^- A+ f3 t7 `* U4 d%------------------------------------------------
% L8 C3 F- r! ~( y%求整體節點位移列陣
5 a$ K9 O- ]' I' i: L( Auu=inv(KX)*FX;+ j" `3 F9 F" k  a( ?$ ~
u=[0;0;uu];
0 a; W) `, u3 y6 [6 @ii=1:2:2*NJ;
9 p0 s! B' n0 }4 X. U( ^4 w7 G) mv=u(ii);                        %各節點撓度
! d+ P, O$ U( `0 s4 V# `xta=u(ii+1);                        %各節點轉角
& N$ j$ Z) b2 ]2 d! g7 D9 D- L' m. L%------------------------------------------------
5 y2 ~, ]1 K/ M0 N- K% 后處理,計算單元應變應力  ]) _! B: g6 A: V- a" W& x/ w3 K/ ?% Q
B=subs(B,x,'L');4 A' b( v2 b/ q, N2 Y
B=subs(B,'L',10/Ne);
# B! D4 ^7 I) [Strain=zeros(1,Ne);                %單元應變,并初始化6 P0 ]- w  p: K4 F8 q, v# n' z$ S
Stress=zeros(1,Ne);                %單元應力,并初始化
4 W4 v+ h- z! ?- R1 Dfor ii=1:Ne; F  \# m* J! J! w
    Strain(1,ii)=B*u(2*ii-1:2*ii+2,1);
& R- k: b, l& {5 S    Stress(1,ii)=E*Strain(1,ii);
+ k! Q4 x) x: n0 E" Hend
6 Q" _& E! @3 k" H7 g) g8 Q%--------------------------------------------------
, U" c% E& ^' U- g0 o+ {% 以下程序為屏幕輸出處理
( J; J- g5 G1 xdisp(' ');: K) N' \1 e4 W: i
disp(' Input:1-print Node displacement ');
. G7 o3 h  M6 u% }' n% B8 odisp(' Input:2-print strain ');
$ k; a. d3 Q1 n5 rdisp(' Input:3-print stress ');& W8 A/ l5 P& A- \4 A
disp(' Input:q-quit or use ''ctrl+c'' Shortcut to end the program');) i8 c; s) w9 D9 c

: s  C4 V! ^  i) l; K0 e6 lwhile 1+ E6 S) W% U4 O* q/ n
    ii=input('Please input1 or 2 or 3:','s');* p. D7 v% E, A2 {+ [3 E- D
        switch(ii)* r1 c; _1 f/ s/ b4 H# `4 [
            case '1'4 n2 K, R  p  }4 o
                disp('Node displacement');3 c- [& t, ^! Q7 W: ]# u) u
                disp(u');* Y+ z4 @( t3 p
            case '2'
8 x3 ]  H% i: N5 `                disp('element strain');9 K& e" F( P  l  g
                disp(Strain);
% Z8 V2 E- J& b& C            case '3'
( N; L: P8 G: U! _3 t% t                disp('element stress');% a( D! @, o' d: Z
                disp(Stress);
  L  [" ]) Q; ]% i1 d  r& q            case 'q'
' R" v5 ?. S5 o" R                disp('End of program');- M1 z. j5 R! Z2 k' ?
                break;; f  D9 G/ v2 V7 m3 U. s
            otherwise/ F* ?" H0 F2 F& j+ T
                disp('error!Please input again');
' \8 {1 D! i  v8 p# K$ S                continue;6 M0 F& h+ `( }
        end
' d+ u4 `5 F; U1 h7 ]1 eend
! J- x9 |& @7 f. \6 }

6 q7 _7 J" b$ d
! q8 D/ j+ F& E" j
回復

使用道具 舉報

2#
 樓主| 發表于 2013-3-24 13:44:32 | 只看該作者
那個地方怎么變成笑臉了,冒號+右括號=笑臉,改了下,應該是下面,把英文括號改為中文,就好了吧
1 b" H, U5 ?% U  l, @6 HK(1,:)=0;        K(:,1)=0;        %可以省略
% m9 ~" l8 i/ b2 z3 Q5 L$ OK(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
/ R# g% {, {! p+ D: c7 D# p謝謝你啊,太感謝你了

+ P$ z8 p6 k' z8 Q: M& _( q這謝啥呢?
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;
6 P/ V, _, `/ E* S8 @v=x.^j;  是干啥的
您需要登錄后才可以回帖 登錄 | 注冊會員

本版積分規則

Archiver|手機版|小黑屋|機械社區 ( 京ICP備10217105號-1,京ICP證050210號,浙公網安備33038202004372號 )

GMT+8, 2025-7-25 15:15 , Processed in 0.071478 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.5 Licensed

© 2001-2025 Discuz! Team.

快速回復 返回頂部 返回列表