CN104850713A - 机械结构随机振动动态应力高精度计算方法 - Google Patents

Abstract

本发明公开了一种机械结构随机振动动态应力高精度计算方法，将随机激励转化为确定性激励，在完全考虑模态耦合效应的同时，提高计算效率；采用模态应力系数计算，避免从位移求导计算应力时精度的损失；进行静力修正，考虑模态截取时忽略的高阶模态的影响，提高随机动态应力的计算精度。本发明在计算时包含了全部参振模态之间的互相关项，得到的是精确解。对于大型复杂工程结构，本发明易于操作实践，其计算效率与传统算法相比有较大提高，并且由于应力模态矩阵和高阶模态静力修正项的引入，提高了计算精度。

Description

机械结构随机振动动态应力高精度计算方法

技术领域

本发明涉及结构力学分析技术领域，具体为一种机械结构随机振动动态应力高精度计算方法。

背景技术

结构受到随机振动载荷作用时，尤其当载荷频率覆盖了结构多阶固有频率时，会导致结构处于较高的动应力水平，以至于发生疲劳破坏。因此，动态应力是结构动强度校核和振动疲劳寿命预估所需的原始数据。只有有了准确的动态应力响应，才能更准确的预估结构疲劳寿命，这直接关系到结构系统的安全可靠性。

结构受迫振动方程描述系统受激励下的动力学问题，是结构动力学行为的时间域描述，在数学上是关于时间变量的二阶微分方程。求解该方程，可以得到位移、速度、加速度等响应。北京航空航天大学的郭兴旺和北京科技大学的邹家祥在1996年学术期刊《振动与冲击》第15卷第2期《对机械振动系统的六种动态响应分析方法的评述》一文中，较为全面地概述了系统动态响应分析的常用方法，指明了其特点和适用范围。其中工程应用较多的求解方法是直接积分法和模态叠加法等。

直接积分法就是直接对运动方程进行积分，将连续时间域离散成n个间距为Δt的时间点，若已经求得了t时刻之前所有时刻的解，再根据设定的前后时刻位移、速度、加速度关系，代入运动方程，继续求t+Δt时刻的解。根据设定的前后时刻位移、速度、加速度关系不同，又分为中心差分法、Wilson-θ法、Newmark法等。直接积分法的计算精度与时间步长有关，随着时间步长的增加，误差逐渐增大，会出现周期延长或振幅衰减的情况。直接积分法计算量巨大，效率极低，适用于求解相对较短时间内的响应。

模态叠加法以系统的模态向量为基底来描述物理量，是一种建立在坐标变换基础上的求解方法。对于多自由度系统，先将系统运动方程通过线性变换解耦为n个单自由度的方程，确定模态坐标响应，再叠加得到物理坐标响应。这种方法只适用于线性系统，在实际计算中，往往只取部分低阶振型进行叠加，忽略了高阶模态的作用。

随机振动中，在求得模态坐标响应的统计量和振型矩阵之后，计算系统响应的统计量时，需要考虑振型组合问题。常见的阵型组合方法有：SRSS法、ABS法、CQC法。SRSS法认为各阶模态响应之间是相互独立的，ABS法则假设所有模态响应之间精确相关，CQC法能得到较好的计算结果，它更准确的考虑了各阶模态的耦合效应，但给出的模态相关系数计算复杂，计算量大，效率低。大连理工大学林家浩在2004年科学出版社出版的其与张亚辉编著的《随机振动的虚拟激励法》一书中提出的虚拟激励算法，将平稳随机振动分析转化为简谐振动分析，计算步骤简化，效率有数量级的提高，但对于振动疲劳寿命预估中关心的随机动应力的计算阐述不多。

