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

 找回密碼
 注冊會員

QQ登錄

只需一步,快速開始

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

分析一段matlab有限元程序

[復(fù)制鏈接]
1#
發(fā)表于 2013-3-24 13:41:49 | 只看該作者 |倒序瀏覽 |閱讀模式
這是關(guān)于梁單元的,可能大家覺得很簡單。。。
4 ?2 }4 O4 \+ X, H今天翻電腦里的東西時發(fā)現(xiàn)的,是我大三時有限元時的作業(yè),記得當(dāng)時花了很多時間研究,學(xué)了不少東西,一個簡單的作業(yè)可以涉及各方面的知識。3 h8 ?( i/ t) L) P8 H9 m
畢業(yè)半年了,雖然記不清彈性矩陣、剛度矩陣等,但是只要一看書(好像現(xiàn)在的心境不下來)應(yīng)該會,學(xué)校交的學(xué)習(xí)方法和理解問題的能力,有些東西不必記,除非每天從事,那該叫熟能生巧。6 K% ?' O5 Y/ T  U/ [
記得當(dāng)時編了兩個程序,梁和三維實體的,我分享一下梁的程序吧,起碼有亮點和創(chuàng)意,可以自己選擇劃分個數(shù),在matlab方面花了不少功夫。
# D& M" l+ O6 m非常感謝當(dāng)時的有限元授課老師“韓清凱”,雖然老師已經(jīng)到了大工。上課的教材為韓清凱 的《彈性力學(xué)及有限元法基礎(chǔ)教程》,書上有個梁單元的例子,在82頁。# o/ r6 Z% n* R: d2 h- A
--------------------------------------------------------------------------------, k$ E5 \0 d9 u0 I
梁程序如下,名字就不寫了,用昵稱吧“太平島”,這是常用的網(wǎng)絡(luò)昵稱
7 u: X: w; v' R1 d, ]: q- _! ~- L" R" H0 _% F
# F) D" H2 T7 T  ^  `, a
%------------------------ BEAM2  ----------------
8 f0 P! e- p. D( g( t  Adisp('========================================');
! n& `- b/ `; L  C* zdisp('            PROGRAM BEAM2               ');9 z6 |( O, S1 K+ ~  A& s+ B
disp('        Beam Bending Analysis           ');
! A6 s; X9 ^+ udisp('        The programmer:太平島           ');
4 X1 ^* N1 Z" P* t# F' l/ O! ndisp('========================================');
) U3 k" R- v% ~' Fdisp(' ');                        %輸出空行7 H! a4 v- o& G6 ~% P
warning off all                        %關(guān)閉所有警告提示(求積分時,警告信息比較多)
( N' h1 S6 P, KNe=input('Please input the number of elements:');        %Ne為劃分的單元數(shù)
/ v6 `1 d' p& \NJ=Ne+1;                        %NJ為節(jié)點數(shù)9 S: e9 `. N0 o
x1=0;
+ Y* T, d& A% w& Y+ y" C1 Vx2=sym('L');
! S: R. O# x5 q% H2 F' |0 Cx=sym('x');                        8 d' [$ b; z( J4 e. H8 R
j=0:3;9 h1 P! |% Z$ S& d6 W- k# h5 Z. c+ `" Q
v=x.^j;1 a, y0 {/ s! o1 T# I4 E4 ?7 H8 \
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];
( K. S2 D; ]5 z$ d8 @- A( y7 DN=v*inv(A);                        %形函數(shù)
& }" G( l! P3 N1 t%-----------------------------------------------$ F( H8 h5 `3 S5 E
% 求單元剛度矩陣
* s* Q& D& u8 J/ h) lE=4.0e11;3 _% ?% A& q$ G. I9 b/ [, u8 c
I=5.2e-7;                        %I=bh^3/12=5.2e-7;
( c3 _1 n% I; ^0 |EI=4.0e11*5.2e-7;
! K3 A" q" G# J, \, k6 ^7 |1 wB=diff(N,2);
) i" ?2 X9 V( O0 t1 Z9 ?k=B'*B;
* E. t3 M0 X' t" [Kee=EI*int(k,x,0,'L');( @3 e4 G! t1 o+ X, c
Ke=sym(zeros(4,4,Ne));        %用三維矩陣存放單元剛度矩陣,每頁存一個,并初始化0
  {3 ^3 n; z' @/ y! ]& B- iKe(:,:,1)=subs(Kee,'L',(10/Ne));
) @  i( h4 S8 @" E4 U; x, [$ C2 sT=eye(4);                % 因為梁于x軸的夾角為0,所以坐標(biāo)變換矩陣為T=eye(4)
% J* _9 H  J# Y* ^% `Ke(:,:,1)=T*Ke(:,:,1)*T';
: \& P! O  s' f6 ?2 _4 afor ii=2:Ne) _2 [; V) C3 E. {; D8 A; A2 X/ W
    Ke(:,:,ii)=Ke(:,:,1);8 U) v, Z# P5 T: j
end
/ `# B8 e& ?+ a& J) o( ?( SKe=double(Ke);* W2 V( }. b' r: M: b, h
%------------------------------------------------6 ^/ \+ Z8 `7 B  `
% 由矩陣裝換法求整體剛度矩陣# g* \9 Y7 [, Q$ B2 T! W) T
% 自由度Ndof=2*NJ
+ c: P; G: n1 i( \# [  ^Ndof=2*NJ;4 X5 D/ U; r/ S" j# h7 E
K=zeros(Ndof);$ T- `2 {; X) @/ R4 Z
for ii=1:Ne
* ]9 H/ A  o6 L- O    G=[zeros(4,2*(ii-1)),eye(4,Ndof-2*(ii-1))];; v; u) v  I1 h, w/ ]# `
    KK=G'*Ke(:,:,Ne)*G;
