一种基于矩阵摄动理论的复模态随机特征值直接方差求解方法

技术领域

本发明适用于复模态结构系统的特征值分析，用以求解结构系统在经受各种扰动的情况 下，其复模态特征值的统计学性质和变化范围，可为非对称结构系统的复模态特征值分析技 术提供指导。

背景技术

矩阵摄动方法作为一种能够进行快速灵敏度分析和快速结构重分析的实用工具，已经在 基础理论和工程应用中受到了广泛的关注并取得了长足的发展。X.W.YANG和S.H.CHEN将 Padeapproximation应用到矩阵摄动理论中，求出了特征向量和特征值变化量的表达式。 KaminskiM和SoleckaM将结构特征值和特征向量用混沌多项式(PCE)方法进行展开，研究 了线性随机系统的受迫振动响应分析。DebrajG应用基于广义摄动理论的随机有限元方法， 开展了桁架结构的可靠性优化研究。郑兆昌等人基于摄动方法对多自由度系统复模态理论开 展了初步研究。关于结构实模态特征值的统计特性，QiuZ.P和QiuH.C提出了直接方差分析 方法(DVA方法)，无需已知或假定结构参数的相关系数矩阵，通过矩阵摄动理论和概率理论 便能直接计算得到实模态结构的随机特征值的方差。

在实际工程中存在着许多非对称系统。例如船舶结构和大型液体燃料火箭中的燃料箱都 存在显著的流固耦合问题，对其研究中采用流体压力作为流体变量，将引起不对称的矩阵方 程。在研究飞机旋翼、转轴等旋转机构时，由科氏加速度引起的科里奥利力与一个反对称矩 阵相联系，而这些系统都是非对称的。不仅如此，对于具有任意阻尼的系统、气动弹性颤振 系统、非保守力作用下的动力系统和阻尼陀螺等控制系统，其有关的系数矩阵不仅不再是实 对称矩阵，而是复非对称矩阵。非对称结构系统的含义是指质量矩阵、刚度矩阵和阻尼矩阵 等结构矩阵中至少有一个是不对称的。由于系统的不对称性，原系统与转置系统不再相同， 实模态摄动理论将不再适用于该类型的问题，但仍然可以用复模态摄动理论对这类问题进行 分析研究。

目前国内外学者对于矩阵摄动理论的研究和应用大多集中在系统的结构参数为实对称 矩阵时的结构实模态矩阵摄动法，在结构复模态理论的摄动方法、结构复特征值及其统计特 性的研究方面还较少。

发明内容

本发明要解决的技术问题是：克服现有技术不足，提供一种计算结构复模态特征值变化 范围的直接方差求解方法，在进行结构复模态特征值分析时，无需已知或者假定结构参数的 相关系数矩阵，更容易快速准确地得到结构系统的复模态特征值的变化范围，因此极大方便 了在大型结构重分析和结构快速灵敏度分析等领域的工程应用。

本发明采用的技术方案为：适用于非对称结构系统，所述非对称结构系统的含义是指结 构系统中的质量矩阵、刚度矩阵和阻尼矩阵中至少有一个是不对称的，首先根据矩阵摄动理 论，推导了结构在刚度矩阵、阻尼矩阵和质量矩阵等参数发生变化时，结构振动的复模态特 征值和相应特征向量的一阶摄动量，然后基于复模态结构的特征值和特征向量的一阶摄动 量，结合概率理论，从而建立计算结构复模态特征值变化范围的直接方差求解方法，其实现 步骤如下：

第一步：根据矩阵摄动理论，在非对称结构经过微扰，并且其质量矩阵、阻尼矩阵和刚 度矩阵发生变化之后，通过比较特征方程等式两边ε的同次幂系数，并且考虑复模态的正交 关系式，得到非对称结构系统的特征值和特征向量的一阶摄动量的表达式；这一步是建立复 模态特征值方差求解算法的基础，后续的推导过程均是以此展开；

第二步：基于第一步中建立的复模态结构特征值和相应特征向量的一阶摄动量的表达 式，将特征方程中的各结构参数(参数矩阵)，包括刚度矩阵、质量矩阵和阻尼矩阵以及特 征值、特征向量都分成确定性部分和随机扰动部分，结合概率理论，在对复特征值平方(sⁱ)²求期望的基础上，进一步得到复特征值sⁱ的方差Var(sⁱ)的表达式，从而建立起复模态特征值 变化范围(方差)的直接求解算法。

所述第一步具体实现如下：

(11)确定非对称系统结构振动的基本方程{y_j}^T[M(s_i+s_j)+C]{xi}＝δ_ij，其中{y_j}为 左特征向量，{x_i}为右特征向量，M为质量矩阵，C为阻尼矩阵，s_i和s_j为特征方程中的不 同阶数的特征根；

(12)使用复模态的矩阵摄动方法，确定复模态特征值的一阶摄动量 以及复模态特征向量的一阶摄动量 $\left\{u_1^i\right\}=\underset{s=1}{\overset{2N}{\Sigma }}{h}_{is}^1\left\{u_0^s\right\}=-\underset{s=1,i=1,s\ne i}{\overset{2N}{\Sigma }}\left[\frac{1}{s_0^i-s_0^s}{\left\{v_0^s\right\}}^T(B_1+s_0^i{A}_1)\right\{u_0^i\left\}\right]\left\{u_0^s\right\}-\left[\frac{1}{2}{\left\{v_0^i\right\}}^T{A}_1\left\{u_0^i\right\}\right]\left\{u_0^i\right\},$其中s₀表示未 经扰动的初始系统的复特征值，{v₀}和{u₀}分别表示初始系统的左、右状态特征向量，和 分别代表初始系统的复特征值和相应特征向量的一阶摄动量，和分别表示初始 系统的左、右特征向量，$A_1=\left(\begin{array}{cc}0& M_1\\ M_1& C_1\end{array}\right),B_1=\left(\begin{array}{cc}-M_1& 0\\ 0& K_1\end{array}\right),$M₁,C₁,K₁分别表示质量矩阵、 阻尼矩阵和刚度矩阵的一阶扰动量，式中的上角标i和s表示各参数的第i阶和第s阶。

