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

 找回密碼
 注冊會(huì)員

QQ登錄

只需一步,快速開始

搜索
查看: 4219|回復(fù): 7

分析一段matlab有限元程序

[復(fù)制鏈接]
1#
發(fā)表于 2013-3-24 13:41:49 | 只看該作者 |倒序?yàn)g覽 |閱讀模式
這是關(guān)于梁單元的,可能大家覺得很簡單。。。) P! K2 m" F+ [: T* K" B
今天翻電腦里的東西時(shí)發(fā)現(xiàn)的,是我大三時(shí)有限元時(shí)的作業(yè),記得當(dāng)時(shí)花了很多時(shí)間研究,學(xué)了不少東西,一個(gè)簡單的作業(yè)可以涉及各方面的知識(shí)。$ V4 i3 B) [+ K# y3 i' o! j* s
畢業(yè)半年了,雖然記不清彈性矩陣、剛度矩陣等,但是只要一看書(好像現(xiàn)在的心境不下來)應(yīng)該會(huì),學(xué)校交的學(xué)習(xí)方法和理解問題的能力,有些東西不必記,除非每天從事,那該叫熟能生巧。6 u# K- B2 H7 H0 g
記得當(dāng)時(shí)編了兩個(gè)程序,梁和三維實(shí)體的,我分享一下梁的程序吧,起碼有亮點(diǎn)和創(chuàng)意,可以自己選擇劃分個(gè)數(shù),在matlab方面花了不少功夫。
' m: c; L. J7 V3 K1 P/ _) v1 \非常感謝當(dāng)時(shí)的有限元授課老師“韓清凱”,雖然老師已經(jīng)到了大工。上課的教材為韓清凱 的《彈性力學(xué)及有限元法基礎(chǔ)教程》,書上有個(gè)梁單元的例子,在82頁。
% w2 R' z/ F2 P3 f# t--------------------------------------------------------------------------------
( Q0 H- j3 r' Q梁程序如下,名字就不寫了,用昵稱吧“太平島”,這是常用的網(wǎng)絡(luò)昵稱
4 C, ?* V. B" a. a% p
) I/ F5 I0 d7 l( Q9 [- g; S# t4 U: v# T  ?$ k3 n
%------------------------ BEAM2  ----------------
" E) J/ b% S  l* I9 v, Xdisp('========================================');
- @  n) [5 J/ I+ f/ S) Ndisp('            PROGRAM BEAM2               ');: s1 I" J5 ~8 F: [; }4 ]  i* E
disp('        Beam Bending Analysis           ');
! s: v. h" P6 \6 L& S7 M5 G9 ldisp('        The programmer:太平島           ');! u( E( h( Y5 D3 X* I2 q! `% B8 s- {: c
disp('========================================');
2 _' B" ]- [8 f- F$ P& Xdisp(' ');                        %輸出空行
. i6 R4 ]. p& T& ^: F4 x6 zwarning off all                        %關(guān)閉所有警告提示(求積分時(shí),警告信息比較多)
& t2 D+ V5 O+ i; ^- u7 ONe=input('Please input the number of elements:');        %Ne為劃分的單元數(shù)
7 e0 s& v. |4 H' h. ]" n0 U) q5 ENJ=Ne+1;                        %NJ為節(jié)點(diǎn)數(shù)
9 Q; N% r% ~+ n" N: v, D& w) G# ex1=0;
: i5 ~2 p3 H9 V# G9 }x2=sym('L');' E" u2 t; R) A) G
x=sym('x');                        
) l4 u* q% Z7 w7 Z2 {# |j=0:3;- t" i0 W3 @& h7 A( J
v=x.^j;
, j9 \! g8 W/ m+ f3 i2 ]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];
$ J3 Y. }. k& JN=v*inv(A);                        %形函數(shù); S7 f# ~- d" q- ~
%-----------------------------------------------
! o* G, f$ B8 p) Q% 求單元?jiǎng)偠染仃?font class="jammer">5 w4 i+ t, V+ K+ s

