?? exam7_1.m
字號(hào):
function exam7_1
% 本程序?yàn)榈谄哒碌牡谝粋€(gè)算例，采用薄板矩形單元計(jì)算矩形板的彎曲
% 輸入?yún)?shù)：
% 無
% 定義板的尺寸
a = 1 ; % 板x方向?qū)挾? b = 1 ; % 板y方向長度
h = 0.01 ; % 板z方向厚度
% 定義網(wǎng)格密度
m = 16 ; % x方向單元數(shù)目
n = 16 ; % y方向單元數(shù)目
% 定義材料性質(zhì)
E = 2.1e11 ; % 彈性模量
mu = 0.3 ; % 泊松比
% 荷載參數(shù)
load = '均布力' ; % 也可以是 集中力
%load = '集中力' ;
p = 1e6 ; % 載荷大小
% 邊界條件
%bc = '簡支' ; % 也可以是 固支
bc = '固支' ;
% 輸出文件名稱
file_out = 'exam7_1.mat' ;
% 檢查輸出文件是否存在
file_prompt = 0 ;
if exist( file_out ) ~= 0 & file_prompt == 1
answer = input( sprintf( '文件 %s 已經(jīng)存在，是否覆蓋? ( Yes / [No] ): ', file_out ), 's' ) ;
if length( answer ) == 0
answer = 'no' ;
end
answer = lower( answer ) ;
if answer ~= 'y' | answer ~= 'yes'
disp( sprintf( '請(qǐng)使用另外的文件名，或備份已有的文件' ) ) ;
disp( sprintf( '程序終止' ) );
return ;
end
end
FemModel(a,b,h,m,n,E,mu,load,p,bc) ; % 定義有限元模型
SolveModel ; % 求解有限元模型
SaveResults( file_out ) ; % 保存計(jì)算結(jié)果
DisplayResults(a,b,h,m,n,E,mu,load,p,bc) ; % 把計(jì)算結(jié)果顯示出來
return
function FemModel(a,b,h,m,n,E,mu,load,p,bc)
% 定義有限元模型
% 輸入?yún)?shù)：
% a --------- 板x方向?qū)挾?% b --------- 板y方向長度
% h --------- 板z方向厚度
% m --------- x方向單元數(shù)目
% n --------- y方向單元數(shù)目
% E --------- 彈性模量
% mu --------- 泊松比
% load ------- 載荷類型
% p --------- 載荷大小
% bc --------- 邊界條件
% 返回值：
% 無
global gNode gElement gMaterial gBC1 gDF gNF gK
% 計(jì)算結(jié)點(diǎn)坐標(biāo)
dx = a / m / 2 ; % 取板的1/4進(jìn)行分析
dy = b / n / 2 ;
gNode = zeros( (m+1)*(n+1), 2 ) ;
for i=1:n+1
for j=1:m+1
gNode( (i-1)*(m+1)+j, : ) = [dx*(j-1),dy*(i-1)] ;
end
end
% 定義單元
gElement = zeros( n*n, 5 ) ;
for i=1:n
for j=1:m
gElement( (i-1)*m+j, 1:4) = [ (i-1)*(m+1)+j, ...
(i-1)*(m+1)+j+1, ...
i*(m+1)+j+1,...
i*(m+1)+j ] ;
end
end
gElement( :, 5 ) = 1 ;
% 定義材料
gMaterial = [ E, mu, h] ;
% 確定邊界條件
if bc == '簡支'
gBC1 = zeros( (n+1)*2+m*2+1, 3 ) ;
for i=1:m+1
gBC1( i, : ) = [i, 2, 0] ; % x=0的邊界上繞x軸的轉(zhuǎn)角等于零
end
for i=1:n+1
gBC1( (m+1)+i, : ) = [(i-1)*(m+1)+1, 3, 0] ; % y=0的邊界上繞y軸的轉(zhuǎn)角等于零
end
for i=1:n+1
gBC1( m+n+2+i, : ) = [i*(m+1), 1, 0 ] ; % x=a/2的邊界上撓度為零
end
for i=1:m
gBC1( (n+1)*2+m+1+i, : ) = [n*(m+1)+i, 1, 0] ; % y=b/2的邊界上撓度為零
end
elseif bc == '固支'
gBC1 = zeros( m*4+(n+1)*3+n, 3 ) ;
for i=1:m
gBC1( i, : ) = [i, 2, 0] ; % x=0的邊界上繞x軸的轉(zhuǎn)角等于零
end
for i=1:n
gBC1( m+i, : ) = [(i-1)*(m+1)+1, 3, 0] ; % y=0的邊界上繞y軸的轉(zhuǎn)角等于零
end
for i=1:n+1
gBC1( m+n+i, : ) = [i*(m+1), 1, 0 ] ; % x=a/2的邊界上撓度為零
end
for i=1:n+1
gBC1( m+n+n+1+i, : ) = [i*(m+1), 2, 0 ] ; % x=a/2的邊界上繞x軸的轉(zhuǎn)角為零
end
for i=1:n+1
gBC1( m+n+(n+1)*2+i, : ) = [i*(m+1), 3, 0 ] ; % x=a/2的邊界上繞y軸的轉(zhuǎn)角為零
end
for i=1:m
gBC1( m+n+(n+1)*3+i, : ) = [n*(m+1)+i, 1, 0] ; % y=b/2的邊界上撓度為零
end
for i=1:m
gBC1( m*2+n+(n+1)*3+i, : ) = [n*(m+1)+i, 2, 0] ; % y=b/2的邊界上繞x軸的轉(zhuǎn)角為零
end
for i=1:m
gBC1( m*3+n+(n+1)*3+i, : ) = [n*(m+1)+i, 3, 0] ; % y=b/2的邊界上繞y軸的轉(zhuǎn)角為零
end
end
% 確定載荷
if load == '均布力'
gNF = [] ;
gDF = gElement ;
gDF(:,5) = p ;
elseif load == '集中力'
gDF = [] ;
gNF = [1, 1, p] ;
end
return
function SolveModel
% 求解有限元模型
% 輸入?yún)?shù)：
% 無
% 返回值：
% 無
% 說明：
% 該函數(shù)求解有限元模型，過程如下
% 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 gDF gNF gK gDelta gElementMoment gNodeMoment
% 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
disp( sprintf( '計(jì)算單元?jiǎng)偠染仃嚥⒓桑?dāng)前單元: %d', ie ) ) ;
k = StiffnessMatrix( ie ) ;
AssembleStiffnessMatrix( ie, k ) ;
end
% step3. 計(jì)算分布?jí)毫Φ牡刃Ч?jié)點(diǎn)力，并集成到整體節(jié)點(diǎn)力向量中
[df_number, dummy] = size( gDF ) ;
for idf = 1:1:df_number
edf = EquivalentDistPressure( gDF(idf,1:4), gDF(idf,5) ) ;
node = gDF( idf, 1:4 ) ;
f( node*3-2 ) = f( node*3-2 ) + edf( 1:3:10 ) ;
f( node*3-1 ) = f( node*3-1 ) + edf( 2:3:11 ) ;
f( node*3 ) = f( node*3 ) + edf( 3:3:12 ) ;
end
% step4. 計(jì)算分布?jí)毫Φ牡刃Ч?jié)點(diǎn)力，并集成到整體節(jié)點(diǎn)力向量中
[nf_number,dummy] = size( gNF ) ;
for inf = 1:nf_number
node = gNF( inf, 1 ) ;
dim = gNF( inf, 2 ) ;
f( (node-1)*3+dim ) = f( (node-1)*3+dim ) + gNF( inf, 3 ) ;
end
% step5. 處理第一類約束條件，修改剛度矩陣和節(jié)點(diǎn)力向量。采用乘大數(shù)法
[bc1_number,dummy] = size( gBC1 ) ;
for ibc=1:1:bc1_number
n = gBC1(ibc, 1 ) ;
d = gBC1(ibc, 2 ) ;
m = (n-1)*3 + d ;
f(m) = gBC1(ibc, 3)* gK(m,m) * 1e10 ;
gK(m,m) = gK(m,m) * 1e10 ;
end
% step5. 求解方程組，得到節(jié)點(diǎn)位移向量
gDelta = gK \ f ;
% step6. 計(jì)算單元彎矩
gElementMoment = zeros( element_number, 4, 3 ) ;
for ie=1:element_number
disp( sprintf( '計(jì)算單元彎矩，當(dāng)前單元: %d', ie ) ) ;
em = ElementMoment( ie ) ;
gElementMoment( ie, :, : ) = em ;
end
% step7. 計(jì)算節(jié)點(diǎn)彎矩(采用繞節(jié)點(diǎn)平均)
gNodeMoment = zeros( node_number, 3 ) ;
nodenum = zeros( node_number,1 ) ;
for i=1:node_number
disp( sprintf( '計(jì)算節(jié)點(diǎn)彎矩，當(dāng)前結(jié)點(diǎn): %d', i ) ) ;
for ie=1:element_number
index = find( gElement( ie, 1:4 ) == i ) ;
if ~isempty(index)
nodenum(i) = nodenum(i) + 1 ;
gNodeMoment(i,1) = gNodeMoment(i,1) + gElementMoment( ie, index, 1 ) ;
gNodeMoment(i,2) = gNodeMoment(i,2) + gElementMoment( ie, index, 2 ) ;
gNodeMoment(i,3) = gNodeMoment(i,3) + gElementMoment( ie, index, 3 ) ;
end
end
end
gNodeMoment(:,1) = gNodeMoment(:,1) ./ nodenum ;
gNodeMoment(:,2) = gNodeMoment(:,2) ./ nodenum ;
gNodeMoment(:,3) = gNodeMoment(:,3) ./ nodenum ;
return
function k = StiffnessMatrix( ie )
% 計(jì)算平面應(yīng)變等參數(shù)單元的剛度矩陣
% 輸入?yún)?shù)：
% ie -- 單元號(hào)
% 返回值：
% k -- 單元?jiǎng)偠染仃?
global gNode gElement gMaterial
k = zeros( 12, 12 ) ;
E = gMaterial( gElement( ie, 5 ), 1 ) ;
mu = gMaterial( gElement( ie, 5 ), 2 ) ;
h = gMaterial( gElement( ie, 5 ), 3 ) ;
D = E * h^3 / 12 / ( 1 - mu^2 ) ;
x = gNode( gElement( ie, : ), 1 ) ;
y = gNode( gElement( ie, : ), 2 ) ;
a = abs( x(3) - x(1) ) / 2 ;
b = abs( y(3) - y(1) ) / 2 ;
ab2 = a^2/b^2 ;
ba2 = b^2/a^2 ;
H = D/60/a/b ;
xi = [ -1, 1, 1, -1 ] ;
eta = [ -1, -1, 1, 1 ] ;
for i=1:4
for j=1:4
xi0 = xi(i)*xi(j) ;
eta0 = eta(i)*eta(j) ;
a11 = 3*H*(15*(ba2*xi0+ab2*eta0) ...
+(14-4*mu+5*ba2+5*ab2)*xi0*eta0) ;
a12 = -3*H*b*((2+3*mu+5*ab2)*xi0*eta(i) ...
+ 15*ab2*eta(i)+5*mu*xi0*eta(j)) ;
a13 = 3*H*a*((2+3*mu+5*ba2)*xi(i)*eta0 ...
+ 15*ba2*xi(i)+5*mu*xi(j)*eta0) ;
a21 = -3*H*b*((2+3*mu+5*ab2)*xi0*eta(j) ...
+ 15*ab2*eta(j)+5*mu*xi0*eta(i)) ;
a22 = H*b^2*(2*(1-mu)*xi0*(3+5*eta0)+...
+ 5*ab2*(3+xi0)*(3+eta0)) ;
a23 = -15*H*mu*a*b*(xi(i)+xi(j))*(eta(i)+eta(j)) ;
a31 = 3*H*a*((2+3*mu+5*ba2)*xi(j)*eta0 ...
+ 15*ba2*xi(j)+5*mu*xi(i)*eta0) ;
a32 = -15*H*mu*a*b*(xi(i)+xi(j))*(eta(i)+eta(j)) ;
a33 = H*a^2*(2*(1-mu)*eta0*(3+5*xi0)+...
+ 5*ba2*(3+xi0)*(3+eta0)) ;
k( i*3-2:i*3, j*3-2:j*3 ) = [ a11 a12 a13;
a21 a22 a23;
a31 a32 a33 ] ;
end
end
return
function AssembleStiffnessMatrix( ie, k )
% 把單元?jiǎng)偠染仃嚰傻秸w剛度矩陣
% 輸入?yún)?shù):
% ie --- 單元號(hào)
% k --- 單元?jiǎng)偠染仃?% 返回值:
% 無
global gElement gK
for i=1:1:4
for j=1:1:4
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 edf = EquivalentDistPressure( node, press )
% 計(jì)算分布?jí)毫Φ牡刃Ч?jié)點(diǎn)力
% 輸入?yún)?shù):
% node ---------- 結(jié)點(diǎn)號(hào)
% press --------- 跟結(jié)點(diǎn)號(hào)對(duì)應(yīng)的壓力值
% 返回值:
% edf ----------- 等效節(jié)點(diǎn)力向量
global gNode
edf = zeros( 12, 1 ) ;
x = gNode( node, 1 ) ;
y = gNode( node, 2 ) ;
a = abs( x(3)-x(1) ) / 2 ;
b = abs( y(3)-y(1) ) / 2 ;
xi = [-1; 1; 1; -1] ;
eta = [-1;-1; 1; 1] ;
for i=1:4
edf(i*3-2:i*3) = [press*a*b; -press*a*b^2/3*eta(i); press*a^2*b/3*xi(i)] ;
end
return
function em = ElementMoment( ie )
% 計(jì)算單元的結(jié)點(diǎn)彎矩
% 輸入?yún)?shù):
% ie ----- 單元號(hào)
% 返回值:
% em ----- 單元結(jié)點(diǎn)彎矩
global gElement gDelta gMaterial
E = gMaterial( gElement( ie, 5 ), 1 ) ;
mu = gMaterial( gElement( ie, 5 ), 2 ) ;
h = gMaterial( gElement( ie, 5 ), 3 ) ;
em = zeros( 4, 3 ) ;
node = gElement( ie, 1:4 ) ;
delta = zeros( 12, 1 ) ;
delta(1:3:10) = gDelta( (node-1)*3+1 ) ;
delta(2:3:11) = gDelta( (node-1)*3+2 ) ;
delta(3:3:12) = gDelta( (node-1)*3+3 ) ;
xi = [ -1, 1, 1, -1 ] ;
eta = [ -1, -1, 1, 1 ] ;
D = E/(1-mu^2)*[1 mu 0
mu 1 0
0 0 (1-mu)/2 ] ;
for i=1:4
B = MatrixB( ie, xi(i), eta(i) ) ;
em( i, : ) = transpose(h^3/12*D*B*delta) ;
end
return
function B = MatrixB( ie, x, e )
% 計(jì)算單元的應(yīng)變矩陣B
% 輸入?yún)?shù):
% ie ----- 單元號(hào)
% 返回值:
% B ------ 單元應(yīng)變矩陣
global gNode gElement
B = zeros( 3, 12 ) ;
xy = gNode( gElement(ie, 1:4), : ) ;
a = abs(xy(3,1)-xy(1,1))/2 ;
b = abs(xy(3,2)-xy(1,2))/2 ;
xi = [ -1, 1, 1, -1 ] ;
eta = [ -1, -1, 1, 1 ] ;
for i=1:4
x0 = xi(i)*x ;
e0 = eta(i)*e ;
B(:,(i-1)*3+1:(i-1)*3+3) = [ 3*b/a*x0*(1+e0) 0 b*xi(i)*(1+3*x0)*(1+e0)
3*a/b*e0*(1+x0) -a*eta(i)*(1+x0)*(1+3*e0) 0
xi(i)*eta(i)*(3*x^2+3*e^2-4) -b*xi(i)*(3*e^2+2*e0-1) a*eta(i)*(3*x^2+2*x0-1) ] ;
end
B = B / a / b / 4 ;
return
function SaveResults( file_out )
% 保存計(jì)算結(jié)果
% 輸入?yún)?shù)：
% 無
% 返回值：
% 無
global gNode gElement gMaterial gBC1 gDF gDelta gNF gK
save( file_out, 'gNode', 'gElement', 'gMaterial', 'gBC1', ...
'gDF', 'gDelta', 'gNF', 'gK' ) ;
return
function DisplayResults(a,b,h,m,n,E,mu,load,p,bc)
% 顯示計(jì)算結(jié)果，并與解析解結(jié)果比較
% 輸入?yún)?shù)：
% a --------- 板x方向?qū)挾?% b --------- 板y方向長度
% h --------- 板z方向厚度
% m --------- x方向單元數(shù)目
% n --------- y方向單元數(shù)目
% E --------- 彈性模量
% mu --------- 泊松比
% load ------- 載荷類型
% p --------- 載荷大小
% bc --------- 邊界條件
% 返回值：
% 無
global gDelta gNodeMoment
wmax = gDelta( 1, 1 ) ; % 中心撓度
D = E*h^3/12/(1-mu^2) ; % 抗彎剛度
mx = gNodeMoment( 1, 1 ) ; % 中心內(nèi)力矩(Mx)
my = gNodeMoment( 1, 2 ) ; % 中心內(nèi)力矩(My)
mxy = gNodeMoment( 1, 3 ) ; % 中心內(nèi)力矩(Mxy)
fprintf( '板的尺寸a×b×h = %g×%g×%g\n', a, b, h ) ;
fprintf( '彈性模量 = %g\n', E ) ;
fprintf( '泊松比 = %g\n', mu ) ;
fprintf( '載荷類型 = %s\n', load ) ;
fprintf( '載荷大小 = %g\n', p ) ;
fprintf( '邊界條件 = %s\n', bc ) ;
fprintf( '網(wǎng)格密度 m×n = %d×%d\n', m, n ) ;
if load == '均布力'
alpha = wmax * D / p / a^4 ;
beta = mx/p/a^2 ;
beta1 = my/p/a^2 ;
elseif load == '集中力'
alpha = wmax * D / p / 4 / a^2 ;
beta = mx/p/4 ;
beta1 = my/p/3 ;
end
fprintf( '撓度系數(shù) alpha = %.7f\n', alpha ) ;
fprintf( '彎矩系數(shù) beta = %.7f\n', beta ) ;
fprintf( '彎矩系數(shù) beta1 = %.7f\n', beta1 ) ;
% fid = fopen( 'exam7_1.out', 'a' ) ;
% fprintf( fid, '%.1f %15.7f %15.7f\n', b/a, alpha, beta ) ;
% fclose( fid ) ;
return?? 快捷鍵說明
復(fù)制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號(hào)
Ctrl + =
減小字號(hào)
Ctrl + -