目前，在研究飞机结构振动疲劳中的应力响应时，需要考虑以下两个影响因素：随机特性及其引起的多模态耦合：飞机结构处于复杂振动环境下，载荷激励和结构响应通常无法确切表示，只能用随机变量来表示其统计意义；当飞机处于随机振动载荷作用下，由于振动载荷形式具有宽频的特点，会导致结构的多阶模态被同时激起，结构的动态响应会更加复杂。2010年9月，在沈阳飞机设计所承办的中国航空结构动力学专业组第十七届学术交流会上，讨论发现现有结构随机动应力工程计算方法存在以下几方面问题：1，完全考虑模态耦合时，相关系数计算复杂，计算量大；2，当随机振动载荷的激励频带较宽时,通过位移结果求导计算应力时，损失了一阶精度，很难反映应力的剧烈变化,造成结构动应力的计算误差偏大，对于振动疲劳寿命分析往往不能满足要求；3，实际的工程结构自由度数很大，计算中不可能把所有模态都予以考虑，而仅能采用有限的低阶模态，忽略了高阶模态的影响，会产生较大的误差。

发明内容

为解决现有机械结构随机振动动态应力计算方法计算效率低，误差大等问题，本发明提出了一种机械结构随机振动动态应力高精度计算方法，在将随机激励转化为确定性激励的基础上，采用模态应力系数计算，避免从位移求导计算应力时精度的损失；并进行静力修正，考虑模态截取时忽略的高阶模态的影响。

本发明的技术方案为：

所述一种机械结构随机振动动态应力高精度计算方法，其特征在于：包括以下步骤：

步骤1：建立机械结构的有限元模型；对有限元模型施加虚拟激励，进行整体静力分析，得到机械结构整体等效静力响应所述虚拟激励其中S_f(ω)为机械结构受到的随机激励的功率谱密度，ω为随机激励的载荷频率，j为虚数符号，t为时间；

步骤2：对机械结构的有限元模型进行模态分析，得到机械结构的前N_d阶固有频率ω_i，以及对应的模态位移向量{φ}_i，i＝1,2,…,N_d；其中第N_d阶固有频率大于随机激励功率谱频率最大值；

步骤3：根据公式[Φ^σ]＝EB[Φ]得到模态应力矩阵[Φ^σ]，其中E为弹性矩阵，B为应变位移矩阵，[Φ]为模态位移矩阵，模态位移矩阵由模态位移向量{φ}_i组成；

步骤4：采用模态位移矩阵[Φ]对机械结构有限元模型受虚拟激励的运动方程

$\left[M\right]\left\{\stackrel{\·\·}{\tilde{y}}\right\}+\left[C\right]\left\{\stackrel{\·}{\tilde{y}}\right\}+\left[K\right]\left\{\tilde{y}\right]=\left\{\tilde{f}\left(t\right)\right\}$

进行模态转换，得到转换后的运动方程

$\left[\stackrel{\‾}{M}\right]\left\{\stackrel{\·\·}{\tilde{q}}\right\}+\left[\stackrel{\‾}{C}\right]\left\{\stackrel{\·}{\tilde{q}}\right\}+\left[\stackrel{\‾}{K}\left]\right\{\tilde{q}\right\}=\left\{\tilde{P}\left(t\right)\right\}$

其中，[M]为结构质量矩阵，[C]为结构阻尼矩阵，[K]为结构刚度矩阵，为虚拟位移向量， $\left[\stackrel{\‾}{M}\right]={\left[\mathrm{\Φ}\right]}^T\left[M\right]\left[\mathrm{\Φ}\right],\left[\stackrel{\‾}{C}\right]={\left[\mathrm{\Φ}\right]}^T\left[C\right]\left[\mathrm{\Φ}\right],\left[\stackrel{\‾}{K}\right]={\left[\mathrm{\Φ}\right]}^T\left[K\right]\left[\mathrm{\Φ}\right],$ 为虚拟模态坐标向量， $\left\{\tilde{y}\right\}=\left[\mathrm{\Φ}\right]\left\{\tilde{q}\right\}={\mathrm{\Σ}}_{i=1}^{N_d}{\tilde{q}}_i{\left\{\mathrm{\φ}\right\}}_i,$ 为虚拟模态坐标，为虚拟广义激励向量， $\left\{\tilde{P}\left(t\right)\right\}={\left[\mathrm{\Φ}\right]}^T\left\{\tilde{f}\left(t\right)\right\};$

