理论分享 | 随机振动分析（2）：随机振动响应的理论基础与工程实现

上篇文章《理论分享 | 随机振动分析（1）：非确定性振动的基本原理与数学描述》提到机械系统对不同频率信号具有天然的滤波特性，因此引入了功率谱密度函数，用以描述随机振动能量在频域中的分布：该函数反映了不同频率下的振动均方值的大小。本文将进一步结合相关理论，介绍自相关函数与功率谱密度的定义、随机振动响应的工程计算实现方法，探讨其背后的物理机理，并与谐响应分析、响应谱分析进行对比。

1、自相关函数与功率谱密度

1.1 自相关函数

相关分析是分析随机信号的重要手段。对于给定的随机信号x ( t ) x(t)x(t)，其自相关函数定义为：

其中，τ \tauτ表示时间偏移。上式表示t tt时刻信号x ( t ) x(t)x(t)与经过τ \tauτ时间偏移信号x ( t + τ ) x(t+\tau)x(t+τ)乘积的时间平均，因此，自相关函数描述了信号本身在不同时刻的相关程度，是关于时间偏移量τ \tauτ的函数。

上式中，若时间偏移τ = 0 \tau=0τ=0，则有：

即R x ( 0 ) R_x(0)Rx​(0)是随机信号的均方值。此外，自相关函数的取值还反映了随机信号的时域特征：

• 若∣ R x ( τ ) ∣ |R_x(\tau)|∣Rx​(τ)∣在某些τ \tauτ处取值较大，说明信号中存在与该τ \tauτ值相对应的频率成分，τ \tauτ值近似反映了该频率的周期；
• 若∣ R x ( τ ) ∣ |R_x(\tau)|∣Rx​(τ)∣在τ = 0 \tau=0τ=0之外迅速衰减至零，说明不同时刻信号几乎不相关，信号变化剧烈，该信号频域能量分布较宽；
• 特别地，若对于任意τ ≠ 0 \tau \neq 0τ=0均有R x ( τ ) = 0 R_x(\tau) = 0Rx​(τ)=0，则信号为理想白噪声，其能量在全频域内为常数；

以下，分别给出了窄带随机信号、宽带随机信号、白噪声信号、白噪声与周期信号叠加的时域信号及其自相关函数的示意图：

根据上图可知：对于窄带随机信号，其自相关函数R x ( τ ) R_x(\tau)Rx​(τ)随时间间隔τ \tauτ的增大而相对缓慢衰减，这意味着在一定时间间隔内，信号仍保持一定的相关性，反映出信号变化相对平缓，能量集中在一定的低频范围内；相比之下，宽带随机信号的自相关函数除τ = 0 \tau = 0τ=0点附近外迅速衰减至零，表明不同时刻的信号值几乎互不相关，信号变化剧烈，频域能量分布较宽。对于理想的白噪声信号，其自相关函数仅在τ = 0 \tau = 0τ=0时有值，其余位置均为0；当白噪声叠加周期信号后，其自相关函数亦呈现周期性，且基本不随τ \tauτ增加而衰减。可见，自相关函数刻画了随机信号振动能量的时域分布特征。

1.2 功率谱密度的定义与计算

正如前面所说的，自相关函数刻画了随机信号振动能量的时域分布特征，可以将其近似等价为一个振动信号。按照傅里叶变换的思想，则自相关函数的傅里叶变换刻画了随机信号振动能量的频域分布特征，满足以下关系：

其中，S x ( ω ) S_x(\omega)Sx​(ω)称为随机信号的功率谱密度函数，简称功率谱。令上述傅里叶逆变换等式中τ = 0 \tau=0τ=0，结合自相关函数的定义，可以推导出如下等式：

上述等式右侧表示随机信号在时域中计算的平均功率：$ x^2(t) $ 反映了 $ t $ 时刻的瞬时功率（或振动强度），对时间积分得到总能量，再取极限将其转换为平均功率。等式左侧则对应频域描述：根据积分意义（忽略因子 $ 1/2\pi，其源于圆频率与频率的转换）， ，其源于圆频率与频率的转换），，其源于圆频率与频率的转换），S_x(\omega) d\omega $ 表示随机信号在微小频率区间 $ [\omega, \omega + d\omega] $ 内的功率贡献，因此 $ S_x(\omega) $ 被称为功率谱密度，其单位为（物理量²/Hz）。由此可见，该等式揭示了一个重要的物理内涵：随机信号的平均功率既可从时域计算，也可从频域计算，两者是等价的。

