分析一段matlab有限元程序

yinzengguang

這是關(guān)于梁單元的,，可能大家覺得很簡單。,。。
今天翻電腦里的東西時(shí)發(fā)現(xiàn)的,，是我大三時(shí)有限元時(shí)的作業(yè),，記得當(dāng)時(shí)花了很多時(shí)間研究，學(xué)了不少東西,，一個(gè)簡單的作業(yè)可以涉及各方面的知識(shí),。
: X( T$ v+ l! c$ k: v$ K畢業(yè)半年了，雖然記不清彈性矩陣,、剛度矩陣等,，但是只要一看書（好像現(xiàn)在的心境不下來）應(yīng)該會(huì)，學(xué)校交的學(xué)習(xí)方法和理解問題的能力,，有些東西不必記,，除非每天從事,，那該叫熟能生巧。
記得當(dāng)時(shí)編了兩個(gè)程序，梁和三維實(shí)體的,，我分享一下梁的程序吧，起碼有亮點(diǎn)和創(chuàng)意,，可以自己選擇劃分個(gè)數(shù),，在matlab方面花了不少功夫。
; k' W. l: T/ L/ X* L非常感謝當(dāng)時(shí)的有限元授課老師“韓清凱”,，雖然老師已經(jīng)到了大工,。上課的教材為韓清凱 的《彈性力學(xué)及有限元法基礎(chǔ)教程》，書上有個(gè)梁單元的例子,，在82頁,。
/ T8 o/ x; N* D--------------------------------------------------------------------------------
梁程序如下，名字就不寫了,，用昵稱吧“太平島”,，這是常用的網(wǎng)絡(luò)昵稱


%------------------------ BEAM2  ----------------
3 t* v6 i+ Q9 p2 A1 Cdisp('========================================');
disp('            PROGRAM BEAM2               ');
disp('        Beam Bending Analysis           ');
disp('        The programmer:太平島           ');
: \0 S9 h; U  d  a8 I$ tdisp('========================================');
) X2 k4 N$ t8 T' h/ tdisp(' ');                        %輸出空行
warning off all                        %關(guān)閉所有警告提示(求積分時(shí),，警告信息比較多)
+ l0 L( [9 p0 s# E3 INe=input('Please input the number of elements:');        %Ne為劃分的單元數(shù)
. m8 n2 |: m* O, \- K0 k4 Y3 bNJ=Ne+1;                        %NJ為節(jié)點(diǎn)數(shù)
# |) M+ {- u* h# r- k' ~. Px1=0;
x2=sym('L');
x=sym('x');                        
j=0:3;
; I" N( S) ]" n+ G" H/ @  Cv=x.^j;
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];
, Q+ X/ J8 x$ A" ~5 IN=v*inv(A);                        %形函數(shù)
%-----------------------------------------------
% 求單元?jiǎng)偠染仃?font class="jammer">/ ~9 a5 Y$ E9 e; T. [( D5 f+ |9 P
E=4.0e11;
. Z$ F, z+ d8 U6 G" H, c* \I=5.2e-7;                        %I=bh^3/12=5.2e-7；
EI=4.0e11*5.2e-7;
B=diff(N,2);
9 R+ |6 |( n; ]4 F& sk=B'*B;
Kee=EI*int(k,x,0,'L');
2 C' R* ?: f2 T! Z( sKe=sym(zeros(4,4,Ne));        %用三維矩陣存放單元?jiǎng)偠染仃?每頁存一個(gè),，并初始化0
Ke(:,:,1)=subs(Kee,'L',(10/Ne));
4 ~% D/ h+ e. vT=eye(4);                % 因?yàn)榱河趚軸的夾角為0,，所以坐標(biāo)變換矩陣為T=eye(4)
Ke(:,:,1)=T*Ke(:,:,1)*T';
* `* _# F9 c9 g3 M; p" M0 U5 zfor ii=2:Ne
    Ke(:,:,ii)=Ke(:,:,1);
9 Z/ S  w1 i5 I! g8 J$ ]end
; S$ P4 f! ?) g7 A& D& T. t4 A+ WKe=double(Ke);
* B, c* s% ^( b, _%------------------------------------------------
% b4 Z4 v; `& o% @* `) u% 由矩陣裝換法求整體剛度矩陣
% 自由度Ndof=2*NJ
: H" r5 o: I, I, lNdof=2*NJ;
K=zeros(Ndof);
/ r" x: u5 b/ X2 |" J6 t5 _/ o7 _6 }for ii=1:Ne
    G=[zeros(4,2*(ii-1)),eye(4,Ndof-2*(ii-1))];
    KK=G'*Ke(:,:,Ne)*G;
+ w  ?3 W9 D1 ?1 @. f    K=K+KK;
end  
%------------------------------------------------
& L6 Y1 z, n. ]& O% l/ X. J% 約束條件,對(duì)整體剛度矩陣進(jìn)行修正，以便計(jì)算
, w7 _6 Z% s( UF=zeros(Ndof,1);
F(Ndof-1)=-100;
% u=F*inv(K),u=[v1,xta1,v2,xta2,v3,xta3,v4,xta4]',v1=0,xta1=0
K(1,=0;        K(:,1)=0;        %可以省略
K(2,=0;        K(:,2)=0;        %可以省略
KX=K(3:Ndof,3:Ndof);
FX=F(3:Ndof,1);
2 B7 }* _! S8 H  R& o& A%------------------------------------------------
& K) f9 I4 _4 {8 K5 i%求整體節(jié)點(diǎn)位移列陣
uu=inv(KX)*FX;
u=[0;0;uu];
ii=1:2:2*NJ;
v=u(ii);                        %各節(jié)點(diǎn)撓度
* s* @0 t  u! Jxta=u(ii+1);                        %各節(jié)點(diǎn)轉(zhuǎn)角
( D" e1 x! U# B2 i%------------------------------------------------
9 p) ]! j! k- `6 E% 后處理,，計(jì)算單元應(yīng)變應(yīng)力
B=subs(B,x,'L');
B=subs(B,'L',10/Ne);
Strain=zeros(1,Ne);                %單元應(yīng)變,，并初始化
Stress=zeros(1,Ne);                %單元應(yīng)力，并初始化
for ii=1:Ne
& v7 V" T, p2 d& {    Strain(1,ii)=B*u(2*ii-1:2*ii+2,1);
( a$ d$ N8 F$ n5 y( h1 a# R# V5 b    Stress(1,ii)=E*Strain(1,ii);
end
5 o) V. g$ k1 a- \2 h! m%--------------------------------------------------
% 以下程序?yàn)槠聊惠敵鎏幚?br /> 9 p' L% B+ m2 _disp(' ');
* K7 h, e, T* x6 zdisp(' Input:1-print Node displacement ');
disp(' Input:2-print strain ');
disp(' Input:3-print stress ');
disp(' Input:q-quit or use ''ctrl+c'' Shortcut to end the program');

while 1
3 l3 V3 R# e% h; c$ |; Z    ii=input('Please input1 or 2 or 3:','s');
& l, h4 K  T/ v        switch(ii)
4 l  u% u8 g8 q+ T8 N8 j            case '1'
                disp('Node displacement');
                disp(u');
            case '2'
                disp('element strain');
8 I* M7 b  x9 M- d+ o- d                disp(Strain);
            case '3'
* T$ ^" N- W" j3 W  j0 U+ y8 K7 Y                disp('element stress');
7 x! D$ ]% M1 p" _! S                disp(Stress);
, t/ _; O& Q# [# y" N            case 'q'
* i; _# |! M/ D# K* |) l( F                disp('End of program');
                break;
8 N# {7 t3 |1 h$ o4 S7 E; |            otherwise
/ ?7 u# N, Z, L, c                disp('error!Please input again');
                continue;
        end
end
) m! B3 J, X4 }( @3 M$ f

yinzengguang

那個(gè)地方怎么變成笑臉了,，冒號(hào)+右括號(hào)=笑臉,，改了下，應(yīng)該是下面,，把英文括號(hào)改為中文,，就好了吧
K(1,:）=0;        K(:,1)=0;        %可以省略
K(2,:）=0;        K(:,2)=0;        %可以省略

孤若流星

沒看懂

jiaweicz

謝謝你啊，太感謝你了

yinzengguang

jiaweicz 發(fā)表于 2013-3-24 16:53
7 p* M. M" _+ T! d1 \+ \謝謝你啊,，太感謝你了

1 B# V+ M  D9 |  ~9 o這謝啥呢,？

懷念昔日熊貓

我以前編過一個(gè)5體展開的多體姿態(tài)動(dòng)力學(xué)計(jì)算程序。,。,。可惜早就忘記了,，sigh

d調(diào)旳華孋

這程序也運(yùn)行不了啊

d調(diào)旳華孋

j=0:3;
9 w& m* p6 k7 Z" U; Q" u5 p  F- ev=x.^j;  是干啥的