将转换后的运动方程分解为N_d个互相独立的单自由度简谐振动方程

${\stackrel{\·\·}{\tilde{q}}}_i+2{\mathrm{\ξ}}_i{\mathrm{\ω}}_i{\stackrel{\·}{\tilde{q}}}_i+{{\mathrm{\ω}}_i}^2{\tilde{q}}_i=\widetilde{P_i}$

其中ξ_i为第i阶模态阻尼比，为第i阶广义虚拟激励量，求解单自由度简谐振动方程得到虚拟模态坐标以及模态等效静力解

${\tilde{q}}_{\mathrm{si}}=\frac{\widetilde{P_i}}{{{\mathrm{\ω}}_i}^2}$

其中H_i为第i阶频响函数，H_i＝(ω_i ²-ω²+2jξ_iω_iω)^-1；

步骤5：根据公式

${\left\{\widetilde{\mathrm{\σ}}\right\}}_d=\underset{i=1}{\overset{N_d}{\mathrm{\Σ}}}{\left\{{\mathrm{\φ}}^{\mathrm{\σ}}\right\}}_i[{\tilde{q}}_i-{\tilde{q}}_{\mathrm{si}}]$

得到虚拟动态应力响应其中{φ^σ}_i为模态应力向量；将虚拟动态应力响应与步骤1得到的机械结构整体等效静力响应相加得到虚拟应力响应量继而得到机械结构随机动态应力功率谱密度

有益效果

本发明将随机激励转化为确定性激励，完全考虑模态耦合效应的同时，提高计算效率；采用模态应力系数计算，避免从位移求导计算应力时精度的损失；进行静力修正，考虑模态截取时忽略的高阶模态的影响，提高随机动态应力的计算精度。本发明在计算时包含了全部参振模态之间的互相关项，得到的是精确解。对于大型复杂工程结构，本发明易于操作实践，其计算效率与传统算法相比有较大提高，并且由于应力模态矩阵和高阶模态静力修正项的引入，提高了计算精度。

附图说明

图1：本发明的流程图。

具体实施方式

下面根据具体实施例描述本发明：

首先给出本方法的推导过程：

多自由度离散结构受随机激励时的运动方程

$\left[M\right]\left\{\stackrel{\·\·}{y}\right\}+\left[C\right]\left\{\stackrel{\·}{y}\right\}+\left[K\right]\left\{y\right\}=\left\{f\left(t\right)\right\}---\left(1\right)$

式中：

[M]——结构质量矩阵；

[C]——结构阻尼矩阵；

[K]——结构刚度矩阵；

{y}——离散结构加速度向量、速度向量、位移向量。

随机激励{f(t)}功率谱密度S_f(ω)已知。利用S_f(ω)构造虚拟激励

$\tilde{f}\left(t\right)=\sqrt{S_f\left(\mathrm{\ω}\right)}{e}^{\mathrm{j\ωt}}---\left(2\right)$

将随机载荷转化为确定性简谐载荷。

得到关于虚拟量的运动方程

$\left[M\right]\left\{\stackrel{\·\·}{\tilde{y}}\right\}+\left[C\right]\left\{\stackrel{\·}{\tilde{y}}\right\}+\left[K\right]\left\{\tilde{y}\right\}=\left\{\tilde{f}\left(t\right)\right\}---\left(3\right)$

其模态位移矩阵为[Φ]，可得其模态应力矩阵为

[Φ^σ]＝EB[Φ]      (4)

式中：

[Φ^σ]——模态应力矩阵；

E——弹性矩阵；

B——应变位移矩阵；

[Φ]——模态位移矩阵。

代入上式，并用[Φ]^T左乘各项，得

