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

 找回密碼
 注冊會員

QQ登錄

只需一步,快速開始

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

分析一段matlab有限元程序

[復(fù)制鏈接]
1#
發(fā)表于 2013-3-24 13:41:49 | 只看該作者 |倒序瀏覽 |閱讀模式
這是關(guān)于梁單元的,可能大家覺得很簡單。。。
3 ?: V3 h! s7 Y3 L5 r- j今天翻電腦里的東西時發(fā)現(xiàn)的,是我大三時有限元時的作業(yè),記得當(dāng)時花了很多時間研究,學(xué)了不少東西,一個簡單的作業(yè)可以涉及各方面的知識。2 w. j- Z5 m! v7 q+ J5 ]
畢業(yè)半年了,雖然記不清彈性矩陣、剛度矩陣等,但是只要一看書(好像現(xiàn)在的心境不下來)應(yīng)該會,學(xué)校交的學(xué)習(xí)方法和理解問題的能力,有些東西不必記,除非每天從事,那該叫熟能生巧。# E7 ^0 B* j6 ]5 N
記得當(dāng)時編了兩個程序,梁和三維實體的,我分享一下梁的程序吧,起碼有亮點和創(chuàng)意,可以自己選擇劃分個數(shù),在matlab方面花了不少功夫。% o  Z- d6 t5 v' @; @- r' ?1 d
非常感謝當(dāng)時的有限元授課老師“韓清凱”,雖然老師已經(jīng)到了大工。上課的教材為韓清凱 的《彈性力學(xué)及有限元法基礎(chǔ)教程》,書上有個梁單元的例子,在82頁。; z4 X4 C' {9 a6 m* e) {. |( E
--------------------------------------------------------------------------------; ~8 m/ c4 X1 K$ @3 ]1 z: ]
梁程序如下,名字就不寫了,用昵稱吧“太平島”,這是常用的網(wǎng)絡(luò)昵稱
; N  F* }, k2 H) D7 {2 @6 W% q  M0 a( N

3 R3 o! P* f: x) u* m%------------------------ BEAM2  ----------------
# l& k" @4 S6 V. ~( Odisp('========================================');
+ ^! `: A( h3 U1 |disp('            PROGRAM BEAM2               ');8 M7 v0 w+ ~4 K
disp('        Beam Bending Analysis           ');
) d! E1 d( t% g& z+ e8 fdisp('        The programmer:太平島           ');
/ K! k0 u0 G9 edisp('========================================');0 ^; v  ]4 i- J; [6 f  Q' Q, A& r
disp(' ');                        %輸出空行
: ~, o" _" R/ n3 Awarning off all                        %關(guān)閉所有警告提示(求積分時,警告信息比較多)
# H; q3 p0 M+ z8 m5 C& L) m) K+ ?Ne=input('Please input the number of elements:');        %Ne為劃分的單元數(shù)
! B  N; N) s" @8 D  iNJ=Ne+1;                        %NJ為節(jié)點數(shù)
8 Z' @( O6 `1 ?* Z: yx1=0;3 t" t* W& Q" x0 a+ H2 ]
x2=sym('L');9 T- D' |- |2 n& f3 K+ W0 p
x=sym('x');                        3 V" N6 R+ j: O& o3 F  g$ D+ D* f
j=0:3;
/ V4 o% z# c9 G" Ov=x.^j;
+ j& p' o6 [+ g/ O5 J9 q* D; WA=[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];( N9 N% O0 R6 o! O! Q, T
N=v*inv(A);                        %形函數(shù)" B3 _5 l! v9 W# R  i
%-----------------------------------------------
1 \# w- q% P4 Z4 k/ y% 求單元剛度矩陣
) ~7 H/ P: J$ m* uE=4.0e11;
9 T  d; s- Y9 u% }I=5.2e-7;                        %I=bh^3/12=5.2e-7;8 ]0 D: g" _0 v2 X4 i) b+ }
EI=4.0e11*5.2e-7;
/ o) T$ m, P- ^/ L) R$ }B=diff(N,2);* M' ]9 I* C  U% g
k=B'*B;
  U9 d3 n6 o) a4 b: iKee=EI*int(k,x,0,'L');- u: S3 w; w5 S  q' o+ h8 g+ T/ V