# N$ f9 [9 {. Z4 n$ l    K=K+KK;' T! ]9 N; C0 a/ ]8 C/ O; @
end  
9 k! |. J) O- w6 k%------------------------------------------------
( s! s) ~' V' p0 d+ O8 E% 約束條件,對整體剛度矩陣進行修正,以便計算
, c$ U  E) Y& iF=zeros(Ndof,1);! i( N5 U( {1 k3 j
F(Ndof-1)=-100;
  ~% z4 v, O+ D% t* z% u=F*inv(K),u=[v1,xta1,v2,xta2,v3,xta3,v4,xta4]',v1=0,xta1=0; w1 U+ h1 u# T0 t( p% t! p
K(1,=0;        K(:,1)=0;        %可以省略
( e# [  Z6 V  L- w) \/ \K(2,=0;        K(:,2)=0;        %可以省略
9 l9 M/ e; z, V4 o; |0 c0 ]# ZKX=K(3:Ndof,3:Ndof);7 ^+ G8 X1 R2 h% d( P+ P
FX=F(3:Ndof,1);
' {3 V4 o; u  |) }%------------------------------------------------# M; j; U* ^' q, L$ ~$ G
%求整體節(jié)點位移列陣
! ^, l, Z5 @& l5 e* g3 zuu=inv(KX)*FX;
5 ~& [3 w' B+ H7 d. yu=[0;0;uu];% K7 O  E2 Q( y+ M7 c- b9 \' K
ii=1:2:2*NJ;+ y2 b$ O* w- B0 a' e! r
v=u(ii);                        %各節(jié)點撓度+ C, k5 b$ d; |" E
xta=u(ii+1);                        %各節(jié)點轉(zhuǎn)角  K9 Z+ ?) R6 p2 B+ J' u
%------------------------------------------------
6 F) b* }* e- t1 V4 t0 w2 {9 i% 后處理,計算單元應(yīng)變應(yīng)力. q3 P! K9 y% q
B=subs(B,x,'L');
3 Y- ]. I' `! |3 M( o6 U) hB=subs(B,'L',10/Ne);0 v$ t7 |: h0 j3 Q
Strain=zeros(1,Ne);                %單元應(yīng)變,并初始化# z- `' E$ ^* r
Stress=zeros(1,Ne);                %單元應(yīng)力,并初始化
9 J3 q  @& @- p1 B" V  A7 L; f* Ifor ii=1:Ne; m. F9 Y* Y5 r" }5 d0 P3 p; }
    Strain(1,ii)=B*u(2*ii-1:2*ii+2,1);
