文章目录1、自相关函数与功率谱密度1.1 自相关函数1.2 功率谱密度的定义与计算2、随机振动的工程计算实现2.1 模态叠加思想2.2 工程案例实现及仿真结果对比3、谐响应分析、响应谱分析与随机振动分析的区别与联系4、总结最后上篇文章《理论分享 | 随机振动分析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 \tau0τ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 \tau0τ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 \tau0τ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), nEig4, 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、总结介绍了自相关函数和功率谱密度函数的定义及其物理意义自相关函数描述随机振动的时域相关特性功率谱密度反映其频域能量分布两者互为傅里叶变换对阐述了随机振动分析中的模态叠加思想及其物理本质其核心计算思路是基于功率谱密度将频率轴离散化为微元在每个微元上利用模态叠加法求解结构的响应功率谱再通过积分获得响应的统计特征同时编写了随机振动分析程序并通过工程案例与商业软件结果进行了对比验证最后截至本文共对工程中常用的几种动力学分析方法模态分析、谐响应分析、响应谱分析及随机振动分析的作用、求解原理以及相互之间的区别与联系进行了介绍瞬态响应因物理意义较为直观明确故未作展开并分享了个人的一些理解希望能为大家提供一些借鉴与参考。在撰写过程中参考了诸多相关资料由于资料来源较多这里就不列举了在此一并进行感谢由于个人水平有限文章内容难免存在不足之处如有错误或疏漏欢迎大家指正。有需要上述完整代码的读者请后台私信作者-完-公众号已同步更新仿真优化工匠欢迎联系