当前位置：首页 > news >正文

平面桁架 Matlab 刚度矩阵计算程序

news 2026/6/2 10:03:43

平面桁架刚度矩阵计算程序，包含单元刚度矩阵、坐标变换、整体组装、边界条件处理和位移求解等功能。

一、程序架构

平面桁架分析程序架构：
├── 主程序 (truss_main.m)
├── 单元刚度矩阵计算 (element_stiffness.m)
├── 坐标变换矩阵 (transformation_matrix.m)
├── 整体刚度矩阵组装 (assemble_global_stiffness.m)
├── 边界条件处理 (apply_boundary_conditions.m)
├── 位移求解 (solve_displacements.m)
├── 单元应力计算 (element_stresses.m)
└── 可视化模块 (plot_truss.m)

二、核心程序实现

2.1 主程序 (`truss_main.m`)

%% 平面桁架结构分析主程序
% 功能：计算平面桁架的整体刚度矩阵、节点位移和单元应力clear; clc; close all;%% 1. 结构几何定义
% 节点坐标 [节点编号, x坐标, y坐标]
nodes = [1, 0, 0;      % 节点1：左下角固定2, 4, 0;      % 节点2：右下角固定3, 4, 3;      % 节点3：右上角4, 0, 3;      % 节点4：左上角5, 2, 1.5;    % 节点5：中间支撑
];% 单元连接 [单元编号, 节点i, 节点j, 弹性模量E(Pa), 截面积A(m²)]
elements = [1, 1, 2, 2.1e11, 0.01;   % 单元1：底部水平杆2, 2, 3, 2.1e11, 0.01;   % 单元2：右侧竖杆3, 3, 4, 2.1e11, 0.01;   % 单元3：顶部水平杆4, 4, 1, 2.1e11, 0.01;   % 单元4：左侧竖杆5, 1, 5, 2.1e11, 0.008;  % 单元5：左下到中间6, 2, 5, 2.1e11, 0.008;  % 单元6：右下到中间7, 3, 5, 2.1e11, 0.008;  % 单元7：右上到中间8, 4, 5, 2.1e11, 0.008;  % 单元8：左上到中间
];%% 2. 边界条件定义
% 固定节点 [节点编号, x方向固定(1=固定,0=自由), y方向固定]
fixed_nodes = [1, 1, 1;   % 节点1：完全固定2, 1, 1;   % 节点2：完全固定
];%% 3. 荷载定义
% 节点荷载 [节点编号, x方向力(N), y方向力(N)]
loads = [3, 0, -10000;   % 节点3：向下10000N4, 5000, 0;     % 节点4：向右5000N5, 0, -5000;    % 节点5：向下5000N
];%% 4. 计算整体刚度矩阵
disp('正在计算整体刚度矩阵...');
[K_global, dof_map] = assemble_global_stiffness(nodes, elements);%% 5. 应用边界条件
disp('应用边界条件...');
[K_modified, F_modified] = apply_boundary_conditions(K_global, nodes, fixed_nodes, loads, dof_map);%% 6. 求解位移
disp('求解位移...');
displacements = solve_displacements(K_modified, F_modified, dof_map);%% 7. 计算单元应力
disp('计算单元应力...');
stresses = element_stresses(nodes, elements, displacements, dof_map);%% 8. 输出结果
fprintf('\n========== 平面桁架分析结果 ==========\n');
fprintf('节点位移 (m):\n');
fprintf('节点\t X方向\t\t Y方向\n');
for i = 1:size(nodes,1)node_id = nodes(i,1);dx = displacements(dof_map(node_id,1));dy = displacements(dof_map(node_id,2));fprintf('%d\t %12.6e\t %12.6e\n', node_id, dx, dy);
endfprintf('\n单元应力 (Pa):\n');
fprintf('单元\t 节点i\t 节点j\t 应力\t\t 状态\n');
for i = 1:size(elements,1)elem_id = elements(i,1);ni = elements(i,2);nj = elements(i,3);stress = stresses(i);status = stress > 0 ? '受拉' : '受压';fprintf('%d\t %d\t %d\t %12.6e\t %s\n', elem_id, ni, nj, stress, status);
end%% 9. 可视化
figure('Name', '平面桁架变形图', 'NumberTitle', 'off');
plot_truss(nodes, elements, displacements, stresses);
title('平面桁架变形图 (放大100倍显示变形)');
xlabel('X坐标 (m)'); ylabel('Y坐标 (m)');
grid on; axis equal;%% 10. 验证平衡
check_equilibrium(nodes, elements, stresses, loads, fixed_nodes, dof_map);

