?? frame.m
字號:
function exam3_1
% 本程序為第三章的第一個算例，采用平面梁單元計算平面剛架的變形和內力
% 輸入參數： 無
% 輸出結果： 節點位移和單元節點力
PlaneFrameModel ; % 定義有限元模型
SolveModel ; % 求解有限元模型
DisplayResults ; % 顯示計算結果
return ;
function PlaneFrameModel
% 定義平面桿系的有限元模型
% 輸入參數：
% 無
% 返回值：
% 無
% 說明：
% 該函數定義平面桿系的有限元模型數據：
% gNode ------- 節點定義
% gElement ---- 單元定義
% gMaterial --- 材料定義，包括彈性模量，梁的截面積和梁的抗彎慣性矩
% gBC1 -------- 約束條件
% gNF --------- 集中力
% gDF --------- 分布力
global gNode gElement gMaterial gBC1 gNF gDF
% 節點坐標
% x y
gNode = [0.0, 0.0 % 節點 1
0.0, 4.0 % 節點 2
3.0, 0.0 % 節點 3
3.0, 4.0 % 節點 4
4.5, 4.0 % 節點 5
6.0, 0.0 % 節點 6
6.0, 4.0 ] ; % 節點 7
% 單元定義
% 節點1 節點2 材料號
gElement = [1, 2, 1 % 單元 1
2, 4, 1 % 單元 2
3, 4, 1 % 單元 3
4, 5, 1 % 單元 4
5, 7, 1 % 單元 5
6, 7, 1] ; % 單元 6
% 材料性質
% 彈性模量 抗彎慣性矩 截面積
gMaterial = [2.1e11, 2.0e-4, 1.0e-2] ; % 材料 1
% 第一類約束條件
% 節點號 自由度號 約束值
gBC1 = [ 1, 1, 0.0
1, 2, 0.0
1, 3, 0.0
3, 1, 0.0
3, 2, 0.0
3, 3, 0.0
6, 1, 0.0
6, 2, 0.0
6, 3, 0.0] ;
% 集中力
% 節點號 自由度號 集中力值
gNF = [ 5, 2, -80e3] ;
% 分布載荷（線性分布）
% 單元號 節點1載荷值 節點2載荷值 自由度號
gDF = [ 1 -30e3 0 2
2 -15e3 -15e3 2 ] ;
return
function SolveModel
% 求解有限元模型
% 輸入參數：
% 無
% 返回值：
% 無
% 說明：
% 該函數求解有限元模型，過程如下
% 1. 計算單元剛度矩陣，集成整體剛度矩陣
% 2. 計算單元的等效節點力，集成整體節點力向量
% 3. 處理約束條件，修改整體剛度矩陣和節點力向量
% 4. 求解方程組，得到整體節點位移向量
global gNode gElement gMaterial gBC1 gNF gDF gK gDelta
% step1. 定義整體剛度矩陣和節點力向量
[node_number,dummy] = size( gNode ) ;
gK = sparse( node_number * 3, node_number * 3 ) ;
f = sparse( node_number * 3, 1 ) ;
% step2. 計算單元剛度矩陣，并集成到整體剛度矩陣中
[element_number,dummy] = size( gElement ) ;
for ie=1:1:element_number
k = StiffnessMatrix( ie, 1 ) ;
AssembleStiffnessMatrix( ie, k ) ;
end
% step3. 把集中力直接集成到整體節點力向量中
[nf_number, dummy] = size( gNF ) ;
for inf=1:1:nf_number
n = gNF( inf, 1 ) ;
d = gNF( inf, 2 ) ;
f( (n-1)*3 + d ) = gNF( inf, 3 ) ;
end
% step4. 計算分布力的等效節點力，不集成到整體節點力向量中
[df_number, dummy] = size( gDF ) ;
for idf = 1:1:df_number
enf = EquivalentNodeForce( gDF(idf,1), gDF(idf, 2), gDF( idf, 3), gDF( idf, 4 ) ) ;
i = gElement( gDF(idf,1), 1 ) ;
j = gElement( gDF(idf,1), 2 ) ;
f( (i-1)*3+1 : (i-1)*3+3 ) = f( (i-1)*3+1 : (i-1)*3+3 ) + enf( 1:3 ) ;
f( (j-1)*3+1 : (j-1)*3+3 ) = f( (j-1)*3+1 : (j-1)*3+3 ) + enf( 4:6 ) ;
end
% step5. 處理約束條件，修改剛度矩陣和節點力向量。采用乘大數法
[bc_number,dummy] = size( gBC1 ) ;
for ibc=1:1:bc_number
n = gBC1(ibc, 1 ) ;
d = gBC1(ibc, 2 ) ;
m = (n-1)*3 + d ;
f(m) = gBC1(ibc, 3)* gK(m,m) * 1e15 ;
gK(m,m) = gK(m,m) * 1e15 ;
end
% step 6. 求解方程組，得到節點位移向量
gDelta = gK \ f ;
return
function DisplayResults
% 顯示計算結果
% 輸入參數：
% 無
% 返回值：
% 無
global gNode gElement gMaterial gBC1 gNF gDF gK gDelta
fprintf( '節點位移\n' ) ;
fprintf( ' 節點號 x方向位移 y方向位移 轉角\n' ) ;
[node_number,dummy] = size( gNode ) ;
for i=1:node_number
fprintf( '%6d %16.8e %16.8e %16.8e\n',...
i, gDelta((i-1)*3+1), gDelta((i-1)*3+2), gDelta((i-1)*3+3) ) ;
end
fprintf( '\n\n節點力\n' ) ;
fprintf( ' 軸力 剪力 彎矩\n' ) ;
[element_number, dummy] = size( gElement ) ;
for ie = 1:element_number
enf = ElementNodeForce( ie ) ;
fprintf( '單元號%6d 節點號%6d %16.8e %16.8e %16.8e\n', ...
ie, gElement(ie,1), enf(1), enf(2), enf(3) ) ;
fprintf( ' 節點號%6d %16.8e %16.8e %16.8e\n', ...
gElement(ie,2), enf(4), enf(5), enf(6) ) ;
end
return
function k = StiffnessMatrix( ie, icoord )
% 計算單元剛度矩陣
% 輸入參數:
% ie ------- 單元號
% icoord -- 坐標系參數，可以是下面兩個之一
% 1 ---- 整體坐標系
% 2 ---- 局部坐標系
% 返回值:
% k ---- 根據icoord的值，相應坐標系下的剛度矩陣
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] ;
if icoord == 1
T = TransformMatrix( ie ) ;
k = T*k*transpose(T) ;
end
return
function AssembleStiffnessMatrix( ie, k )
% 把單元剛度矩陣集成到整體剛度矩陣
% 輸入參數:
% ie --- 單元號
% k --- 單元剛度矩陣
% 返回值:
% 無
global gElement gK
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) + k(m,n) ;
end
end
end
end
return
function enf = EquivalentNodeForce( ie, p1, p2, idof )
% 計算線性分布荷載的等效節點力
% 輸入參數:
% ie ----- 單元號
% p1 ----- 第一個節點上的分布力集度值
% p2 ----- 第二個節點上的分布力集度值
% idof --- 分布力的種類，它可以是下面幾種
% 1 --- 分布軸向力
% 2 --- 分布橫向力
% 3 --- 分布彎矩
% 返回值:
% enf ----- 整體坐標系下等效節點力向量
global gElement gNode
enf = zeros( 6, 1 ) ; % 定義 6x1 的等效節點力向量
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 ) ;
switch idof
case 1 % 分布軸向力
enf( 1 ) = (2*p1+p2)*L/6 ;
enf( 4 ) = (p1+2*p2)*L/6 ;
case 2 % 分布橫向力
enf( 2 ) = (7*p1+3*p2)*L/20 ;
enf( 3 ) = (3*p1+2*p2)*L^2/60 ;
enf( 5 ) = (3*p1+7*p2)*L/20 ;
enf( 6 ) = -(2*p1+3*p2)*L^2/60 ;
case 3 % 分布彎矩
enf( 2 ) = -(p1+p2)/2 ;
enf( 3 ) = (p1-p2)*L/12 ;
enf( 5 ) = (p1+p2)/2 ;
enf( 6 ) = -(p1-p2)*L/12 ;
otherwise
disp( sprintf( '分布力的種類錯誤，單元號:%d',ie ) ) ;
end
T = TransformMatrix( ie ) ; % 計算單元的轉換矩陣
enf = T * enf ; % 把等效節點力轉換到整體坐標下
return
function T = TransformMatrix( ie )
% 計算單元的坐標轉換矩陣( 局部坐標 -> 整體坐標 )
% 輸入參數
% ie ----- 節點號
% 返回值
% T ------- 從局部坐標到整體坐標的坐標轉換矩陣
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 enf = ElementNodeForce( ie )
% 計算單元的節點力
% 輸入參數
% ie ----- 節點號
% 返回值
% enf ----- 單元局部坐標系下的節點力
global gElement gNode gDelta gDF
i = gElement( ie, 1 ) ;
j = gElement( ie, 2 ) ;
de = zeros( 6, 1 ) ;
de( 1:3 ) = gDelta( (i-1)*3+1:(i-1)*3+3 ) ;
de( 4:6 ) = gDelta( (j-1)*3+1:(j-1)*3+3 ) ;
k = StiffnessMatrix( ie, 1 ) ;
enf = k * de ;
[df_number, dummy] = size( gDF ) ;
for idf = 1:1:df_number
if ie == gDF( idf, 1 )
enf = enf - EquivalentNodeForce( gDF(idf,1), ...
gDF(idf, 2), gDF( idf, 3), gDF( idf, 4 ) ) ;
break ;
end
end
T = TransformMatrix( ie ) ;
enf = transpose( T ) * enf ;
return?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -