代码输出结果中含有较多NaN,给出修改意见:
function dengcan8exam4_2( m, n )
clc;
% 本程序为第四章的第一个算例,采用矩形单元计算纯弯梁的变形
% exam4_1(m,n)
% 输入参数:
% m ------ x方向单元数目
% n ------ y方向单元数目
% 定义全局变量
% gNode ------ 节点坐标
% gElement --- 单元定义
% gMaterial -- 材料性质
% gBC1 ------- 第一类约束条件
% gDF -------- 分布力
% gK --------- 整体刚度矩阵
% gDelta ----- 整体节点坐标
global gNode gElement gMaterial gBC1 gDF gK gDelta
if nargin < 1
m = 2 ;
n = 2 ;
elseif nargin < 2
n = 4 ;
end
FemModel(m, n) ; % 定义有限元模型
SolveModel ; % 求解有限元模型
DisplayResults ; % 显示计算结果
return ;
function FemModel(m, n)
% 定义有限元模型
% 输入参数:
% m --- x方向单元数目
% n --- y方向单元数目
% 返回值:
% 无
% 说明:
% 该函数定义平面杆系的有限元模型数据:
% gNode ------- 节点定义
% gElement ---- 单元定义
% gMaterial --- 材料定义,包括弹性模量,梁的截面积和梁的抗弯惯性矩
% gBC --------- 约束条件
% gDF --------- 分布力
global gNode gElement gMaterial gBC1 gDF
length = 0.04 ; % 计算部分的长度(x方向)
height = 0.04 ; % 计算部分的高度(y方向)
dx = length/(2*m) ; % 矩形单元的宽度
dy = height/(2*n) ; % 矩形单元的高度
p = -1e11 ; % 弯曲应力10MPa
% 节点坐标
gNode = zeros( (2*m+1)*(2*n+1), 2 ) ;
for i=1:2*n+1
for j=1:2*m+1
k = (i-1)*(2*m+1)+j ; % 节点号
xk = (j-1)*dx; % 节点的x坐标
yk = (i-1)*dy ; % 节点的y坐标
gNode(k,:) = [xk, yk] ;
end
end
% 单元定义
gElement = zeros( m*n, 8 ) ;
for i=1:n
for j=1:m
k = (i-1)*m + j; % 单元号
col = 2*m + 1; % gNode的x方向节点数(列数)
% 角节点(对应gNode的(2m+1)列,每单元占2个间隔)
% 单元角节点编号(左下n1、右下n2、右上n3、左上n4)
n1 = 2*(i-1) * col + 2*(j-1) + 1; % 左下
n2 = 2*(i-1) * col + (2*(j-1) + 2) + 1; % 右下
n4 = (2*(i-1) + 2) * col + 2*(j-1) + 1; % 左上
n3 = (2*(i-1) + 2) * col + (2*(j-1) + 2) + 1;% 右上
% 边中点(角节点中间的节点)
n5 = n1 + 1; % 1-2边中点
n6 = n2 + col; % 2-3边中点
n7 = n3 - 1; % 3-4边中点
n8 = n4 - col; % 4-1边中点
gElement(k,:) = [n1, n2, n3, n4, n5, n6, n7, n8];
end
end
% 材料性质
% 弹性模量 泊松比 厚度
gMaterial = [2.2e11, 0.25, 0.01] ; % 材料 1
% 第一类约束条件
col = 2*m + 1; % x方向节点总数(列数,含中点)
row = 2*n + 1; % y方向节点总数(行数,含中点)
left_nodes = (0:row-1)*col +1; % 生成左端节点编号向量,长度=row
under_nodes= (0:col-1)+1; % 生成下端节点编号向量,长度=row-1
% 定义gBC1的大小
gBC1 = zeros(col+1, 3);
for j = 1 :col
gBC1(j, :) = [under_nodes(j), 2, 0.0];
end
gBC1(col+1,:) = [1, 1, 0.0];
% 分布载荷(线性分布)
gDF = zeros( n, 5 ) ;
for i=1:n
k =(n-1)*m+i;
gDF(i,:) = [ k, 3, (i-1)*p/n, i*p/n, 2] ;
end
return
function SolveModel
% 求解有限元模型
% 输入参数:
% 无
% 返回值:
% 无
% 说明:
% 该函数求解有限元模型,过程如下
% 1. 计算单元刚度矩阵,集成整体刚度矩阵
% 2. 计算单元的等效节点力,集成整体节点力向量
% 3. 处理约束条件,修改整体刚度矩阵和节点力向量
% 4. 求解方程组,得到整体节点位移向量
global gNode gElement gMaterial gBC1 gDF gK gDelta
% step1. 定义整体刚度矩阵和节点力向量
[node_number,dummy] = size( gNode ) ;
gK = zeros( node_number * 2, node_number * 2 ) ;
f = zeros( node_number * 2, 1 ) ;
% step2. 计算单元刚度矩阵,并集成到整体刚度矩阵中
[element_number,dummy] = size( gElement ) ;
for ie=1:1:element_number
k = StiffnessMatrix( ie );
AssembleStiffnessMatrix( ie, k ) ;
end
% step3. 计算分布力的等效节点力,并集成到整体节点力向量中(四节点适配)
[df_number, dummy] = size(gDF);
for idf = 1:df_number
% 1. 提取当前荷载的参数,调用等效节点力函数(enf为8×1向量)
ie = gDF(idf, 1); % 单元号
iedge = gDF(idf, 2); % 边号(1-4)
p1 = gDF(idf, 3); % 起始力集度
p2 = gDF(idf, 4); % 终止力集度
idof = gDF(idf, 5); % 力方向(1=x,2=y)
enf = EquivalentNodeForce(ie, iedge, p1, p2, idof); % 四节点enf:8×1(节点1-4的x/y)
% 2. 获取当前单元的8个节点全局编号
node_ids = gElement(ie, 1:8); % node_ids(1)=节点1,node_ids(2)=节点2,node_ids(3)=节点3,node_ids(4)=节点4
% 3. 将enf的8维力按节点分配到整体力向量f
for i = 1:8
node_global = node_ids(i); % 当前节点的全局编号
% 提取enf中当前节点的x/y力(i=1→enf(1:2),i=2→enf(3:4),以此类推)
enf_node = enf(2*(i-1)+1 : 2*i); % 每个节点对应enf的2个元素(x/y)
% 分配到整体力向量f的对应位置(全局节点的x/y自由度)
f(2*(node_global-1)+1 : 2*node_global) = f(2*(node_global-1)+1 : 2*node_global) + enf_node;
end
end
% step4. 处理约束条件
global modK p q;
modK = gK;
[i_bc, ~] = size(gBC1);
for x=1:i_bc
node = gBC1(x,1); % 约束节点编号
dof = gBC1(x,2); % 约束自由度
bc_p = (node-1)*2 + dof; % 约束自由度的全局索引
bc_q = gBC1(x,3); % 约束位移值
% 放大刚度,确保约束生效
modK(bc_p, :) = 0;
modK(:, bc_p) = 0;
modK(bc_p, bc_p) = gK(bc_p, bc_p) * 1e15; % 关键:放大1e15倍
f(bc_p) = modK(bc_p, bc_p) * bc_q; % 修正力向量
end
% step 5. 求解方程组,得到节点位移向量
gDelta = modK \ f ;
return
function k = StiffnessMatrix( ie )
% 计算平面应变等参数单元的刚度矩阵
% 输入参数:
%
% ie -- 单元号
% 返回值:
% k-- 单元刚度矩阵
% 说明:
% 用高斯积分法求解平面等参数单元的刚度矩阵
global gElement gMaterial
k=zeros(16,16);
E = gMaterial( 1,1 ) ;
mu = gMaterial( 1,2 ) ;
h = gMaterial( 1,3 ) ;
%平面应力的弹性常数
A1 = mu;
A2 = (1-mu)/2;
A3 = E/(1-mu^2);
D = A3* [ 1 A1 0
A1 1 0
0 0 A2] ;
% 2 x 2 高斯积分点和权系数
x = [-0.577350269189626, 0.577350269189626] ;
w = [1, 1] ;
for i=1:1:length(x)
for j=1:1:length(x)
B = MatrixB( ie, x(i), x(j) ) ;
J = Jacobi( ie, x(i), x(j) ) ;
k = k + h*w(i)*w(j)*transpose(B)*D*B*det(J) ;
end
end
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 B = MatrixB( ie, xi, eta )
% 计算单元的应变矩阵B
% 输入参数:
% ie --------- 单元号
% xi,eta ----- 局部坐标
% 返回值:
% B --------- 在局部坐标处的应变矩阵B
[N_x,N_y] = N_xy( ie, xi, eta );
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 enf= EquivalentNodeForce(ie,iedge,p1,p2,idof)
% 计算线性分布荷载的等效节点力
% 输入参数:
% ie ----- 单元号
% ieged --- 作用的边号
% p1 ----- 第一个节点上的分布力集度值
% p2 ----- 第二个节点上的分布力集度值
% idof --- 分布力的方向
% 1 --- x方向
% 2 --- y方向
% 返回值:
% enf ----- 等效节点力向量
global gElement gNode gMaterial
enf=zeros(16,1);
h = gMaterial( 3 ) ;
% 获取当前单元8个节点的整体坐标
node = gElement(ie, 1:8); %单元节点编号(按1-2-3-4逆时针顺序)
% 高斯积分参数
g = [-0.577350269189626, 0.577350269189626] ;
w = [1, 1] ;
% 确定边界对应的局部坐标规则
[fixed_coord, var_coord, coord_type] = getBoundaryParam(iedge);
% 积分计算等效节点力
node_force= zeros(8, 1);
for i=1:length(g)
if strcmp(coord_type, 'xi')
xi = fixed_coord; % 固定ξ
eta = g(i); % 变化η
else
xi = g(i); % 变化ξ
eta = fixed_coord; % 固定η
end
%计算当前积分点的形函数(N为1×4向量,对应4个节点)
N = ShapeFunction( xi, eta );
J = Jacobi( ie, xi, eta ) ;
x_xi=J(1,1);
x_eta=J(1,2);
y_xi=J(2,1);
y_eta=J(2,2);
if strcmp(coord_type, 'xi')
% 固定xi,ds由eta变化决定:
ds = sqrt((x_eta)^2 + (y_eta)^2);
else
% 固定eta,ds由xi变化决定:ds = sqrt((dx/dxi)^2 + (dy/dxi)^2)
ds = sqrt((x_xi)^2 + (y_xi)^2);
end
% 线性分布力插值
if strcmp(coord_type, 'xi')
s = (eta + 1) / 2;
else
s = (xi + 1) / 2;
end
p = p1*(1-s) + p2*s;
node_force = node_force + N'* p * h * ds * w(i);
end
for i = 1:8
enf(2*(i-1) + idof) = node_force(i); % 对应节点的idof方向力
end
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, x_eta ;
y_xi, y_eta ];
if det(J) <= 1e-12
warning('单元 %d 在积分点 (%f,%f) 雅克比矩阵接近奇异,det(J)=%e', ie, xi, eta, det(J));
end
return
function [fixed_coord, var_coord, coord_type] = getBoundaryParam(iedge)
% 定义八节点等参元各边的局部坐标规则
% fixed_coord: 固定的局部坐标值(如边2(右边界)固定ξ=1)
% var_coord: 变化的局部坐标(包含3个节点:角点1-中点-角点2)
% coord_type: 标记固定坐标类型('xi'或'eta')
switch iedge
case 1 % 边1:节点1-5-2(下边界,η=-1固定,ξ变化)
fixed_coord = -1;
var_coord = [-1, 0, 1]; % ξ从-1(节点1)→0(节点5)→1(节点2)
coord_type = 'eta'; % 固定的是eta
case 2 % 边2:节点2-6-3(右边界,ξ=1固定,η变化)
fixed_coord = 1;
var_coord = [-1, 0, 1]; % η从-1(节点2)→0(节点6)→1(节点3)
coord_type = 'xi'; % 固定的是xi
case 3 % 边3:节点3-7-4(上边界,η=1固定,ξ变化)
fixed_coord = 1;
var_coord = [1, 0, -1]; % ξ从1(节点3)→0(节点7)→-1(节点4)
coord_type = 'eta'; % 固定的是eta
case 4 % 边4:节点4-8-1(左边界,ξ=-1固定,η变化)
fixed_coord = -1;
var_coord = [1, 0, -1]; % η从1(节点4)→0(节点8)→-1(节点1)
coord_type = 'xi'; % 固定的是xi
otherwise
error('边号错误!iedge必须为1-4(对应节点1-5-2、2-6-3、3-7-4、4-8-1)');
end
return
function N = ShapeFunction( xi, eta )
% 计算形函数的值
% 输入参数:
% ie ---------- 单元号
% xi, eta ----- 单元内局部坐标
% 返回值:
% N ----------- 形函数的值
N5 = ( 1-eta ) * ( 1-xi^2 ) / 2 ;
N6 = ( xi + 1 ) * ( 1 - eta^2 ) / 2 ;
N7 = ( eta + 1 ) * ( 1 - xi^2 ) / 2 ;
N8 = (1- xi ) * ( 1-eta^2 ) / 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 DisplayResults
% 显示计算结果
% 输入参数:
% 无
% 返回值:
% 无
global gNode gDelta
gDelta=full(gDelta);
fprintf( '节点位移\n' ) ;
fprintf( ' 节点号 x方向位移 y方向位移\n' ) ;
[node_number,dummy] = size( gNode ) ;
for i=1:node_number
fprintf( '%6d %16.8e %16.8e\n',...
i, gDelta((i-1)*2+1), gDelta((i-1)*2+2) ) ;
end
return