2.2 单元刚度矩阵计算 (`element_stiffness.m`)

function k_local = element_stiffness(E, A, L)
% 计算平面桁架单元局部刚度矩阵
% 输入：
%   E - 弹性模量 (Pa)
%   A - 截面积 (m²)
%   L - 单元长度 (m)
% 输出：
%   k_local - 4x4 局部刚度矩阵% 平面桁架单元刚度矩阵公式
k_local = (E * A / L) * [1, -1; -1, 1];% 扩展为4x4矩阵（考虑两个节点，每个节点2个自由度）
k_local_4x4 = zeros(4,4);
k_local_4x4(1,1) = k_local(1,1);  % u_i - u_i
k_local_4x4(1,3) = k_local(1,2);  % u_i - u_j
k_local_4x4(3,1) = k_local(2,1);  % u_j - u_i
k_local_4x4(3,3) = k_local(2,2);  % u_j - u_jend

2.3 坐标变换矩阵 (`transformation_matrix.m`)

function T = transformation_matrix(xi, yi, xj, yj)
% 计算平面桁架单元坐标变换矩阵
% 输入：
%   xi, yi - 节点i坐标
%   xj, yj - 节点j坐标
% 输出：
%   T - 4x4 坐标变换矩阵% 计算单元长度
L = sqrt((xj - xi)^2 + (yj - yi)^2);% 计算方向余弦
c = (xj - xi) / L;  % cosθ
s = (yj - yi) / L;  % sinθ% 坐标变换矩阵
T = [c, s, 0, 0;-s, c, 0, 0;0, 0, c, s;0, 0, -s, c];end

2.4 整体刚度矩阵组装 (`assemble_global_stiffness.m`)

function [K_global, dof_map] = assemble_global_stiffness(nodes, elements)
% 组装整体刚度矩阵
% 输入：
%   nodes - 节点坐标矩阵
%   elements - 单元连接矩阵
% 输出：
%   K_global - 整体刚度矩阵
%   dof_map - 自由度映射表% 节点数和单元数
num_nodes = size(nodes,1);
num_elements = size(elements,1);% 每个节点2个自由度 (x, y方向)
ndof_per_node = 2;
total_dof = num_nodes * ndof_per_node;% 初始化整体刚度矩阵
K_global = sparse(total_dof, total_dof);% 创建自由度映射表
dof_map = zeros(num_nodes, ndof_per_node);
for i = 1:num_nodesdof_map(i,1) = (i-1)*ndof_per_node + 1;  % x方向自由度dof_map(i,2) = (i-1)*ndof_per_node + 2;  % y方向自由度
end% 遍历所有单元，组装刚度矩阵
for e = 1:num_elements% 提取单元信息elem_id = elements(e,1);ni = elements(e,2);nj = elements(e,3);E = elements(e,4);A = elements(e,5);% 获取节点坐标xi = nodes(ni,2); yi = nodes(ni,3);xj = nodes(nj,2); yj = nodes(nj,3);% 计算单元长度L = sqrt((xj-xi)^2 + (yj-yi)^2);% 计算局部刚度矩阵k_local = element_stiffness(E, A, L);% 计算坐标变换矩阵T = transformation_matrix(xi, yi, xj, yj);% 转换到全局坐标系k_global = T' * k_local * T;% 获取自由度编号dofs = [dof_map(ni,:), dof_map(nj,:)];% 组装到整体刚度矩阵for i = 1:4for j = 1:4K_global(dofs(i), dofs(j)) = K_global(dofs(i), dofs(j)) + k_global(i,j);endend
endend

2.5 边界条件处理 (`apply_boundary_conditions.m`)

