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

機(jī)械社區(qū)

標(biāo)題: 分析一段matlab有限元程序 [打印本頁]
作者: yinzengguang    時(shí)間: 2013-3-24 13:41
標(biāo)題: 分析一段matlab有限元程序
這是關(guān)于梁單元的，可能大家覺得很簡單。。。
2 X8 O9 r9 `7 S- t" b7 }% R今天翻電腦里的東西時(shí)發(fā)現(xiàn)的，是我大三時(shí)有限元時(shí)的作業(yè)，記得當(dāng)時(shí)花了很多時(shí)間研究，學(xué)了不少東西，一個(gè)簡單的作業(yè)可以涉及各方面的知識。1 i) G. n6 c" H7 _  F+ G
畢業(yè)半年了，雖然記不清彈性矩陣、剛度矩陣等，但是只要一看書（好像現(xiàn)在的心境不下來）應(yīng)該會，學(xué)校交的學(xué)習(xí)方法和理解問題的能力，有些東西不必記，除非每天從事，那該叫熟能生巧。
6 X& X! j) @+ ]6 X  R! ?記得當(dāng)時(shí)編了兩個(gè)程序，梁和三維實(shí)體的，我分享一下梁的程序吧，起碼有亮點(diǎn)和創(chuàng)意，可以自己選擇劃分個(gè)數(shù)，在matlab方面花了不少功夫。
" l$ ^3 v  G  g$ ^: x7 U6 b/ ~8 F非常感謝當(dāng)時(shí)的有限元授課老師“韓清凱”，雖然老師已經(jīng)到了大工。上課的教材為韓清凱 的《彈性力學(xué)及有限元法基礎(chǔ)教程》，書上有個(gè)梁單元的例子，在82頁。& E6 c/ \# I, m0 O" S& t
--------------------------------------------------------------------------------7 R) z* U8 a  w' v. \
梁程序如下，名字就不寫了，用昵稱吧“太平島”，這是常用的網(wǎng)絡(luò)昵稱( r/ [, ^# x! t) K) D0 s' C
$ A! S& d. M2 P
4 y+ m6 w6 z, [" p
%------------------------ BEAM2  ----------------$ j. l6 ~8 d6 g- C8 _5 N3 a9 R& B
disp('========================================');
: h, W6 \4 u' y; x3 ~1 }) ldisp('            PROGRAM BEAM2               ');
7 b3 m& v1 [2 {6 |! Jdisp('        Beam Bending Analysis           ');, W5 ?0 u! N0 w! ?) O
disp('        The programmer:太平島           ');8 g0 r! H* N# ^5 a6 G
disp('========================================');6 Y- i+ H9 D* F0 b7 W. x$ z3 K& [
disp(' ');                        %輸出空行
, S7 x" G8 d* Q  q1 x* Owarning off all                        %關(guān)閉所有警告提示(求積分時(shí)，警告信息比較多)
6 T. r; E# j! m3 x' ?Ne=input('Please input the number of elements:');        %Ne為劃分的單元數(shù)
. {1 L4 h6 `! \$ u+ tNJ=Ne+1;                        %NJ為節(jié)點(diǎn)數(shù)
6 N: k& q" v5 b3 C7 p7 fx1=0;% O. o3 K* p; m$ P
x2=sym('L');8 N/ Q3 R, M& M! \2 R6 N
x=sym('x');                        
) _9 \2 z& s: g, D7 n5 Ij=0:3;2 Y. ^4 h& S, a3 t7 X7 I
v=x.^j;
, u$ `" W" a% p  k8 sA=[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];
3 S0 C, K* M( cN=v*inv(A);                        %形函數(shù)
" G4 O! x+ L; K%-----------------------------------------------
/ ^4 O7 b5 N& \. S$ `2 A2 A5 ~0 ?% 求單元?jiǎng)偠染仃?font class="jammer">8 [2 `: }, {5 I1 p; I
E=4.0e11;, d8 h4 r4 g& U: Y' M! t$ b9 n! X
I=5.2e-7;                        %I=bh^3/12=5.2e-7；
' o3 s* G0 \( `+ V" mEI=4.0e11*5.2e-7;* V* @& L2 x# O0 f2 i( j& C
B=diff(N,2);: D7 i9 V& w, V0 W4 y$ I( F; k
k=B'*B;# |& l7 w: e' W/ x, O( c
Kee=EI*int(k,x,0,'L');
# P7 Q7 {( J: A: H. k' YKe=sym(zeros(4,4,Ne));        %用三維矩陣存放單元?jiǎng)偠染仃?每頁存一個(gè)，并初始化0- j4 O3 x+ U1 K7 ]; W' ^
Ke(:,:,1)=subs(Kee,'L',(10/Ne));4 A5 C" j. |+ U9 |
T=eye(4);                % 因?yàn)榱河趚軸的夾角為0，所以坐標(biāo)變換矩陣為T=eye(4). ]+ c: H+ P4 D. m
Ke(:,:,1)=T*Ke(:,:,1)*T';2 g, H# p. d- x3 j2 H# n" k
for ii=2:Ne5 ]% P" D. o/ p$ {% V
    Ke(:,:,ii)=Ke(:,:,1);( |0 U# u* s4 W1 P! `5 G) u( e
end 3 ]' a6 [  {- [/ e4 @' P
Ke=double(Ke);
1 `! m$ n  X/ k+ C8 G%------------------------------------------------6 h/ l3 ?8 p2 t: w! P
% 由矩陣裝換法求整體剛度矩陣
& u* U: E( v: g) ^% 自由度Ndof=2*NJ
# ?6 Z' U. i  `) A8 E( W( z8 lNdof=2*NJ;
' D$ b4 r* ^* _0 q2 T8 p" @; R0 @K=zeros(Ndof);
2 r9 n+ M( S) ofor ii=1:Ne2 P% R* g- Z" z' j$ E9 s
    G=[zeros(4,2*(ii-1)),eye(4,Ndof-2*(ii-1))];/ f& V0 G, p" `5 V# h2 I# x* h5 S
    KK=G'*Ke(:,:,Ne)*G;
/ S2 m% O. q7 k9 I7 u. f. N    K=K+KK;
4 R" F, Y. P3 Z' Q- K" {end  ' l, N# F( M+ d# B1 @
%------------------------------------------------8 r4 T' R- W3 n5 U
% 約束條件,對整體剛度矩陣進(jìn)行修正，以便計(jì)算
) y: v  V3 E  h0 N% E* wF=zeros(Ndof,1);
5 w: w* L8 G# r8 CF(Ndof-1)=-100;) V( A5 j, t9 W/ t
% u=F*inv(K),u=[v1,xta1,v2,xta2,v3,xta3,v4,xta4]',v1=0,xta1=0
* q- G0 N% u9 m& Q' s$ WK(1,=0;        K(:,1)=0;        %可以省略
# |/ |9 y9 R4 U- u+ d! I" T9 }K(2,=0;        K(:,2)=0;        %可以省略) ]+ R% P6 |; k. g3 f
KX=K(3:Ndof,3:Ndof);8 [$ B' E' S( A
FX=F(3:Ndof,1);
  s( P  C$ Z; O. \# B, d. n  x%------------------------------------------------8 M. D# L! q. l4 o
%求整體節(jié)點(diǎn)位移列陣/ v, a9 C+ d+ m1 A9 I0 W/ `
uu=inv(KX)*FX;
5 ~0 g7 B, K) d" |( m& _$ `8 H4 ru=[0;0;uu];
) v7 z/ }7 u6 h- {, Cii=1:2:2*NJ;
) A; S1 _  L# _% u5 x3 zv=u(ii);                        %各節(jié)點(diǎn)撓度
" ~! i- W' }$ o) ?+ G* E0 Yxta=u(ii+1);                        %各節(jié)點(diǎn)轉(zhuǎn)角/ |' H! J1 @/ y* A. S
%------------------------------------------------! P1 K% R4 ?) Z8 Q: z, [' @, L) |
% 后處理，計(jì)算單元應(yīng)變應(yīng)力
3 F& q; P, i- ~+ _& v/ O& ?B=subs(B,x,'L');
5 D9 X: \7 M! N  `* ]; w$ vB=subs(B,'L',10/Ne);* `) ~8 e& Y: _- y  T/ ?# [: y4 _
Strain=zeros(1,Ne);                %單元應(yīng)變，并初始化( d) ]/ n* ]/ I$ P' ~
Stress=zeros(1,Ne);                %單元應(yīng)力，并初始化
) c& y8 r. s$ d! u  e$ J' Rfor ii=1:Ne+ o/ n; y+ V0 Q7 H: N
    Strain(1,ii)=B*u(2*ii-1:2*ii+2,1);! c6 D; S) i( A2 Q+ z. u
    Stress(1,ii)=E*Strain(1,ii);
9 b# U0 f: M# {end6 f- {" D' n1 a" y7 m* ~
%--------------------------------------------------3 m8 ]% |9 W, j
% 以下程序?yàn)槠聊惠敵鎏幚?font class="jammer">  m( f9 I8 J3 P, H0 U7 q8 n
disp(' ');
9 D# ~; v7 l+ b& D/ _$ udisp(' Input:1-print Node displacement ');- C0 u% k: ?/ C8 x
disp(' Input:2-print strain ');9 }1 ]3 l0 u) s  O9 K1 X
disp(' Input:3-print stress ');
# H0 l" I/ n- }2 }; E# K( I6 ldisp(' Input:q-quit or use ''ctrl+c'' Shortcut to end the program');: o9 o0 }) k% s/ G& _2 M