Ke=sym(zeros(4,4,Ne));        %用三維矩陣存放單元剛度矩陣,每頁存一個,并初始化02 Z4 @0 r; n1 |" b8 ~0 Q. ]
Ke(:,:,1)=subs(Kee,'L',(10/Ne));. u, t" C3 n1 G9 a* r' E1 i' Y3 v# S
T=eye(4);                % 因為梁于x軸的夾角為0,所以坐標(biāo)變換矩陣為T=eye(4)8 N8 g) S  h/ T* E8 F6 q
Ke(:,:,1)=T*Ke(:,:,1)*T';- L2 j8 o8 K* v6 ^8 f5 H
for ii=2:Ne
* S8 s7 {+ q, L1 A- x5 z    Ke(:,:,ii)=Ke(:,:,1);
9 S8 j2 E2 |5 I. b0 O' cend - f3 o) x* v! V' Y+ F9 X2 H8 R; {
Ke=double(Ke);" L9 R; F7 Q( q/ C9 B
%------------------------------------------------
9 j9 H+ j% U; t4 v' g! o9 l% 由矩陣裝換法求整體剛度矩陣. C7 E/ j2 b' W" k8 c8 y: X
% 自由度Ndof=2*NJ% b% ^3 A2 N! {: ^9 i
Ndof=2*NJ;6 A- O, d: j5 Q1 U
K=zeros(Ndof);
6 K7 y& ?7 T8 u3 n- j. qfor ii=1:Ne
+ ]; a# g* L5 Y, T! @    G=[zeros(4,2*(ii-1)),eye(4,Ndof-2*(ii-1))];0 L1 l8 U$ E* m+ v& y
    KK=G'*Ke(:,:,Ne)*G;
$ }) u6 p3 y! P  }    K=K+KK;* T5 u* Z' b0 \
end  
$ `1 [1 _( z; P8 d& n%------------------------------------------------  i% w" M& P/ b' y- T2 q, j
% 約束條件,對整體剛度矩陣進行修正,以便計算
* t" [& N7 K. b: H; qF=zeros(Ndof,1);# j! |9 C4 m/ L4 h2 P: H9 M
F(Ndof-1)=-100;5 |/ @+ D9 [- v' J/ B, F+ V: ]7 y
% u=F*inv(K),u=[v1,xta1,v2,xta2,v3,xta3,v4,xta4]',v1=0,xta1=0! u3 n8 e) j7 z+ r4 N$ H# Y) f
K(1,=0;        K(:,1)=0;        %可以省略$ ~- n7 H$ |* r
K(2,=0;        K(:,2)=0;        %可以省略7 F) M" z  w% H* Y2 `
KX=K(3:Ndof,3:Ndof);
1 g1 L- k% B. @7 B- L+ ^FX=F(3:Ndof,1);9 P7 c7 n* V* ?) C  B
%------------------------------------------------
/ k4 w8 o' p2 S%求整體節(jié)點位移列陣
# L0 e' k; u0 `$ l  Xuu=inv(KX)*FX;
+ i7 i2 \* S- s4 v+ y$ T3 H6 Su=[0;0;uu];, W+ E! k) `& S' P" f6 H2 U, L
ii=1:2:2*NJ;
/ }% L/ ^5 o( o2 Fv=u(ii);                        %各節(jié)點撓度, P$ u9 u* @4 b3 R) l: ^! v
xta=u(ii+1);                        %各節(jié)點轉(zhuǎn)角: Z0 k% }; ?, N$ |4 }8 T. v
%------------------------------------------------
: a4 P% i+ E6 ]( D% 后處理,計算單元應(yīng)變應(yīng)力0 B* C0 L4 D9 |1 B& j6 \& o# e+ o
B=subs(B,x,'L');
$ O5 |0 Y8 x5 `! \: rB=subs(B,'L',10/Ne);; M7 w/ F4 v, T/ P9 Z
Strain=zeros(1,Ne);                %單元應(yīng)變,并初始化
5 X+ W4 w7 A* x* {+ \1 t% s5 b# o6 `Stress=zeros(1,Ne);                %單元應(yīng)力,并初始化4 E( X* j( V7 B7 u0 p
for ii=1:Ne
- s" E$ _+ t# @% h7 W+ W2 k    Strain(1,ii)=B*u(2*ii-1:2*ii+2,1);
" A' D0 }( r4 e# ?9 b1 C" a: z    Stress(1,ii)=E*Strain(1,ii);" T9 x+ e7 `1 h& P
end
' x4 Y" |. Q: [' g  {3 V6 C%--------------------------------------------------
1 G9 A' k& G) t# A8 U5 [$ y% 以下程序為屏幕輸出處理" s& c$ U' Y# `' G# }
disp(' ');) ~: Z8 t! f% x# b# w# _% W
disp(' Input:1-print Node displacement ');
9 }$ a# ?" O' ]disp(' Input:2-print strain ');
2 n- K! |" T  h* L5 @$ R- [- adisp(' Input:3-print stress ');
  b' n& j# ^' h0 E" \: Y: z+ Hdisp(' Input:q-quit or use ''ctrl+c'' Shortcut to end the program');
