工程计算与分析课程报告-Matlab

如下图所示的单自由度系统: m=20kg, c=30 Ns/m, and k=3000 N/m.

1、应用数学、机械等理论知识推导出以下系统的数学模型：

其中：

x坐标如上图所示，取其稳态点为零点。

M为重物质量

K为弹簧弹性系数

C为阻尼系数

F(t)为施加在重物上的外激励力量

2、计算当外激励为周期性激励：F(t)=40cos(10t)(N),而所有的初值条件为零，即 x(0)=dx(0)/dt=0, 应用MATLAB计算系统的动态响应x(t)=?并绘制图形。（计算时间域为[0 5]秒）。c=30 Ns/m, and k=3000 N/m

% 定义系统参数  
m = 20.0;  % 质量，单位：kg  
c = 30.0;  % 阻尼系数，单位：Ns/m  
k = 3000.0;  % 弹簧刚度，单位：N/m  
% 定义外部激励函数  
F = @(t) 40*cos(10*t);  
% 定义一阶方程系统的匿名函数  
system_equation = @(t, y) [y(2); (F(t) - c*y(2) - k*y(1)) / m];  
% 初始条件  
y0 = [0; 0];  % 初始位移和初始速度  
% 时间范围  
tspan = [0 5];  
% 使用ode45求解  
[t, y] = ode45(system_equation, tspan, y0);  
% 提取位移  
x = y(:, 1);  
figure;  
plot(t, x);  
xlabel('T');  
title('动态响应');  

3、当外激励为脉冲激励时, 计算系统的稳态响应x(t)=?并绘制图形。（计算时间域为[0 5]秒）。 IMPLUSE

c=30 Ns/m, and k=3000 N/m
即：考虑具有粘性阻尼的弹簧 - 质量系统在 t = 0 时受到一个单位冲量作用:

% 定义系统参数  
m = 20.0;  % 质量，单位：kg  
c = 30.0;  % 阻尼系数，单位：Ns/m  
k = 3000.0;  % 弹簧刚度，单位：N/m  
  
% 定义ODE的匿名函数  
ode_system = @(t, y) [y(2); (-c*y(2) - k*y(1)) / m];  
  
% 初始条件（假设脉冲导致了一个初始速度）  
y0 = [0; 1];  % 初始位移为0，初始速度为1 m/s  
  
% 时间范围  
tspan = [0 5];  
  
% 使用ode45求解  
[t, y] = ode45(ode_system, tspan, y0);  
% 提取位移  
x = y(:, 1);  
% 绘制图形  
figure;  
plot(t, x);  
xlabel('T');  
title('脉冲激励响应');  

4、当外激励为阶跃激励时, 计算系统的稳态响应x(t)=?并绘制图形。（计算时间域为[0 5]秒）。 STEP

即：F(t)=20 N

c=30 Ns/m, and k=3000 N/m

% 参数  
m = 20.0;  % 质量，单位：kg  
c = 30.0;  % 阻尼系数，单位：Ns/m  
k = 3000.0;  % 弹簧刚度，单位：N/m  
F = 20;    % 阶跃激励的大小，单位：N  
% 定义ODE的匿名函数，考虑阶跃激励  
ode_system = @(t, y) [y(2); (F * (t > 0) - c*y(2) - k*y(1)) / m];  
% 初始条件  
y0 = [0; 0];  % 初始位移和速度都为0  
% 时间范围  
tspan = [0 5];  
% 使用ode45求解  
[t, y] = ode45(ode_system, tspan, y0);  
  
% 提取位移  
x = y(:, 1);  
  
% 绘制图形  
figure;  
plot(t, x);  
xlabel('T');  
title('阶跃激励响应');  

5、当外激励为阶跃激励（即：F(t)=20 N）时, 改变不同弹簧弹性系数k、阻尼系数c的值（见下述条件），计算系统的稳态响应x(t)=?并绘制图形，分析结果。（计算时间域为[0 5]秒）。

条件：

c=10；15 ；20 ；25；30；35；40；Ns/m,

 k=1000；2000；3000；4000  N/m

% 参数定义  
F = 20; % 阶跃激励的大小，单位：N  
m = 20.0; % 质量，单位：kg  
tspan = [0 5]; % 时间范围  
y0 = [0; 0]; % 初始条件：位移和速度都为0  
  
% 定义ODE函数（不需要单独的文件，直接作为匿名函数）  
ode_system = @(t, y, F, m, c, k) [y(2); (F - c*y(2) - k*y(1)) / m];  
  
% 初始化绘图参数  
figure; % 创建一个新的图形窗口  
hold on; % 保持当前图形，以便在同一窗口上绘制多个曲线  
grid on; % 打开网格  
xlabel('Time (s)'); % x轴标签  
ylabel('Displacement (m)'); % y轴标签  
title('Response of the System to a Step Input'); % 图形标题  
  
% 绘制不同k和c值的响应  
colors = 'rbgcykm'; % 颜色向量，用于区分不同曲线  
k_values = [1000, 2000, 3000, 4000]; % 弹簧弹性系数  
c_values = [10, 15, 20, 25, 30, 35, 40]; % 阻尼系数  
  
subplot(221)
k = 1000;
for i_c = 1:length(c_values)
    c = c_values(i_c);
    % 使用ode45求解ODE
    [t, y] = ode45(@(t, y) ode_system(t, y, F, m, c, k), tspan, y0);
    x = y(:, 1); % 提取位移
    % 绘制位移响应曲线
    plot(t, x, 'LineWidth', 2, 'Color', colors(mod(i_k + i_c, length(colors)) + 1)); % 使用不同的颜色和线宽
    hold on;
end
title("K=1000")
legend('c=10','c=15','c=20','c=25','c=30','c=35','c=40')

subplot(222)
k = 2000;
for i_c = 1:length(c_values)
    c = c_values(i_c);
    % 使用ode45求解ODE
    [t, y] = ode45(@(t, y) ode_system(t, y, F, m, c, k), tspan, y0);
    x = y(:, 1); % 提取位移
    % 绘制位移响应曲线
    plot(t, x, 'LineWidth', 2, 'Color', colors(mod(i_k + i_c, length(colors)) + 1)); % 使用不同的颜色和线宽
    hold on;
end
title("K=2000")
legend('c=10','c=15','c=20','c=25','c=30','c=35','c=40')

subplot(223)
k = 3000;
for i_c = 1:length(c_values)
    c = c_values(i_c);
    % 使用ode45求解ODE
    [t, y] = ode45(@(t, y) ode_system(t, y, F, m, c, k), tspan, y0);
    x = y(:, 1); % 提取位移
    % 绘制位移响应曲线
    plot(t, x, 'LineWidth', 2, 'Color', colors(mod(i_k + i_c, length(colors)) + 1)); % 使用不同的颜色和线宽
    hold on;
end
title("K=3000")
legend('c=10','c=15','c=20','c=25','c=30','c=35','c=40')

subplot(224)
k = 4000;
for i_c = 1:length(c_values)
    c = c_values(i_c);
    % 使用ode45求解ODE
    [t, y] = ode45(@(t, y) ode_system(t, y, F, m, c, k), tspan, y0);
    x = y(:, 1); % 提取位移
    % 绘制位移响应曲线
    plot(t, x, 'LineWidth', 2, 'Color', colors(mod(i_k + i_c, length(colors)) + 1)); % 使用不同的颜色和线宽
        hold on;
end
title("K=4000")
legend('c=10','c=15','c=20','c=25','c=30','c=35','c=40')