简支梁在荷载作用下的变形计算
用于计算简支梁在各种荷载作用下的变形。程序包括解析解和有限元数值解两种方法，并提供了可视化功能。
% 简支梁荷载作用下的变形计算
clear; clc; close all;
%% 1. 参数设置
fprintf('设置梁参数和荷载条件...\n');
% 梁的几何和材料参数
L = 10; % 梁长度 (m)
E = 200e9; % 弹性模量 (Pa) - 钢的典型值
I = 8.33e-6; % 截面惯性矩 (m^4) - 矩形截面 0.1m x 0.2m
rho = 7850; % 材料密度 (kg/m^3)
A = 0.02; % 横截面积 (m^2) - 0.1m x 0.2m
% 荷载参数
load_type = 'uniform'; % 荷载类型: 'point', 'uniform', 'triangular'
P = 10000; % 集中荷载大小 (N) - 用于点荷载
w = 1000; % 分布荷载集度 (N/m) - 用于均布和三角形荷载
a = 3; % 点荷载作用位置 (距左端距离, m)
% 有限元参数
n_elements = 20; % 单元数量
%% 2. 解析解计算
fprintf('计算解析解...\n');
% 创建位置向量
x = linspace(0, L, 1000);
% 根据荷载类型计算变形
switch load_type
case 'point'
% 集中荷载作用下的简支梁变形
deflection = zeros(size(x));
for i = 1:length(x)
if x(i) <= a
deflection(i) = (P * b * x(i) * (L^2 - x(i)^2 - b^2)) / (6 * E * I * L);
else
deflection(i) = (P * b * (L^2 - x(i)^2 - b^2) * x(i) + P * (x(i) - a)^3) / (6 * E * I * L);
end
end
title_str = '集中荷载作用下的简支梁变形';
case 'uniform'
% 均布荷载作用下的简支梁变形
deflection = (w * x .* (L^3 - 2 * L * x.^2 + x.^3)) / (24 * E * I);
title_str = '均布荷载作用下的简支梁变形';
case 'triangular'
% 三角形分布荷载作用下的简支梁变形
% 假设三角形荷载从左侧最大w到右侧为0
deflection = (w * x .* (7 * L^4 - 10 * L^2 * x.^2 + 3 * x.^4)) / (360 * E * I * L);
title_str = '三角形分布荷载作用下的简支梁变形';
otherwise
error('不支持的荷载类型。请选择 ''point'', ''uniform'', 或 ''triangular''');
end
% 计算转角和弯矩
slope = gradient(deflection, x(2)-x(1));
moment = -E * I * gradient(slope, x(2)-x(1));
shear = gradient(moment, x(2)-x(1));
%% 3. 有限元数值解
fprintf('进行有限元分析...\n');
% 调用有限元函数
[fe_deflection, fe_slope, fe_moment, fe_shear, nodes] = beam_fem_analysis(L, E, I, rho, A, load_type, P, w, a, n_elements);
%% 4. 结果可视化
fprintf('绘制结果...\n');
% 绘制变形图
figure('Position', [100, 100, 1200, 800]);
subplot(2,2,1);
plot(x, deflection * 1000, 'b-', 'LineWidth', 2);
hold on;
plot(nodes, fe_deflection * 1000, 'ro--', 'LineWidth', 1, 'MarkerSize', 4);
xlabel('位置 (m)');
ylabel('挠度 (mm)');
title('梁的挠度曲线');
legend('解析解', '有限元解', 'Location', 'Best');
grid on;
% 绘制转角图
subplot(2,2,2);
plot(x, slope * 1000, 'b-', 'LineWidth', 2);
hold on;
plot(nodes, fe_slope * 1000, 'ro--', 'LineWidth', 1, 'MarkerSize', 4);
xlabel('位置 (m)');
ylabel('转角 (mrad)');
title('梁的转角曲线');
legend('解析解', '有限元解', 'Location', 'Best');
grid on;
% 绘制弯矩图
subplot(2,2,3);
plot(x, moment / 1000, 'b-', 'LineWidth', 2);
hold on;
plot(nodes, fe_moment / 1000, 'ro--', 'LineWidth', 1, 'MarkerSize', 4);
xlabel('位置 (m)');
ylabel('弯矩 (kN·m)');
title('梁的弯矩图');
legend('解析解', '有限元解', 'Location', 'Best');
grid on;
% 绘制剪力图
subplot(2,2,4);
plot(x, shear / 1000, 'b-', 'LineWidth', 2);
hold on;
plot(nodes, fe_shear / 1000, 'ro--', 'LineWidth', 1, 'MarkerSize', 4);
xlabel('位置 (m)');
ylabel('剪力 (kN)');
title('梁的剪力图');
legend('解析解', '有限元解', 'Location', 'Best');
grid on;
sgtitle(sprintf('%s (L = %.1f m, E = %.0f GPa, I = %.2e m^4)', ...
title_str, L, E/1e9, I), 'FontSize', 14, 'FontWeight', 'bold');
% 绘制梁的变形形状
figure;
plot_beam_deformation(nodes, fe_deflection, L, title_str);
%% 5. 结果输出
fprintf('\n===== 计算结果摘要 =====\n');
fprintf('最大挠度: %.4f mm\n', max(deflection) * 1000);
fprintf('最大转角: %.4f mrad\n', max(abs(slope)) * 1000);
fprintf('最大弯矩: %.4f kN·m\n', max(abs(moment)) / 1000);
fprintf('最大剪力: %.4f kN\n', max(abs(shear)) / 1000);
fprintf('\n有限元分析结果:\n');
fprintf('最大挠度: %.4f mm\n', max(fe_deflection) * 1000);
fprintf('最大转角: %.4f mrad\n', max(abs(fe_slope)) * 1000);
fprintf('最大弯矩: %.4f kN·m\n', max(abs(fe_moment)) / 1000);
fprintf('最大剪力: %.4f kN\n', max(abs(fe_shear)) / 1000);
%% 6. 有限元分析函数
function [deflection, slope, moment, shear, nodes] = beam_fem_analysis(L, E, I, rho, A, load_type, P, w, a, n_elements)
% 简支梁有限元分析
% 节点和单元定义
n_nodes = n_elements + 1;
nodes = linspace(0, L, n_nodes);
h = L / n_elements; % 单元长度
% 初始化全局刚度矩阵和荷载向量
K = zeros(2 * n_nodes);
F = zeros(2 * n_nodes, 1);
% 组装全局刚度矩阵
for i = 1:n_elements
% 单元刚度矩阵
ke = (E * I / h^3) * [
12, 6*h, -12, 6*h;
6*h, 4*h^2, -6*h, 2*h^2;
-12, -6*h, 12, -6*h;
6*h, 2*h^2, -6*h, 4*h^2
];
% 组装到全局矩阵
dofs = [2*i-1, 2*i, 2*i+1, 2*i+2];
K(dofs, dofs) = K(dofs, dofs) + ke;
end
% 组装荷载向量
switch load_type
case 'point'
% 找到最接近荷载位置的节点
[~, node_idx] = min(abs(nodes - a));
F(2*node_idx-1) = P; % 施加集中力
case 'uniform'
% 均布荷载等效节点荷载
for i = 1:n_elements
fe = w * h * [1/2; h/12; 1/2; -h/12];
dofs = [2*i-1, 2*i, 2*i+1, 2*i+2];
F(dofs) = F(dofs) + fe;
end
case 'triangular'
% 三角形分布荷载等效节点荷载
for i = 1:n_elements
x1 = nodes(i);
x2 = nodes(i+1);
w1 = w * (1 - x1/L); % 左侧荷载强度
w2 = w * (1 - x2/L); % 右侧荷载强度
fe = h * [ (3*w1 + w2)/12;
h*(2*w1 + w2)/30;
(w1 + 3*w2)/12;
-h*(w1 + 2*w2)/30 ];
dofs = [2*i-1, 2*i, 2*i+1, 2*i+2];
F(dofs) = F(dofs) + fe;
end
end
% 应用边界条件 (简支梁: 左端挠度和转角为0，右端挠度为0)
fixed_dofs = [1, 2, 2*n_nodes-1]; % 左端挠度和转角，右端挠度
free_dofs = setdiff(1:2*n_nodes, fixed_dofs);
% 求解方程组
U = zeros(2*n_nodes, 1);
U(free_dofs) = K(free_dofs, free_dofs) \ F(free_dofs);
% 提取结果
deflection = U(1:2:end); % 挠度
slope = U(2:2:end); % 转角
% 计算弯矩和剪力
moment = zeros(n_nodes, 1);
shear = zeros(n_nodes, 1);
for i = 1:n_elements
% 单元节点位移
ue = U([2*i-1, 2*i, 2*i+1, 2*i+2]);
% 单元弯矩和剪力
for xi = [0, h] % 在单元两端计算
N1 = 1 - 3*(xi/h)^2 + 2*(xi/h)^3;
N2 = xi * (1 - xi/h)^2;
N3 = 3*(xi/h)^2 - 2*(xi/h)^3;
N4 = xi * ((xi/h)^2 - xi/h);
N = [N1, N2, N3, N4];
dN1 = -6*xi/h^2 + 6*xi^2/h^3;
dN2 = 1 - 4*xi/h + 3*xi^2/h^2;
dN3 = 6*xi/h^2 - 6*xi^2/h^3;
dN4 = -2*xi/h + 3*xi^2/h^2;
dN = [dN1, dN2, dN3, dN4];
d2N1 = -6/h^2 + 12*xi/h^3;
d2N2 = -4/h + 6*xi/h^2;
d2N3 = 6/h^2 - 12*xi/h^3;
d2N4 = -2/h + 6*xi/h^2;
d2N = [d2N1, d2N2, d2N3, d2N4];
% 弯矩 M = EI * d²w/dx²
M = E * I * d2N * ue;
% 剪力 V = dM/dx = EI * d³w/dx³
d3N1 = 12/h^3;
d3N2 = 6/h^2;
d3N3 = -12/h^3;
d3N4 = 6/h^2;
d3N = [d3N1, d3N2, d3N3, d3N4];
V = E * I * d3N * ue;
% 存储结果
node_idx = i + (xi > 0);
moment(node_idx) = M;
shear(node_idx) = V;
end
end
end
%% 7. 梁变形形状绘制函数
function plot_beam_deformation(nodes, deflection, L, title_str)
% 绘制梁的变形形状
% 放大变形以便可视化
scale_factor = 50;
deflected_shape = nodes + deflection * scale_factor;
% 绘制原始梁和变形后的梁
figure;
plot(nodes, zeros(size(nodes)), 'k-', 'LineWidth', 3);
hold on;
plot(deflected_shape, zeros(size(nodes)), 'r-', 'LineWidth', 2);
% 绘制连接线
for i = 1:length(nodes)
plot([nodes(i), deflected_shape(i)], [0, 0], 'b--', 'LineWidth', 0.5);
end
% 设置图形属性
axis equal;
xlim([-0.1*L, 1.1*L]);
ylim([-0.1*L, 0.1*L]);
xlabel('位置 (m)');
title(sprintf('%s\n变形形状 (放大 %d 倍)', title_str, scale_factor));
legend('原始梁', '变形后梁', 'Location', 'Best');
grid on;
% 添加支座符号
hold on;
plot(0, 0, '^', 'MarkerSize', 15, 'MarkerFaceColor', 'g');
plot(L, 0, '^', 'MarkerSize', 15, 'MarkerFaceColor', 'g');
% 添加荷载示意
if contains(title_str, '集中荷载')
text(a, -0.02*L, '↓ P', 'FontSize', 12, 'FontWeight', 'bold', ...
'HorizontalAlignment', 'center');
elseif contains(title_str, '均布荷载')
for x = linspace(0, L, 10)
plot([x, x], [0, -0.02*L], 'k-', 'LineWidth', 1.5);
end
text(L/2, -0.04*L, 'w (N/m)', 'FontSize', 12, 'FontWeight', 'bold', ...
'HorizontalAlignment', 'center');
elseif contains(title_str, '三角形')
for x = linspace(0, L, 10)
h = 0.02 * L * (1 - x/L);
plot([x, x], [0, -h], 'k-', 'LineWidth', 1.5);
end
text(L/2, -0.04*L, 'w (N/m)', 'FontSize', 12, 'FontWeight', 'bold', ...
'HorizontalAlignment', 'center');
end
end
%% 8. 交互式参数调整函数（可选）
function interactive_beam_analysis()
% 交互式简支梁分析界面
fig = figure('Position', [100, 100, 800, 600], 'Name', '简支梁变形分析');
% 创建控件
uicontrol('Style', 'text', 'Position', [20, 560, 100, 20], 'String', '梁长度 (m):');
length_edit = uicontrol('Style', 'edit', 'Position', [120, 560, 60, 20], 'String', '10');
uicontrol('Style', 'text', 'Position', [20, 530, 100, 20], 'String', '弹性模量 (GPa):');
E_edit = uicontrol('Style', 'edit', 'Position', [120, 530, 60, 20], 'String', '200');
uicontrol('Style', 'text', 'Position', [20, 500, 100, 20], 'String', '惯性矩 (m^4):');
I_edit = uicontrol('Style', 'edit', 'Position', [120, 500, 60, 20], 'String', '8.33e-6');
uicontrol('Style', 'text', 'Position', [20, 470, 100, 20], 'String', '荷载类型:');
load_popup = uicontrol('Style', 'popupmenu', 'Position', [120, 470, 100, 20], ...
'String', {'集中荷载', '均布荷载', '三角形荷载'});
uicontrol('Style', 'text', 'Position', [20, 440, 100, 20], 'String', '荷载大小:');
load_edit = uicontrol('Style', 'edit', 'Position', [120, 440, 60, 20], 'String', '10000');
uicontrol('Style', 'text', 'Position', [20, 410, 100, 20], 'String', '荷载位置 (m):');
pos_edit = uicontrol('Style', 'edit', 'Position', [120, 410, 60, 20], 'String', '3');
uicontrol('Style', 'pushbutton', 'Position', [20, 370, 160, 30], 'String', '计算', ...
'Callback', @calculate_callback);
% 回调函数
function calculate_callback(~, ~)
% 获取参数值
L = str2double(get(length_edit, 'String'));
E = str2double(get(E_edit, 'String')) * 1e9;
I = str2double(get(I_edit, 'String'));
load_val = str2double(get(load_edit, 'String'));
a = str2double(get(pos_edit, 'String'));
% 获取荷载类型
load_types = {'point', 'uniform', 'triangular'};
load_idx = get(load_popup, 'Value');
load_type = load_types{load_idx};
% 执行计算
main_beam_analysis(L, E, I, load_type, load_val, a);
end
end
%% 9. 主分析函数（供交互式界面调用）
function main_beam_analysis(L, E, I, load_type, load_val, a)
% 简化的主分析函数
% 设置默认参数
rho = 7850;
A = 0.02;
n_elements = 20;
% 根据荷载类型分配参数
if strcmp(load_type, 'point')
P = load_val;
w = 1000; % 默认值，不会使用
else
P = 10000; % 默认值，不会使用
w = load_val;
end
% 创建位置向量
x = linspace(0, L, 1000);
% 根据荷载类型计算变形
switch load_type
case 'point'
b = L - a; % 荷载到右端的距离
deflection = zeros(size(x));
for i = 1:length(x)
if x(i) <= a
deflection(i) = (P * b * x(i) * (L^2 - x(i)^2 - b^2)) / (6 * E * I * L);
else
deflection(i) = (P * b * (L^2 - x(i)^2 - b^2) * x(i) + P * (x(i) - a)^3) / (6 * E * I * L);
end
end
title_str = '集中荷载作用下的简支梁变形';
case 'uniform'
deflection = (w * x .* (L^3 - 2 * L * x.^2 + x.^3)) / (24 * E * I);
title_str = '均布荷载作用下的简支梁变形';
case 'triangular'
deflection = (w * x .* (7 * L^4 - 10 * L^2 * x.^2 + 3 * x.^4)) / (360 * E * I * L);
title_str = '三角形分布荷载作用下的简支梁变形';
end
% 有限元分析
[fe_deflection, fe_slope, fe_moment, fe_shear, nodes] = beam_fem_analysis(...
L, E, I, rho, A, load_type, P, w, a, n_elements);
% 绘制结果
figure('Position', [100, 100, 1000, 800]);
subplot(2,2,1);
plot(x, deflection * 1000, 'b-', 'LineWidth', 2);
hold on;
plot(nodes, fe_deflection * 1000, 'ro--', 'LineWidth', 1, 'MarkerSize', 4);
xlabel('位置 (m)');
ylabel('挠度 (mm)');
title('梁的挠度曲线');
legend('解析解', '有限元解', 'Location', 'Best');
grid on;
% 计算转角和弯矩
slope = gradient(deflection, x(2)-x(1));
moment = -E * I * gradient(slope, x(2)-x(1));
subplot(2,2,2);
plot(x, slope * 1000, 'b-', 'LineWidth', 2);
hold on;
plot(nodes, fe_slope * 1000, 'ro--', 'LineWidth', 1, 'MarkerSize', 4);
xlabel('位置 (m)');
ylabel('转角 (mrad)');
title('梁的转角曲线');
legend('解析解', '有限元解', 'Location', 'Best');
grid on;
subplot(2,2,3);
plot(x, moment / 1000, 'b-', 'LineWidth', 2);
hold on;
plot(nodes, fe_moment / 1000, 'ro--', 'LineWidth', 1, 'MarkerSize', 4);
xlabel('位置 (m)');
ylabel('弯矩 (kN·m)');
title('梁的弯矩图');
legend('解析解', '有限元解', 'Location', 'Best');
grid on;
subplot(2,2,4);
shear = gradient(moment, x(2)-x(1));
plot(x, shear / 1000, 'b-', 'LineWidth', 2);
hold on;
plot(nodes, fe_shear / 1000, 'ro--', 'LineWidth', 1, 'MarkerSize', 4);
xlabel('位置 (m)');
ylabel('剪力 (kN)');
title('梁的剪力图');
legend('解析解', '有限元解', 'Location', 'Best');
grid on;
sgtitle(sprintf('%s (L = %.1f m, E = %.0f GPa, I = %.2e m^4)', ...
title_str, L, E/1e9, I), 'FontSize', 14, 'FontWeight', 'bold');
% 输出结果
fprintf('\n最大挠度: %.4f mm\n', max(deflection) * 1000);
fprintf('有限元最大挠度: %.4f mm\n', max(fe_deflection) * 1000);
end
% 运行交互式界面（取消注释以使用）
% interactive_beam_analysis();
程序说明
这个MATLAB程序提供了简支梁在荷载作用下的变形计算，包括以下功能：
1. 解析解计算
- 支持三种常见荷载类型：集中荷载、均布荷载和三角形分布荷载
- 计算梁的挠度、转角、弯矩和剪力
- 基于梁弯曲理论的基本公式
2. 有限元数值解
- 使用梁单元进行有限元分析
- 可以处理各种荷载条件和边界条件
- 提供与解析解的对比
3. 可视化功能
- 绘制梁的挠度曲线、转角曲线、弯矩图和剪力图
- 显示梁的变形形状（放大效果）
- 比较解析解和有限元解
4. 交互式界面（可选）
- 提供图形用户界面，方便参数调整
- 实时计算和显示结果
使用方法
- 直接运行程序：程序将使用默认参数进行计算并显示结果
- 修改参数：在代码开头部分修改梁的参数和荷载条件
- 使用交互式界面：取消最后一行注释，运行
interactive_beam_analysis()函数
理论基础
简支梁弯曲基本方程
简支梁的弯曲行为由以下微分方程描述：
d²/dx² (EI d²w/dx²) = q(x)
其中：
- w(x) 是梁的挠度
- E 是弹性模量
- I 是截面惯性矩
- q(x) 是分布荷载
参考代码 简支梁荷载作用下的变形计算 www.youwenfan.com/contentcnj/63960.html
有限元方法
程序使用欧拉-伯努利梁单元进行有限元分析，每个节点有两个自由度：挠度和转角。单元刚度矩阵为：
k = (EI/L³) * [12 6L -12 6L
6L 4L² -6L 2L²
-12 -6L 12 -6L
6L 2L² -6L 4L²]
扩展功能
程序可以进一步扩展以包括：
- 更多荷载类型：梯形荷载、任意分布荷载等
- 多种边界条件：固支梁、悬臂梁等
- 材料非线性：考虑塑性变形
- 几何非线性：大变形分析
- 动力分析：考虑振动和动力响应
这个程序提供了一个完整的简支梁变形分析工具，既可以用于教学演示，也可以用于工程计算。
浙公网安备 33010602011771号