这一关系本质上与帕塞瓦尔定理相通。帕塞瓦尔定理指出：信号在时域的总能量等于其在频域各频率分量能量之和，即：

将上述等式两边同时除以时间长度 $ T $（对于功率有限信号取 $ T \to \infty $ 的极限），即可得到平均功率的表达式。因此，帕塞瓦尔定理对应的是能量，而功率谱密度描述的是功率，两者相差一个时间平均因子 $ 1/T $。

基于上述关系，功率谱密度还可通过另一种方式计算：先对信号做傅里叶变换得到 $ X(\omega) $，再取其幅值平方并除以时间长度，即：

基于上述两种计算方式，随机信号（以前述宽带信号为例）的功率谱密度可直接计算得到，如下图所示：

一般来说，通常不需要自己处理随机振动信号获取PSD数据，可以通过查阅手册即可获取相应的数据。根据前述定义可知，功率谱密度的单位为（物理量²/Hz），其中“/Hz”表示每赫兹带宽，而非横轴频率值本身——横轴的频率f ff是自变量，反映的是信号能量随频率的分布。这一概念可与质量线密度ρ \rhoρ进行类比：线密度的单位为 kg/m，表示不同位置处单位长度上的质量；若要求得从位置x 1 x_1x1​到x 2 x_2x2​一段铁丝的总质量，需对线密度进行积分∫ x 1 x 2 ρ d x \int_{x_1}^{x_2}\rho dx∫x1​x2​​ρdx。类似地，若要计算随机信号在某一频带[ f 1 , f 2 ] [f_1, f_2][f1​,f2​]内的总能量，也需对功率谱密度在该频带内进行积分∫ f 1 f 2 S x ( f ) d f \int_{f_1}^{f_2}S_x(f)df∫f1​f2​​Sx​(f)df，这正是随机振动分析进行数值计算的基本逻辑。

2、随机振动的工程计算实现

2.1 模态叠加思想

正如前文所述，功率谱密度函数描述了随机振动能量在不同频率下的分布特征。对于线性结构而言，随机振动响应的核心计算思想可以理解为：将输入的随机激励分解为一系列具有不同频率和不同振动水平的激励，分别计算各频率下的结构响应，再进行能量叠加。具体而言，设输入功率谱密度函数为S x ( ω ) S_x(\omega)Sx​(ω)，考虑其在微小频率区间[ ω 1 , ω 1 + d ω ] [\omega_1, \omega_1 + d\omega][ω1​,ω1​+dω]内所包含的能量。在此频段内，可将结构视为受一频率为ω 1 \omega_1ω1​、振动强度为S x ( ω 1 ) S_x(\omega_1)Sx​(ω1​)的激励作用，对应的响应幅值由频响函数H ( ω ) H(\omega)H(ω)确定。依据上篇文章中介绍的随机振动理论，对于单输入单输出系统，响应功率谱密度与激励功率谱密度之间存在如下简单关系：

对于一个多自由度系统的动力学行为可通过如下二阶线性微分方程组来描述：

由于功率谱密度是频域分布，因此，自然而然地需要将上述时域方程转换到频域方程中：

频响函数矩阵H ( ω ) = U ( ω ) F ( ω ) H(\omega)=\frac{U(\omega)}{F(\omega)}H(ω)=F(ω)U(ω)​，可以表示为：

假设输入为激励力功率谱密度S F F ( ω ) S_{FF}(\omega)SFF​(ω)，按照响应功率谱的定义可知：

其中，S F F ( ω ) S_{FF}(\omega)SFF​(ω)为激励力的功率谱密度矩阵，S U U ( ω ) S_{UU}(\omega)SUU​(ω)为位移响应的功率谱密度矩阵，上标∗ ^*∗表示共轭转置。需要指出的是，频响函数矩阵H ( ω ) \mathbf{H}(\omega)H(ω)通常为满秩的非对角矩阵，直接通过求逆方法计算效率低下，尤其在矩阵规模较大时更为明显；而且，随机振动分析还涉及多个频率点，上述公式的工程可行性受到严重制约，仅仅理论可行。（类似问题也出现在谐响应分析的直接法中，不过在谐响应分析中可通过求解线性方程组替代直接求逆，以提升计算效率。）

为解决上述计算效率问题，可以通过模态叠加法进行随机振动问题的求解，其基本流程如下：

• 进行模态分析、坐标变换和模态解耦，从而获得各阶模态坐标的频率响应函数：
  
• 建立模态响应谱与模态力谱之间的关系：
  