2 K* `. I  ^4 @0 \; z, a' O- q3 ^3 j3 X% w, `2 g( ~0 V0 E
while 13 f. d/ A6 u8 e
    ii=input('Please input1 or 2 or 3:','s');2 f* j5 z3 V* N% [/ z  a. J: h* ?
        switch(ii)
2 p7 x& D: x8 l            case '1'
" t) r8 p6 W$ o* N9 X+ t" Y                disp('Node displacement');
" K0 g( a9 P! c2 r                disp(u');
- W, ~' K; B" A: L            case '2'# h3 v& H4 M. h; e& K9 c& `' B5 v) l
                disp('element strain');
8 L6 l& P8 L$ B1 `4 s                disp(Strain);
9 `2 b( E# g6 I$ X- u            case '3', ?1 n& g9 t; @3 d
                disp('element stress');) Q$ t2 t7 ]6 S" F+ [) ]" b
                disp(Stress);
; S% p, _& y7 k& ~* r2 y            case 'q'8 _0 l: N4 B$ ?
                disp('End of program');
; l5 y: u  N% @$ D% p  W5 P' j/ A! b                break;
+ N' W0 _0 v8 V            otherwise8 V* ~/ N) h: ~0 x# ?" b
                disp('error!Please input again');, C0 Z/ t7 L/ M+ L* ?( C# I6 Q
                continue;( e& R) D5 v4 q. B
        end
9 E  A. F$ t1 u* zend
( _/ H6 X6 ~1 [, w
/ }: g+ |* K- a+ ~+ V5 _2 E- v4 ?

! Y4 ?/ A. ?. t
回復(fù)

使用道具 舉報

2#
 樓主| 發(fā)表于 2013-3-24 13:44:32 | 只看該作者
那個地方怎么變成笑臉了,冒號+右括號=笑臉,改了下,應(yīng)該是下面,把英文括號改為中文,就好了吧
6 l7 x5 ]6 m/ }* lK(1,:)=0;        K(:,1)=0;        %可以省略4 {- x  s* V9 }- h3 S8 k$ W
K(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 ( `2 E+ o8 h2 _) Z" A' Z
謝謝你啊,太感謝你了

+ p* [' m4 y0 U4 S) X這謝啥呢?
6#
發(fā)表于 2013-3-24 21:21:19 | 只看該作者
我以前編過一個5體展開的多體姿態(tài)動力學(xué)計算程序。。。可惜早就忘記了,sigh

點評

是啊,很多東西不用就忘記啦  發(fā)表于 2013-3-26 12:34
7#
發(fā)表于 2013-11-7 20:39:02 | 只看該作者
這程序也運行不了啊
8#
發(fā)表于 2013-11-7 20:44:34 | 只看該作者
j=0:3;
0 |( Z9 ?" P; |9 b: q& h- O) nv=x.^j;  是干啥的
您需要登錄后才可以回帖 登錄 | 注冊會員

本版積分規(guī)則

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

GMT+8, 2025-7-24 18:00 , Processed in 0.079684 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.5 Licensed

© 2001-2025 Discuz! Team.

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