对单自由度系统的强制振动进行数值方法研究的实现(Matlab代码实现)

 💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势:🌞🌞🌞博客内容尽量做到思维缜密,逻辑清晰,为了方便读者。

⛳️座右铭:行百里者,半于九十。

📋📋📋本文目录如下:🎁🎁🎁

目录

💥1 概述

单自由度系统强制振动的数值方法研究

一、引言

二、理论基础与数学模型

1. 单自由度系统动力学方程

2. 强制振动的解析解

3. 共振现象

三、数值求解方法

1. 有限差分法

2. Runge-Kutta法(四阶)

四、方法对比与误差分析

五、数值仿真案例(MATLAB实现)

案例1:阻尼系统受迫振动

案例2:共振条件下的有限差分法

六、研究文档标准框架

七、结论

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现


💥1 概述

单自由度系统强制振动的数值方法研究

摘要
本文研究了受谐波力激励的阻尼单自由度(SDOF)系统的精确解,并对其与Matlab内置函数ode45、中心差分法、Newmark法和四阶Runge-Kutta法提供的数值解进行了比较。这些数值解的实现基于S. Rao的著作。为了实现这些数值方法,我们编写了RK4.m函数(使用四阶Runge-Kutta法数值求解阻尼系统的运动方程)、Newmark.m函数(使用Newmark法数值求解阻尼系统的运动方程)、CentDiff.m函数(使用中心差分法数值求解阻尼系统的运动方程)以及Matlab LiveScript Documentation.mlx(用于文档说明)。

内容
本文对受谐波力激励的阻尼单自由度系统进行了深入研究,旨在探索数值方法对系统振动的准确性和稳定性。我们首先介绍了系统的精确解,并将其作为基准与其他数值解进行比较。通过编写并实现RK4.m、Newmark.m和CentDiff.m等函数,我们成功地使用四阶Runge-Kutta法、Newmark法和中心差分法对系统的运动方程进行了数值求解。这些数值方法的实现基于文献【2】.

内容

RK4.m函数,使用四阶Runge-Kutta法数值求解阻尼系统的运动方程
Newmark.m函数,使用Newmark法数值求解阻尼系统的运动方程
CentDiff.m函数,使用中心差分法数值求解阻尼系统的运动方程
Matlab LiveScript Documentation.mlx,用于文档说明


一、引言

单自由度系统是振动理论的基础模型,其动力学特性为多自由度系统和连续体分析提供了重要参考。强制振动指系统在周期性外力作用下的持续振动,其核心问题包括共振机制、稳态响应特性及数值求解方法。本文将从数学模型构建、常用数值方法(有限差分法、Runge-Kutta法)的原理与实现、误差分析及仿真案例等方面展开系统性研究。


二、理论基础与数学模型
1. 单自由度系统动力学方程

系统由质量块mm、弹簧刚度kk和阻尼器cc组成,受外力F(t)F(t)作用时,运动方程为:

其中:

  • u(t)为位移响应
  • ω为激励频率
  • 无阻尼固有频率
2. 强制振动的解析解

方程解分为瞬态响应和稳态响应:

3. 共振现象

当ω≈ωn时,系统进入共振区,振幅急剧增大(无阻尼时理论幅值趋向无穷大)。实际工程中需通过阻尼设计抑制共振危害。


三、数值求解方法
1. 有限差分法
  • 原理:用差分近似代替微分,将连续方程离散化。例如,中心差分法:

    代入原方程可递推求解各时间步位移。

  • 适用条件:时间步长需满足Δt<2/ωn(稳定性条件),否则解会发散。

  • 误差来源:截断误差(O(Δt2))、舍入误差累积,尤其在长时间仿真中显著。

2. Runge-Kutta法(四阶)
  • 原理:将二阶微分方程降阶为两个一阶方程,通过加权平均斜率计算下一步状态:

    四阶RK法的递推公式为:

    其中k1,k2,k3,k4为不同斜率估计值。

  • 优势:自启动、较高精度(局部截断误差O(Δt5)),适用于非线性或变刚度系统。

  • 适用性:对刚性方程需减小步长或采用隐式变体(如RKF45)。