所述第二步具体实现如下：

(21)将特征方程中的各结构参数矩阵包括刚度矩阵、质量矩阵和阻尼矩阵、特征值以 及特征向量均分成确定性部分和随机扰动部分，其中下角标d表示各参数的确定性部分，下 角标r表示各参数的随机扰动部分，ε表示一个小参数 K＝K_d+εK_r,M＝M_d+εM_r,C＝C_d+εC_r,A＝A_d+εA_r,B＝B_d+εB_r,$\left\{y^i\right\}=\left\{y_d^i\right\}+ϵ\left\{y_r^i\right\},\left\{x^i\right\}=\left\{x_d^i\right\}+ϵ\left\{x_r^i\right\}.$

通过以上的表述，为下一步对各结构参数矩阵的期望和方差运算作准备；

(22)结合第一步的求解方法，进行特征值的求解，得到特征值的随机部分和特征向 量的随机部分

$s_r^i=-{\left\{y_d^i\right\}}^T({\left(s_d^i\right)}^2{M}_r+s_d^i{C}_r+K_r)\left\{x_d^i\right\},$

$\left\{u_r^i\right\}=\underset{s=1}{\overset{2N}{\Sigma }}{h}_{is}^1\left\{u_d^s\right\}=-\underset{s=1,i=1,s\ne i}{\overset{2N}{\Sigma }}\left[\frac{1}{s_d^i-s_d^s}{\left\{v_d^s\right\}}^T(B_r+s_d^i{A}_r)\right\{u_d^i\left\}\right]\left\{u_d^s\right\}-\left[\frac{1}{2}{\left\{v_d^i\right\}}^T{A}_r\left\{u_d^i\right\}\right]\left\{u_d^i\right\}.$

(23)由上一步获得的特征值的随机部分和特征向量的随机部分的表达式，结合 概率理论，得到非对称系统复特征值扰动量的方差：

$\left(\begin{array}{c}Var\left(s_r^i\right)=E\left[{\left(s_r^i\right)}^2\right]=E\left[{(-{\left\{y_d^i\right\}}^T({\left(s_d^i\right)}^2{M}_r+s_d^i{C}_r+K_r\left)\right\{x_d^i\left\}\right)}^2\right]\\ =E\left[{\left({\left\{y_d^i\right\}}^T\right({\left(s_d^i\right)}^2{M}_r+s_d^i{C}_r+K_r\left)\right\{x_d^i\left\}\right)}^2\right]\\ ={\left(s_d^i\right)}^{4}E\left[{\left({\left\{y_d^i\right\}}^T{M}_r\right\{x_d^i\left\}\right)}^2\right]+{\left(s_d^i\right)}^{2}E\left[{\left({\left\{y_d^i\right\}}^T{C}_r\right\{x_d^i\left\}\right)}^2\right]+E\left[{\left({\left\{y_d^i\right\}}^T{K}_r\right\{x_d^i\left\}\right)}^2\right]\end{array}\right)$

(24)根据上一步获得的复特征值扰动量的方差，进一步整理，得到复特征值sⁱ的方 差Var(sⁱ)的表达式：

$\begin{array}{c}Var\left(s^i\right)={ϵ}^{2}Var\left(s_r^i\right)\\ ={ϵ}^2\left(\begin{array}{c}({\Theta }^2-4{\Lambda }^2)E({{\tilde{M}}_a}^2-{{\tilde{M}}_b}^2)+\Theta E({{\tilde{C}}_a}^2-{{\tilde{C}}_b}^2)+\\ E({{\tilde{K}}_a}^2-{{\tilde{K}}_b}^2)-4\Lambda E\left({\tilde{M}}_a{\tilde{M}}_b\right)-4\Lambda E\left({\tilde{C}}_a{\tilde{C}}_b\right)\end{array}\right)+\\ {\mathrm{iϵ}}^2\left(\begin{array}{c}2({\Theta }^2-4{\Lambda }^2)E\left({\tilde{M}}_a{\tilde{M}}_b\right)+2\Theta E\left({\tilde{C}}_a{\tilde{C}}_b\right)+2E\left({\tilde{K}}_a{\tilde{K}}_b\right)\\ +4\Lambda \Theta E({{\tilde{M}}_a}^2-{{\tilde{M}}_b}^2)+2\Lambda E({{\tilde{C}}_a}^2-{{\tilde{C}}_b}^2)\end{array}\right)\end{array},$

