📄 exam6_1.m
字号:
B = zeros( 3, 16 ) ;
for i=1:1:8
B(1:3,(2*i-1):2*i) = [ N_x(i) 0
0 N_y(i)
N_y(i), N_x(i)];
end
return
function [N_x, N_y] = N_xy( ie, xi, eta )
% 计算形函数对整体坐标的导数
% 输入参数：
% ie --------- 单元号
% xi,eta ----- 局部坐标
% 返回值：
% N_x ------- 在局部坐标处的形函数对x坐标的导数
% N_y ------- 在局部坐标处的形函数对y坐标的导数
J = Jacobi( ie, xi, eta ) ;
[N_xi,N_eta] = N_xieta( ie, xi, eta ) ;
A=inv(J)*[N_xi;N_eta] ;
N_x = A(1,:) ;
N_y = A(2,:) ;
return
function [N_xi, N_eta] = N_xieta( ie, xi, eta )
% 计算形函数对局部坐标的导数
% 输入参数：
% ie --------- 单元号
% xi,eta ----- 局部坐标
% 返回值：
% N_xi ------- 在局部坐标处的形函数对xi坐标的导数
% N_eta ------- 在局部坐标处的形函数对eta坐标的导数
x = [ -1, 1, 1, -1 ] ;
e = [ -1, -1, 1, 1 ] ;
N_xi = zeros( 1, 8 ) ;
N_eta = zeros( 1, 8 ) ;
N_xi( 5 ) = xi*(eta-1) ;
N_eta( 5 ) = 0.5*(xi^2-1) ;
N_xi( 6 ) = 0.5*(1-eta^2) ;
N_eta( 6 ) = -eta*(xi+1) ;
N_xi( 7 ) = -xi*(eta+1) ;
N_eta( 7 ) = 0.5*(1-xi^2) ;
N_xi( 8 ) = 0.5*(eta^2-1) ;
N_eta( 8 ) = eta*(xi-1) ;
N_xi(1) = x(1)*(1+e(1)*eta)/4 - 0.5*( N_xi(5) + N_xi(8) );
N_eta(1) = e(1)*(1+x(1)*xi)/4 - 0.5*( N_eta(5) + N_eta(8) ) ;
N_xi(2) = x(2)*(1+e(2)*eta)/4 - 0.5*( N_xi(5) + N_xi(6) );
N_eta(2) = e(2)*(1+x(2)*xi)/4 - 0.5*( N_eta(5) + N_eta(6) ) ;
N_xi(3) = x(3)*(1+e(3)*eta)/4 - 0.5*( N_xi(6) + N_xi(7) );
N_eta(3) = e(3)*(1+x(3)*xi)/4 - 0.5*( N_eta(6) + N_eta(7) ) ;
N_xi(4) = x(4)*(1+e(4)*eta)/4 - 0.5*( N_xi(7) + N_xi(8) );
N_eta(4) = e(4)*(1+x(4)*xi)/4 - 0.5*( N_eta(7) + N_eta(8) ) ;
return
function J = Jacobi( ie, xi, eta )
% 计算雅克比矩阵
% 输入参数：
% ie --------- 单元号
% xi,eta ----- 局部坐标
% 返回值：
% J ------- 在局部坐标(xi,eta)处的雅克比矩阵
global gNode gElement
x = gNode(gElement(ie,1:8),1) ;
y = gNode(gElement(ie,1:8),2) ;
[N_xi,N_eta] = N_xieta( ie, xi, eta ) ;
x_xi = N_xi * x ;
x_eta = N_eta * x ;
y_xi = N_xi * y ;
y_eta = N_eta * y ;
J = [ x_xi, y_xi; x_eta, y_eta ];
return
function AssembleStiffnessMatrix( ie, k )
% 把单元刚度矩阵集成到整体刚度矩阵
% 输入参数:
% ie --- 单元号
% k --- 单元刚度矩阵
% 返回值:
% 无
global gElement gK
for i=1:1:8
for j=1:1:8
for p=1:1:2
for q=1:1:2
m = (i-1)*2+p ;
n = (j-1)*2+q ;
M = (gElement(ie,i)-1)*2+p ;
N = (gElement(ie,j)-1)*2+q ;
gK(M,N) = gK(M,N) + k(m,n) ;
end
end
end
end
return
function edf = EquivalentDistPressure( node, pressure )
% 计算分布压力的等效节点力
% 输入参数:
% node ---------- 结点号
% pressure ------ 跟结点号对应的压力值
% 返回值:
% edf ----------- 等效节点力向量
global gNode
x = gNode( node, 1 ) ;
y = gNode( node, 2 ) ;
xi = x(1) ; xj = x(2) ; xk = x(3) ;
yi = y(1) ; yj = y(2) ; yk = y(3) ;
sigi = pressure(1) ;
sigj = pressure(2) ;
sigk = pressure(3) ;
X1 = -10*xi-2*xj+12*xk ;
Y1 = -10*yi-2*yj+12*yk ;
X2 = xi - xj ;
Y2 = yi - yj ;
X3 = -6*xi-2*xj+8*xk ;
Y3 = -6*yi-2*yj+8*yk ;
X4 = 2*xi+10*xj-12*xk ;
Y4 = 2*yi+10*yj-12*yk ;
X5 = 2*xi+6*xj-8*xk ;
Y5 = 2*yi+6*yj-8*yk ;
X6 = -16*xi+16*xj ;
Y6 = -16*yi+16*yj ;
XY = [ Y1 Y2 Y3
-X1 -X2 -X3
Y2 Y4 Y5
-X2 -X4 -X5
Y3 Y5 Y6
-X3 -X5 -X6 ] ;
sig = [sigi; sigj; sigk ] ;
edf = 1/30*XY*sig ;
return
function evf = EquivalentVolumeForce( ie )
% 计算离心力的等效节点力
% 输入参数:
% ie ---------- 单元号
% 返回值:
% evf --------- 等效节点力向量
global gNode gElement gMaterial
evf = zeros( 16, 1 ) ;
omega = 2094.0 ; % 旋转角速度
ro = gMaterial( gElement( ie, 9 ), 3 ) ;
x = gNode( gElement(ie,1:8), 1 ) ;
y = gNode( gElement(ie,1:8), 2 ) ;
xi = [-0.774596669241483,0.0, 0.774596669241483] ;
w = [ 0.555555555555556,0.888888888888889,0.555555555555556] ;
for i=1:length(xi)
for j=1:length(xi)
J = Jacobi( ie, xi(i), xi(j) ) ;
detJ = det(J);
N = ShapeFunction( xi(i), xi(j) ) ;
evf(1:2:15) = evf(1:2:15) + N'*N*x*detJ*w(i)*w(j) ;
evf(2:2:16) = evf(2:2:16) + N'*N*y*detJ*w(i)*w(j) ;
end
end
evf = evf * ro * omega^2 ;
return
function N = ShapeFunction( xi, eta )
% 计算形函数的值
% 输入参数:
% ie ---------- 单元号
% xi, eta ----- 单元内局部坐标
% 返回值:
% N ----------- 形函数的值
N5 = ( eta - 1 ) * ( xi^2 - 1 ) / 2 ;
N6 = ( xi + 1 ) * ( 1 - eta^2 ) / 2 ;
N7 = ( eta + 1 ) * ( 1 - xi^2 ) / 2 ;
N8 = ( xi - 1 ) * ( eta^2 - 1 ) / 2 ;
N1 = ( 1 - xi ) * ( 1 - eta ) / 4 - 0.5 * ( N8 + N5 ) ;
N2 = ( 1 + xi ) * ( 1 - eta ) / 4 - 0.5 * ( N5 + N6 ) ;
N3 = ( 1 + xi ) * ( 1 + eta ) / 4 - 0.5 * ( N6 + N7 ) ;
N4 = ( 1 - xi ) * ( 1 + eta ) / 4 - 0.5 * ( N7 + N8 ) ;
N = [ N1 N2 N3 N4 N5 N6 N7 N8 ];
return
function SaveResults( file_out )
% 保存计算结果
% 输入参数：
% 无
% 返回值：
% 无
global gNode gElement gMaterial gBC1 gBC3 gDF gK gDelta gDelta1 gElementStress
save( file_out, 'gNode', 'gElement', 'gMaterial', 'gBC1', 'gBC3', ...
'gDF', 'gDelta', 'gDelta1', 'gElementStress' ) ;
return
function DisplayResults
% 显示计算结果，并与解析解结果比较
% 输入参数：
% 无
% 返回值：
% 无
% 说明:
% 1。受内压的轴对称圆筒平面应力问题的解析解是(参考文献[1])
% 径向位移: ur = -C1/E*(1+mu)/r + 2*C2/E*(1-mu)*r
% C1 = -a^2*b^2*p/(b^2-a^2)
% C2 = p*a^2/2/(b^2-a^2)
% 径向应力: sr = a^2*p/(b^2-a^2)*(1-b^2/r^2)
% 周向应力: st = a^2*p/(b^2-a^2)*(1+b^2/r^2)
% 2。以匀角速度w旋转的圆筒平面应力问题的解析解(参考文献[2])
% 径向位移: ur = 1/E*((1-mu)*C3*r - (1+mu)*C4/r - (1-mu^2)/8*ro*w^2*r^3 )
% C3 = (3+mu)/8*ro*w^2*(b^2+a^2)
% C4 = -(3+mu)/8*ro*w^2*a^2*b^2
% 径向应力: sr = (3+mu)/8*ro*w^2*(b^2+a^2-a^2*b^2/r^2-r^2)
% 周向应力: st = (3+mu)/8*ro*w^2*(b^2+a^2+a^2*b^2/r^2-(1+3*mu)/(3+mu)*r^2)
%
% 参数说明
% E ----- 弹性模量
% mu ----- 泊松比
% a ----- 内半径
% b ----- 外半径
% r ----- 径向坐标
% p ----- 内压
% ro ----- 密度
% w ----- 旋转角速度(rad/s)
%
% 参考文献
% [1] 钱伟长 叶开元, 1956, 弹性力学, pp.230-231
% [2] Timoshenko S. P., Goodier J. N., 1970 Theory of elasticity. New York: McGraw-Hill, pp.80-83
global gNode gElement gMaterial gDelta gDelta1 gElementStress
E = gMaterial( 1, 1 ) ;
mu = gMaterial( 1, 2 ) ;
ro = gMaterial( 1, 3 ) ;
a = 0.05 ;
b = 0.10 ;
p = 120e6 ;
w = 2094 ;
r = 0.05:0.005:0.1 ;
% 计算受内压圆筒的位移和应力的解析解
C1 = -a^2*b^2*p/(b^2-a^2) ;
C2 = p*a^2/2/(b^2-a^2) ;
ur1 = -C1/E*(1+mu)./r + 2*C2/E*(1-mu)*r ;
sr1 = a^2*p/(b^2-a^2)*(1-b^2./r.^2) ;
st1 = a^2*p/(b^2-a^2)*(1+b^2./r.^2) ;
% 计算旋转圆盘的位移和应力的解析解
C3 = (3+mu)/8*ro*w^2*(b^2+a^2) ;
C4 = -(3+mu)/8*ro*w^2*a^2*b^2 ;
ur2 = 1/E*((1-mu)*C3*r - (1+mu)*C4./r - (1-mu^2)/8*ro*w^2*r.^3 ) ;
sr2 = (3+mu)/8*ro*w^2*(b^2+a^2-a^2*b^2./r.^2-r.^2) ;
st2 = (3+mu)/8*ro*w^2*(b^2+a^2+a^2*b^2./r.^2-(1+3*mu)/(3+mu)*r.^2) ;
% 读取有限元的计算结果
node = [3, 5, 8, 10, 13, 15, 18, 20, 23, 25, 28];
node_num = zeros( size(node) ) ;
ur1_fem = gDelta( node*2-1, 1 ) ;
ur2_fem = gDelta( node*2-1, 2 ) ;
[element_number,dummy] = size( gElement ) ;
sr1_fem = zeros( size( sr1 ) ) ;
st1_fem = zeros( size( st1 ) ) ;
sr2_fem = zeros( size( sr2 ) ) ;
st2_fem = zeros( size( st2 ) ) ;
for n=1:length(node)
for ie=1:element_number
index = find( gElement( ie, 1:8 ) == node( n ) ) ;
if ~isempty(index)
sr1_fem( n ) = sr1_fem( n ) + gElementStress( ie, index, 1, 1 ) ;
st1_fem( n ) = st1_fem( n ) + gElementStress( ie, index, 2, 1 ) ;
sr2_fem( n ) = sr2_fem( n ) + gElementStress( ie, index, 1, 2 ) ;
st2_fem( n ) = st2_fem( n ) + gElementStress( ie, index, 2, 2 ) ;
node_num(n) = node_num(n) + 1 ;
end
end
end
sr1_fem = sr1_fem ./ node_num ;
st1_fem = st1_fem ./ node_num ;
sr2_fem = sr2_fem ./ node_num ;
st2_fem = st2_fem ./ node_num ;
fprintf( '\n\n\n' ) ;
fprintf( ' 表1 径向位移 (mm)\n' ) ;
fprintf( '==================================================================\n' ) ;
fprintf( '结点号 坐标 承受内压 承受离心力\n' ) ;
fprintf( ' (mm) 有限元解 解析解 有限元解 解析解\n' ) ;
fprintf( '------------------------------------------------------------------\n' ) ;
for i=1:length(node)
fprintf( '%4d %4.0f %.6f %.6f %.6f %.6f\n', ...
node(i), r(i)*1000, ur1_fem(i)*1000, ur1(i)*1000, ur2_fem(i)*1000, ur2(i)*1000 ) ;
end
fprintf( '==================================================================\n\n\n' ) ;
fprintf( ' 表2 径向应力 (MPa)\n' ) ;
fprintf( '==================================================================\n' ) ;
fprintf( '结点号 坐标 承受内压 承受离心力\n' ) ;
fprintf( ' (mm) 有限元解 解析解 有限元解 解析解\n' ) ;
fprintf( '------------------------------------------------------------------\n' ) ;
for i=1:length(node)
fprintf( '%4d %4.0f %8.2f %8.2f %8.2f %8.2f\n', ...
node(i), r(i)*1000, sr1_fem(i)/1e6, sr1(i)/1e6, sr2_fem(i)/1e6, sr2(i)/1e6 ) ;
end
fprintf( '==================================================================\n\n\n' ) ;
fprintf( ' 表3 周向应力 (MPa)\n' ) ;
fprintf( '==================================================================\n' ) ;
fprintf( '结点号 坐标 承受内压 承受离心力\n' ) ;
fprintf( ' (mm) 有限元解 解析解 有限元解 解析解\n' ) ;
fprintf( '------------------------------------------------------------------\n' ) ;
for i=1:length(node)
fprintf( '%4d %4.0f %8.2f %8.2f %8.2f %8.2f\n', ...
node(i), r(i)*1000, st1_fem(i)/1e6, st1(i)/1e6, st2_fem(i)/1e6, st2(i)/1e6 ) ;
end
fprintf( '==================================================================\n\n\n' ) ;
return⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -