国产精品乱码一区-性开放网站-少妇又紧又爽视频-西西大胆午夜人体视频-国产极品一区-欧美成人tv-四虎av在线-国产无遮挡无码视频免费软件-中文字幕亚洲乱码熟女一区二区-日产精品一区二区三区在线观看-亚洲国产亚综合在线区-五月婷婷综合色-亚洲日本视频在线观看-97精品人人妻人人-久久久久久一区二区三区四区别墅-www.免费av-波多野结衣绝顶大高潮-日本在线a一区视频高清视频-强美女免费网站在线视频-亚洲永久免费

<noframes id="4mu6g"></noframes>

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

標(biāo)題: 分析一段matlab有限元程序 [打印本頁]

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

%------------------------ BEAM2  ----------------
disp('========================================');
disp('          PROGRAM BEAM2             ');
disp('       Beam Bending Analysis          ');
disp('       The programmer:太平島          ');
disp('========================================');
disp(' ');                      %輸出空行
warning off all                      %關(guān)閉所有警告提示(求積分時，警告信息比較多)
Ne=input('Please input the number of elements:');       %Ne為劃分的單元數(shù)
NJ=Ne+1;                      %NJ為節(jié)點(diǎn)數(shù)
x1=0;
x2=sym('L');
x=sym('x');
j=0:3;
v=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];
N=v*inv(A);                      %形函數(shù)
%-----------------------------------------------
% 求單元剛度矩陣
E=4.0e11;
I=5.2e-7;                      %I=bh^3/12=5.2e-7；
EI=4.0e11*5.2e-7;
B=diff(N,2);
k=B'*B;
Kee=EI*int(k,x,0,'L');
Ke=sym(zeros(4,4,Ne));       %用三維矩陣存放單元剛度矩陣,每頁存一個，并初始化0
Ke(:,:,1)=subs(Kee,'L',(10/Ne));
T=eye(4);             % 因?yàn)榱河趚軸的夾角為0，所以坐標(biāo)變換矩陣為T=eye(4)
Ke(:,:,1)=T*Ke(:,:,1)*T';
for ii=2:Ne
Ke(:,:,ii)=Ke(:,:,1);
end
Ke=double(Ke);
%------------------------------------------------
% 由矩陣裝換法求整體剛度矩陣
% 自由度Ndof=2*NJ
Ndof=2*NJ;
K=zeros(Ndof);
for ii=1:Ne
G=[zeros(4,2*(ii-1)),eye(4,Ndof-2*(ii-1))];
KK=G'*Ke(:,:,Ne)*G;
K=K+KK;
end
%------------------------------------------------
% 約束條件,對整體剛度矩陣進(jìn)行修正，以便計(jì)算
F=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);
%------------------------------------------------
%求整體節(jié)點(diǎn)位移列陣
uu=inv(KX)*FX;
u=[0;0;uu];
ii=1:2:2*NJ;
v=u(ii);                      %各節(jié)點(diǎn)撓度
xta=u(ii+1);                      %各節(jié)點(diǎn)轉(zhuǎn)角
%------------------------------------------------
% 后處理，計(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
Strain(1,ii)=B*u(2*ii-1:2*ii+2,1);
Stress(1,ii)=E*Strain(1,ii);
end
%--------------------------------------------------
% 以下程序?yàn)槠聊惠敵鎏幚?br /> disp(' ');
disp(' 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
ii=input('Please input1 or 2 or 3:','s');
      switch(ii)
         case '1'
            disp('Node displacement');
            disp(u');
         case '2'
            disp('element strain');
            disp(Strain);
         case '3'
            disp('element stress');
            disp(Stress);
         case 'q'
            disp('End of program');
            break;
         otherwise
            disp('error!Please input again');
            continue;
      end
end

作者: yinzengguang 時間: 2013-3-24 13:44
那個地方怎么變成笑臉了，冒號+右括號=笑臉，改了下，應(yīng)該是下面，把英文括號改為中文，就好了吧
K(1,:）=0; K(:,1)=0; %可以省略
K(2,:）=0; K(:,2)=0; %可以省略

作者: 孤若流星 時間: 2013-3-24 14:55
沒看懂

作者: jiaweicz 時間: 2013-3-24 16:53
謝謝你啊，太感謝你了

作者: yinzengguang 時間: 2013-3-24 18:35

jiaweicz 發(fā)表于 2013-3-24 16:53 5 Y: X3 D$ z7 Z
謝謝你啊，太感謝你了

這謝啥呢？

作者: 懷念昔日熊貓 時間: 2013-3-24 21:21
我以前編過一個5體展開的多體姿態(tài)動力學(xué)計(jì)算程序。。。可惜早就忘記了，sigh

作者: d調(diào)旳華孋 時間: 2013-11-7 20:39
這程序也運(yùn)行不了啊

作者: d調(diào)旳華孋 時間: 2013-11-7 20:44
j=0:3;
v=x.^j; 是干啥的

歡迎光臨機(jī)械社區(qū) (http://www.whclglass.com.cn/)