$\left[\stackrel{\‾}{M}\right]\left\{\stackrel{\·\·}{\tilde{q}}\right\}+\left[\stackrel{\‾}{C}\right]\left\{\stackrel{\·}{\tilde{q}}\right\}+\left[\stackrel{\‾}{K}\right]\left\{\tilde{q}\right\}={\left[\mathrm{\Φ}\right]}^T\left\{\tilde{f}\left(t\right)\right\}=\left\{\tilde{P}\left(t\right)\right\}---\left(5\right)$

其中， $\left[\stackrel{\‾}{M}\right]={\left[\mathrm{\Φ}\right]}^T\left[M\right]\left[\mathrm{\Φ}\right],\left[\stackrel{\‾}{C}\right]={\left[\mathrm{\Φ}\right]}^T\left[C\right]\left[\mathrm{\Φ}\right],\left[\stackrel{\‾}{K}\right]={\left[\mathrm{\Φ}\right]}^T\left[K\right]\left[\mathrm{\Φ}\right].$

将运动方程可分解为N_d个互相独立的单自由度简谐振动方程

${\stackrel{\·\·}{\tilde{q}}}_i+2{\mathrm{\ξ}}_i{\mathrm{\ω}}_i{\stackrel{\·}{\tilde{q}}}_i+{{\mathrm{\ω}}_i}^2{\tilde{q}}_i=\widetilde{P_i}---\left(6\right)$

式中：

ξ_i——第i阶模态阻尼比；

ω_i——第i阶固有频率；

——广义虚拟激励量，

得到虚拟模态坐标稳态解

${\tilde{q}}_i=H_i\widetilde{P_i}---\left(7\right)$

式中：

H_i——第i阶频响函数，H_i＝(ω_i ²-ω²+2jξ_iω_iω)^-1。

虚拟应力响应其中{φ^σ}_i为模态应力向量，表征在时对应力响应的贡献。对于线性系统，结构响应可以表示成为模态坐标的叠加形式。参考模态位移的概念，建立模态应力向量，其物理意义是其相应的模态坐标为单位量，其他模态坐标为零时，结构的应力分布向量。借鉴位移可以表示成模态向量与模态位移乘积的叠加形式，将应力表示成模态应力向量与模态位移乘积的叠加形式。这样求得模态坐标，进行叠加后就可直接获取应力，避免了通过位移求导计算应力时损失的精度，这一点在荷载的激励频率范围大，位移响应变化剧烈时尤其重要。

采用模态叠加法进行求解时，考虑所有高阶模态的影响，但高阶模态相应的模态坐标反应的求解采用简化的静力分析法。注意到，高阶模态坐标对应的自振频率较大，即结构中惯性力和阻尼力的影响相对较小，可以忽略。这样，对阶数足够高的模态，其模态坐标可以通过静力法求解。

考虑虚拟应力响应

$\left\{\widetilde{\mathrm{\σ}}\right\}={\mathrm{\Σ}}_{i=1}^{N_d}{\tilde{q}}_i{\left\{{\mathrm{\φ}}^{\mathrm{\σ}}\right\}}_i+{\mathrm{\Σ}}_{i=N_d+1}^N{\tilde{q}}_i{\left\{{\mathrm{\φ}}^{\mathrm{\σ}}\right\}}_i---\left(8\right)$

其中前N_d阶相应于模态截断后，实际采用的低阶项，动力影响明显，而N_d+1项至N项为高阶模态项，其模态反应通过静力法计算：

${\tilde{q}}_i=\frac{{\tilde{P}}_i}{{{\mathrm{\ω}}_i}^2},i=N_d+1,...,N$

代入上式，得到采用静力修正后的应力计算的公式为

$\left\{\widetilde{\mathrm{\σ}}\right\}={\mathrm{\Σ}}_{i=1}^{N_d}{\tilde{q}}_i{\left\{{\mathrm{\φ}}^{\mathrm{\σ}}\right\}}_i+{\mathrm{\Σ}}_{i=N_d+1}^N\frac{{\tilde{P}}_i}{{{\mathrm{\ω}}_i}^2}{\left\{{\mathrm{\φ}}^{\mathrm{\σ}}\right\}}_i---\left(10\right)$

后一项计算采用简化计算方法，推导如下：