• 将模态响应谱转换为物理响应谱，并进行工程简化，输出仿真计算结果：
  
  对上述随机响应的计算过程进行简要总结：首先，通过模态解耦建立模态坐标响应与模态力之间的频响函数关系。仿真中输入的是若干自由度的力谱，即一个力谱矩阵。通常假设各节点力相互独立，因此物理坐标下的力谱矩阵为对角阵。根据模态变换思想，将物理坐标系下的力转换至模态坐标系，得到模态力谱。理论上模态力谱不再是对角形式，但实际计算中往往忽略耦合项，从而可由模态力谱直接求解模态响应功率谱。得到模态响应功率谱后，需将其恢复至物理响应。工程中通常关注响应的自谱，因此可将矩阵运算简化为标量运算，以提高效率。需要指出的是，以上过程是在单一激励频率下完成的。随机振动分析涉及一定频率范围内的激励谱，需采用微元法逐个频点计算响应，最后积分得到总均方响应。如上篇文章所述，均方根值即为关注的1σ值。

2.2 工程案例实现及仿真结果对比

作为理论分享的验证，这里针对随机振动分析进行了编程实现。采用悬臂梁案例，梁长100 mm，截面尺寸8 mm×4 mm，单元类型为SOLID185，材料为铝合金（弹性模量68 GPa，泊松比0.3）。梁一端固支，另一端中点施加沿4 mm厚度方向的力谱，频率范围为1–50 Hz，功率谱密度为常值0.1。输出激励方向下1 σ 1\sigma1σ变形值，ANSYS仿真结果与自编程序结果对比如下：

简要说明：上图给出了沿激励方向的1 σ 1\sigma1σ变形云图（实际出现的位移小于该仿真结果的概率为68.269%，亦可查看3 σ 3\sigma3σ值），两者分布基本一致，幅值误差在2%以内（误差来源暂未确定）。

关键部分Matlab代码如下：

%% 模态求解 nEig = 6; % 求解6阶 matFun = @(x) dK\(M(free, free)*x); warning off;[eivecs(:, :), D(:, :)] = eigs(matFun, length(free), nEig+4, 'la', ... 'Tolerance', 1e-12, 'MaxIterations', 200); % 10阶特征值计算 freqModal = zeros(ndof, nEig); gama_num = zeros(nEig, 1); [mu, ii] = sort(diag(D(:, :)), 'descend'); % 圆频率平方倒数 eivSort = eivecs(:, ii(1:nEig)); gama_num(:) = mu(1:nEig); freqModal(free, :) = eivSort./sqrt(diag(eivSort'*M(free, free)*eivSort)'); % 振型向量关于质量矩阵归一化 lambda = 1./mu(1:nEig); frequency = sqrt(lambda)/2/pi; % 固有频率 omega = 2*pi*frequency; % 圆频率 %% 随机振动响应求解 freq_range = [1, 50]; % 频率分析范围 (Hz) n_freq_points = 100; % 频率点数 freq_points = linspace(freq_range(1), freq_range(2), n_freq_points); zeta = 0.02; % 模态阻尼比 PSD_input = 0.1; % 输入力谱常数PSD（可根据实际情况修改） PSD_response = zeros(ndof, n_freq_points); % 初始化所有自由度的响应PSD RMS_response = zeros(ndof, 1); % 初始化均方根响应 SFF = zeros(ndof); SFF(dof_excitation,dof_excitation) = 0.1; Sff = diag(diag(freqModal'*SFF*freqModal)); for i = 1:n_freq_points % 矩阵计算法，逐个频点计算 w = 2*pi*freq_points(i); H_mod = diag(1./(omega.^2 - w^2 + 2i*zeta.*omega.*w)); % 模态频响对角矩阵（位移/力） S_qq = H_mod * Sff * H_mod'; % 模态响应PSD矩阵 S_yy = freqModal * S_qq * freqModal'; % 物理响应PSD矩阵 PSD_response(:,i) = real(diag(S_yy)); % 取自谱，应该为实数，数值计算原因导致存在虚部，但很小 end for i = 1:n_freq_points % 标量计算法，将矩阵乘法离散化为标量相乘，效率更快，同样逐个频点计算 w = 2*pi*freq_points(i); for j = 1:nEig H_j = 1 / (omega(j)^2 - w^2 + 2i*zeta*omega(j)*w); % 频率响应函数 PSD_response(:, i) = PSD_response(:, i) + freqModal(:,j).^2*abs(H_j)^2*Sff(j,j); % 仅计算自谱 end end for i = 1:ndof % 计算各自由度的RMS响应 RMS_response(i) = sqrt(trapz(freq_points, PSD_response(i, :))); end

此外，ANSYS Workbench平台无法直接施加力谱，通常仅支持加速度谱，但可通过命令流实现力谱加载，具体方式如下：

PSDUNIT,1,FORCE ! 定义PSD类型为力 PSDFRQ,1,0,1,50 ! 定义PSD谱的频率和幅值 PSDVAL,1,0.1,0.1 NSEL,S,LOC,Y,100 ! 选择该面上的中心节点作为加载 NSEL,R,LOC,X,4 NSEL,R,LOC,Z,2 F,all,FZ,1 !在Z方向施加1倍的力谱 ALLSEL PFACT,1,NODE ! 计算模态参与因子

本文在进行随机振动响应计算时，采用了相对简单的恒定幅值力谱。然而在实际工程中，更为常用的是加速度功率谱密度。在处理基础加速度激励下的随机振动问题时，需要区分约束节点位移与自由节点位移，并将基础运动的加速度谱转换为等效载荷进行计算。该转换过程涉及之前文章中介绍过的模态参与因子，其实前文构建模态力谱矩阵时也涉及模态参与，只是表达形式有所不同，分别为Φ i M λ \Phi_i \mathbf{M} \lambdaΦi​Mλ和Φ i F \Phi_i FΦi​F。相关推导相对复杂一些，但求解思路基本一致，这里不再展开。

3、谐响应分析、响应谱分析与随机振动分析的区别与联系

针对工程中常见的谐响应分析、响应谱分析与随机振动分析，以下从作用、求解原理的区别与联系进行简要梳理：

• 谐响应分析用于求解结构在确定性简谐载荷作用下的稳态响应。基于模态叠加原理，对于每一个外界激励频率，各阶模态均以不同的幅值和相位参与响应，叠加后得到该频率下的总响应。不同激励频率下的响应相互独立、无需叠加，描述的是响应随激励频率变化的特性；
• 响应谱分析用于求解结构在已知激励谱（如地震谱）作用下的最大响应。其输入谱记录了不同固有频率的单自由度系统在该激励下可能产生的峰值响应。实际计算时，直接按结构的各阶固有频率从谱曲线中读取对应的响应最大值，再通过组合（如SRSS、CQC）得到总体响应。这里输入谱的横轴“频率”指的是结构的固有频率，而非外界激励的频率；
• 随机振动分析用于求解结构在随机激励下的概率响应。以功率谱密度描述激励在频域上的统计分布。实际计算时，针对每一小段外界激励频率（按照微元法思想，将功率谱密度分为各个频率对应的激励谱值），求解各阶模态的响应功率谱，叠加得到该激励频率下的物理响应功率谱；再通过对频率轴积分，得到响应的统计特征（如位移均方值）。其中功率谱密度的横轴“频率”与谐响应分析中的频率物理意义一致，指外界激励的频率；

总的来说，三者均蕴含模态叠加的思想，但谐响应分析与随机振动分析涉及在外界不同激励频率下对各阶模态响应进行叠加的过程；而响应谱分析则是基于输入谱（对应单自由度系统在复合激励下的峰值响应）直接读取各阶模态的最大值并进行组合，不存在对每个激励频率逐点求解的步骤。

4、总结

• 介绍了自相关函数和功率谱密度函数的定义及其物理意义：自相关函数描述随机振动的时域相关特性，功率谱密度反映其频域能量分布，两者互为傅里叶变换对；
• 阐述了随机振动分析中的模态叠加思想及其物理本质：其核心计算思路是基于功率谱密度，将频率轴离散化为微元，在每个微元上利用模态叠加法求解结构的响应功率谱，再通过积分获得响应的统计特征；同时，编写了随机振动分析程序，并通过工程案例与商业软件结果进行了对比验证；

最后

截至本文，共对工程中常用的几种动力学分析方法：模态分析、谐响应分析、响应谱分析及随机振动分析的作用、求解原理以及相互之间的区别与联系进行了介绍（瞬态响应因物理意义较为直观明确，故未作展开），并分享了个人的一些理解，希望能为大家提供一些借鉴与参考。在撰写过程中参考了诸多相关资料，由于资料来源较多，这里就不列举了，在此一并进行感谢！由于个人水平有限，文章内容难免存在不足之处，如有错误或疏漏，欢迎大家指正。

有需要上述完整代码的读者请后台私信作者！

-完-
公众号已同步更新：仿真优化工匠
欢迎联系！