📄 exam6_2.m
字号:
return
function [N_x, N_y, N_z] = N_xyz( ie, xi, eta, zeta )
% 计算形函数对整体坐标的导数
% 输入参数：
% ie -------------- 单元号
% xi,eta,zeta ----- 局部坐标
% 返回值：
% N_x ------------ 在局部坐标处的形函数对x坐标的导数
% N_y ------------ 在局部坐标处的形函数对y坐标的导数
% N_z ------------ 在局部坐标处的形函数对z坐标的导数
J = Jacobi( ie, xi, eta, zeta ) ;
[N_xi,N_eta,N_zeta] = N_xietazeta( ie, xi, eta, zeta ) ;
N = [N_xi;N_eta;N_zeta] ;
A = inv(J)*N ;
N_x = A(1,:) ;
N_y = A(2,:) ;
N_z = A(3,:) ;
return
function [N_xi, N_eta, N_zeta] = N_xietazeta( ie, xi, eta, zeta )
% 计算形函数对局部坐标的导数
% 输入参数：
% ie -------------- 单元号
% xi,eta,zeta ----- 局部坐标
% 返回值：
% N_xi ---------- 在局部坐标处的形函数对xi坐标的导数
% N_eta ---------- 在局部坐标处的形函数对eta坐标的导数
% N_zeta --------- 在局部坐标处的形函数对zeta坐标的导数
x = [-1, 1, 1, -1, -1, 1, 1, -1, 0, 0, 0, 0, 1, 1, -1, -1, -1, 1, 1, -1];
e = [-1, -1, 1, 1, -1, -1, 1, 1, -1, 1, 1, -1, 0, 0, 0, 0, -1, -1, 1, 1];
z = [-1, -1, -1, -1, 1, 1, 1, 1, -1, -1, 1, 1, -1, 1, 1, -1, 0, 0, 0, 0];
N_xi = zeros( 1, 20 ) ;
N_eta = zeros( 1, 20 ) ;
N_zeta = zeros( 1, 20 ) ;
N_xi(1:8) = x(1:8).*(1+e(1:8)*eta).*(1+z(1:8)*zeta)/8 ;
N_eta(1:8) = (1+x(1:8)*xi).*e(1:8).*(1+z(1:8)*zeta)/8 ;
N_zeta(1:8) = (1+x(1:8)*xi).*(1+e(1:8)*eta).*z(1:8)/8 ;
N_xi(9:12) = -2*xi*(1+e(9:12)*eta).*(1+z(9:12)*zeta)/4 ;
N_eta(9:12) = ( 1-xi^2)*e(9:12).*(1+z(9:12)*zeta)/4 ;
N_zeta(9:12) = ( 1-xi^2)*(1+e(9:12)*eta).*z(9:12)/4 ;
N_xi(13:16) = ( 1-eta^2)*x(13:16).*(1+z(13:16)*zeta)/4 ;
N_eta(13:16) = -2*eta*(1+x(13:16)*xi).*(1+z(13:16)*zeta)/4 ;
N_zeta(13:16) = ( 1-eta^2)*(1+x(13:16)*xi).*z(13:16)/4 ;
N_xi(17:20) = ( 1-zeta^2)*x(17:20).*(1+e(17:20)*eta)/4 ;
N_eta(17:20) = ( 1-zeta^2)*(1+x(17:20)*xi).*e(17:20)/4 ;
N_zeta(17:20) = -2*zeta*(1+x(17:20)*xi).*(1+e(17:20)*eta)/4 ;
N_xi(1) = N_xi(1) - 0.5*(N_xi(9) +N_xi(16)+N_xi(17)) ;
N_xi(2) = N_xi(2) - 0.5*(N_xi(9) +N_xi(13)+N_xi(18)) ;
N_xi(3) = N_xi(3) - 0.5*(N_xi(10)+N_xi(13)+N_xi(19)) ;
N_xi(4) = N_xi(4) - 0.5*(N_xi(10)+N_xi(16)+N_xi(20)) ;
N_xi(5) = N_xi(5) - 0.5*(N_xi(12)+N_xi(15)+N_xi(17)) ;
N_xi(6) = N_xi(6) - 0.5*(N_xi(12)+N_xi(14)+N_xi(18)) ;
N_xi(7) = N_xi(7) - 0.5*(N_xi(11)+N_xi(14)+N_xi(19)) ;
N_xi(8) = N_xi(8) - 0.5*(N_xi(11)+N_xi(15)+N_xi(20)) ;
N_eta(1) = N_eta(1) - 0.5*(N_eta(9) +N_eta(16)+N_eta(17)) ;
N_eta(2) = N_eta(2) - 0.5*(N_eta(9) +N_eta(13)+N_eta(18)) ;
N_eta(3) = N_eta(3) - 0.5*(N_eta(10)+N_eta(13)+N_eta(19)) ;
N_eta(4) = N_eta(4) - 0.5*(N_eta(10)+N_eta(16)+N_eta(20)) ;
N_eta(5) = N_eta(5) - 0.5*(N_eta(12)+N_eta(15)+N_eta(17)) ;
N_eta(6) = N_eta(6) - 0.5*(N_eta(12)+N_eta(14)+N_eta(18)) ;
N_eta(7) = N_eta(7) - 0.5*(N_eta(11)+N_eta(14)+N_eta(19)) ;
N_eta(8) = N_eta(8) - 0.5*(N_eta(11)+N_eta(15)+N_eta(20)) ;
N_zeta(1) = N_zeta(1) - 0.5*(N_zeta(9) +N_zeta(16)+N_zeta(17)) ;
N_zeta(2) = N_zeta(2) - 0.5*(N_zeta(9) +N_zeta(13)+N_zeta(18)) ;
N_zeta(3) = N_zeta(3) - 0.5*(N_zeta(10)+N_zeta(13)+N_zeta(19)) ;
N_zeta(4) = N_zeta(4) - 0.5*(N_zeta(10)+N_zeta(16)+N_zeta(20)) ;
N_zeta(5) = N_zeta(5) - 0.5*(N_zeta(12)+N_zeta(15)+N_zeta(17)) ;
N_zeta(6) = N_zeta(6) - 0.5*(N_zeta(12)+N_zeta(14)+N_zeta(18)) ;
N_zeta(7) = N_zeta(7) - 0.5*(N_zeta(11)+N_zeta(14)+N_zeta(19)) ;
N_zeta(8) = N_zeta(8) - 0.5*(N_zeta(11)+N_zeta(15)+N_zeta(20)) ;
return
function J = Jacobi( ie, xi, eta, zeta )
% 计算雅克比矩阵
% 输入参数：
% ie -------------- 单元号
% xi,eta,zeta ----- 局部坐标
% 返回值：
% J ------------- 在局部坐标(xi,eta)处的雅克比矩阵
global gNode gElement
x = gNode(gElement(ie,1:20),1) ;
y = gNode(gElement(ie,1:20),2) ;
z = gNode(gElement(ie,1:20),3) ;
[N_xi, N_eta, N_zeta] = N_xietazeta( ie, xi, eta, zeta ) ;
x_xi = N_xi * x ;
x_eta = N_eta * x ;
x_zeta = N_zeta * x ;
y_xi = N_xi * y ;
y_eta = N_eta * y ;
y_zeta = N_zeta* y;
z_xi = N_xi * z ;
z_eta = N_eta * z ;
z_zeta = N_zeta * z ;
J = [ x_xi, y_xi, z_xi
x_eta, y_eta, z_eta
x_zeta, y_zeta, z_zeta ];
return
function AssembleStiffnessMatrix( ie, k )
% 把单元刚度矩阵集成到整体刚度矩阵
% 输入参数:
% ie --- 单元号
% k --- 单元刚度矩阵
% 返回值:
% 无
global gElement gK
for i=1:1:20
for j=1:1:20
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, pressure )
% 计算分布压力的等效节点力
% 输入参数:
% node ---------- 结点号
% pressure ------ 跟结点号对应的压力值
% 返回值:
% edf ----------- 等效节点力向量
global gNode
x = gNode( node, 1 ) ;
y = gNode( node, 2 ) ;
z = gNode( node, 3 ) ;
g = [-0.774596669241483, 0.0,0.774596669241483] ;
w = [ 0.555555555555556,0.888888888888889,0.555555555555556] ;
edf = zeros( 8*3, 1 ) ;
for i=1:length(g)
for j=1:length(g)
xi = g(i) ;
eta = g(j) ;
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 ] ;
N5_xi = (eta-1)*xi ;
N6_xi = 1/2*(1-eta^2) ;
N7_xi = -(eta+1)*xi ;
N8_xi = 1/2*(eta^2-1) ;
N1_xi = -1/4*(1-eta) - 1/2 * ( N5_xi + N8_xi) ;
N2_xi = 1/4*(1-eta) - 1/2 * (N5_xi + N6_xi) ;
N3_xi = 1/4*(eta+1) - 1/2 * (N6_xi + N7_xi) ;
N4_xi = -1/4*(eta+1) - 1/2 * (N7_xi + N8_xi ) ;
N_xi = [ N1_xi N2_xi N3_xi N4_xi N5_xi N6_xi N7_xi N8_xi ] ;
N5_eta = 1/2*(xi^2-1) ;
N6_eta = -eta*(1+xi) ;
N7_eta = 1/2*(1-xi^2) ;
N8_eta = eta*(xi-1) ;
N1_eta = -1/4*(1-xi) - 1/2 * ( N5_eta + N8_eta) ;
N2_eta = -1/4*(1+xi) - 1/2 * (N5_eta + N6_eta) ;
N3_eta = 1/4*(1+xi) - 1/2 * (N6_eta + N7_eta) ;
N4_eta = 1/4*(1-xi) - 1/2 * (N7_eta + N8_eta) ;
N_eta = [ N1_eta N2_eta N3_eta N4_eta N5_eta N6_eta N7_eta N8_eta ] ;
x_xi = N_xi * x ;
y_xi = N_xi * y ;
z_xi = N_xi * z ;
x_eta = N_eta * x ;
y_eta = N_eta * y ;
z_eta = N_eta * z ;
px = y_xi * z_eta - y_eta * z_xi ;
py = z_xi * x_eta - z_eta * x_xi ;
pz = x_xi * y_eta - x_eta * y_xi ;
sig = N * transpose( pressure ) ;
edf = edf + w(i)*w(j)*sig*[N1*eye(3);N2*eye(3);N3*eye(3);N4*eye(3);N5*eye(3);N6*eye(3);N7*eye(3);N8*eye(3)]*[px;py;pz] ;
end
end
return
function SaveResults( file_out )
% 保存计算结果
% 输入参数：
% 无
% 返回值：
% 无
global gNode gElement gMaterial gBC1 gDF gK gDelta gElementStress
save( file_out, 'gNode', 'gElement', 'gMaterial', 'gBC1', ...
'gDF', 'gDelta', 'gElementStress' ) ;
return
function DisplayResults
% 显示计算结果
% 输入参数：
% 无
% 返回值：
% 无
% 读取正应力沿截面高度的变化
global gElement gNode gMaterial gElementStress
[element_number,dummy] = size( gElement ) ;
node = [ 112 426 445 463 482 500 519 537 285 ] ;
node_num = zeros( size(node) ) ;
sigz = zeros( size( node ) ) ;
for n=1:length(node)
for ie=1:element_number
index = find( gElement(ie,1:20) == node(n) ) ;
if ~isempty( index )
sigz( n ) = sigz( n ) + gElementStress( ie, index, 3 ) ;
node_num(n) = node_num(n) + 1;
end
end
end
sigz = sigz./node_num ;
y = gNode( node, 2 ) ;
sigz1 = 16e6*(0-gNode(node(1),3))/pi * y ;
figure ;
plot( sigz1, y, '-', sigz, y, 'o' ) ;
xlabel( 'sigma-z' ) ;
ylabel( 'y' ) ;
% 读取剪应力沿截面宽度的变化
node = [285 366 365 364 363 362 361 360 137] ;
node_num = zeros( size(node) ) ;
taoyz = zeros( size( node ) ) ;
for n=1:length(node)
for ie=1:element_number
index = find( gElement(ie,1:20) == node(n) ) ;
if ~isempty( index )
taoyz( n ) = taoyz( n ) + gElementStress( ie, index, 4 ) ;
node_num(n) = node_num(n) + 1;
end
end
end
taoyz = taoyz./node_num ;
x = gNode( node, 1 ) ;
mu = gMaterial( 1, 2 ) ;
taoyz1 = -(3+2*mu)/8/(1+mu)*16e6/pi*(1-(1-2*mu)/(3+2*mu)*x.^2) ;
figure ;
plot( x, taoyz1, '-', x, taoyz, 'o' ) ;
xlabel( 'x' ) ;
ylabel( 'tao-xz' ) ;
set( gca, 'ylim', [min(taoyz), 0]*1.2 ) ;
fprintf( '\n\n\n 有限元与解析解的应力比较(MPa)\n' ) ;
fprintf( '==========================================================\n' ) ;
fprintf( ' 正应力(sigma-x) 剪应力(tao-yz)\n' ) ;
fprintf( ' 位置(y) 有限元 解析解 位置(x) 有限元 解析解\n' ) ;
fprintf( '----------------------------------------------------------\n' ) ;
for n=1:length(y)
fprintf( ' %6.3f %5.2f %5.2f %6.3f %5.2f %5.2f\n', ...
y(n), sigz(n)/1e6, sigz1(n)/1e6, ...
x(n), taoyz(n)/1e6, taoyz1(n)/1e6 ) ;
end
fprintf( '==========================================================\n' ) ;
return ⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -