亚洲欧美日韩国产一区二区精品_亚洲国产精品一区二区动图_级婬片A片手机免费播放_亚洲国产成人Av毛片大全,男女爱爱好爽好疼视频免费,中文日韩AV在线,无码视频免费,欧美在线观看成人高清视频,在线播放免费人成毛片,成 人 网 站 在 线 视 频A片 ,亚洲AV成人精品一区二区三区

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

標(biāo)題: 分析一段matlab有限元程序 [打印本頁]
作者: yinzengguang    時(shí)間: 2013-3-24 13:41
標(biāo)題: 分析一段matlab有限元程序
這是關(guān)于梁單元的,，可能大家覺得很簡單,。,。。8 m& Q0 S) N+ c! q6 I0 Q' x+ P8 p
今天翻電腦里的東西時(shí)發(fā)現(xiàn)的,，是我大三時(shí)有限元時(shí)的作業(yè)，記得當(dāng)時(shí)花了很多時(shí)間研究,，學(xué)了不少東西,，一個(gè)簡單的作業(yè)可以涉及各方面的知識(shí)。3 ^  `7 P7 X, R
畢業(yè)半年了,，雖然記不清彈性矩陣,、剛度矩陣等，但是只要一看書（好像現(xiàn)在的心境不下來）應(yīng)該會(huì),，學(xué)校交的學(xué)習(xí)方法和理解問題的能力,，有些東西不必記，除非每天從事,，那該叫熟能生巧,。
8 B/ f0 l8 I8 K- I* S8 t記得當(dāng)時(shí)編了兩個(gè)程序,，梁和三維實(shí)體的，我分享一下梁的程序吧,，起碼有亮點(diǎn)和創(chuàng)意,，可以自己選擇劃分個(gè)數(shù)，在matlab方面花了不少功夫,。
; ^2 y& ~. x; O1 U# \非常感謝當(dāng)時(shí)的有限元授課老師“韓清凱”,，雖然老師已經(jīng)到了大工。上課的教材為韓清凱 的《彈性力學(xué)及有限元法基礎(chǔ)教程》,，書上有個(gè)梁單元的例子,，在82頁。
' J. ~$ X3 ~# e$ E' l: P) S5 n% \--------------------------------------------------------------------------------
. `' r' M& x2 }梁程序如下,，名字就不寫了,，用昵稱吧“太平島”，這是常用的網(wǎng)絡(luò)昵稱* _9 a2 b; \2 Y% r( E

2 E, m7 P7 \, R& K  d
: F4 H5 W4 C& z% b7 b%------------------------ BEAM2  ----------------
* h* I6 N0 h6 gdisp('========================================');
( t: L3 ~! v5 V) {  B! o0 u0 b7 [* ?disp('            PROGRAM BEAM2               ');
1 a$ z- X( x0 y: g9 {8 |; A* [disp('        Beam Bending Analysis           ');. ?% I( c* A( O+ B5 N% ]
disp('        The programmer:太平島           ');  c8 h. E7 n0 h* P$ G/ j
disp('========================================');
0 V. x) G7 Y  G$ Edisp(' ');                        %輸出空行
. _" t1 s- H' t" x4 E$ x) _+ Iwarning off all                        %關(guān)閉所有警告提示(求積分時(shí),，警告信息比較多)' K( c) F) @8 X) j) z1 ]4 W
Ne=input('Please input the number of elements:');        %Ne為劃分的單元數(shù)0 m( g) O# ?1 h) }9 r8 }
NJ=Ne+1;                        %NJ為節(jié)點(diǎn)數(shù)1 D$ S4 f& N/ y+ ?9 d/ T
x1=0;
, O, s% y& J$ u1 K* s0 Rx2=sym('L');
0 r# H2 ?$ A, [2 Nx=sym('x');                        
. @( g$ h6 [; q0 H5 X6 T6 dj=0:3;
: O6 V$ E! v2 Av=x.^j;& P0 {  B: G2 H9 U6 n3 P+ m
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];
( d& ?: ~0 M+ y" U/ M" I. `N=v*inv(A);                        %形函數(shù)
! p& T% i4 T0 J4 T( h# m# A%-----------------------------------------------/ p& F% i  r# s! Z$ s9 T$ a
% 求單元?jiǎng)偠染仃?font class="jammer">8 K* [: U/ ?3 ]% m$ v% q# C
E=4.0e11;
+ G) _0 t+ y* S4 Q7 r& n+ w$ hI=5.2e-7;                        %I=bh^3/12=5.2e-7,；
5 Q) {0 a: i  b. n3 ~( |2 ~EI=4.0e11*5.2e-7;
$ N) i& k; I, A) ^B=diff(N,2);
( {7 o) j3 j0 o3 h4 {: q- ik=B'*B;
+ S& f6 n7 a/ x" _! a4 AKee=EI*int(k,x,0,'L');
& ?9 Y; p, ~  vKe=sym(zeros(4,4,Ne));        %用三維矩陣存放單元?jiǎng)偠染仃?每頁存一個(gè)，并初始化0
7 H5 }1 o: }, f. c! |Ke(:,:,1)=subs(Kee,'L',(10/Ne));* t: u2 e) C8 E% \
T=eye(4);                % 因?yàn)榱河趚軸的夾角為0,，所以坐標(biāo)變換矩陣為T=eye(4)
+ Y2 N3 N- h; `" \/ rKe(:,:,1)=T*Ke(:,:,1)*T';* S" h! c8 T2 M5 T
for ii=2:Ne6 l4 V9 R; ?9 F" a/ T( D
    Ke(:,:,ii)=Ke(:,:,1);6 o7 h7 Q1 S3 d7 a( m' m5 ^( v
end
; E9 }* P* C6 l' ^. hKe=double(Ke);
: X( F6 `" s% {%------------------------------------------------% U8 R; C" n7 ?; `
% 由矩陣裝換法求整體剛度矩陣
" R! v# v9 y# L3 ^9 D" z3 t' f+ ^% 自由度Ndof=2*NJ7 I+ m7 B  J0 P
Ndof=2*NJ;2 c+ d. s& W+ n  M) Y
K=zeros(Ndof);
$ z8 a* S2 S/ g( B) Mfor ii=1:Ne8 z$ u: S! g* i3 L& K" V# y
    G=[zeros(4,2*(ii-1)),eye(4,Ndof-2*(ii-1))];
0 W6 e+ G% E! z) t8 M    KK=G'*Ke(:,:,Ne)*G;9 \3 I- o. C, d
    K=K+KK;
' L0 u3 G" ?# s% d( c/ b2 kend  # e, j( \0 Z5 v- P$ G8 q
%------------------------------------------------/ p. r- O# B0 c" s' z1 M5 f
% 約束條件,對(duì)整體剛度矩陣進(jìn)行修正,，以便計(jì)算
# W, f9 x3 j% w0 y5 I' J* I" \5 t2 V" d  Z7 dF=zeros(Ndof,1);  M. ^3 W9 w/ _1 F& c4 M. {
F(Ndof-1)=-100;( V! w( A' b3 F: Z  }: Y# z
% u=F*inv(K),u=[v1,xta1,v2,xta2,v3,xta3,v4,xta4]',v1=0,xta1=0
* M, h) X& [  v7 EK(1,=0;        K(:,1)=0;        %可以省略
9 S2 W9 h( a0 }* x# wK(2,=0;        K(:,2)=0;        %可以省略* |4 U4 ^. `& g2 S+ {8 t
KX=K(3:Ndof,3:Ndof);4 O5 {& S% T$ y! A, p7 q+ F
FX=F(3:Ndof,1);3 m- ^  Y: E4 ~3 Z* a! s! ^6 v5 E
%------------------------------------------------
7 e$ Q1 o" y6 l+ R8 w6 e%求整體節(jié)點(diǎn)位移列陣
- A) d0 |& v$ f! |- [uu=inv(KX)*FX;
. @6 _. m3 k! [2 Vu=[0;0;uu];, p$ t* h& b$ [) l1 z9 E
ii=1:2:2*NJ;
) e5 d$ G0 I4 a( f) i3 uv=u(ii);                        %各節(jié)點(diǎn)撓度7 _4 _0 J7 t$ N1 P; V* j& }
xta=u(ii+1);                        %各節(jié)點(diǎn)轉(zhuǎn)角1 G- C( X, G' d" Z; i6 H  w' G; _3 J
%------------------------------------------------2 P+ f& X; F+ P& l" `7 y' d) x
% 后處理，計(jì)算單元應(yīng)變應(yīng)力, u7 y6 W4 A+ q( L5 ?0 N5 h
B=subs(B,x,'L');$ C! P$ p. D& L9 n5 m
B=subs(B,'L',10/Ne);2 m3 {9 R2 f  q. U" Y
Strain=zeros(1,Ne);                %單元應(yīng)變,，并初始化
/ u3 U8 a" {# M, t! r3 z' @Stress=zeros(1,Ne);                %單元應(yīng)力,，并初始化
7 T7 j. W! I# y! d2 G7 u3 Yfor ii=1:Ne
" S; x  p6 l  J0 B, D! B    Strain(1,ii)=B*u(2*ii-1:2*ii+2,1);
7 _4 z3 {/ V- \9 ~    Stress(1,ii)=E*Strain(1,ii);
% c: }2 l9 D$ `& @5 U4 Jend2 p, ~9 T0 \. y! R3 E
%--------------------------------------------------
, h; h/ f2 B/ n7 J! |2 O% 以下程序?yàn)槠聊惠敵鎏幚?font class="jammer">5 F* l. o6 ?6 ^: l0 Q5 }: a* ^: D+ f
disp(' ');5 p* F2 D" H. N& g  d
disp(' Input:1-print Node displacement ');. [" o7 Q3 |* Y" f' G8 O- x+ }
disp(' Input:2-print strain ');# N! j) H: z( ^& @) n
disp(' Input:3-print stress ');
6 w6 f; h# t* n" r. X  s6 G; tdisp(' Input:q-quit or use ''ctrl+c'' Shortcut to end the program');, f& S0 ^: c- E, E8 E0 K5 i( B
9 R, y, \7 Y4 j4 i' X; S
while 13 x: C" ]4 J, N6 d; ?. X' x' R6 E
    ii=input('Please input1 or 2 or 3:','s');
0 O; C8 L% ?8 p% v# D3 C. k        switch(ii)! b0 T: u! @8 K# ?  R
            case '1'' ]3 I3 F. E( P6 y8 Z
                disp('Node displacement');
  z8 S  y+ X. H( |& b( L3 B3 h                disp(u');/ k4 k2 g* x* h* V) l5 S
            case '2'1 X, [: a2 |/ r  l5 ?8 z6 H
                disp('element strain');3 _+ X. v! U: w$ J9 E1 x
                disp(Strain);
$ ^! H4 T8 W6 ]# _% c1 d8 k3 ]            case '3'3 s! M3 ^6 V9 a2 A
                disp('element stress');
% F: S3 r1 m6 y. r                disp(Stress);4 F* {, @, z* ]) I* F+ I1 S- ?
            case 'q'
" z5 l3 U2 n$ Z( ?( ~8 k/ ?                disp('End of program');
4 V0 S$ H% c3 U! |% M                break;
% z% ~/ y. E3 L+ J" @, c+ l) S            otherwise4 H3 q) T1 Z  s, u6 |
                disp('error!Please input again');! ]  `7 ^& ]4 x$ S. T! @2 R
                continue;" H% C9 v& F1 I+ G$ Y
        end
7 T' ^. P# g1 Q. Qend" h# a. O* \8 `: D
" N5 `) `& a" [+ _% }

' o3 i% L+ u$ \. ]# z- |+ R. x8 r
作者: yinzengguang    時(shí)間: 2013-3-24 13:44
那個(gè)地方怎么變成笑臉了，冒號(hào)+右括號(hào)=笑臉,，改了下,，應(yīng)該是下面，把英文括號(hào)改為中文,，就好了吧7 k% X& l- _! P7 q  |
K(1,:）=0;        K(:,1)=0;        %可以省略% e6 E* j4 Z' i
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 3 F5 H3 L1 [- |4 Q
謝謝你啊，太感謝你了

' ]0 y6 `- ]$ |+ \6 W/ h, E這謝啥呢,？
作者: 懷念昔日熊貓    時(shí)間: 2013-3-24 21:21
我以前編過一個(gè)5體展開的多體姿態(tài)動(dòng)力學(xué)計(jì)算程序,。。,�,？上г缇屯浟耍瑂igh
作者: d調(diào)旳華孋    時(shí)間: 2013-11-7 20:39
這程序也運(yùn)行不了啊
作者: d調(diào)旳華孋    時(shí)間: 2013-11-7 20:44
j=0:3;
9 I, H% h9 V% x& vv=x.^j;  是干啥的



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