MATLAB实现梁单元刚度矩阵的计算和总体刚度矩阵的组装
1. 理论核心：梁单元刚度矩阵
1.1 局部坐标系下的梁单元刚度矩阵
对于平面梁单元（每个节点有3个自由度：u, v, θ），在局部坐标系下的刚度矩阵为：
局部坐标系下的刚度矩阵公式：
k_local = [EA/L, 0, 0, -EA/L, 0, 0;
0, 12EI/L³, 6EI/L², 0, -12EI/L³, 6EI/L²;
0, 6EI/L², 4EI/L, 0, -6EI/L², 2EI/L;
-EA/L, 0, 0, EA/L, 0, 0;
0, -12EI/L³, -6EI/L², 0, 12EI/L³, -6EI/L²;
0, 6EI/L², 2EI/L, 0, -6EI/L², 4EI/L]
其中：
- E = 弹性模量
- A = 横截面积
- I = 截面惯性矩
- L = 单元长度
1.2 坐标变换：局部→全局
这是最关键的一步。实际结构中梁单元方向各异，需要将局部坐标系的刚度矩阵转换到全局坐标系：
坐标转换公式：
k_global = Tᵀ * k_local * T
其中变换矩阵 T 为：
T = [cosβ, sinβ, 0, 0, 0, 0;
-sinβ, cosβ, 0, 0, 0, 0;
0, 0, 1, 0, 0, 0;
0, 0, 0, cosβ, sinβ, 0;
0, 0, 0, -sinβ, cosβ, 0;
0, 0, 0, 0, 0, 1]
β 是梁单元与全局X轴的夹角。
2. MATLAB实现代码
2.1 主程序框架
% 有限元分析：梁结构总体刚度矩阵组装
clear; clc; close all;
%% 1. 输入参数
E = 2.1e11; % 弹性模量 (Pa, 钢)
A = 0.01; % 截面积 (m²)
I = 8.33e-6; % 惯性矩 (m⁴)
% 节点坐标 (x, y) [单位: m]
nodes = [0, 0; % 节点1
2, 0; % 节点2
4, 0; % 节点3
2, 2]; % 节点4
% 单元连接关系 [起始节点, 终止节点]
elements = [1, 2; % 单元1
2, 3; % 单元2
2, 4; % 单元3
3, 4]; % 单元4
num_nodes = size(nodes, 1); % 节点总数
num_elements = size(elements, 1); % 单元总数
dof_per_node = 3; % 每个节点的自由度
total_dof = num_nodes * dof_per_node; % 总自由度
%% 2. 初始化总体刚度矩阵
K_global = zeros(total_dof, total_dof);
%% 3. 循环处理每个单元
for e = 1:num_elements
% 3.1 获取单元信息
node1 = elements(e, 1);
node2 = elements(e, 2);
% 节点坐标
x1 = nodes(node1, 1); y1 = nodes(node1, 2);
x2 = nodes(node2, 1); y2 = nodes(node2, 2);
% 3.2 计算单元长度和角度
L = sqrt((x2 - x1)^2 + (y2 - y1)^2); % 单元长度
cos_beta = (x2 - x1) / L; % 方向余弦
sin_beta = (y2 - y1) / L; % 方向正弦
% 3.3 计算局部刚度矩阵
k_local = beamLocalStiffness(E, A, I, L);
% 3.4 计算变换矩阵
T = transformationMatrix(cos_beta, sin_beta);
% 3.5 转换到全局坐标系
k_global_element = T' * k_local * T;
% 3.6 组装到总体刚度矩阵
K_global = assembleGlobalStiffness(K_global, k_global_element, ...
node1, node2, dof_per_node);
end
%% 4. 输出结果
fprintf('=== 梁结构有限元分析结果 ===\n');
fprintf('节点数: %d\n', num_nodes);
fprintf('单元数: %d\n', num_elements);
fprintf('总自由度: %d\n', total_dof);
fprintf('总体刚度矩阵维度: %d × %d\n', size(K_global));
% 显示总体刚度矩阵（部分）
disp('总体刚度矩阵（前12×12部分）:');
disp(K_global(1:min(12,total_dof), 1:min(12,total_dof)));
2.2 核心函数：梁单元局部刚度矩阵
function k_local = beamLocalStiffness(E, A, I, L)
% 计算局部坐标系下的梁单元刚度矩阵
% 输入: E-弹性模量, A-截面积, I-惯性矩, L-单元长度
% 输出: 6×6局部刚度矩阵
% 预计算常用项
EA_L = E * A / L;
EI_L3 = 12 * E * I / L^3;
EI_L2 = 6 * E * I / L^2;
EI_L = 4 * E * I / L;
EI_L_half = 2 * E * I / L;
% 组装局部刚度矩阵
k_local = zeros(6, 6);
% 第一行
k_local(1,1) = EA_L; k_local(1,4) = -EA_L;
% 第二行
k_local(2,2) = EI_L3; k_local(2,3) = EI_L2;
k_local(2,5) = -EI_L3; k_local(2,6) = EI_L2;
% 第三行
k_local(3,2) = EI_L2; k_local(3,3) = EI_L;
k_local(3,5) = -EI_L2; k_local(3,6) = EI_L_half;
% 第四行 (对称于第一行)
k_local(4,1) = -EA_L; k_local(4,4) = EA_L;
% 第五行 (对称于第二行)
k_local(5,2) = -EI_L3; k_local(5,3) = -EI_L2;
k_local(5,5) = EI_L3; k_local(5,6) = -EI_L2;
% 第六行 (对称于第三行)
k_local(6,2) = EI_L2; k_local(6,3) = EI_L_half;
k_local(6,5) = -EI_L2; k_local(6,6) = EI_L;
% 利用对称性填充上三角部分
for i = 1:6
for j = i+1:6
k_local(i,j) = k_local(j,i);
end
end
end
2.3 核心函数：坐标变换矩阵
function T = transformationMatrix(cos_beta, sin_beta)
% 计算局部到全局坐标的变换矩阵
% 输入: cos_beta-方向余弦, sin_beta-方向正弦
% 输出: 6×6变换矩阵
T = zeros(6, 6);
% 变换矩阵的各个子块
T_sub = [cos_beta, sin_beta, 0;
-sin_beta, cos_beta, 0;
0, 0, 1];
% 填充变换矩阵
T(1:3, 1:3) = T_sub;
T(4:6, 4:6) = T_sub;
end
2.4 核心函数：总体刚度矩阵组装
function K_global = assembleGlobalStiffness(K_global, k_element, ...
node1, node2, dof_per_node)
% 将单元刚度矩阵组装到总体刚度矩阵
% 输入: K_global-当前总体矩阵, k_element-单元矩阵
% node1, node2-节点编号, dof_per_node-每节点自由度
% 输出: 更新后的总体刚度矩阵
% 计算在总体矩阵中的位置索引
dof_index1 = (node1-1)*dof_per_node + (1:dof_per_node);
dof_index2 = (node2-1)*dof_per_node + (1:dof_per_node);
% 所有相关自由度索引
dof_indices = [dof_index1, dof_index2];
% 组装过程 (直接加和法)
for i = 1:length(dof_indices)
row_global = dof_indices(i);
for j = 1:length(dof_indices)
col_global = dof_indices(j);
% 将单元矩阵的贡献加到总体矩阵
K_global(row_global, col_global) = ...
K_global(row_global, col_global) + k_element(i, j);
end
end
end
3. 结果验证与调试方法
3.1 验证总体刚度矩阵的特性
组装完成后，务必检查以下数学特性：
%% 验证总体刚度矩阵的特性
fprintf('\n=== 矩阵特性验证 ===\n');
% 1. 对称性检查
sym_error = max(max(abs(K_global - K_global')));
fprintf('对称性误差: %.2e (应接近0)\n', sym_error);
% 2. 半正定性检查（无约束时应有零特征值）
eigenvalues = eig(K_global);
fprintf('最小特征值: %.2e\n', min(eigenvalues));
fprintf('零特征值个数: %d (刚体模态)\n', sum(abs(eigenvalues) < 1e-9));
% 3. 对角线元素检查（应全为正，边界条件处理前）
diag_elements = diag(K_global);
fprintf('对角线元素最小值: %.2e\n', min(diag_elements));
fprintf('对角线元素最大值: %.2e\n', max(diag_elements));
3.2 可视化函数（可选但推荐）
function plotBeamStructure(nodes, elements, K_global)
% 可视化梁结构和刚度矩阵模式
figure('Position', [100, 100, 1200, 500]);
% 子图1：结构几何
subplot(1,2,1);
hold on; grid on; axis equal;
title('梁结构几何', 'FontSize', 12);
% 绘制单元
for e = 1:size(elements, 1)
node1 = elements(e, 1);
node2 = elements(e, 2);
plot([nodes(node1,1), nodes(node2,1)], ...
[nodes(node1,2), nodes(node2,2)], ...
'b-o', 'LineWidth', 2, 'MarkerSize', 8, ...
'MarkerFaceColor', 'r');
end
% 标注节点编号
for i = 1:size(nodes, 1)
text(nodes(i,1)+0.1, nodes(i,2)+0.1, ...
sprintf('%d', i), 'FontSize', 12, 'Color', 'k');
end
xlabel('X (m)'); ylabel('Y (m)');
% 子图2：刚度矩阵稀疏模式
subplot(1,2,2);
spy(K_global, 'b.', 10);
title('总体刚度矩阵稀疏模式', 'FontSize', 12);
xlabel('自由度编号'); ylabel('自由度编号');
% 添加矩阵信息
matrix_info = sprintf('维度: %d×%d\n非零元素: %d\n稀疏度: %.2f%%', ...
size(K_global,1), size(K_global,2), ...
nnz(K_global), 100*nnz(K_global)/numel(K_global));
text(0.05, 0.95, matrix_info, 'Units', 'normalized', ...
'FontSize', 10, 'BackgroundColor', 'w');
end
参考代码 利用MATLAB计算梁单元刚度矩阵，并组装成总体刚度矩阵 www.3dddown.com/cnb/54607.html
4. 使用示例与扩展建议
4.1 运行示例
% 在主程序末尾调用可视化函数
plotBeamStructure(nodes, elements, K_global);
% 保存结果
save('beam_analysis_results.mat', 'K_global', 'nodes', 'elements', 'E', 'A', 'I');
4.2 扩展功能建议
- 添加边界条件处理：在组装后引入约束，消除刚体位移
- 载荷向量组装：添加集中力、分布载荷的处理
- 求解位移和内力：扩展求解功能
U = K\F - 支持三维梁单元：扩展为每节点6自由度的空间梁
- 材料非线性：引入非线性本构关系
- 动态分析：添加质量矩阵形成特征值问题
4.3 常见问题排查
- 矩阵奇异性：施加约束前总体刚度矩阵必然奇异（有刚体模态）
- 数值精度：使用双精度运算，避免大数吃小数
- 内存优化：对于大规模问题，使用稀疏矩阵存储
K_sparse = sparse(K_global)
浙公网安备 33010602011771号