小型模型的有限元分析
2D 弹性有限元：
- 手算单元刚度矩阵
- 手写组装成稀疏全局刚度
- 用
K(r,r)\F(r)（反斜杠）解方程 - 恢复应力并画云图
例子用 带圆孔的平板单轴拉伸（你想要“小型模型”，这个正好，而且很经典）。
1）问题：带圆孔的平板 · 平面应力 · 小变形
尺寸（示意）：
- 板：宽
W，高H - 中心圆孔半径
r0 - 左端固定（
u=v=0） - 右端施加均匀拉应力
σ0（用 nodal 力等效）
你不想用搜索工具，所以我不拟合材料参数、不做优化，只做确定性的 静力 FE solve。
2）主脚本（单文件可跑）
%% ========== 小型 FE：带孔平板拉伸 ==========
clear; clc; close all;
%% ---- 几何/材料/网格参数 ----
L = 200; % 板半宽 (总宽 2L)
H = 120; % 板半高 (总高 2H)
R = 40; % 圆孔半径
nelx = 60; nely = 36; % 网格密度（小型即可跑）
E = 210e3; % MPa (钢约 210GPa→210e3MPa)
nu = 0.3; % 泊松
t = 10; % 厚度 (mm)
sig0= 100; % 名义拉应力 (MPa)
planeStress = true; % true=平面应力；false=平面应变
fprintf('=== 手写小型FE：带孔平板 ===\n');
%% ---- 1) 结构化网格 + 挖圆孔 ----
% 节点：按单元中心是否在孔外决定是否保留
nx = nelx+1; ny = nely+1;
x = linspace(-L, L, nx);
y = linspace(-H, H, ny);
[X,Y] = meshgrid(x, y); % Y=row=y, X=col=x
node = [X(:), Y(:)]; % N x 2
Nnodes = size(node,1);
% 挖孔：删掉孔内节点（并留一圈“合理”点）
hole = sqrt(X.^2+Y.^2) <= R;
hole = hole(:);
keep = ~hole;
% 重新编号
oldID = zeros(Nnodes,1); oldID(keep) = 1:sum(keep);
node = node(keep,:);
N = size(node,1);
fprintf('节点数: %d (原始 %d，挖掉 %d)\n', N, Nnodes, sum(hole));
%% ---- 2) 连接成 4-node四边形单元（手编拓扑） ----
elem = zeros((nelx)*(nely), 4);
eidx = 0;
for j = 1:nely
for i = 1:nelx
n1 = (j-1)*nx + i;
n2 = n1 + 1;
n3 = n2 + nx;
n4 = n1 + nx;
quad = [n1,n2,n3,n4];
if all(keep(quad))
eidx = eidx+1;
elem(eidx,:) = oldID(quad);
end
end
end
elem = elem(1:eidx,:);
Ne = size(elem,1);
fprintf('单元数: %d\n', Ne);
%% ---- 3) 材料矩阵 C ----
if planeStress
C = E/(1-nu^2)*[1 nu 0; nu 1 0; 0 0 (1-nu)/2];
else
lam = E*nu/((1+nu)*(1-2*nu));
mu = E/(2*(1+nu));
C = [lam+2*mu lam 0;
lam lam+2*mu 0;
0 0 mu];
end
%% ---- 4) 组装 K,F ----
% 每个节点2个DOF: [u1 v1 u2 v2 ...]
freeDOF = false(2*N,1);
K = sparse(2*N,2*N);
F = zeros(2*N,1);
for e = 1:Ne
conn = elem(e,:);
xe = node(conn,:);
[Ke, We] = quad4_stiffness(xe, C, t);
dof = [2*conn-1; 2*conn];
dof = dof(:);
K(dof,dof) = K(dof,dof) + Ke;
end
%% ---- 5) 边界条件（固定左边，拉右边） ----
tolBC = 1e-6;
for i = 1:N
xi = node(i,1); yi = node(i,2);
% 固定左边缘
if abs(xi + L) < tolBC
freeDOF(2*i-1) = false; % u=0
freeDOF(2*i) = false; % v=0
else
freeDOF(2*i-1) = true;
freeDOF(2*i) = true;
end
end
% 右边缘：施加等效拉力（分布力 → 节点力）
for i = 1:N
xi = node(i,1); yi = node(i,2);
if abs(xi - L) < tolBC
% 这里用最简单的“按邻接边长度分配”
% 找邻居节点（同一右边缘）来近似edge segment长度
dy_right = 0;
for j = 1:N
if j~=i && abs(node(j,1)-L)<tolBC
dy_right = max(dy_right, abs(node(j,2)-yi));
end
end
if dy_right < tolBC, dy_right = y(2)-y(1); end
Le = dy_right; % 该节点代表的右边缘段长度
F(2*i-1) = F(2*i-1) + 0; % x方向无牵引（或你自己加）
F(2*i) = F(2*i) + 0;
% 若要纯 σx=sig0 的均匀拉：
F(2*i-1) = F(2*i-1) + sig0 * t * Le; % 向右拉（正x）
end
end
% 强制钉死左固定（反斜杠法，无搜索）
u = zeros(2*N,1);
allDOF = (1:2*N)';
fixedDOF = allDOF(~freeDOF);
freeSet = allDOF(freeDOF);
Kff = K(freeSet,freeSet);
Ff = F(freeSet) - K(freeSet,fixedDOF)*u(fixedDOF);
u(freeSet) = Kff \ Ff;
%% ---- 6) 恢复应力（单元内高斯点） ----
sx = zeros(Ne,1); sy = sx; sxy = sx; VonMises = sx;
for e = 1:Ne
conn = elem(e,:);
xe = node(conn,:);
ue = u([2*conn-1; 2*conn]); ue = ue(:);
[S] = quad4_stress(xe, ue, C);
sx(e)=S(1); sy(e)=S(2); sxy(e)=S(3);
VonMises(e)=sqrt( S(1)^2 - S(1)*S(2) + S(2)^2 + 3*S(3)^2 );
end
%% ---- 7) 可视化 ----
figure('Color','w','Pos',[100 80 1050 440]);
% 网格
subplot(1,3,1); hold on;
for e=1:Ne
patch(node(elem(e,:),1), node(elem(e,:),2),'w','EdgeColor',[0.6 0.6 0.6]);
end
plot(node(:,1),node(:,2),'.k','MarkerS',3);
title('网格（带孔）'); axis equal; grid on;
% 位移放大
subplot(1,3,2); hold on;
fact = 5e3 / max(abs(u)+eps)*t; % 适量放大
udisp = [u(1:2:end-1), u(2:2:end)];
patch('XData',node(elem(:,1),1)+fact*udisp(elem(:,1),1), ...
'YData',node(elem(:,1),2)+fact*udisp(elem(:,1),2), ...
'Vertices',node+fact*udisp,'Faces',elem,'FaceVertexCData',sqrt(sum(udisp.^2,2)), ...
'FaceColor','interp','EdgeColor','k','LineW',0.2);
colorbar; title('|u| 位移场'); axis equal; grid on;
% Von Mises 云图（分片常数→画单元色块）
subplot(1,3,3); hold on;
col = VonMises;
for e=1:Ne
fill(node(elem(e,:),1), node(elem(e,:),2), 'flat', ...
'CData',col(e)*[1 1 1], 'EdgeColor',0.6*[1 1 1],'LineW',0.2);
end
% 用 scatter trick 给个 colorbar
scatter(node(:,1),node(:,2),0,VonMises(elem2node_mean(elem)),'filled');
colormap(jet); colorbar;
title('Von Mises 应力 (单元常数)'); axis equal; grid on;
sgtitle('手写小型FE：带孔平板拉伸');
fprintf('Done. 最大 |u|≈ %.3f mm\n', max(abs(udisp(:))));
参考代码 小型模型的有限元分析 www.youwenfan.com/contentcsv/113038.html
3）两个“力学内核”函数
3.1 四节点四边形单元刚度 quad4_stiffness.m
function [Ke, We] = quad4_stiffness(xe, C, t)
% xe: 4x2 节点坐标
% C: 3x3 本构矩阵（平面应力/应变）
% t: 厚度
% 返回 Ke: 8x8
% 2x2 Gauss
g = 1/sqrt(3);
gp = [-g g; -g g; g -g; g g];
w = [1 1; 1 1];
Ke = zeros(8);
for q = 1:4
xi = gp(q,1); eta = gp(q,2);
% 形函数及其导对(xi,eta)
N = 0.25*[ (1-xi)*(1-eta)
(1+xi)*(1-eta)
(1+xi)*(1+eta)
(1-xi)*(1+eta) ];
dNdxi = 0.25*[ -(1-eta) (1-eta) (1+eta) -(1+eta)
-(1-xi) -(1+xi) (1+xi) (1-xi) ];
dN = reshape(dNdxi,2,4);
% Jacobian
J = dN * xe; % 2x2
detJ = det(J);
invJ = inv(J);
dNdx = invJ * dN; % 2x4 → 每行是 dN/dx 或 dN/dy
% B matrix (3x8)
B = zeros(3,8);
for a = 1:4
B(1,2*a-1) = dNdx(1,a);
B(2,2*a) = dNdx(2,a);
B(3,2*a-1) = dNdx(2,a);
B(3,2*a) = dNdx(1,a);
end
Ke = Ke + B'*C*B * (detJ) * t * w(q,1)*w(q,2);
end
end
3.2 应力恢复 quad4_stress.m
function S = quad4_stress(xe, ue, C)
% 在单元中心 (0,0) 估计应力
xi=0; eta=0;
dNdxi = 0.25*[ -(1-eta) (1-eta) (1+eta) -(1+eta)
-(1-xi) -(1+xi) (1+xi) (1-xi) ];
dN = reshape(dNdxi,2,4);
J = dN*xe;
invJ = inv(J);
dNdx = invJ*dN;
B = zeros(3,8);
for a = 1:4
B(1,2*a-1)=dNdx(1,a); B(2,2*a)=dNdx(2,a);
B(3,2*a-1)=dNdx(2,a); B(3,2*a)=dNdx(1,a);
end
S = C*B*ue;
end
3.3（可选）辅助：把单元值摊到节点作平滑显示
function nval = elem2node_mean(elem, eval)
% 把单元常数 eval(Ne,1) 平均到节点
N = max(elem(:));
nval = zeros(N,1); cnt = zeros(N,1);
for e = 1:size(elem,1)
for a = 1:4
n = elem(e,a);
nval(n) = nval(n)+eval(e);
cnt(n) = cnt(n)+1;
end
end
nval = nval./cnt;
end