其中$\left(\begin{array}{c}{\left\{y_{d_1}\right\}}^T{M}_r\left\{x_{d_1}\right\}-{\left\{y_{d_2}\right\}}^T{M}_r\left\{x_{d_2}\right\}={\tilde{M}}_a\\ {\left\{y_{d_1}\right\}}^T\left)M_r\right\{x_{d_2}\}+{\left\{y_{d_2}\right\}}^T{M}_r\{x_{d_1}\}={\tilde{M}}_b\\ {\left\{y_{d_1}\right\}}^T{C}_r\left\{x_{d_1}\right\}-{\left\{y_{d_2}\right\}}^T{C}_r\left\{x_{d_2}\right\}={\tilde{C}}_a\\ {\left\{y_{d_1}\right\}}^T\left)C_r\right\{x_{d_2}\}+{\left\{y_{d_2}\right\}}^T{C}_r\{x_{d_1}\}={\tilde{C}}_b\\ {\left\{y_{d_1}\right\}}^T{K}_r\left\{x_{d_1}\right\}-{\left\{y_{d_2}\right\}}^T{K}_r\left\{x_{d_2}\right\}={\tilde{K}}_a\\ {\left\{y_{d_1}\right\}}^T\left)K_r\right\{x_{d_2}\}+{\left\{y_{d_2}\right\}}^T{K}_r\{x_{d_1}\}={\tilde{K}}_b\end{array}\right),\left(\begin{array}{c}{s}_{d_1}^2-s_{d_2}^2=\Theta \\ s_{d_1}{s}_{d_2}=\Lambda \end{array}\right);$下角标d₁和d₂分别表示各参数的 实数部分和复数部分。

本发明与现有技术相比的优点在于：

(1)本发明克服了在非对称结构系统中，实模态矩阵摄动理论不再适用的困难，采用 复模态矩阵摄动理论对该类特征值问题进行分析研究；

(2)与以往的复模态特征值分析方法相比，本发明无需假定或者已知各结构参数的相 关系数矩阵，就可以快速准确地得到结构系统的复模态特征值的变化范围，应用更加广泛。

附图说明

图1为本发明方法实现流程图；

图2为本发明方法的实施例示意图。

具体实施方式

本发明提出了一种基于矩阵摄动理论的复模态随机特征值直接方差求解方法，其具体实 施步骤是：

第一步：根据矩阵摄动理论，在非对称结构经过微扰，并且其质量矩阵、阻尼矩阵和刚 度矩阵发生变化之后，通过比较特征方程等式两边ε的同次幂系数，并且考虑复模态的正交 关系式，得到非对称结构系统的特征值和特征向量的一阶摄动量的表达式；这一步是建立复 模态特征值方差求解算法的基础，后续的推导过程均是以此展开，下面给出具体的过程：

(1)确定非对称系统结构振动的基本方程

具有N个自由度的线性系统的振动方程为：

$M\stackrel{··}{q}+C\stackrel{·}{q}+Kq=Q\left(t\right),$

对于非对称系统，其自由振动方程为：

$M\stackrel{··}{q}+C\stackrel{·}{q}+Kq=0.$

令q＝{x}e^st，将其代入上式，得到相应的右特征值问题为：

(Ms²+Cs+K){x}＝0,

相应的伴随特征值问题为(Ms²+Cs+K)^T{y}＝0，转置上式得到：

{y}^T(Ms²+Cs+K)＝0.

将向量{x}和{y}分别称为右特征向量和左特征向量。

引入状态向量：

$\left\{u\right\}=\left(\begin{array}{c}sx\\ x\end{array}\right)=\left[T\right]\left\{x\right\},$

其中[T]为状态变换矩阵，并且有：

$\left[T\right]=\left(\begin{array}{c}sI\\ I\end{array}\right),$

相似地引入状态向量：

$\left\{v\right\}=\left(\begin{array}{c}sy\\ y\end{array}\right)=\left[T\right]\left\{y\right\},$

上面的{u}和{v}为对应复模态向量{x}和{y}的状态特征向量。从而得到：

(As+B){u}＝0,

以及：

{v}^T(As+B)＝0,

其中：

$A=\left(\begin{array}{cc}0& M\\ M& C\end{array}\right),B=\left(\begin{array}{cc}-M& 0\\ 0& K\end{array}\right).$

特征问题(As+B){u}＝0和其伴随特征问题{v}^T(As+B)＝0有相同的特征值，其特征方程 为：

det(As+B)＝0.

上式特征方程为2N次的代数方程，在复数域中有2N个特征根s_i(i＝1,2,...,2N)，对于每个 s_i，其左、右模态向量{v_i}和{u_i}应满足：

(As_i+B){u_i}＝0,

和

{v_i}^T(As_i+B)＝0.

根据正交关系有：

{v_j}^TA{u_i}＝δ_ij,

{v_j}^TB{u_i}＝-s_iδ_ij.

综合以上各式得到：

{y_j}^T[T_j]^TA[T_i]{x_i}＝δ_ij，

{y_j}^T[T_j]^TB[T_i]{x_i}＝-s_iδ_ij

将$A=\left(\begin{array}{cc}0& M\\ M& C\end{array}\right),B=\left(\begin{array}{cc}-M& 0\\ 0& K\end{array}\right)$代入上式可得：

{y_j}^T[M(s_i+s_j)+C]{x_i}＝δ_ij,

{y_j}^T[-Ms_is_j+K]{x_i}＝-s_iδ_ij.

(2)确定复模态特征值和特征向量的一阶摄动量

结构参数的变化是通过描述系统的质量、阻尼和刚度矩阵的改变而实现的，故设结构在 经过微扰之后的质量矩阵、阻尼矩阵和刚度矩阵分别为：

$\left(\begin{array}{c}M=M_0+ϵM_1\\ C=C_0+{\mathrm{ϵC}}_1\\ K=K_0+ϵK_1\end{array}\right),$

由上式结合$A=\left(\begin{array}{cc}0& M\\ M& C\end{array}\right),B=\left(\begin{array}{cc}-M& 0\\ 0& K\end{array}\right),$可得：