对等效静力反应，直接采用整体平衡方程，满足

${\left\{\widetilde{\mathrm{\σ}}\right\}}_s=\left\{\tilde{f}\left(t\right)\right\}---\left(11\right)$

再采用模态叠加

${\left\{\widetilde{\mathrm{\σ}}\right\}}_s={\mathrm{\Σ}}_{i=1}^N{\tilde{q}}_{\mathrm{si}}{\left\{{\mathrm{\φ}}^{\mathrm{\σ}}\right\}}_i---\left(12\right)$

可得模态反应

$K_i{\tilde{q}}_{\mathrm{si}}=\widetilde{P_i}$

则等效静力解为

${\left\{\widetilde{\mathrm{\σ}}\right\}}_s={\mathrm{\Σ}}_{i=1}^N\frac{{\tilde{P}}_i}{K_i}{\left\{{\mathrm{\φ}}^{\mathrm{\σ}}\right\}}_i---\left(13\right)$

因此根据上式可把N_d+1项至N项的求和用前N_d项表示如下

${\mathrm{\Σ}}_{i=N_d+1}^N\frac{{\tilde{P}}_i}{K_i}{\left\{{\mathrm{\φ}}^{\mathrm{\σ}}\right\}}_i={\left\{\widetilde{\mathrm{\σ}}\right\}}_s-{\mathrm{\Σ}}_{i=1}^{N_d}\frac{{\tilde{P}}_i}{K_i}{\left\{{\mathrm{\φ}}^{\mathrm{\σ}}\right\}}_i,K_i={{\mathrm{\ω}}_i}^2---\left(14\right)$

这样就得到可以计算的静力修正方法公式

$\begin{array}{c}\left\{\widetilde{\mathrm{\σ}}\right\}=\underset{i=1}{\overset{N_d}{\mathrm{\Σ}}}{\tilde{q}}_i{\left\{{\mathrm{\φ}}^{\mathrm{\σ}}\right\}}_i+{\left\{\widetilde{\mathrm{\σ}}\right\}}_s-\underset{i=1}{\overset{N_d}{\mathrm{\Σ}}}\frac{\widetilde{P_i}}{K_i}{\left\{{\mathrm{\φ}}^{\mathrm{\σ}}\right\}}_i\\ ={\left\{\widetilde{\mathrm{\σ}}\right\}}_s+{\mathrm{\Σ}}_{i=1}^{N_d}{\left\{{\mathrm{\φ}}^{\mathrm{\σ}}\right\}}_i[{\tilde{q}}_i-\frac{{\tilde{P}}_i}{{{\mathrm{\ω}}_i}^2}]\end{array}---\left(15\right)$

应力响应功率谱矩阵

$\left[S_{\mathrm{\σ}}\left(\mathrm{\ω}\right)\right]={\left\{\widetilde{\mathrm{\σ}}\right\}}^{*}{\left[\widetilde{\mathrm{\σ}}\right]}^T.---\left(16\right)$

根据上述推导过程，本实施例的机械结构随机振动动态应力高精度计算方法包括以下步骤：

步骤1：建立机械结构的有限元模型；对有限元模型施加虚拟激励，进行整体静力分析，得到机械结构整体等效静力响应所述虚拟激励其中S_f(ω)为机械结构受到的随机激励的功率谱密度，ω为随机激励的载荷频率，j为虚数符号，t为时间。

步骤2：对机械结构的有限元模型进行模态分析，得到机械结构的前N_d阶固有频率ω_i，以及对应的模态位移向量{φ}_i，i＝1,2,…,N_d；其中第N_d阶固有频率大于随机激励功率谱频率最大值。

步骤3：根据公式[Φ^σ]＝EB[Φ]得到模态应力矩阵[Φ^σ]，其中E为弹性矩阵，B为应变位移矩阵，[Φ]为模态位移矩阵，模态位移矩阵由模态位移向量{φ}_i组成。

步骤4：采用模态位移矩阵[Φ]对机械结构有限元模型受虚拟激励的运动方程

