动刚度，在之前的文章，一篇文章搞清楚——橡胶的粘弹性，已经大致讲清楚了，如果有人给了你个件儿，让你给他讲清楚这个件的动刚度是什么，你需要给他展示哪些图片，以及动刚度中，储能刚度、损耗刚度与损耗因子之间的关系，如何通过正弦加载推导出来，等等，这篇文章不再重复。

本文着重讲，对于一个构件，动刚度应如何获取。

动刚度 是指隔振器在受到周期性交变力作用时，其力振幅与位移振幅的复数比值。它不是一个简单的标量，而是一个与频率和振幅强相关的复数量。

因此，动刚度在展示的时候，一般都是横坐标是频率，纵坐标是动刚度的一幅图。

于是我们现在的目的，就是拿着所测量得到的时域上的力-位移曲线，获取其在频域内的动刚度曲线。

1.试验获取动刚度曲线

这是获取橡胶隔振器真实动刚度的标准途径。通常在液压伺服疲劳试验机或专用的动态材料试验机上完成。按照之前文章我们讨论的动刚度曲线结果图，直接测量其不同频率下的动刚度即可。

基本测试流程：

1. 试样安装：将橡胶隔振器按照实际工作状态（包括预压缩量）安装在试验机上。
2. 施加激励：对试件施加正弦波位移（或力）激励。
• x(t) = X * sin(ωt) （位移控制模式更常见）
3. 测量响应：精确测量产生的作用力响应。
• F(t) = F * sin(ωt + δ)，其中 δ 是力相对于位移的相位角。
4. 数据处理：
• 幅值比：|K_d(ω)| = F / X （这是动刚度的幅值）
• 储能刚度：K'(ω) = |K_d(ω)| * cosδ
• 耗能刚度：K''(ω) = |K_d(ω)| * sinδ
• 损耗因子：η = tanδ = K''(ω) / K'(ω)

关键测试条件：

• 频率扫描：在关心的频率范围（如1-200Hz或更高）内，固定振幅，测量动刚度随频率的变化。这是分析频率依赖性的关键。
• 振幅扫描：在固定频率下，改变激励振幅，测量动刚度随振幅的变化。用于评估非线性程度。
• 预载荷：必须模拟隔振器在实际工作中的静态压缩量（预载），因为预载对动刚度影响极大。
• 温度控制：橡胶性能对温度敏感，测试应在可控温度下进行，或注明测试温度。

输出结果：一系列曲线图，如K'-频率曲线、η-频率曲线。

2.任意力-位移响应的动刚度

液压伺服疲劳试验机测试是获取动刚度的常规方法，但有时候我们想要知道构件在工作过程中，承受各种载荷下，的瞬时动刚度，如何求解？

三种方法，椭圆拟合、傅里叶变换法、切线刚度法。

其中，椭圆拟合法只适用于正弦激励，傅里叶变换法最合理，切线刚度法无法计算动刚度，只是应力应变曲线瞬时的切线斜率。

2.1 椭圆拟合法