$\left(\begin{array}{c}A=A_0+{\mathrm{ϵA}}_1\\ B=B_0+{\mathrm{ϵB}}_1\end{array}\right).$

式中ε是一个小参数，而ε＝0对应的系统称为原系统。M₀、C₀和K₀是原系统的质量、阻尼 和刚度矩阵。εM₁、εC₁和εK₁表示各矩阵相应微小的变化，而且满足：

$\left(\begin{array}{c}M\to M_0\\ C\to C_0\\ K\to K_0\\ A\to A_0\\ B\to B_0\end{array}\right),s.t.\left(\begin{array}{c}ϵM_1\to 0\\ ϵC_1\to 0\\ ϵK_1\to 0\\ ϵA_1\to 0\\ ϵB_1\to 0\end{array}\right).$

本发明中讨论原系统的特征值是特征方程为单根的情形。

根据矩阵摄动理论，将特征值和特征向量按小参数ε展开为如下的幂级数形式：

$s^i=s_0^i+{\mathrm{ϵs}}_1^i+{ϵ}^2{s}_2^i+...,$

$\left\{u^i\right\}=\left\{u_0^i\right\}+ϵ\left\{u_1^i\right\}+{ϵ}^2\left\{u_2^i\right\}+...,$

$\left\{v^i\right\}=\left\{v_0^i\right\}+ϵ\left\{v_1^i\right\}+{ϵ}^2\left\{v_2^i\right\}+...,$

$\left\{x^i\right\}=\left\{x_0^i\right\}+ϵ\left\{x_1^i\right\}+{ϵ}^2\left\{x_2^i\right\}+...,$

$\left\{y^i\right\}=\left\{y_0^i\right\}+ϵ\left\{y_1^i\right\}+{ϵ}^2\left\{y_2^i\right\}+....$

综合以上各式得到：

$\left(\right(A_0+{\mathrm{ϵA}}_1\left)\right(s_0^i+{\mathrm{ϵs}}_1^i+{ϵ}^2{s}_2^i+...)+(B_0+{\mathrm{ϵB}}_1\left)\right)\left(\right\{u_0^i\}+ϵ\{u_1^i\}+{ϵ}^2\{u_2^i\}+...)=0,$

${\left(\right\{v_0^i\}+ϵ\{v_1^i\}+{ϵ}^2\{v_2^i\}+...)}^T\left(\right(A_0+{\mathrm{ϵA}}_1\left)\right(s_0^i+{\mathrm{ϵs}}_1^i+{ϵ}^2{s}_2^i+...)+(B_0+{\mathrm{ϵB}}_1\left)\right)=0.$

将$\left(\right(A_0+{\mathrm{ϵA}}_1\left)\right(s_0^i+{\mathrm{ϵs}}_1^i+{ϵ}^2{s}_2^i+...)+(B_0+{\mathrm{ϵB}}_1\left)\right)\left(\right\{u_0^i\}+ϵ\{u_1^i\}+{ϵ}^2\{u_2^i\}+...)=0$展开并略去O(ε³)项 之后，比较ε的同次幂系数可得：

${ϵ}^0:(A_0{s}_0^i+B_0)\left\{u_0^i\right\}=0,$

${ϵ}^1:(B_0+s_0^i{A}_0)\left\{u_1^i\right\}+B_1\left\{u_0^i\right\}+s_0^i{A}_1\left\{u_o^i\right\}+s_1^i{A}_0\left\{u_o^i\right\}=0,$

${ϵ}^2:(B_0+s_0^i{A}_0)\left\{u_2^i\right\}+(B_1+s_0^i{A}_1)\left\{u_1^i\right\}+s_1^i{A}_0\left\{u_1^i\right\}+s_1^i{A}_1\left\{u_0^i\right\}+s_2^i{A}_0\left\{u_0^i\right\}=0.$

同理，将${\left(\right\{v_0^i\}+ϵ\{v_1^i\}+{ϵ}^2\{v_2^i\}+...)}^T\left(\right(A_0+{\mathrm{ϵA}}_1\left)\right(s_0^i+{\mathrm{ϵs}}_1^i+{ϵ}^2{s}_2^i+...)+(B_0+{\mathrm{ϵB}}_1\left)\right)=0$展开并略去 O(ε³)项之后，比较ε的同次幂系数可得：

${ϵ}^0:{(A_0{s}_0^i+B_0)}^T\left\{v_0^i\right\}=0,$

${ϵ}^1:{(B_0+s_0^i{A}_0)}^T\left\{v_1^i\right\}+{(B_1+s_0^i{A}_1+s_1^i{A}_0)}^T\left\{v_o^i\right\}=0,$

${ϵ}^2:{(B_0+s_0^i{A}_0)}^T\left\{v_2^i\right\}+{(B_1+s_0^i{A}_1)}^T\left\{v_1^i\right\}+s_1^i{A_0}^T\left\{v_1^i\right\}+s_1^i{A_1}^T\left\{v_0^i\right\}+s_2^i{A_0}^T\left\{v_0^i\right\}=0.$

由此可以求出特征值的一阶摄动量和二阶摄动量以及左、右特征向量的一阶摄动量和 二阶摄动量。

将按照原系统的右特征向量展开如下：

$\left\{u_1^i\right\}=\underset{s=1}{\overset{2N}{\Sigma }}{h}_{is}^1\left\{u_0^s\right\}.$

将上式代入ε¹对应的方程有：

$(B_0+s_0^i{A}_0)\underset{s=1}{\overset{2N}{\Sigma }}{h}_{is}^1\left\{u_0^s\right\}+(B_1+s_0^i{A}_1)\left\{u_0^i\right\}+s_1^i{A}_0\left\{u_0^i\right\}=0,$

