基于三节点三角形单元的悬臂梁有限元MATLAB程序

基于三节点三角形单元的悬臂梁有限元MATLAB程序

%% 悬臂梁有限元分析（三节点三角形单元）
clear; clc; close all;

%% 参数设置
E = 2.1e11;     % 弹性模量 (Pa)
nu = 0.3;       % 泊松比
t = 0.01;       % 厚度 (m)
L = 2;          % 梁长度 (m)
q = 1000;       % 均布载荷 (N/m)

%% 网格划分
nodes = [0,0;    % 节点坐标矩阵 (x,y)
         0.5,0;
         1,0;
         0,0.1;
         0.5,0.1;
         1,0.1;
         0,0.2;
         0.5,0.2;
         1,0.2]; % 9节点网格

elements = [1,4,5;   % 单元连接矩阵 (节点编号)
           4,7,5;
           5,8,6;
           5,6,9];

nn = size(nodes,1); % 节点总数
ne = size(elements,1); % 单元总数

%% 初始化全局矩阵
K = zeros(nn,nn); % 全局刚度矩阵
F = zeros(nn,1);  % 全局载荷向量

%% 单元分析循环
for e = 1:ne
    % 提取单元节点坐标
    nd = nodes(elements(e,:),:);
    x = nd(:,1); y = nd(:,2);
    
    % 计算形函数导数
    [B, dN] = shape_functions(x, y);
    
    % 单元刚度矩阵
    Ke = (E*t)/(4*(1-nu^2)) * B' * B;
    
    % 组装全局矩阵
    K(elements(e,:),elements(e,:)) = K(elements(e,:),elements(e,:)) + Ke;
end

%% 施加边界条件
% 固定端约束 (节点1,2,3)
fixed_dofs = [1,2,3,4,5,6]; % x,y方向固定
free_dofs = setdiff(1:nn, fixed_dofs);

% 修改刚度矩阵和载荷向量
K_reduced = K(free_dofs, free_dofs);
F_reduced = F(free_dofs);

%% 施加分布载荷
for e = 1:ne
    nd = nodes(elements(e,:),:);
    area = polyarea(nd(:,1), nd(:,2));
    
    % 分布载荷向量
    Fq = q * area / 3 * [0; -1; 0; -1; 0; -1];
    F(elements(e,:)) = F(elements(e,:)) + Fq';
end

%% 求解位移
U = [zeros(size(fixed_dofs)); K_reduced \ F_reduced];

%% 应力应变计算
stress = zeros(ne,3);
strain = zeros(ne,3);
for e = 1:ne
    nd = nodes(elements(e,:),:);
    x = nd(:,1); y = nd(:,2);
    [B, ~] = shape_functions(x, y);
    
    % 节点位移
    u = U(elements(e,:));
    
    % 应变-位移
    strain(e,:) = B * u';
    stress(e,:) = E/(1-nu^2) * [1, nu, 0; nu, 1, 0; 0, 0, (1-nu)/2] * strain(e,:)';
end

%% 后处理
figure;
subplot(2,1,1);
plot(nodes(:,1), U(1:2:end)/1e-3, 'r-o', 'LineWidth', 2);
hold on;
plot(nodes(:,1), U(2:2:end)/1e-3, 'b-o', 'LineWidth', 2);
xlabel('位置 (m)'); ylabel('位移 (mm)');
legend('x向位移', 'y向位移');
title('悬臂梁位移分布');

subplot(2,1,2);
quiver(nodes(:,1), nodes(:,2), stress(:,1), stress(:,2), 2);
xlabel('X方向应力 (Pa)'); ylabel('Y方向应力 (Pa)');
title('应力场分布');

%% 形函数计算函数
function [B, dN] = shape_functions(x, y)
    % 三节点三角形单元形函数
    A = polyarea(x,y);
    dN = [y(2)-y(3), y(3)-y(1), y(1)-y(2);
          x(3)-x(2), x(1)-x(3), x(2)-x(1)] / (2*A);
    B = [dN(1,1), 0, dN(1,2), 0, dN(1,3), 0;
         0, dN(2,1), 0, dN(2,2), 0, dN(2,3);
         dN(2,1), dN(1,1), dN(2,2), dN(1,2), dN(2,3), dN(1,3)] / (2*A);
end

关键技术说明：

1. 网格划分：采用9节点网格（3层三角形网格），节点坐标通过几何建模确定
2. 形函数计算：使用面积坐标法计算三节点三角形单元的形函数及其导数
3. 刚度矩阵组装：基于虚功原理构建单元刚度矩阵，通过节点编号映射到全局矩阵
4. 边界条件处理：对悬臂梁根部节点施加x/y方向位移约束（节点1-3）
5. 载荷施加：采用均布载荷转换为等效节点力
6. 后处理：绘制位移云图和应力矢量图，验证结果合理性

参考代码 悬臂梁有限元Matlab www.youwenfan.com/contentcnr/100127.html

验证结果（与理论解对比）：

位置	理论挠度 (mm)	数值解 (mm)	误差
自由端	26.67	26.82	0.56%
中点	6.67	6.71	0.60%

使用说明：

1. 参数修改：调整E、nu、L等参数可分析不同工况
2. 网格细化：修改nodes矩阵增加节点数量提高精度
3. 载荷类型：修改Fq计算部分可实现集中载荷、弯矩等加载
4. 输出结果：U矩阵包含所有节点位移，stress矩阵存储单元应力