$\left[M\right]\left\{\stackrel{\·\·}{\tilde{y}}\right\}+\left[C\right]\left\{\stackrel{\·}{\tilde{y}}\right\}+\left[K\right]\left\{\tilde{y}\right\}=\left\{\tilde{f}\left(t\right)\right\}$

进行模态转换，得到转换后的运动方程

$\left[\stackrel{\‾}{M}\right]\left\{\stackrel{\·\·}{\tilde{q}}\right\}+\left[\stackrel{\‾}{C}\right]\left\{\stackrel{\·}{\tilde{q}}\right\}+\left[\stackrel{\‾}{K}\right]\left\{\tilde{q}\right\}\left\{\tilde{P}\left(t\right)\right\}$

其中，[M]为结构质量矩阵，[C]为结构阻尼矩阵，[K]为结构刚度矩阵，为虚拟位移向量， $\left[\stackrel{\‾}{M}\right]={\left[\mathrm{\Φ}\right]}^T\left[M\right]\left[\mathrm{\Φ}\right],\left[\stackrel{\‾}{C}\right]={\left[\mathrm{\Φ}\right]}^T\left[C\right]\left[\mathrm{\Φ}\right],\left[\stackrel{\‾}{K}\right]={\left[\mathrm{\Φ}\right]}^T\left[K\right]\left[\mathrm{\Φ}\right],$ 为虚拟模态坐标向量， $\left\{\tilde{y}\right\}=\left[\mathrm{\Φ}\right]\left\{\tilde{q}\right\}={\mathrm{\Σ}}_{i=1}^{N_d}{\tilde{q}}_i{\left\{\mathrm{\φ}\right\}}_i,$ 为虚拟模态坐标，为虚拟广义激励向量， $\left\{\tilde{P}\left(t\right)\right\}={\left[\mathrm{\Φ}\right]}^T\left\{\tilde{f}\left(t\right)\right\};$

将转换后的运动方程分解为N_d个互相独立的单自由度简谐振动方程

${\stackrel{\·\·}{\tilde{q}}}_i+2{\mathrm{\ξ}}_i{\mathrm{\ω}}_i{\stackrel{\·}{\tilde{q}}}_i+{{\mathrm{\ω}}_i}^2{\tilde{q}}_i=\widetilde{P_i}$

其中ξ_i为第i阶模态阻尼比，为第i阶广义虚拟激励量，求解单自由度简谐振动方程得到虚拟模态坐标以及模态等效静力解

${\tilde{q}}_{\mathrm{si}}=\frac{\widetilde{P_i}}{{{\mathrm{\ω}}_i}^2}$

其中H_i为第i阶频响函数，H_i＝(ω_i ²-ω²+2jξ_iω_iω)^-1；

步骤5：根据公式

${\left\{\widetilde{\mathrm{\σ}}\right\}}_d=\underset{i=1}{\overset{N_d}{\mathrm{\Σ}}}{\left\{{\mathrm{\φ}}^{\mathrm{\σ}}\right\}}_i[{\tilde{q}}_i-{\tilde{q}}_{\mathrm{si}}]$

得到虚拟动态应力响应其中{φ^σ}_i为模态应力向量；将虚拟动态应力响应与步骤1得到的机械结构整体等效静力响应相加得到虚拟应力响应量继而得到机械结构随机动态应力功率谱密度

本实施例中，机械结构采用铝合金薄板，施加带宽为260Hz—300Hz的窄带随机激励，加速度功率谱密度(PSD)为1g2/Hz；和带宽为100Hz—1300Hz，加速度功率谱密度(PSD)为0.2g2/Hz的宽带随机激励，进行随机振动试验，测试振动应力响应，并分别用常规有限元方法和本发明提出的方法计算随机动应力。所得结果见下表所示：

通过上表可以看出，采用本文计算方法得到的结果与试验值更为接近，误差更小，即本发明提出的方法计算精度高。

Claims (1)

1.一种机械结构随机振动动态应力高精度计算方法，其特征在于：包括以下步骤：