用左特征向量的转置左乘上式，有：

${\left\{v_0^s\right\}}^T(B_0+s_0^i{A}_0)\underset{s=1}{\overset{2N}{\Sigma }}{h}_{is}^1\left\{u_0^s\right\}+{\left\{v_0^s\right\}}^T\{B_1+s_0^i{A}_1\}\left\{u_0^i\right\}+{\left\{v_0^s\right\}}^T{A}_0\left\{u_0^i\right\}{s}_1^i=0,$

由复模态正交关系式，可以将上式化为：

$h_{is}^1(s_0^i-s_0^s)+{\left\{v_0^s\right\}}^T(B_1+s_0^i{A}_1)\left\{u_0^i\right\}+s_1^i{\delta }_{is}=0.$

当s＝i时，由上式可得：

$s_1^i=-{\left\{v_0^i\right\}}^T(B_1+s_0^i{A}_1)\left\{u_0^i\right\},$

上式又可写成：

$s_1^i=-{\left\{y_0^i\right\}}^T({\left(s_0^i\right)}^2{M}_1+s_0^i{C}_1+K_1)\left\{x_0^i\right\},$

当s≠i时，δ_is＝0，由式$h_{is}^1(s_0^i-s_0^s)+{\left\{v_0^s\right\}}^T(B_1+s_0^i{A}_1)\left\{u_0^i\right\}+s_1^i{\delta }_{is}=0$可得：

$h_{is}^1=-\frac{1}{s_0^i-s_0^s}{\left\{v_0^s\right\}}^T(B_1+s_0^i{A}_1)\left\{u_0^i\right\},(i\ne s,i,s=1,2,...,2N).$

同理，此式又可以写成：

$h_{is}^1=-\frac{1}{s_0^i-s_0^s}{\left\{y_0^s\right\}}^T({\left(s_0^i\right)}^2{M}_1+s_0^i{C}_1+K_1)\left\{x_0^i\right\},(i\ne s,i,s=1,2,...,2N).$

另一方面，当s＝i，系数则由模态的正则化条件来确定。

特征向量的正则化条件为：

{vⁱ}^TA{uⁱ}＝1,

因此可得：

${\left(\right\{v_0^i\}+ϵ\{v_1^i\}+{ϵ}^2\{v_2^i\}+...)}^T(A_0+{\mathrm{ϵA}}_1)\left(\right\{u_0^i\}+ϵ\{u_1^i\}+{ϵ}^2\{u_2^i\}+...)=1$

展开上式并略去O(ε³)项，比较ε同次幂系数得：

${ϵ}^0:{\left\{v_0^i\right\}}^T{A}_0\left\{u_0^i\right\}=1,$

${ϵ}^1:{\left\{v_0^i\right\}}^T{A}_0\left\{u_1^i\right\}+{\left\{v_0^i\right\}}^T{A}_1\left\{u_0^i\right\}+{\left\{v_1^i\right\}}^T{A}_0\left\{u_0^i\right\}=0,$

${ϵ}^2:{\left\{v_0^i\right\}}^T{A}_0\left\{u_2^i\right\}+{\left\{v_0^i\right\}}^T{A}_0\left\{u_1^i\right\}+{\left\{v_1^i\right\}}^T{A}_0\left\{u_1^i\right\}+{\left\{v_1^i\right\}}^T{A}_1\left\{u_0^i\right\}+{\left\{v_2^i\right\}}^T{A}_0\left\{u_0^i\right\}=0.$

用左乘并且根据正交关系式，得到：

$h_{ii}^1={\left\{v_0^i\right\}}^T{A}_0\left\{u_1^i\right\}.$

与相似，将按照原模态展开如下：

${\left\{v_1^i\right\}}^T=\underset{s=1}{\overset{2N}{\Sigma }}{d}_{is}^1{\left\{v_0^s\right\}}^T,$

用右乘上式，并且根据正交关系式，得到：

$d_{ii}^1={\left\{v_1^i\right\}}^T{A}_0{u}_0^i.$

将和代入ε¹对应的方程，得到：

$h_{ii}^1+d_{ii}^1=-{\left\{v_0^i\right\}}^T{A}_1\left\{u_0^i\right\}.$

不妨取则有：

$h_{ii}^1=d_{ii}^1=-\frac{1}{2}{\left\{v_0^i\right\}}^T{A}_1\left\{u_0^i\right\}.$

因此有：

$\left\{u_1^i\right\}=\underset{s=1}{\overset{2N}{\Sigma }}{h}_{is}^1\left\{u_0^s\right\}=-\underset{s=1,i=1,s\ne i}{\overset{2N}{\Sigma }}\left[\frac{1}{s_0^i-s_0^s}{\left\{v_0^s\right\}}^T(B_1+s_0^i{A}_1)\right\{u_0^i\left\}\right]\left\{u_0^s\right\}-\left[\frac{1}{2}{\left\{v_0^i\right\}}^T{A}_1\left\{u_0^i\right\}\right]\left\{u_0^i\right\}.$

第二步：基于第一步中建立的复模态结构特征值和相应特征向量的一阶摄动量的表达 式，将特征方程中的各结构参数(参数矩阵)，包括刚度矩阵、质量矩阵和阻尼矩阵以及特 征值、特征向量都分成确定性部分和随机扰动部分，结合概率理论，在对复特征值平方(sⁱ)²求期望的基础上，进一步得到复特征值sⁱ的方差Var(sⁱ)的表达式，从而建立起复模态特征值 变化范围(方差)的直接求解算法。具体实施步骤如下：