E=4.0e11;7 |" k; m  y; h; G
I=5.2e-7;                        %I=bh^3/12=5.2e-7;
: i. z: `8 S7 ^4 N+ t9 ~" z+ `EI=4.0e11*5.2e-7;( V+ |/ ^& j  r+ |9 U6 t1 _$ K
B=diff(N,2);
% N1 t( i8 v  H: w- N7 f. ~k=B'*B;
% T0 L; W* a5 F9 i- J* b1 G. xKee=EI*int(k,x,0,'L');
6 m& u- d7 r3 y3 y7 {( n% I, TKe=sym(zeros(4,4,Ne));        %用三維矩陣存放單元?jiǎng)偠染仃?每頁存一個(gè),并初始化0
  u! y, L9 D, h7 HKe(:,:,1)=subs(Kee,'L',(10/Ne));
" B4 v3 u$ d: OT=eye(4);                % 因?yàn)榱河趚軸的夾角為0,所以坐標(biāo)變換矩陣為T=eye(4)' |; ~- \1 p0 v* q$ R0 ?( U
Ke(:,:,1)=T*Ke(:,:,1)*T';/ R* t" g& l0 Z# |( z4 }
for ii=2:Ne! }! @0 Z, i; Y3 @5 _. u
    Ke(:,:,ii)=Ke(:,:,1);7 I' q, A: W9 ?
end
) ?# ~" x. |; H  S% [( XKe=double(Ke);& S2 F) `8 P- {" W
%------------------------------------------------/ q2 @) Z: T. P2 T5 r6 t8 P
% 由矩陣裝換法求整體剛度矩陣
8 s4 x' t5 A6 Q. |, j% 自由度Ndof=2*NJ2 L) o  ~& V/ g) T% d% _1 u% N
Ndof=2*NJ;
3 S+ ]- e& E$ P+ l2 x4 \% i$ oK=zeros(Ndof);
) j6 v: B) F! N+ `& i) Pfor ii=1:Ne
( P4 y# s( f4 Z% g7 H9 }1 O6 e, |2 C    G=[zeros(4,2*(ii-1)),eye(4,Ndof-2*(ii-1))];* y/ H6 @/ C7 f1 T1 [/ e6 d" a4 V
    KK=G'*Ke(:,:,Ne)*G;4 {1 k# i2 v) M& j1 I; N1 H4 \
    K=K+KK;
1 W5 B1 E) i+ A# F& Vend  6 [0 N2 A0 u6 k. ?& y$ h% A) Z  `
%------------------------------------------------
. ?% f9 a- t; Q% 約束條件,對整體剛度矩陣進(jìn)行修正,以便計(jì)算
+ e( N- m# I, A0 z7 L2 q( x; k3 \F=zeros(Ndof,1);. \# P; [$ r" q8 D  V0 E
F(Ndof-1)=-100;
1 [- W9 Z2 |  ~1 a3 Z8 ~2 C% u=F*inv(K),u=[v1,xta1,v2,xta2,v3,xta3,v4,xta4]',v1=0,xta1=0) Y# [2 @. V( ?$ b- K% J8 Z
K(1,=0;        K(:,1)=0;        %可以省略
  i8 v8 \1 ?  {1 g, D+ DK(2,=0;        K(:,2)=0;        %可以省略
9 `( e5 T4 Z; q+ m( mKX=K(3:Ndof,3:Ndof);
7 H0 n) O, w" x6 NFX=F(3:Ndof,1);
! ~0 u  ]1 \" g' q9 j* ~%------------------------------------------------
* Y; p1 h" V$ U+ c7 O, C" B%求整體節(jié)點(diǎn)位移列陣5 `; J2 X. c: M  _! _
uu=inv(KX)*FX;/ Z( N  u4 [& R+ }; v* M
u=[0;0;uu];
& X8 V+ r) P$ k9 f. w4 Lii=1:2:2*NJ;* y+ p: Z+ v! o2 j
v=u(ii);                        %各節(jié)點(diǎn)撓度
" c. I. M' L  @0 f( Cxta=u(ii+1);                        %各節(jié)點(diǎn)轉(zhuǎn)角" d  a9 A* T& a( E
%------------------------------------------------  ^$ e& l: N" l4 s4 p
% 后處理,計(jì)算單元應(yīng)變應(yīng)力
1 U& Y/ C( T* C: H+ kB=subs(B,x,'L');" L. d! X; F: g" |
B=subs(B,'L',10/Ne);2 d, G. k; w" R( B* Q- B
Strain=zeros(1,Ne);                %單元應(yīng)變,并初始化- `- V! R9 I0 e9 Q
Stress=zeros(1,Ne);                %單元應(yīng)力,并初始化
/ x4 s: A4 i% Rfor ii=1:Ne0 q  c; V, V/ c
    Strain(1,ii)=B*u(2*ii-1:2*ii+2,1);
& l0 Z2 t% N  l# I# h$ _9 L/ d    Stress(1,ii)=E*Strain(1,ii);: B1 M: d" B5 L" A5 h
end
1 c0 ^2 {: y( B8 l; D%--------------------------------------------------
0 U. x! v$ v( ^" a1 r$ e* U2 u% 以下程序?yàn)槠聊惠敵鎏幚?br /> " @( q' p- \& r& K/ j5 Odisp(' ');9 Z2 y6 ?; c& R' i
disp(' Input:1-print Node displacement ');% G/ A6 v7 z# y4 Q4 J+ Y" ?- X. `
disp(' Input:2-print strain ');& H/ ]8 U! U0 l9 s' R
disp(' Input:3-print stress ');* n  @- h+ s1 `/ w# B" C
disp(' Input:q-quit or use ''ctrl+c'' Shortcut to end the program');
& ^: P! M! r: Z, p% i3 l8 ?5 S
) R  u3 e4 R$ \- |4 C' `; y7 u/ Swhile 1/ x# L. B9 g* Z8 {* Y7 Z. m# K' n
    ii=input('Please input1 or 2 or 3:','s');
4 n6 d, }  f  [4 [8 I2 u1 G1 {, `$ x        switch(ii)0 {# m; f$ @% ?  E% [/ s' \, P
            case '1'0 P! n  x, U% T' ]4 H. E+ \/ g+ J
                disp('Node displacement');% z# i6 S1 H# R9 }3 X7 K# a
                disp(u');! e2 P' L- K5 w& j; J
            case '2'; V. |1 R( A! `7 i* z+ q( J; h$ D- s
                disp('element strain');" K5 V3 k6 z, w$ I8 G  v
                disp(Strain);) C& g' `$ f% x
            case '3'- |# w8 |' A9 n1 T5 ?
                disp('element stress');- b) [4 a% V6 e  ~* }4 {/ `
                disp(Stress);# ]* N6 U; |' o! z& [; @
            case 'q'8 ?" W/ F. r- W  Q% W* r
                disp('End of program');
8 `8 g# b) C5 J, M/ F% O                break;
. F9 X% x& w8 @  n( s" C            otherwise! R. x: E* B& D0 a) s9 C
                disp('error!Please input again');( F* E% D' E2 f3 r, x8 n( G
                continue;
- X0 n; O  M3 j  Q        end- F4 O4 K6 f: a+ }0 M+ M
end. O8 P% l7 Q0 x( v5 E% C( {7 r8 I# h, f

( `; I' i7 g: l3 {9 k
. S; z% q6 R) ]; F7 H4 q& A
回復(fù)

使用道具 舉報(bào)

2#
 樓主| 發(fā)表于 2013-3-24 13:44:32 | 只看該作者
那個(gè)地方怎么變成笑臉了,冒號+右括號=笑臉,改了下,應(yīng)該是下面,把英文括號改為中文,就好了吧
& J: z3 r# {& m% R- Q" eK(1,:)=0;        K(:,1)=0;        %可以省略
7 u# w. A1 a& s  YK(2,:)=0;        K(:,2)=0;        %可以省略
3#
發(fā)表于 2013-3-24 14:55:52 | 只看該作者
沒看懂
4#
發(fā)表于 2013-3-24 16:53:09 | 只看該作者
謝謝你啊,太感謝你了
5#
 樓主| 發(fā)表于 2013-3-24 18:35:59 | 只看該作者
jiaweicz 發(fā)表于 2013-3-24 16:53 ) ?* k% s6 {+ L- H, ?
謝謝你啊,太感謝你了
* m: R  E+ S5 E. l
這謝啥呢?
6#
發(fā)表于 2013-3-24 21:21:19 | 只看該作者
我以前編過一個(gè)5體展開的多體姿態(tài)動(dòng)力學(xué)計(jì)算程序。。。可惜早就忘記了,sigh

點(diǎn)評

是啊,很多東西不用就忘記啦  發(fā)表于 2013-3-26 12:34
7#
發(fā)表于 2013-11-7 20:39:02 | 只看該作者
這程序也運(yùn)行不了啊
8#
發(fā)表于 2013-11-7 20:44:34 | 只看該作者
j=0:3;" B5 u/ u1 z( M3 B! h
v=x.^j;  是干啥的
您需要登錄后才可以回帖 登錄 | 注冊會(huì)員

本版積分規(guī)則

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

GMT+8, 2025-7-25 13:57 , Processed in 0.095008 second(s), 24 queries , Gzip On.

Powered by Discuz! X3.5 Licensed

© 2001-2025 Discuz! Team.

快速回復(fù) 返回頂部 返回列表