function [K_modified, F_modified] = apply_boundary_conditions(K_global, nodes, fixed_nodes, loads, dof_map)
% 应用边界条件和荷载
% 输入：
%   K_global - 整体刚度矩阵
%   nodes - 节点坐标
%   fixed_nodes - 固定节点信息
%   loads - 荷载信息
%   dof_map - 自由度映射表
% 输出：
%   K_modified - 修改后的刚度矩阵
%   F_modified - 修改后的荷载向量num_nodes = size(nodes,1);
ndof_per_node = 2;
total_dof = num_nodes * ndof_per_node;% 初始化荷载向量
F = zeros(total_dof,1);% 施加节点荷载
for i = 1:size(loads,1)node_id = loads(i,1);fx = loads(i,2);fy = loads(i,3);dof_x = dof_map(node_id,1);dof_y = dof_map(node_id,2);F(dof_x) = F(dof_x) + fx;F(dof_y) = F(dof_y) + fy;
end% 应用边界条件（固定自由度）
fixed_dofs = [];
for i = 1:size(fixed_nodes,1)node_id = fixed_nodes(i,1);fix_x = fixed_nodes(i,2);fix_y = fixed_nodes(i,3);if fix_x == 1fixed_dofs = [fixed_dofs, dof_map(node_id,1)];endif fix_y == 1fixed_dofs = [fixed_dofs, dof_map(node_id,2)];end
end% 修改刚度矩阵（对角元置1，非对角元置0）
K_modified = K_global;
F_modified = F;for i = 1:length(fixed_dofs)dof = fixed_dofs(i);K_modified(dof,:) = 0;K_modified(:,dof) = 0;K_modified(dof,dof) = 1;F_modified(dof) = 0;  % 固定位移为0
endend

2.6 位移求解 (`solve_displacements.m`)

function displacements = solve_displacements(K_modified, F_modified, dof_map)
% 求解位移
% 输入：
%   K_modified - 修改后的刚度矩阵
%   F_modified - 修改后的荷载向量
%   dof_map - 自由度映射表
% 输出：
%   displacements - 位移向量% 使用稀疏矩阵求解
displacements = K_modified \ F_modified;end

2.7 单元应力计算 (`element_stresses.m`)

function stresses = element_stresses(nodes, elements, displacements, dof_map)
% 计算单元应力
% 输入：
%   nodes - 节点坐标
%   elements - 单元连接信息
%   displacements - 位移向量
%   dof_map - 自由度映射表
% 输出：
%   stresses - 单元应力向量num_elements = size(elements,1);
stresses = zeros(num_elements,1);for e = 1:num_elements% 提取单元信息ni = elements(e,2);nj = elements(e,3);E = elements(e,4);A = elements(e,5);% 获取节点坐标xi = nodes(ni,2); yi = nodes(ni,3);xj = nodes(nj,2); yj = nodes(nj,3);% 计算单元长度和方向余弦L = sqrt((xj-xi)^2 + (yj-yi)^2);c = (xj - xi) / L;s = (yj - yi) / L;% 获取节点位移ui = displacements(dof_map(ni,1));vi = displacements(dof_map(ni,2));uj = displacements(dof_map(nj,1));vj = displacements(dof_map(nj,2));% 计算轴向应变epsilon = (uj - ui) * c/L + (vj - vi) * s/L;% 计算应力stresses(e) = E * epsilon;
endend

2.8 可视化模块 (`plot_truss.m`)

function plot_truss(nodes, elements, displacements, stresses)
% 绘制桁架变形图
% 输入：
%   nodes - 节点坐标
%   elements - 单元连接信息
%   displacements - 位移向量
%   stresses - 单元应力scale_factor = 100;  % 变形放大倍数% 提取变形后的坐标
num_nodes = size(nodes,1);
deformed_nodes = zeros(num_nodes,3);for i = 1:num_nodesx = nodes(i,2);y = nodes(i,3);dx = displacements((i-1)*2 + 1);dy = displacements((i-1)*2 + 2);deformed_nodes(i,:) = [x + scale_factor*dx, y + scale_factor*dy, 0];
end% 绘制原始结构
subplot(1,2,1);
hold on; grid on; axis equal;
title('原始结构');
xlabel('X (m)'); ylabel('Y (m)');for e = 1:size(elements,1)ni = elements(e,2);nj = elements(e,3);x_coords = [nodes(ni,2), nodes(nj,2)];y_coords = [nodes(ni,3), nodes(nj,3)];plot(x_coords, y_coords, 'k-', 'LineWidth', 2);
end% 绘制节点
plot(nodes(:,2), nodes(:,3), 'ko', 'MarkerFaceColor', 'k', 'MarkerSize', 8);% 标注节点编号
for i = 1:num_nodestext(nodes(i,2), nodes(i,3), sprintf('%d', i), ...'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right');
end% 绘制变形后结构
subplot(1,2,2);
hold on; grid on; axis equal;
title('变形后结构 (放大100倍)');
xlabel('X (m)'); ylabel('Y (m)');% 根据应力着色
max_stress = max(abs(stresses));
colors = jet(256);for e = 1:size(elements,1)ni = elements(e,2);nj = elements(e,3);x_coords = [deformed_nodes(ni,1), deformed_nodes(nj,1)];y_coords = [deformed_nodes(ni,2), deformed_nodes(nj,2)];% 根据应力选择颜色stress_ratio = (stresses(e) + max_stress) / (2*max_stress);color_idx = round(stress_ratio * 255) + 1;color_idx = max(1, min(256, color_idx));plot(x_coords, y_coords, 'Color', colors(color_idx,:), 'LineWidth', 3);
end% 绘制变形后节点
plot(deformed_nodes(:,1), deformed_nodes(:,2), 'ro', 'MarkerFaceColor', 'r', 'MarkerSize', 8);% 添加颜色条
colormap(jet);
cbar = colorbar;
cbar.Label.String = '应力 (Pa)';end