(1)首先，把刚度矩阵K、质量矩阵M、阻尼矩阵C、复特征值sⁱ，左特征向量{y}和 右特征向量{x}表示为：

K＝K_d+εK_r,

M＝M_d+εM_r,

C＝C_d+εC_r,

A＝A_d+εA_r,

B＝B_d+εB_r,

$s^i=s_d^i+{\mathrm{ϵs}}_r^i,$

$\left\{y^i\right\}=\left\{y_d^i\right\}+ϵ\left\{y_r^i\right\},$

$\left\{x^i\right\}=\left\{x_d^i\right\}+ϵ\left\{x_r^i\right\}.$

不妨对以上各式做出假设，设ε均为一个小参数，K_d、M_d、C_d、A_d、B_d、为相应的矩阵或向量中的确定性部分，K_r、M_r、C_r、A_r、B_r、为相应 的矩阵或向量的随机部分，并且满足均值为零。

(2)对以上各式取数学期望得到：

E[K]＝E[K_d]+εE[K_r]＝K_d,

E[M]＝E[M_d]+εE[M_r]＝M_d,

E[C]＝E[C_d]+εE[C_r]＝C_d,

E[A]＝E[A_d]+εE[A_r]＝A_d,

E[B]＝E[B_d]+εE[B_r]＝B_d,

$E\left[s^i\right]=E\left[s_d^i\right]+ϵE\left[s_r^i\right]=s_d^i,$

$E\left[\right\{y^i\left\}\right]=E\left[\right\{y_d^i\left\}\right]+ϵE\left[\right\{y_r^i\left\}\right]=\left\{y_d^i\right\},$

$E\left[\right\{x^i\left\}\right]=E\left[\right\{x_d^i\left\}\right]+ϵE\left[\right\{x_r^i\left\}\right]=\left\{x_d^i\right\}.$

对左右两边取平方得：

${\left(s^i\right)}^2={\left(s_d^i\right)}^2+2{\mathrm{ϵs}}_d^i{s}_r^i+{ϵ}^2{\left(s_r^i\right)}^2,$

对上式求数学期望：

$E\left[{\left(s^i\right)}^2\right]=E\left[{\left(s_d^i\right)}^2\right]+{ϵ}^{2}E\left[{\left(s_r^i\right)}^2\right],$

由概率理论可以得到复特征值sⁱ的方差满足：

Var(sⁱ)＝E[(sⁱ)²]-(E[sⁱ])²,

联立以上两式得到：

$Var\left(s^i\right)={ϵ}^{2}E\left[{\left(s_r^i\right)}^2\right].$

采用本发明第一步中的求解方法可以容易得到如下结果：

$s_r^i=-{\left\{y_d^i\right\}}^T({\left(s_d^i\right)}^2{M}_r+s_d^i{C}_r+K_r)\left\{x_d^i\right\},$

$\left\{u_r^i\right\}=\underset{s=1}{\overset{2N}{\Sigma }}{h}_{is}^1\left\{u_d^s\right\}=-\underset{s=1,i=1,s\ne i}{\overset{2N}{\Sigma }}\left[\frac{1}{s_d^i-s_d^s}{\left\{v_d^s\right\}}^T(B_r+s_d^i{A}_r)\right\{u_d^i\left\}\right]\left\{u_d^s\right\}-\left[\frac{1}{2}{\left\{v_d^i\right\}}^T{A}_r\left\{u_d^i\right\}\right]\left\{u_d^i\right\}.$

(3)下面推导sⁱ的方差，由可知sⁱ的数学期望是

$\begin{array}{c}Var\left(s^i\right)=Var(s_d^i+{\mathrm{ϵs}}_r^i)=Var\left(s_d^i\right)+Var\left({\mathrm{ϵs}}_r^i\right)=0+Var\left({\mathrm{ϵs}}_r^i\right)\\ ={ϵ}^{2}Var\left(s_r^i\right)={ϵ}^{2}E\left[{(s_r^i-E(s_r^i\left)\right)}^2\right]={ϵ}^{2}E\left[{\left(s_r^i\right)}^2\right]\end{array}.$

因此有：

$Var\left(s_r^i\right)=E\left[{\left(s_r^i\right)}^2\right],$

将$s_r^i=-{\left\{y_d^i\right\}}^T({\left(s_d^i\right)}^2{M}_r+s_d^i{C}_r+K_r)\left\{x_d^i\right\}$代入上式得：

$\left(\begin{array}{c}Var\left(s_r^i\right)=E\left[{\left(s_r^i\right)}^2\right]=E\left[{(-{\left\{y_d^i\right\}}^T({\left(s_d^i\right)}^2{M}_r+s_d^i{C}_r+K_r\left)\right\{x_d^i\left\}\right)}^2\right]\\ =E\left[{(-{\left\{y_d^i\right\}}^T({\left(s_d^i\right)}^2{M}_r+s_d^i{C}_r+K_r\left)\right\{x_d^i\left\}\right)}^2\right]\end{array}\right)$

考虑到E[M_r]＝E[C_r]＝E[K_r]＝0，故上式中交叉项的数学期望为零，可以化为：