function[K_prime, K_double_prime, loss_factor] = calculate_dynamic_stiffness_ellipse(force, displacement)% 方法1：椭圆拟合法% 适用于近似正弦激励的情况% 1. 绘制滞回环并移除可能的偏移    F = force - mean(force);    X = displacement - mean(displacement);% 2. 椭圆拟合（使用最小二乘法）% 椭圆方程：A*x^2 + B*x*y + C*y^2 + D*x + E*y = 1    D = [X.^2, X.*F, F.^2, X, F];    params = D \ ones(length(X), 1);    A = params(1); B = params(2); C = params(3);% 3. 计算椭圆几何参数% 长轴角度    theta = 0.5 * atan2(B, A - C);% 4. 旋转到主轴坐标系    R = [cos(theta), -sin(theta); sin(theta), cos(theta)];    coords = R * [X'; F'];    X_rot = coords(1, :)';    F_rot = coords(2, :)';% 5. 计算动态参数% 长轴斜率 = 储能刚度 K'    K_prime = (max(F_rot) - min(F_rot)) / (max(X_rot) - min(X_rot));% 椭圆面积 = π * a * b * e，其中e为椭圆度    a = (max(X_rot) - min(X_rot)) / 2;  % 半长轴    b = (max(F_rot) - min(F_rot)) / 2;  % 半短轴% 能量耗散 = 椭圆面积    energy_dissipation = pi * a * b;% 耗能刚度 K'' (近似)% 基于：能量耗散 = π * X0^2 * ω * η * K'，其中η=tanδ    X0 = max(abs(X));  % 位移幅值    K_double_prime = energy_dissipation / (pi * X0^2);% 损耗因子    loss_factor = K_double_prime / K_prime;% 绘图figure('Position', [100, 100, 1200, 400]);    subplot(1,3,1);plot(X, F, 'b.', 'MarkerSize', 8); hold on;% 绘制拟合椭圆    theta_ellipse = linspace(0, 2*pi, 100);    x_ellipse = a * cos(theta_ellipse);    y_ellipse = b * sin(theta_ellipse);% 旋转回原坐标系    coords_ellipse = R' * [x_ellipse; y_ellipse];plot(coords_ellipse(1,:), coords_ellipse(2,:), 'r-', 'LineWidth', 2);    xlabel('位移 (m)'); ylabel('力 (N)');    title('滞回环与椭圆拟合');legend('原始数据', '椭圆拟合', 'Location', 'best');    grid on;    subplot(1,3,2);plot(displacement, force, 'b-');    xlabel('位移 (m)'); ylabel('力 (N)');    title('力-位移时程曲线');    grid on;    subplot(1,3,3);    bar([1,2,3], [K_prime, K_double_prime, loss_factor]);    set(gca, 'XTickLabel', {'储能刚度K''', '耗能刚度K''''', '损耗因子η'});    ylabel('数值');    title('动态参数');    grid on;    fprintf('椭圆拟合法结果：\n');    fprintf('  储能刚度 K'' = %.2f N/m\n', K_prime);    fprintf('  耗能刚度 K'''' = %.2f N/m\n', K_double_prime);    fprintf('  损耗因子 η = %.4f\n', loss_factor);end

2.2 傅里叶变换法

function[K_prime, K_double_prime, loss_factor, freq] = calculate_dynamic_stiffness_fft(force, displacement, fs)% 方法2：傅里叶变换法（最准确）% fs: 采样频率 (Hz)    N = length(force);% 1. 去除直流分量    F = force - mean(force);    X = displacement - mean(displacement);% 2. 计算FFT    fft_F = fft(F);    fft_X = fft(X);% 3. 频率向量    freq = (0:N-1) * fs / N;% 4. 找到激励频率（最大位移对应的频率）    [~, idx] = max(abs(fft_X(2:floor(N/2))));  % 忽略直流    idx = idx + 1;  % 补偿索引    f0 = freq(idx);% 5. 提取基波分量    F1 = fft_F(idx);    X1 = fft_X(idx);% 6. 计算复数刚度    K_complex = F1 / X1;  % 复数刚度% 7. 分解为实部和虚部    K_prime = real(K_complex);      % 储能刚度    K_double_prime = imag(K_complex); % 耗能刚度% 8. 计算其他参数    K_mag = abs(K_complex);        % 动刚度幅值    phase_angle = angle(K_complex); % 相位角 (rad)    loss_factor = tan(phase_angle); % 损耗因子% 9. 绘图figure('Position', [100, 100, 1200, 500]);% 频谱图    subplot(2,3,1);plot(freq(1:floor(N/2)), abs(fft_X(1:floor(N/2))), 'b-', 'LineWidth', 1.5);hold on;plot(freq(idx), abs(X1), 'ro', 'MarkerSize', 10, 'LineWidth', 2);    xlabel('频率 (Hz)'); ylabel('位移幅值');    title('位移频谱');    grid on;legend('频谱', '基波分量');    subplot(2,3,2);plot(freq(1:floor(N/2)), abs(fft_F(1:floor(N/2))), 'b-', 'LineWidth', 1.5);hold on;plot(freq(idx), abs(F1), 'ro', 'MarkerSize', 10, 'LineWidth', 2);    xlabel('频率 (Hz)'); ylabel('力幅值');    title('力频谱');    grid on;% 滞回环    subplot(2,3,3);plot(X, F, 'b.');    xlabel('位移 (m)'); ylabel('力 (N)');    title('滞回环');    grid on;% 复数刚度示意图    subplot(2,3,4);plot([0, real(K_complex)], [0, imag(K_complex)], 'r-', 'LineWidth', 2);hold on;plot(real(K_complex), imag(K_complex), 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r');    xlabel('实部 (K'')'); ylabel('虚部 (K'''')');    title('复数刚度');    axis equal; grid on;% 参数表格    subplot(2,3,5);    axis off;    text(0.1, 0.9, sprintf('动态参数分析结果'), 'FontSize', 12, 'FontWeight', 'bold');    text(0.1, 0.7, sprintf('激励频率: %.2f Hz', f0), 'FontSize', 10);    text(0.1, 0.6, sprintf('储能刚度 K'': %.2f N/m', K_prime), 'FontSize', 10);    text(0.1, 0.5, sprintf('耗能刚度 K'''': %.2f N/m', K_double_prime), 'FontSize', 10);    text(0.1, 0.4, sprintf('动刚度幅值: %.2f N/m', K_mag), 'FontSize', 10);    text(0.1, 0.3, sprintf('相位角: %.2f°', rad2deg(phase_angle)), 'FontSize', 10);    text(0.1, 0.2, sprintf('损耗因子 η: %.4f', loss_factor), 'FontSize', 10);% 相位关系    subplot(2,3,6);    t = (0:N-1)/fs;plot(t, X/max(abs(X)), 'b-', 'LineWidth', 1.5);hold on;plot(t, F/max(abs(F)), 'r-', 'LineWidth', 1.5);    xlabel('时间 (s)'); ylabel('归一化幅值');    title('位移与力相位关系');legend('位移', '力');    grid on;    fprintf('\n傅里叶变换法结果：\n');    fprintf('  激励频率: %.2f Hz\n', f0);    fprintf('  储能刚度 K'': %.2f N/m\n', K_prime);    fprintf('  耗能刚度 K'''': %.2f N/m\n', K_double_prime);    fprintf('  动刚度幅值 |K*|: %.2f N/m\n', K_mag);    fprintf('  相位角: %.2f°\n', rad2deg(phase_angle));    fprintf('  损耗因子 η: %.4f\n', loss_factor);end

2.3 瞬时刚度法

function[instant_stiffness, time] = calculate_instant_stiffness(force, displacement, fs, filter_cutoff)% 方法3：计算瞬时（切线）刚度% filter_cutoff: 低通滤波截止频率 (Hz)，用于平滑if nargin < 4        filter_cutoff = 50;  % 默认截止频率end% 1. 计算导数（差分法）    dt = 1/fs;    dF = gradient(force, dt);    dX = gradient(displacement, dt);% 2. 瞬时刚度 = dF/dX% 注意：当dX接近0时会出现奇异值    instant_stiffness = dF ./ dX;% 3. 去除奇异值（dX太小的情况）    dX_threshold = max(abs(dX)) * 1e-3;    instant_stiffness(abs(dX) < dX_threshold) = NaN;% 4. 低通滤波（可选，用于平滑）if filter_cutoff > 0% 设计低通滤波器        nyquist = fs/2;        Wn = filter_cutoff/nyquist;        [b, a] = butter(4, Wn, 'low');% 滤波前处理NaN值        valid_idx = ~isnan(instant_stiffness);        instant_stiffness_filt = instant_stiffness;        instant_stiffness_filt(valid_idx) = filtfilt(b, a, instant_stiffness(valid_idx));        instant_stiffness = instant_stiffness_filt;end% 5. 时间向量    time = (0:length(force)-1)/fs;% 6. 绘图figure('Position', [100, 100, 1000, 800]);% 力-位移曲线    subplot(3,1,1);plot(displacement, force, 'b-');    xlabel('位移 (m)'); ylabel('力 (N)');    title('力-位移曲线');    grid on;% 位移和力时程    subplot(3,1,2);    yyaxis left;plot(time, displacement, 'b-', 'LineWidth', 1.5);    ylabel('位移 (m)');    yyaxis right;plot(time, force, 'r-', 'LineWidth', 1.5);    ylabel('力 (N)');    xlabel('时间 (s)');    title('位移和力时程');legend('位移', '力', 'Location', 'best');    grid on;% 瞬时刚度    subplot(3,1,3);plot(time, instant_stiffness, 'k-', 'LineWidth', 1.5);    xlabel('时间 (s)'); ylabel('瞬时刚度 (N/m)');    title('瞬时（切线）刚度');    grid on;% 统计信息    valid_stiffness = instant_stiffness(~isnan(instant_stiffness));    fprintf('\n瞬时刚度统计：\n');    fprintf('  平均值: %.2f N/m\n', mean(valid_stiffness));    fprintf('  标准差: %.2f N/m\n', std(valid_stiffness));    fprintf('  最小值: %.2f N/m\n', min(valid_stiffness));    fprintf('  最大值: %.2f N/m\n', max(valid_stiffness));end

2.4 主程序：调用示例

% 主程序示例：橡胶隔振器动态刚度分析clear; close all; clc;% ========== 1. 生成模拟数据（用于测试）==========fs = 1000;          % 采样频率 (Hz)T = 2;              % 总时长 (s)f0 = 10;            % 激励频率 (Hz)t = 0:1/fs:T-1/fs;% 位移信号（正弦波）X0 = 0.005;         % 位移幅值 (m)displacement = X0 * sin(2*pi*f0*t);% 模拟橡胶力-位移关系（包含刚度和阻尼）K_prime = 1e6;      % 储能刚度 (N/m)eta = 0.1;          % 损耗因子K_double_prime = K_prime * eta;  % 耗能刚度% 力响应（包含相位差）phase_angle = atan(eta);  % 相位角force = K_prime * displacement + K_double_prime/(2*pi*f0) * gradient(displacement, 1/fs);% 添加噪声（模拟真实数据）noise_level = 0.02; % 2%噪声force = force + noise_level * max(force) * randn(size(force));displacement = displacement + noise_level * max(displacement) * randn(size(displacement));fprintf('=== 橡胶隔振器动态刚度分析 ===\n');fprintf('模拟参数：\n');fprintf('  储能刚度 K'': %.2e N/m\n', K_prime);fprintf('  损耗因子 η: %.3f\n', eta);fprintf('  激励频率: %.1f Hz\n', f0);fprintf('  位移幅值: %.4f m\n\n', X0);% ========== 2. 调用分析方法 ==========% 方法1：椭圆拟合法[K_prime_ellipse, K_double_prime_ellipse, eta_ellipse] = ...    calculate_dynamic_stiffness_ellipse(force, displacement);% 方法2：傅里叶变换法（推荐）[K_prime_fft, K_double_prime_fft, eta_fft, f_identified] = ...    calculate_dynamic_stiffness_fft(force, displacement, fs);% 方法3：瞬时刚度法[instant_stiffness, time] = calculate_instant_stiffness(force, displacement, fs, 50);% ========== 3. 结果对比 ==========fprintf('\n=== 结果对比 ===\n');fprintf('参数         椭圆法       傅里叶法     理论值\n');fprintf('K'' (N/m)   %9.2e   %9.2e   %9.2e\n', ...    K_prime_ellipse, K_prime_fft, K_prime);fprintf('η           %9.4f     %9.4f     %9.4f\n', ...    eta_ellipse, eta_fft, eta);fprintf('频率(Hz)       -        %9.2f     %9.2f\n', f_identified, f0);% ========== 4. 建议 ==========fprintf('\n=== 工程建议 ===\n');fprintf('1. 对于正弦激励，推荐使用傅里叶变换法（最准确）\n');fprintf('2. 椭圆法适用于滞回环接近椭圆的情况\n');fprintf('3. 瞬时刚度可用于分析非线性程度，但需谨慎解释\n');fprintf('4. 储能刚度K''用于计算系统固有频率\n');fprintf('5. 损耗因子η用于评估减振效果（η越大，阻尼越大）\n');

3.对您具体数据的建议

1. 先用主程序测试，确保代码理解
2. 导入您的实际数据（确保格式一致）
3. 先尝试傅里叶变换法，查看频谱确认主频
4. 观察滞回环形状，判断非线性程度
5. 如果滞回环不规则，可能需要分段分析或考虑振幅依赖性