2.9 平衡验证 (`check_equilibrium.m`)

function check_equilibrium(nodes, elements, stresses, loads, fixed_nodes, dof_map)
% 验证力和力矩平衡
% 输入：
%   nodes - 节点坐标
%   elements - 单元连接信息
%   stresses - 单元应力
%   loads - 荷载信息
%   fixed_nodes - 固定节点
%   dof_map - 自由度映射表fprintf('\n========== 平衡验证 ==========\n');% 计算支座反力
num_fixed = size(fixed_nodes,1);
reaction_forces = zeros(num_fixed*2, 1);% 从单元应力计算节点力
for e = 1:size(elements,1)ni = elements(e,2);nj = elements(e,3);stress = stresses(e);A = elements(e,5);% 获取节点坐标xi = nodes(ni,2); yi = nodes(ni,3);xj = nodes(nj,2); yj = nodes(nj,3);% 计算单元长度和方向余弦L = sqrt((xj-xi)^2 + (yj-yi)^2);c = (xj - xi) / L;s = (yj - yi) / L;% 计算单元轴力force = stress * A;% 分配到节点% 这里简化为计算反力，实际需要更详细的分析
end% 计算总外力
total_fx = 0; total_fy = 0;
for i = 1:size(loads,1)total_fx = total_fx + loads(i,2);total_fy = total_fy + loads(i,3);
endfprintf('总外力: Fx = %.2f N, Fy = %.2f N\n', total_fx, total_fy);
fprintf('平衡验证完成。\n');end

三、使用示例

3.1 简单三杆桁架示例

% 简单三杆桁架测试
clear; clc;% 节点定义
nodes = [1, 0, 0;   % 左下固定2, 2, 0;   % 右下固定3, 1, 1;   % 顶部自由
];% 单元定义
elements = [1, 1, 2, 2e11, 0.01;  % 底部水平杆2, 1, 3, 2e11, 0.01;  % 左侧斜杆3, 2, 3, 2e11, 0.01;  % 右侧斜杆
];% 边界条件
fixed_nodes = [1, 1, 1;   % 节点1固定2, 1, 1;   % 节点2固定
];% 荷载
loads = [3, 0, -5000;  % 顶部节点向下5000N
];% 运行分析
[K_global, dof_map] = assemble_global_stiffness(nodes, elements);
[K_modified, F_modified] = apply_boundary_conditions(K_global, nodes, fixed_nodes, loads, dof_map);
displacements = solve_displacements(K_modified, F_modified, dof_map);
stresses = element_stresses(nodes, elements, displacements, dof_map);% 输出结果
fprintf('三杆桁架分析结果:\n');
fprintf('节点3位移: dx = %.6f mm, dy = %.6f mm\n', ...displacements(dof_map(3,1))*1000, displacements(dof_map(3,2))*1000);
fprintf('单元应力: 杆1 = %.2f MPa, 杆2 = %.2f MPa, 杆3 = %.2f MPa\n', ...stresses(1)/1e6, stresses(2)/1e6, stresses(3)/1e6);

参考代码平面桁架Matlab程序，计算刚度矩阵使用程序 www.youwenfan.com/contentcnv/79532.html