$\left(\begin{array}{c}Var\left(s_r^i\right)=E\left[{\left(s_r^i\right)}^2\right]=E\left[{(-{\left\{y_d^i\right\}}^T({\left(s_d^i\right)}^2{M}_r+s_d^i{C}_r+K_r\left)\right\{x_d^i\left\}\right)}^2\right]\\ =E\left[{\left({\left\{y_d^i\right\}}^T\right({\left(s_d^i\right)}^2{M}_r+s_d^i{C}_r+K_r\left)\right\{x_d^i\left\}\right)}^2\right]\\ ={\left(s_d^i\right)}^{4}E\left[{\left({\left\{y_d^i\right\}}^T{M}_r\right\{x_d^i\left\}\right)}^2\right]+{\left(s_d^i\right)}^{2}E\left[{\left({\left\{y_d^i\right\}}^T{C}_r\right\{x_d^i\left\}\right)}^2\right]+E\left[{\left({\left\{y_d^i\right\}}^T{K}_r\right\{x_d^i\left\}\right)}^2\right]\end{array}\right)$

(4)为了推导的方便，去掉上式中各项的上角标i。需要说明的是，上式中各项中只 有特征值s_d，特征向量{y_d}^T和{x_d}会出现为复数的情况，而各参数矩阵的一阶摄动量如 M_r、C_r、K_r则均为实数矩阵。因此，令复特征值s_d，左特征向量{y_d}^T，右特征向量{x_d} 分别满足：

$\left(\begin{array}{c}{s}_d=s_{d_1}+i{s}_{d_2}\\ {\left\{y_d\right\}}^T={\left\{y_{d_1}\right\}}^T+i{\left\{y_{d_2}\right\}}^T,\\ \left\{x_d\right\}=\left\{x_{d_1}\right\}+i\left\{x_{d_2}\right\}\end{array}\right)$

由上式，可得到：

$\begin{array}{c}Var\left(s_r^i\right)={\left(s_d^i\right)}^{4}E\left[{\left({\left\{y_d^i\right\}}^T{M}_r\right\{x_d^i\left\}\right)}^2\right]+{\left(s_d^i\right)}^{2}E\left[{\left({\left\{y_d^i\right\}}^T{C}_r\right\{x_d^i\left\}\right)}^2\right]+E\left[{\left({\left\{y_d^i\right\}}^T{K}_r\right\{x_d^i\left\}\right)}^2\right].\\ ={(s_{d_1}+{\mathrm{is}}_{d_2})}^{4}E\left[{\left(({\left\{y_{d_1}\right\}}^T+i{\left\{y_{d_2}\right\}}^T)M_r\left(\right\{x_{d_1}\}+i\{x_{d_2}\left\}\right)\right)}^2\right]+\\ {(s_{d_1}+{\mathrm{is}}_{d_2})}^{2}E\left[{\left(({\left\{y_{d_1}\right\}}^T+i{\left\{y_{d_2}\right\}}^T)C_r\left(\right\{x_{d_1}\}+i\{x_{d_2}\left\}\right)\right)}^2\right]+\\ E\left[{\left(({\left\{y_{d_1}\right\}}^T+i{\left\{y_{d_2}\right\}}^T)K_r\left(\right\{x_{d_1}\}+i\{x_{d_2}\left\}\right)\right)}^2\right]\\ ={(s_{d_1}+{\mathrm{is}}_{d_2})}^{4}E\left[{(\left({\left\{y_{d_1}\right\}}^T{M}_r\right\{x_{d_1}\}-{\left\{y_{d_2}\right\}}^T{M}_r\{x_{d_2}\left\}\right)+i\left({\left\{y_{d_1}\right\}}^T\right)M_r\{x_{d_2}\}+{\left\{y_{d_2}\right\}}^T{M}_r\{x_{d_1}\left\}\right)}^2\right]\\ {(s_{d_1}+{\mathrm{is}}_{d_2})}^{2}E\left[{(\left({\left\{y_{d_1}\right\}}^T{C}_r\right\{x_{d_1}\}-{\left\{y_{d_2}\right\}}^T{C}_r\{x_{d_2}\left\}\right)+i\left({\left\{y_{d_1}\right\}}^T\right)C_r\{x_{d_2}\}+{\left\{y_{d_2}\right\}}^T{C}_r\{x_{d_1}\left\}\right)}^2\right]+\\ E\left[{(\left({\left\{y_{d_1}\right\}}^T{K}_r\right\{x_{d_1}\}-{\left\{y_{d_2}\right\}}^T{K}_r\{x_{d_2}\left\}\right)+i\left({\left\{y_{d_1}\right\}}^T\right)K_r\{x_{d_2}\}+{\left\{y_{d_2}\right\}}^T{K}_r\{x_{d_1}\left\}\right)}^2\right]\end{array}.$

令

$\left(\begin{array}{c}{\left\{y_{d_1}\right\}}^T{M}_r\left\{x_{d_1}\right\}-{\left\{y_{d_2}\right\}}^T{M}_r\left\{x_{d_2}\right\}={\tilde{M}}_a\\ {\left\{y_{d_1}\right\}}^T\left)M_r\right\{x_{d_2}\}+{\left\{y_{d_2}\right\}}^T{M}_r\{x_{d_1}\}={\tilde{M}}_b\\ {\left\{y_{d_1}\right\}}^T{C}_r\left\{x_{d_1}\right\}-{\left\{y_{d_2}\right\}}^T{C}_r\left\{x_{d_2}\right\}={\tilde{C}}_a\\ {\left\{y_{d_1}\right\}}^T\left)C_r\right\{x_{d_2}\}+{\left\{y_{d_2}\right\}}^T{C}_r\{x_{d_1}\}={\tilde{C}}_b\\ {\left\{y_{d_1}\right\}}^T{K}_r\left\{x_{d_1}\right\}-{\left\{y_{d_2}\right\}}^T{K}_r\left\{x_{d_2}\right\}={\tilde{K}}_a\\ {\left\{y_{d_1}\right\}}^T\left)K_r\right\{x_{d_2}\}+{\left\{y_{d_2}\right\}}^T{K}_r\{x_{d_1}\}={\tilde{K}}_b\end{array}\right)$