$ N5 G0 q7 _! U( Y% J' b    Stress(1,ii)=E*Strain(1,ii);
; m( Y6 @1 }+ f1 q# \% cend5 ^; z9 _- c, I, _- L/ ?& [
%--------------------------------------------------
! B" Y% F  h5 K" V* \( g- U; O% 以下程序為屏幕輸出處理2 a& M+ T5 d# v$ v( R5 _& D3 u1 A. ?
disp(' ');' h8 w' l- C0 S! j2 e
disp(' Input:1-print Node displacement ');
- w& h3 `% f8 f* K. udisp(' Input:2-print strain ');
5 ~' U& P* K% \% E9 Q- ~disp(' Input:3-print stress ');
5 k/ E/ G# r2 l+ C8 l8 Q" edisp(' Input:q-quit or use ''ctrl+c'' Shortcut to end the program');# F9 v. e* F2 ?: k

9 z, K0 z2 n+ N# Uwhile 1
7 b6 E- D4 \7 C# V0 o* \4 u; g    ii=input('Please input1 or 2 or 3:','s');1 v9 k! ], S9 N! d8 v
        switch(ii)
/ l- v; f; f7 ]7 l            case '1'. |/ I& ]" X7 |: T4 y
                disp('Node displacement');  x! i. k& U* }5 G# Y* q
                disp(u');0 q, M- v. L# c- Z4 B
            case '2'& ]6 {" R' G& s$ ^4 a9 h
                disp('element strain');
) I. \5 f( C5 ]& U( d                disp(Strain);% T* s; T0 [- C; j( T
            case '3'
1 y0 i! O. w/ ~5 e6 Z$ p. u: ]                disp('element stress');' l  y% P# Q6 U! O3 s+ H7 @( k7 R
                disp(Stress);* y$ f" O. K2 S6 H5 u1 j, S2 B
            case 'q'  P" S$ _( z: w8 b% \3 p: L
                disp('End of program');! @3 l- U% ]0 c! V, Y5 @
                break;
$ m+ E. i& w$ S2 ~: x            otherwise( y) \% O5 i7 [, D4 W+ h! E
                disp('error!Please input again');9 M2 ^2 h& A# p" H6 e' v0 l
                continue;& h( l/ K3 q7 U( m
        end$ v3 Z1 r8 D8 ]
end) O7 g5 R& b4 e. m

( B4 Y2 \0 X8 b2 {' W' O
; J2 l; L4 t! a- _, e
回復(fù)

使用道具 舉報

2#
 樓主| 發(fā)表于 2013-3-24 13:44:32 | 只看該作者
那個地方怎么變成笑臉了,冒號+右括號=笑臉,改了下,應(yīng)該是下面,把英文括號改為中文,就好了吧3 s! J$ A8 {6 F
K(1,:)=0;        K(:,1)=0;        %可以省略7 a! J/ |9 H) O, x' d6 w7 l
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
/ r4 K) b( R# ^( v/ u謝謝你啊,太感謝你了

1 o( j( v# j, Z& U$ w5 z這謝啥呢?
6#
發(fā)表于 2013-3-24 21:21:19 | 只看該作者
我以前編過一個5體展開的多體姿態(tài)動力學(xué)計算程序。。??上г缇屯浟耍瑂igh

點評

是啊,很多東西不用就忘記啦  發(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;9 @1 Q) d. Z5 J4 O' r
v=x.^j;  是干啥的
您需要登錄后才可以回帖 登錄 | 注冊會員

本版積分規(guī)則

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

GMT+8, 2025-7-26 18:39 , Processed in 0.076520 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.5 Licensed

© 2001-2025 Discuz! Team.

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