步骤1：建立机械结构的有限元模型；对有限元模型施加虚拟激励，进行整体静力分析，得到机械结构整体等效静力响应所述虚拟激励其中S_f(ω)为机械结构受到的随机激励的功率谱密度，ω为随机激励的载荷频率，j为虚数符号，t为时间；

步骤2：对机械结构的有限元模型进行模态分析，得到机械结构的前N_d阶固有频率ω_i，以及对应的模态位移向量{φ}_i，i＝1,2,…,N_d；其中第N_d阶固有频率大于随机激励功率谱频率最大值；

步骤3：根据公式[Φ^σ]＝EB[Φ]得到模态应力矩阵[Φ^σ]，其中E为弹性矩阵，B为应变位移矩阵，[Φ]为模态位移矩阵，模态位移矩阵由模态位移向量{φ}_i组成；

步骤4：采用模态位移矩阵[Φ]对机械结构有限元模型受虚拟激励的运动方程

$\left[M\right]\left\{\stackrel{\·\·}{\tilde{y}}\right\}+\left[C\right]\left\{\stackrel{\·}{\tilde{y}}\right\}+\left[K\right]\left\{\tilde{y}\right\}=\left\{\tilde{f}\left(t\right)\right\}$

进行模态转换，得到转换后的运动方程

$\left[\stackrel{\‾}{M}\right]\left\{\stackrel{\·\·}{\tilde{q}}\right\}+\left[\stackrel{\‾}{C}\right]\left\{\stackrel{\·}{\tilde{q}}\right\}+\left[\stackrel{\‾}{K}\right]\left\{\tilde{q}\right\}=\left\{\tilde{P}\left(t\right)\right\}$

其中，[M]为结构质量矩阵，[C]为结构阻尼矩阵，[K]为结构刚度矩阵，为虚拟位移向量， $\left[\stackrel{\‾}{M}\right]={\left[\mathrm{\Φ}\right]}^T\left[M\right]\left[\mathrm{\Φ}\right],\left[\stackrel{\‾}{C}\right]={\left[\mathrm{\Φ}\right]}^T\left[C\right]\left[\mathrm{\Φ}\right],\left[\stackrel{\‾}{K}\right]={\left[\mathrm{\Φ}\right]}^T\left[K\right]\left[\mathrm{\Φ}\right],$ 为虚拟模态坐标向量， 为虚拟模态坐标，为虚拟广义激励向量， $\left\{\tilde{P}\left(t\right)\right\}={\left[\mathrm{\Φ}\right]}^T\left\{\tilde{f}\left(t\right)\right\};$

将转换后的运动方程分解为N_d个互相独立的单自由度简谐振动方程

${\stackrel{\·\·}{\tilde{q}}}_i+2{\mathrm{\ξ}}_i{\mathrm{\ω}}_i{\stackrel{\·}{\tilde{q}}}_i+{{\mathrm{\ω}}_i}^2{\tilde{q}}_i={\tilde{P}}_i$

其中ξ_i为第i阶模态阻尼比，为第i阶广义虚拟激励量，求解单自由度简谐振动方程得到虚拟模态坐标以及模态等效静力解

${\tilde{q}}_{\mathrm{si}}=\frac{{\tilde{P}}_i}{{{\mathrm{\ω}}_i}^2}$

其中H_i为第i阶频响函数，H_i＝(ω_i ²-ω²+2jξ_iω_iω)^-1；

步骤5：根据公式

${\left\{\widetilde{\mathrm{\σ}}\right\}}_d=\underset{i=1}{\overset{N_d}{\mathrm{\Σ}}}{\left\{{\mathrm{\φ}}^{\mathrm{\σ}}\right\}}_i[{\tilde{q}}_i-{\tilde{q}}_{\mathrm{si}}]$

得到虚拟动态应力响应其中{φ^σ}_i为模态应力向量；将虚拟动态应力响应与步骤1得到的机械结构整体等效静力响应相加得到虚拟应力响应量继而得到机械结构随机动态应力功率谱密度