四、方法对比与误差分析
方法精度稳定性计算效率适用场景
有限差分法中等条件稳定线性系统、实时仿真
Runge-Kutta法较好(显式)中等非线性系统、高精度需求
  • 误差控制策略
    • 自适应步长调整:根据局部误差估计动态调整ΔtΔt,平衡精度与计算量。
    • 能量守恒检验:验证数值解是否满足机械能守恒(无阻尼自由振动时)。

五、数值仿真案例(MATLAB实现)
案例1:阻尼系统受迫振动
% 参数设置
m = 1; c = 0.2; k = 10; F0 = 1; omega = 3;
tspan = [0 20]; u0 = [0; 0];

% 四阶Runge-Kutta法求解
[t, u] = ode45(@(t,u) [u(2); (F0*sin(omega*t) - c*u(2) - k*u(1))/m], tspan, u0);

% 可视化
plot(t, u(:,1), 'LineWidth', 1.5);
xlabel('Time (s)'); ylabel('Displacement (m)');
title('Forced Vibration Response with RK4');

结果分析:初期瞬态振动逐渐衰减,稳态响应幅值与理论解吻合。

案例2:共振条件下的有限差分法
dt = 0.01; t = 0:dt:10;
u = zeros(size(t)); v = zeros(size(t));
u(1) = 0; v(1) = 0;

for i = 2:length(t)-1
    a = (F0*sin(omega*t(i)) - c*v(i) - k*u(i)) / m;
    v(i+1) = v(i) + a*dt;
    u(i+1) = u(i) + v(i)*dt;
end

注意:需验证步长是否满足稳定性条件,避免数值发散。


六、研究文档标准框架
  1. 引言:研究背景、意义及文献综述。
  2. 理论基础:系统方程推导、解析解与共振机理。
  3. 数值方法:算法原理、实现步骤及代码示例。
  4. 案例分析:仿真设置、结果对比与误差讨论。
  5. 结论:方法优劣总结、工程应用建议及未来研究方向。

七、结论

单自由度系统强制振动的数值研究需综合解析理论与数值方法。有限差分法适合快速实现线性系统仿真,而Runge-Kutta法在高精度和非线性问题中更具优势。未来研究可探索机器学习加速算法或多物理场耦合模型的数值求解。

📚2 运行结果

部分代码:

function [y] = CentDiff(F,M,K,C,dt,x0,v0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GOAl: Numerical solve the equation of motion of a damped system


% INPUT
% F : vector  -- size: [1x N] -- Time series representinf the time history of the load. 
% M : scalar  -- size: [1 x 1] -- Modal mass
% K : scalar  -- size: [1 x 1] -- Modal stifness
% C : scalar  -- size: [1 x 1] -- Modal damping
% dt : scalar  -- size: [1 x 1] -- time step
% x0 : scalar  -- size: [1 x 1] -- initial displacement
% v0 : scalar  -- size: [1 x 1] -- initial velocity

% OUTPUT
% y: time history of the system response to the load

% author: Etienne Cheyet  ----  last updated: 07/11/2015
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
% Initialisation
N = size(F,2);
% preallocation
y = zeros(size(F));

% initial acceleration
a0 = M\(F(1)-C.*v0-K.*x0);
% initialisation of y (first 2 values).
y0 = x0-dt.*v0+dt^2/2*a0;
y(:,1) = x0;
A = (M./dt.^2+C./(2*dt));
B = ((2*M./dt.^2-K).*y(:,1)+(C./(2*dt)-M./dt.^2).*y0+F(:,1));
y(:,2) = A\B;

% For the rest of integration points
for ii=2:N-1,
    A = (M./dt.^2+C./(2*dt));
    B = ((2*M./dt.^2-K).*y(:,ii)+(C./(2*dt)-M./dt.^2).*y(:,ii-1)+F(:,ii)); 
    y(:,ii+1) = A\B;
end

end

🎉3 参考文献

文章中一些内容引自网络,会注明出处或引用为参考文献,难免有未尽之处,如有不妥,请随时联系删除。

[1]林恩强.流体—结构耦合数值方法研究及其应用[D].昆明理工大学,2007.

[2]许媛欣.方柱绕流涡激振动及控制方法的数值研究[D].天津大学,2015.

[3]张予,刘晓薇,等.基于Matlab的单自由度振动系统的数字仿真实验[J].华东交通大学学报, 2002.

🌈4 Matlab代码实现

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

荔枝科研社

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值