?? exam3_1.m
字號(hào):
function exam3_1
% 本程序?yàn)榈谌碌牡谝粋€(gè)算例，采用平面梁?jiǎn)卧?jì)算平面剛架的變形和內(nèi)力
% 輸入?yún)?shù)： 無(wú)
% 輸出結(jié)果： 節(jié)點(diǎn)位移和單元節(jié)點(diǎn)力
PlaneFrameModel ; % 定義有限元模型
SolveModel ; % 求解有限元模型
DisplayResults ; % 顯示計(jì)算結(jié)果
return ;
function PlaneFrameModel
% 定義平面桿系的有限元模型
% 輸入?yún)?shù)：
% 無(wú)
% 返回值：
% 無(wú)
% 說(shuō)明：
% 該函數(shù)定義平面桿系的有限元模型數(shù)據(jù)：
% gNode ------- 節(jié)點(diǎn)定義
% gElement ---- 單元定義
% gMaterial --- 材料定義，包括彈性模量，梁的截面積和梁的抗彎慣性矩
% gBC1 -------- 約束條件
% gNF --------- 集中力
% gDF --------- 分布力
global gNode gElement gMaterial gBC1 gNF gDF
% 節(jié)點(diǎn)坐標(biāo)
% x y
gNode = [0.0, 0.0 % 節(jié)點(diǎn) 1
0.0, 4.0 % 節(jié)點(diǎn) 2
3.0, 0.0 % 節(jié)點(diǎn) 3
3.0, 4.0 % 節(jié)點(diǎn) 4
4.5, 4.0 % 節(jié)點(diǎn) 5
6.0, 0.0 % 節(jié)點(diǎn) 6
6.0, 4.0 ] ; % 節(jié)點(diǎn) 7
% 單元定義
% 節(jié)點(diǎn)1 節(jié)點(diǎn)2 材料號(hào)
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
% 材料性質(zhì)
% 彈性模量 抗彎慣性矩 截面積
gMaterial = [2.1e11, 2.0e-4, 1.0e-2] ; % 材料 1
% 第一類(lèi)約束條件
% 節(jié)點(diǎn)號(hào) 自由度號(hào) 約束值
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] ;
% 集中力
% 節(jié)點(diǎn)號(hào) 自由度號(hào) 集中力值
gNF = [ 5, 2, -80e3] ;
% 分布載荷（線(xiàn)性分布）
% 單元號(hào) 節(jié)點(diǎn)1載荷值 節(jié)點(diǎn)2載荷值 自由度號(hào)
gDF = [ 1 -30e3 0 2
2 -15e3 -15e3 2 ] ;
return
function SolveModel
% 求解有限元模型
% 輸入?yún)?shù)：
% 無(wú)
% 返回值：
% 無(wú)
% 說(shuō)明：
% 該函數(shù)求解有限元模型，過(guò)程如下
% 1. 計(jì)算單元?jiǎng)偠染仃嚕烧w剛度矩陣
% 2. 計(jì)算單元的等效節(jié)點(diǎn)力，集成整體節(jié)點(diǎn)力向量
% 3. 處理約束條件，修改整體剛度矩陣和節(jié)點(diǎn)力向量
% 4. 求解方程組，得到整體節(jié)點(diǎn)位移向量
global gNode gElement gMaterial gBC1 gNF gDF gK gDelta
% step1. 定義整體剛度矩陣和節(jié)點(diǎn)力向量
[node_number,dummy] = size( gNode ) ;
gK = sparse( node_number * 3, node_number * 3 ) ;
f = sparse( node_number * 3, 1 ) ;
% step2. 計(jì)算單元?jiǎng)偠染仃嚕⒓傻秸w剛度矩陣中
[element_number,dummy] = size( gElement ) ;
for ie=1:1:element_number
k = StiffnessMatrix( ie, 1 ) ;
AssembleStiffnessMatrix( ie, k ) ;
end
% step3. 把集中力直接集成到整體節(jié)點(diǎn)力向量中
[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. 計(jì)算分布力的等效節(jié)點(diǎn)力，不集成到整體節(jié)點(diǎn)力向量中
[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. 處理約束條件，修改剛度矩陣和節(jié)點(diǎn)力向量。采用乘大數(shù)法
[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. 求解方程組，得到節(jié)點(diǎn)位移向量
gDelta = gK \ f ;
return
function DisplayResults
% 顯示計(jì)算結(jié)果
% 輸入?yún)?shù)：
% 無(wú)
% 返回值：
% 無(wú)
global gNode gElement gMaterial gBC1 gNF gDF gK gDelta
fprintf( '節(jié)點(diǎn)位移\n' ) ;
fprintf( ' 節(jié)點(diǎn)號(hào) x方向位移 y方向位移 轉(zhuǎn)角\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節(jié)點(diǎn)力\n' ) ;
fprintf( ' 軸力 剪力 彎矩\n' ) ;
[element_number, dummy] = size( gElement ) ;
for ie = 1:element_number
enf = ElementNodeForce( ie ) ;
fprintf( '單元號(hào)%6d 節(jié)點(diǎn)號(hào)%6d %16.8e %16.8e %16.8e\n', ...
ie, gElement(ie,1), enf(1), enf(2), enf(3) ) ;
fprintf( ' 節(jié)點(diǎn)號(hào)%6d %16.8e %16.8e %16.8e\n', ...
gElement(ie,2), enf(4), enf(5), enf(6) ) ;
end
return
function k = StiffnessMatrix( ie, icoord )
% 計(jì)算單元?jiǎng)偠染仃?% 輸入?yún)?shù):
% ie ------- 單元號(hào)
% icoord -- 坐標(biāo)系參數(shù)，可以是下面兩個(gè)之一
% 1 ---- 整體坐標(biāo)系
% 2 ---- 局部坐標(biāo)系
% 返回值:
% k ---- 根據(jù)icoord的值，相應(yīng)坐標(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] ;
if icoord == 1
T = TransformMatrix( ie ) ;
k = T*k*transpose(T) ;
end
return
function AssembleStiffnessMatrix( ie, k )
% 把單元?jiǎng)偠染仃嚰傻秸w剛度矩陣
% 輸入?yún)?shù):
% ie --- 單元號(hào)
% k --- 單元?jiǎng)偠染仃?% 返回值:
% 無(wú)
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 )
% 計(jì)算線(xiàn)性分布荷載的等效節(jié)點(diǎn)力
% 輸入?yún)?shù):
% ie ----- 單元號(hào)
% p1 ----- 第一個(gè)節(jié)點(diǎn)上的分布力集度值
% p2 ----- 第二個(gè)節(jié)點(diǎn)上的分布力集度值
% idof --- 分布力的種類(lèi)，它可以是下面幾種
% 1 --- 分布軸向力
% 2 --- 分布橫向力
% 3 --- 分布彎矩
% 返回值:
% enf ----- 整體坐標(biāo)系下等效節(jié)點(diǎn)力向量
global gElement gNode
enf = zeros( 6, 1 ) ; % 定義 6x1 的等效節(jié)點(diǎn)力向量
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( '分布力的種類(lèi)錯(cuò)誤，單元號(hào):%d',ie ) ) ;
end
T = TransformMatrix( ie ) ; % 計(jì)算單元的轉(zhuǎn)換矩陣
enf = T * enf ; % 把等效節(jié)點(diǎn)力轉(zhuǎn)換到整體坐標(biāo)下
return
function T = TransformMatrix( ie )
% 計(jì)算單元的坐標(biāo)轉(zhuǎn)換矩陣( 局部坐標(biāo) -> 整體坐標(biāo) )
% 輸入?yún)?shù)
% ie ----- 節(jié)點(diǎn)號(hào)
% 返回值
% 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 enf = ElementNodeForce( ie )
% 計(jì)算單元的節(jié)點(diǎn)力
% 輸入?yún)?shù)
% ie ----- 節(jié)點(diǎn)號(hào)
% 返回值
% enf ----- 單元局部坐標(biāo)系下的節(jié)點(diǎn)力
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?? 快捷鍵說(shuō)明
復(fù)制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號(hào)
Ctrl + =
減小字號(hào)
Ctrl + -