?? exam8_1.m
字號:
function exam8_1
% 本程序為第八章的第一個算例，采用平面梁單元計算兩鉸拋物線拱的自由振動特性
% 輸入?yún)?shù)： 無
% 輸出結(jié)果： 前3階振動頻率及其相應(yīng)的振型
% 定義全局變量
% gNode ------ 節(jié)點坐標(biāo)
% gElement --- 單元定義
% gMaterial -- 材料性質(zhì)
% gBC1 ------- 第一類約束條件
% gK --------- 整體剛度矩陣
% gDelta ----- 整體節(jié)點坐標(biāo)
PlaneFrameModel ; % 定義有限元模型
SolveModel ; % 求解有限元模型
DisplayResults ; % 顯示計算結(jié)果
return ;
function PlaneFrameModel
% 定義平面桿系的有限元模型
% 輸入?yún)?shù)：
% 無
% 返回值：
% 無
% 說明：
% 該函數(shù)定義平面桿系的有限元模型數(shù)據(jù)：
% gNode ------- 節(jié)點定義
% gElement ---- 單元定義
% gMaterial --- 材料定義，包括彈性模量，梁的截面積和梁的抗彎慣性矩
% gBC --------- 約束條件
global gNode gElement gMaterial gBC1
% 給定拋物線拱的幾何特征
L = 60 ; % 計算跨徑
f = 7.5 ; % 計算矢高
n = 100 ; % 單元數(shù)目
x = -L/2:L/n:L/2 ; % 結(jié)點的x坐標(biāo)
a = f/L^2*4 ;
y = - a * x.^2 ; % 結(jié)點的y坐標(biāo)
% 節(jié)點坐標(biāo)
gNode = [x' y'] ;
% 單元定義
gElement = zeros( n, 3 ) ;
for i=1:n
gElement( i, : ) = [ i, i+1, 1 ] ;
end
% 材料性質(zhì)
% 彈性模量 抗彎慣性矩 截面積 密度
gMaterial = [2.06e11, 0.03622, 0.0815, 1435.2/0.0815]; % 材料 1
% 第一類約束條件
% 節(jié)點號 自由度號 約束值
gBC1 = [ 1, 1, 0.0
1, 2, 0.0
n+1, 1, 0.0
n+1, 2, 0.0] ;
return
function SolveModel
% 求解有限元模型
% 輸入?yún)?shù)：
% 無
% 返回值：
% 無
% 說明：
% 該函數(shù)求解有限元模型，過程如下
% 1. 計算單元的剛度和質(zhì)量矩陣，集成整體剛度和質(zhì)量矩陣
% 2. 處理約束條件，修改整體剛度矩陣
% 3. 求解特征值問題
global gNode gElement gMaterial gBC1 gK gM gEigValue gEigVector
% step1. 定義整體剛度矩陣和節(jié)點力向量
[node_number,dummy] = size( gNode ) ;
gK = sparse( node_number * 3, node_number * 3 ) ;
gM = sparse( node_number * 3, node_number * 3 ) ;
% step2. 計算單元剛度和質(zhì)量矩陣，并集成到整體剛度和質(zhì)量矩陣中
[element_number,dummy] = size( gElement ) ;
for ie=1:1:element_number
k = StiffnessMatrix( ie ) ;
m = MassMatrix( ie ) ;
AssembleGlobalMatrix( ie, k, m ) ;
end
% step3. 處理第一類約束條件，修改剛度矩陣和質(zhì)量矩陣。（采用劃行劃列法）
[bc1_number,dummy] = size( gBC1 ) ;
w2max = max( diag(gK)./diag(gM) ) ;
for ibc=1:1:bc1_number
n = gBC1(ibc, 1 ) ;
d = gBC1(ibc, 2 ) ;
m = (n-1)*3 + d ;
gK(:,m) = zeros( node_number*3, 1 ) ;
gK(m,:) = zeros( 1, node_number*3 ) ;
gK(m,m) = 1;
gM(:,m) = zeros( node_number*3, 1 ) ;
gM(m,:) = zeros( 1, node_number*3 ) ;
gM(m,m) = gK(m,m)/w2max/1e10 ;
end
% step4. 求解特征值問題
% step4.1為了使剛度矩陣和質(zhì)量矩陣對稱（在計算時可能引入舍入誤差）
for i=1:node_number*3
for j=i:node_number*3
gK(j,i) = gK(i,j) ;
gM(j,i) = gM(i,j) ;
end
end
% step4.2 計算前6階特征值和特征向量
[gEigVector, gEigValue] = eigs(gK, gM, 3, 'SM' ) ;
% step4.3 修改特征向量中受約束的自由度
for ibc=1:1:bc1_number
n = gBC1(ibc, 1 ) ;
d = gBC1(ibc, 2 ) ;
m = (n-1)*3 + d ;
gEigVector(m,:) = gBC1(ibc,3) ;
end
return
function k = StiffnessMatrix( ie )
% 計算單元剛度矩陣
% 輸入?yún)?shù):
% ie ------- 單元號
% 返回值:
% k ---- 整體坐標(biāo)系下的剛度矩陣
global gNode gElement gMaterial
k = zeros( 6, 6 ) ;
E = gMaterial( gElement(ie, 3), 1 ) ;
I = gMaterial( gElement(ie, 3), 2 ) ;
A = gMaterial( gElement(ie, 3), 3 ) ;
xi = gNode( gElement( ie, 1 ), 1 ) ;
yi = gNode( gElement( ie, 1 ), 2 ) ;
xj = gNode( gElement( ie, 2 ), 1 ) ;
yj = gNode( gElement( ie, 2 ), 2 ) ;
L = ( (xj-xi)^2 + (yj-yi)^2 )^(1/2) ;
k = [ E*A/L 0 0 -E*A/L 0 0
0 12*E*I/L^3 6*E*I/L^2 0 -12*E*I/L^3 6*E*I/L^2
0 6*E*I/L^2 4*E*I/L 0 -6*E*I/L^2 2*E*I/L
-E*A/L 0 0 E*A/L 0 0
0 -12*E*I/L^3 -6*E*I/L^2 0 12*E*I/L^3 -6*E*I/L^2
0 6*E*I/L^2 2*E*I/L 0 -6*E*I/L^2 4*E*I/L] ;
T = TransformMatrix( ie ) ;
k = T*k*transpose(T) ;
return
function m = MassMatrix( ie )
% 計算單元質(zhì)量矩陣
% 輸入?yún)?shù):
% ie ------- 單元號
% 返回值:
% m ---- 整體坐標(biāo)系下的質(zhì)量矩陣
global gNode gElement gMaterial
m = zeros( 6, 6 ) ;
E = gMaterial( gElement(ie, 3), 1 ) ;
A = gMaterial( gElement(ie, 3), 3 ) ;
ro = gMaterial( gElement(ie, 3 ), 4 ) ;
xi = gNode( gElement( ie, 1 ), 1 ) ;
yi = gNode( gElement( ie, 1 ), 2 ) ;
xj = gNode( gElement( ie, 2 ), 1 ) ;
yj = gNode( gElement( ie, 2 ), 2 ) ;
L = ( (xj-xi)^2 + (yj-yi)^2 )^(1/2) ;
m = ro*A*L/420*[140 0 0 70 0 0
0 156 22*L 0 54 -13*L
0 22*L 4*L^2 0 13*L -3*L^2
70 0 0 140 0 0
0 54 13*L 0 156 -22*L
0 -13*L -3*L^2 0 -22*L 4*L^2 ] ;
T = TransformMatrix( ie ) ;
m = T*m*transpose(T) ;
return
function AssembleGlobalMatrix( ie, ke, me )
% 把單元剛度和質(zhì)量矩陣集成到整體剛度矩陣
% 輸入?yún)?shù):
% ie --- 單元號
% ke --- 單元剛度矩陣
% me --- 單元質(zhì)量矩陣
% 返回值:
% 無
global gElement gK gM
for i=1:1:2
for j=1:1:2
for p=1:1:3
for q =1:1:3
m = (i-1)*3+p ;
n = (j-1)*3+q ;
M = (gElement(ie,i)-1)*3+p ;
N = (gElement(ie,j)-1)*3+q ;
gK(M,N) = gK(M,N) + ke(m,n) ;
gM(M,N) = gM(M,N) + me(m,n) ;
end
end
end
end
return
function T = TransformMatrix( ie )
% 計算單元的坐標(biāo)轉(zhuǎn)換矩陣( 局部坐標(biāo) -> 整體坐標(biāo) )
% 輸入?yún)?shù)
% ie ----- 節(jié)點號
% 返回值
% T ------- 從局部坐標(biāo)到整體坐標(biāo)的坐標(biāo)轉(zhuǎn)換矩陣
global gElement gNode
xi = gNode( gElement( ie, 1 ), 1 ) ;
yi = gNode( gElement( ie, 1 ), 2 ) ;
xj = gNode( gElement( ie, 2 ), 1 ) ;
yj = gNode( gElement( ie, 2 ), 2 ) ;
L = sqrt( (xj-xi)^2 + (yj-yi)^2 ) ;
c = (xj-xi)/L ;
s = (yj-yi)/L ;
T=[ c -s 0 0 0 0
s c 0 0 0 0
0 0 1 0 0 0
0 0 0 c -s 0
0 0 0 s c 0
0 0 0 0 0 1] ;
return
function DisplayResults
% 顯示計算結(jié)果
% 輸入?yún)?shù)：
% 無
% 返回值：
% 無
global gNode gElement gMaterial gBC1 gEigValue gEigVector
fre_number = length(diag(gEigValue)) ;
% 打印特征向量(振型)
fprintf( '\n\n 表一 特征向量(振型) \n' ) ;
for i=1:fre_number
fprintf( '----------------') ;
end
fprintf( '\n' ) ;
for i=1:fre_number
fprintf( ' %6d ', i ) ;
end
fprintf( '\n' ) ;
for i=1:fre_number
fprintf( '----------------') ;
end
fprintf( '\n' ) ;
[dof,dummy]=size(gEigVector) ;
for i=1:dof
for j=fre_number:-1:1
fprintf( '%15.7e ', gEigVector(i,j) ) ;
end
fprintf( '\n' ) ;
end
for i=1:fre_number
fprintf( '----------------') ;
end
fprintf( '\n' ) ;
% 打印特征值
fprintf( '\n\n\n\n 表二 特征值(頻率)列表 \n' ) ;
fprintf( '----------------------------------------------------------\n') ;
fprintf( ' 階數(shù) 特征值 頻率(Hz) 圓頻率(Hz) \n' ) ;
fprintf( '----------------------------------------------------------\n') ;
for i=fre_number:-1:1
fprintf( '%6d %15.7e %15.7e %15.7e\n', fre_number-i+1, ...
gEigValue(i,i), sqrt(gEigValue(i,i))/2/pi, sqrt(gEigValue(i,i)) ) ;
end
fprintf( '----------------------------------------------------------\n') ;
% 繪制振型圖
for j=fre_number:-1:1
figure ;
x = gNode(:,1) ;
y = gNode(:,2) ;
dx = gEigVector(1:3:length(x)*3, j ) ;
dy = gEigVector(2:3:length(x)*3, j ) ;
factor = max( [max(abs(x))/max(abs(dx)), max(abs(y))/max(abs(dy))] )* 0.05;
plot(x,y,'-', x+factor*dx, y+factor*dy, ':') ;
title( sprintf( '第%d階頻率: %.3f Hz', fre_number-j+1, sqrt(gEigValue(j,j))/2/pi ) ) ;
axis equal;
axis off ;
end
return
?? 快捷鍵說明
復(fù)制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -