简介本资源是一套面向工程力学、计算力学初学者及MATLAB编程实践者的三角形有限元分析完整实现代码包聚焦二维线性三角形单元Tri3的刚度矩阵构建、边界条件处理与系统求解全流程。包内共13个.m文件涵盖几何建模InputData系列、形函数计算Tri3Shape、单元刚度矩阵生成Tri3EleStif、全局矩阵组装LoopCalcGlobalStif、位移/载荷边界施加LoopBoundaryDisp/Load、弹性矩阵定义DMtrElas、主程序调度FemMain及结果可视化PlotResults结构清晰、模块职责明确便于逐层理解有限元核心逻辑。压缩包仅8KB轻量易读全部为可直接运行或调试的MATLAB脚本无外部依赖。目前已有2381人学习下载适合希望从零掌握三角形单元理论推导与MATLAB工程化落地的本科生、研究生及仿真入门工程师。1. 三角形单元不是“画个三角形就完事”MATLAB里真正跑通一个线性位移场有限元求解器需要6个核心函数协同、3类边界处理逻辑、2种刚度矩阵组装路径很多初学者拿到Tri3EleStif.m就以为“三角形单元”已经跑起来了——结果FemMain.m报错Index exceeds matrix dimensions或PlotResults.m画出的位移云图全为零。根本原因在于三角形单元在 MATLAB 中不是单个函数而是一套节点-单元-材料-边界-求解-后处理闭环系统。它强制你面对真实工程建模的底层约束每个三角形单元必须满足协调性C⁰连续、静力等效性虚功原理离散化、以及稀疏矩阵组装时的全局索引映射。这套代码包tri3.rar之所以能稳定复现经典平面应力/应变问题关键在于LoopCalcGlobalStif.m用列主元遍历所有单元并累加刚度LoopBoundaryDisp.m和LoopBoundaryLoad.m分别处理位移约束与载荷施加——二者不可互换顺序否则刚度矩阵会因自由度编号错位而奇异。适合正在用 MATLAB 实现《弹性力学有限元》课设、或需快速验证薄板/地基局部应力分布的结构工程师对只调用 PDE Toolbox 的用户价值有限但对想理解assembleFEMMatrix底层逻辑、调试自定义单元如含初始应力的 Tri3的人这是少有的可逐行 debug 的轻量级参考实现。2. 从形函数推导到刚度矩阵为什么 Tri3Shape.m 的插值系数必须用面积坐标而非笛卡尔坐标2.1 线性形函数的本质是面积坐标的仿射映射三角形单元的三个节点i, j, k构成的任意点 P 的位移 u(x,y) 被假定为线性插值u(x,y) Nᵢ(x,y)·uᵢ Nⱼ(x,y)·uⱼ Nₖ(x,y)·uₖ其中形函数 Nᵢ, Nⱼ, Nₖ 必须满足Nᵢ1 在节点 iNᵢ0 在 j,k且 ΣNₐ 1。若直接用笛卡尔坐标 (x,y) 构造需解三元一次方程组计算冗余且易受坐标系缩放影响。Tri3Shape.m采用面积坐标Barycentric coordinatesfunction [N, dNdx, dNdy] Tri3Shape(xy, xy_el) % xy: [x y] 坐标点 % xy_el: 3×2 矩阵每行是节点坐标 [xi yi; xj yj; xk yk] A 0.5 * det([xy_el(1,:) 1; xy_el(2,:) 1; xy_el(3,:) 1]); % 单元面积 L1 ((xy_el(2,1)-xy_el(3,1))*(xy(1)-xy_el(3,1)) ... (xy_el(2,2)-xy_el(3,2))*(xy(2)-xy_el(3,2))) / (2*A); L2 ((xy_el(3,1)-xy_el(1,1))*(xy(1)-xy_el(3,1)) ... (xy_el(3,2)-xy_el(1,2))*(xy(2)-xy_el(3,2))) / (2*A); L3 1 - L1 - L2; N [L1, L2, L3]; % 形函数对x,y的偏导数用于应变-位移矩阵B dNdx [(xy_el(2,2)-xy_el(3,2))/(2*A), ... (xy_el(3,2)-xy_el(1,2))/(2*A), ... (xy_el(1,2)-xy_el(2,2))/(2*A)]; dNdy [(xy_el(3,1)-xy_el(2,1))/(2*A), ... (xy_el(1,1)-xy_el(3,1))/(2*A), ... (xy_el(2,1)-xy_el(1,1))/(2*A)]; end提示det([... 1])计算有向面积符号决定单元定向逆时针为正。若A 0说明节点顺序错误会导致刚度矩阵负定——这是Tri3EleStif.m运行后位移发散的首要排查点。2.2 刚度矩阵组装Tri3EleStif.m 如何避免数值积分陷阱线性三角形单元的刚度矩阵 Kᵉ ∫∫_Ω Bᵀ D B |J| dξ dη 可解析积分高斯点数1无需数值积分。Tri3EleStif.m直接利用面积坐标性质function Ke Tri3EleStif(xy_el, E, nu, t, plane_type) % xy_el: 3×2 节点坐标 % E, nu: 弹性模量、泊松比t: 厚度plane_type: plane_stress or plane_strain A 0.5 * abs(det([xy_el(1,:) 1; xy_el(2,:) 1; xy_el(3,:) 1])); % 绝对面积 % 材料矩阵 D平面应力/应变切换 if strcmp(plane_type, plane_stress) D (E/(1-nu^2)) * [1, nu, 0; nu, 1, 0; 0, 0, (1-nu)/2]; else D (E/((1nu)*(1-2*nu))) * [1-nu, nu, 0; nu, 1-nu, 0; 0, 0, (1-2*nu)/2]; end % 计算应变-位移矩阵 B6×6因每个节点2个自由度 dNdx zeros(1,3); dNdy zeros(1,3); for i1:3 dNdx(i) (xy_el(mod(i,3)1,2) - xy_el(mod(i1,3)1,2)) / (2*A); dNdy(i) (xy_el(mod(i1,3)1,1) - xy_el(mod(i,3)1,1)) / (2*A); end B zeros(3,6); B(1,1:2:end) dNdx; % ε_x Σ dN_i/dx * u_i B(2,2:2:end) dNdy; % ε_y Σ dN_i/dy * v_i B(3,1:2:end) dNdy; % γ_xy Σ dN_i/dy * u_i Σ dN_i/dx * v_i B(3,2:2:end) dNdx; % 刚度矩阵 Ke t * A * B * D * B Ke t * A * B * D * B; end注意mod(i,3)1确保循环取节点索引1→2, 2→3, 3→1这是保证dNdx,dNdy符号一致的关键。若手写时误用i1导致越界B矩阵将含 NaN后续linsolve直接崩溃。2.3 全局刚度矩阵组装LoopCalcGlobalStif.m 的稀疏索引策略LoopCalcGlobalStif.m不用full(K)而用sparse因 1000 个单元对应 2000 自由度时稠密矩阵占内存约 300MB而稀疏矩阵仅约 8MB。其核心是预分配I,J,S向量function K_global LoopCalcGlobalStif(node_conn, xy, E, nu, t, plane_type) % node_conn: nelem×3 矩阵每行是单元节点编号 [i j k] nelem size(node_conn, 1); nnode size(xy, 1); ndof 2 * nnode; % 每个节点2个自由度 % 预分配稀疏矩阵三元组 max_nnz nelem * 36; % 每个3节点单元贡献6×6子矩阵 → 36非零元 I zeros(max_nnz, 1); J zeros(max_nnz, 1); S zeros(max_nnz, 1); idx 0; for e 1:nelem nodes node_conn(e, :); xy_el xy(nodes, :); % 提取当前单元节点坐标 Ke Tri3EleStif(xy_el, E, nu, t, plane_type); % 映射局部自由度到全局节点i → 全局自由度 [2*i-1, 2*i] dofs zeros(6,1); dofs(1:2:end) 2*nodes-1; % ux 分量 dofs(2:2:end) 2*nodes; % uy 分量 % 展开 Ke 到全局索引 for ii 1:6 for jj 1:6 idx idx 1; I(idx) dofs(ii); J(idx) dofs(jj); S(idx) Ke(ii,jj); end end end K_global sparse(I(1:idx), J(1:idx), S(1:idx), ndof, ndof); end关键参数说明dofs构造必须严格按[ux_i, uy_i, ux_j, uy_j, ux_k, uy_k]顺序否则Ke的行列无法正确映射。sparse(I,J,S,ndof,ndof)中ndof必须是总自由度数若漏算约束自由度K_global维度错误将导致linsolve报错Matrix dimensions must agree。3. 边界条件与载荷施加LoopBoundaryDisp.m 与 LoopBoundaryLoad.m 的执行时序不可逆3.1 位移边界条件为什么必须先修改刚度矩阵再施加载荷LoopBoundaryDisp.m的核心任务是将指定自由度设为已知位移如固定支座 u0这需两步操作置零行/列将 K_global 对应行、列全置零仅保留对角元为1修正右端项将已知位移乘以原刚度矩阵该行减去右端项function [K_mod, F_mod] LoopBoundaryDisp(K, F, bc_dof, bc_val) % bc_dof: 已知位移的自由度编号向量如 [1,2,5] 表示 ux1,uy1,ux20 % bc_val: 对应位移值向量如 [0,0,0] K_mod K; F_mod F; for i 1:length(bc_dof) dof bc_dof(i); K_mod(dof, :) 0; % 清零整行 K_mod(:, dof) 0; % 清零整列 K_mod(dof, dof) 1; % 对角元1 F_mod(dof) bc_val(i); % 右端项已知位移 % 修正其他行F_j - K_jdof * bc_val(i) F_mod F_mod - K(dof, :) * (bc_val(i) - F_mod(dof)); end end注意F_mod F_mod - K(dof, :) * (bc_val(i) - F_mod(dof))这行是关键——它将已知位移对其他自由度的影响从右端项中扣除。若跳过此步求解后非约束自由度位移将包含虚假刚体位移。3.2 力边界条件LoopBoundaryLoad.m 的节点力分配原则集中力必须分配到直接受力节点分布力需按形函数等效节点力。LoopBoundaryLoad.m处理两类输入load_nodes: 节点编号向量如[10,15]load_vals: 对应节点的[Fx, Fy]向量如[0, -1000]表示 y 方向 -1000Nfunction F LoopBoundaryLoad(F, load_nodes, load_vals) % load_vals: n×2 矩阵每行是 [Fx, Fy] for node load_nodes(i) for i 1:length(load_nodes) node load_nodes(i); F(2*node-1) F(2*node-1) load_vals(i, 1); % ux 分量 F(2*node) F(2*node) load_vals(i, 2); % uy 分量 end end提示若载荷作用于单元边如均布载荷 q100N/m需先计算等效节点力对边 ijq 等效为[0, -q*L/2, 0, -q*L/2]L 为边长再调用LoopBoundaryLoad。InputData*.m中load_nodes和load_vals的维度必须严格匹配否则F索引越界。3.3 执行时序FemMain.m 中不可颠倒的三步链完整求解流程在FemMain.m中固化为K_global LoopCalcGlobalStif(...)→ 组装未约束刚度矩阵[K_mod, F_mod] LoopBoundaryDisp(K_global, F_initial, bc_dof, bc_val)→ 施加位移约束F_final LoopBoundaryLoad(F_mod, load_nodes, load_vals)→ 施加力载荷警告若将步骤2和3颠倒LoopBoundaryLoad会向已被置零的行添加载荷导致F_final(dof)≠bc_val(i)求解后约束点位移偏离设定值。这是新手最常犯的错误调试时应打印F_mod和F_final前10行验证。4. 输入数据分层设计InputData*.m 如何通过文件名后缀控制网格密度与精度平衡4.1 文件命名隐含的网格规模协议InputData8.m,InputData16.m,InputData64.m,InputData256.m并非随意编号而是对应正方形域划分的单元边数文件名单元边数总单元数节点数无重合典型用途InputData8.m812881快速验证算法逻辑秒级求解InputData16.m16512289教学演示应力梯度初步显现InputData64.m6481924225工程级精度捕捉局部应力集中InputData256.m25613107266049高精度需求需稀疏求解器优化InputData*.m内部结构统一function [xy, node_conn, bc_dof, bc_val, load_nodes, load_vals] InputData64() % 正方形域 [0,1]×[0,1]64×64 网格 → 128×128 个三角形单元 nx 64; ny 64; [x, y] meshgrid(linspace(0,1,nx1), linspace(0,1,ny1)); xy [x(:), y(:)]; % 节点坐标 % 生成三角形单元连接表Delaunay 三角剖分简化版 node_conn []; for i 1:ny for j 1:nx n1 (i-1)*(nx1) j; n2 (i-1)*(nx1) j1; n3 i*(nx1) j; n4 i*(nx1) j1; % 划分两个三角形n1-n2-n3 和 n2-n4-n3 node_conn(end1,:) [n1, n2, n3]; node_conn(end1,:) [n2, n4, n3]; end end % 左侧固定x0 所有节点 uxuy0 left_nodes find(xy(:,1)0); bc_dof [2*left_nodes-1; 2*left_nodes]; bc_val zeros(size(bc_dof)); % 右侧中点施加向下力 right_mid_node find(xy(:,1)1 abs(xy(:,2)-0.5)1e-6); load_nodes right_mid_node; load_vals [0, -1000]; end4.2 材料与几何参数的集中管理DMtrElas.m 的模块化设计DMtrElas.m封装材料属性支持多材料区域function D_mat DMtrElas(E, nu, plane_type, mat_id) % mat_id: 材料编号用于区分不同区域如复合材料层 switch mat_id case 1 E_val E(1); nu_val nu(1); case 2 E_val E(2); nu_val nu(2); otherwise E_val E(1); nu_val nu(1); end if strcmp(plane_type, plane_stress) D_mat (E_val/(1-nu_val^2)) * [1, nu_val, 0; nu_val, 1, 0; 0, 0, (1-nu_val)/2]; else D_mat (E_val/((1nu_val)*(1-2*nu_val))) * ... [1-nu_val, nu_val, 0; nu_val, 1-nu_val, 0; 0, 0, (1-2*nu_val)/2]; end end技巧若要模拟带孔板hole in plate可在InputData*.m中用inpolygon排除圆内节点并调整node_conn删除包含这些节点的单元——Tri3BMtr.m边界矩阵计算会自动适配新拓扑无需修改核心求解器。5. 后处理可视化与应力精度验证PlotResults.m 的等效应力云图与理论解比对5.1 位移场绘制用 trisurf 实现单元中心插值PlotResults.m不直接画节点位移而是计算每个单元质心处的位移更符合物理意义function PlotResults(xy, node_conn, U, E, nu, plane_type, title_str) % U: 2*nnode×1 位移向量 [ux1;uy1;ux2;uy2;...] nnode size(xy,1); Ux U(1:2:end); Uy U(2:2:end); % 计算每个单元质心坐标及位移 n_elem size(node_conn, 1); centroid_x zeros(n_elem, 1); centroid_y zeros(n_elem, 1); U_centroid zeros(n_elem, 1); for e 1:n_elem nodes node_conn(e, :); centroid_x(e) mean(xy(nodes, 1)); centroid_y(e) mean(xy(nodes, 2)); % 线性插值U_centroid N(centroid) * [Ux(nodes); Uy(nodes)] xc centroid_x(e); yc centroid_y(e); [N,~,~] Tri3Shape([xc,yc], xy(nodes,:)); U_centroid(e) N * [Ux(nodes); Uy(nodes)]; end % 绘制位移云图 figure(Name, [Displacement Magnitude - title_str]); trisurf(node_conn, centroid_x, centroid_y, U_centroid, ... FaceColor,interp,EdgeColor,none); colormap(jet); colorbar; xlabel(x); ylabel(y); zlabel(Displacement); title([Displacement Magnitude (mm) - title_str]); end参数说明FaceColor,interp启用面内颜色插值使云图平滑trisurf的X,Y,Z,C参数中C是标量场此处为位移幅值node_conn定义三角面片连接关系。5.2 等效应力von Mises计算与理论验证Tri3BMtr.m计算每个单元的应力张量PlotResults.m进而生成 von Mises 应力% 在 PlotResults.m 中追加 % 计算每个单元的 von Mises 应力 sigma_vm zeros(n_elem, 1); for e 1:n_elem nodes node_conn(e, :); U_e [Ux(nodes); Uy(nodes)]; % 局部位移向量 xy_el xy(nodes, :); [N, dNdx, dNdy] Tri3Shape(mean(xy_el,1), xy_el); % 在质心处计算形函数导数 B zeros(3,6); B(1,1:2:end) dNdx; B(2,2:2:end) dNdy; B(3,1:2:end) dNdy; B(3,2:2:end) dNdx; D DMtrElas(E, nu, plane_type, 1); epsilon B * U_e; sigma D * epsilon; sigma_vm(e) sqrt( sigma(1)^2 sigma(2)^2 - sigma(1)*sigma(2) 3*sigma(3)^2 ); end % 绘制应力云图 figure(Name, [von Mises Stress - title_str]); trisurf(node_conn, centroid_x, centroid_y, sigma_vm, ... FaceColor,interp,EdgeColor,none); colormap(parula); colorbar; xlabel(x); ylabel(y); zlabel(\sigma_{vm} (Pa)); title([von Mises Stress (Pa) - title_str]);验证技巧对悬臂梁受端部剪力问题理论最大弯曲应力 σ_max 6FL/t²hF载荷L长度t厚度h高度。运行InputData64.m后用max(sigma_vm)与理论值比对相对误差应 5%网格足够密时。若误差 10%检查Tri3EleStif.m中plane_type是否误设为plane_strain悬臂梁属平面应力。5.3 快速定位收敛性问题三行命令诊断网格质量当sigma_vm出现尖锐跳变或负值大概率是网格畸变。用以下命令检查% 在 FemMain.m 求解后插入 % 计算每个单元的最小角弧度和长宽比 min_angles zeros(n_elem,1); aspect_ratios zeros(n_elem,1); for e 1:n_elem nodes node_conn(e,:); coords xy(nodes,:); % 计算三边长 a norm(coords(2,:)-coords(1,:)); b norm(coords(3,:)-coords(2,:)); c norm(coords(1,:)-coords(3,:)); % 最小角用余弦定理 cosA (b^2 c^2 - a^2)/(2*b*c); cosB (a^2 c^2 - b^2)/(2*a*c); cosC (a^2 b^2 - c^2)/(2*a*b); min_angles(e) min([acos(cosA), acos(cosB), acos(cosC)]); % 长宽比最长边 / 最短边 aspect_ratios(e) max([a,b,c]) / min([a,b,c]); end fprintf(Min angle: %.2f deg, Max aspect ratio: %.2f\n, min(min_angles)*180/pi, max(aspect_ratios)); % 若 min angle 15° 或 aspect ratio 10需重新生成网格实操建议对复杂边界用delaunay生成初始网格后调用pdetool的refine功能或adaptmesh局部加密比手动改InputData*.m更可靠。本文还有配套的精品资源点击获取