, ]7 h) q4 @( h5 W& Qwhile 18 g0 a. F* k/ J" F+ g1 b/ a) O
    ii=input('Please input1 or 2 or 3:','s');
1 r8 ?( \; R9 Q        switch(ii)$ |9 W6 J) s, V! U0 l+ T, ^3 d, J
            case '1'$ ~! s; ?" V* v
                disp('Node displacement');
5 [/ A' E% B: i0 i, k                disp(u');
6 B7 H* H  k, z2 Z: q7 m            case '2'; n. I- I! [( v$ G
                disp('element strain');
- A4 X- U3 z  e# s                disp(Strain);- i, n( K  T* b, u/ J
            case '3'4 d" e* w. t  A9 p% H
                disp('element stress');8 R3 A" [( `0 K' {/ H
                disp(Stress);
2 i9 n. w4 T: r2 [- ]            case 'q'
6 P% s. [  m2 {# C) d; p                disp('End of program');" N( s5 T, c' k" p- p* e2 ]3 k$ j
                break;. _, A$ ?* |$ S$ Y% P9 D
            otherwise+ p) M9 G' O4 G6 ?( _
                disp('error!Please input again');& B2 ~* B" J6 D# s
                continue;
# F7 D$ t& T) m) k6 |        end
4 g* O3 F$ k, G- d% W1 ?end
3 ?& o$ G7 h; Z# U3 k: W- D4 a4 R! N) Q/ q/ t- P( _/ [, r& N

1 q* T* ~1 A  N0 L5 j! A7 G- P
作者: yinzengguang    時(shí)間: 2013-3-24 13:44
那個(gè)地方怎么變成笑臉了，冒號+右括號=笑臉，改了下，應(yīng)該是下面，把英文括號改為中文，就好了吧' b: A3 J- M/ [. |. `8 R- C$ T
K(1,:）=0;        K(:,1)=0;        %可以省略+ c) e* V' n+ ~- Q/ s6 q0 w. p
K(2,:）=0;        K(:,2)=0;        %可以省略
作者: 孤若流星    時(shí)間: 2013-3-24 14:55
沒看懂
作者: jiaweicz    時(shí)間: 2013-3-24 16:53
謝謝你啊，太感謝你了
作者: yinzengguang    時(shí)間: 2013-3-24 18:35
jiaweicz 發(fā)表于 2013-3-24 16:53 ( L5 x* r+ }' ]& Q: Q- Z
謝謝你啊，太感謝你了

3 B$ W9 Y( i$ ~9 h- i- t這謝啥呢？
作者: 懷念昔日熊貓    時(shí)間: 2013-3-24 21:21
我以前編過一個(gè)5體展開的多體姿態(tài)動(dòng)力學(xué)計(jì)算程序。。。可惜早就忘記了，sigh
作者: d調(diào)旳華孋    時(shí)間: 2013-11-7 20:39
這程序也運(yùn)行不了啊
作者: d調(diào)旳華孋    時(shí)間: 2013-11-7 20:44
j=0:3;
& O6 T. g: `& S" f, Rv=x.^j;  是干啥的



歡迎光臨 機(jī)械社區(qū) (http://www.ytsybjq.com/) Powered by Discuz! X3.5