并且

$\left(\begin{array}{c}{s}_{d_1}^2-s_{d_2}^2=\Theta \\ s_{d_1}{s}_{d_2}=\Lambda \end{array}\right),$

联立并整理以上三式，得到：

$\left(\begin{array}{c}Var\left(s_r^i\right)=\left(\begin{array}{c}({\Theta }^2-4{\Lambda }^2)E({{\tilde{M}}_a}^2-{{\tilde{M}}_b}^2)+\Theta E({{\tilde{C}}_a}^2-{{\tilde{C}}_b}^2)+\\ E({{\tilde{K}}_a}^2-{{\tilde{K}}_b}^2)-4\Lambda E\left({\tilde{M}}_a{\tilde{M}}_b\right)-4\Lambda E\left({\tilde{C}}_a{\tilde{C}}_b\right)\end{array}\right)+\\ i\left(\begin{array}{c}2({\Theta }^2-4{\Lambda }^2)E\left({\tilde{M}}_a{\tilde{M}}_b\right)+2\Theta E\left({\tilde{C}}_a{\tilde{C}}_b\right)+2E\left({\tilde{K}}_a{\tilde{K}}_b\right)\\ +4\Lambda \Theta E({{\tilde{M}}_a}^2-{{\tilde{M}}_b}^2)+2\Lambda E({{\tilde{C}}_a}^2-{{\tilde{C}}_b}^2)\end{array}\right)\end{array}\right)$

又由$\begin{array}{c}Var\left(s^i\right)=Var(s_d^i+{\mathrm{ϵs}}_r^i)=Var\left(s_d^i\right)+Var\left({\mathrm{ϵs}}_r^i\right)=0+Var\left({\mathrm{ϵs}}_r^i\right)\\ ={ϵ}^{2}Var\left(s_r^i\right)={ϵ}^{2}E\left[{(s_r^i-E(s_r^i\left)\right)}^2\right]={ϵ}^{2}E\left[{\left(s_r^i\right)}^2\right]\end{array}.$可得：

$\left(\begin{array}{c}Var\left(s^i\right)={ϵ}^{2}Var\left(s_r^i\right)\\ ={ϵ}^2\left(\begin{array}{c}({\Theta }^2-4{\Lambda }^2)E({{\tilde{M}}_a}^2-{{\tilde{M}}_b}^2)+\Theta E({{\tilde{C}}_a}^2-{{\tilde{C}}_b}^2)+\\ E({{\tilde{K}}_a}^2-{{\tilde{K}}_b}^2)-4\Lambda E\left({\tilde{M}}_a{\tilde{M}}_b\right)-4\Lambda E\left({\tilde{C}}_a{\tilde{C}}_b\right)\end{array}\right)+\\ {\mathrm{iϵ}}^2\left(\begin{array}{c}2({\Theta }^2-4{\Lambda }^2)E\left({\tilde{M}}_a{\tilde{M}}_b\right)+2\Theta E\left({\tilde{C}}_a{\tilde{C}}_b\right)+2E\left({\tilde{K}}_a{\tilde{K}}_b\right)\\ +4\Lambda \Theta E({{\tilde{M}}_a}^2-{{\tilde{M}}_b}^2)+2\Lambda E({{\tilde{C}}_a}^2-{{\tilde{C}}_b}^2)\end{array}\right)\end{array}\right)$

综上，可以通过复模态的矩阵摄动理论，结合概率统计方法，直接得到复模态特征值的 统计学性质(复模态变化范围)，提高了相关方法的适用范围。

实施例：

为了更充分的了解该发明的特点及其对工程实际的适用性，本发明以图2的结构系统为 例进行复模态的随机特征值分析验证。图2中c₁,c₂,c₃分别代表系统中三个阻尼器的阻尼系 数，k代表系统中弹簧的刚性系数，m表示滑块的质量，x₁,x₂分别表示系统中两个滑块的位 置坐标。

考虑两自由度振动系统，满足c＝1，k＝9，m＝1，其中阻尼系数c₁＝c₂＝c₃＝c；利用 达朗伯原理容易建立系统的运动微分方程：

$\left(\begin{array}{cc}m& 0\\ 0& m\end{array}\right)\left(\begin{array}{c}{\stackrel{··}{x}}_1\\ {\stackrel{··}{x}}_2\end{array}\right)+\left(\begin{array}{cc}{c}_1+c_2& -c_2\\ -c_2& c_2+c_3\end{array}\right)\left(\begin{array}{c}{\stackrel{·}{x}}_1\\ {\stackrel{·}{x}}_2\end{array}\right)+\left(\begin{array}{cc}2k& -k\\ -k& 2k\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)=0.$

系统的状态向量{u}、矩阵A和矩阵B为：

$\left\{u\right\}=\left(\begin{array}{c}{\stackrel{·}{x}}_1\\ {\stackrel{·}{x}}_2\\ x_1\\ x_2\end{array}\right),A=\left(\begin{array}{cccc}0& 0& m& 0\\ 0& 0& 0& m\\ m& 0& 2c& -c\\ 0& m& -c& 2c\end{array}\right),B=\left(\begin{array}{cccc}-m& 0& 0& 0\\ 0& -m& 0& 0\\ 0& 0& 2k& -k\\ 0& 0& -k& 2k\end{array}\right),$