CN106021711A - 一种面向密集频率结构振动特征值的随机摄动方法 - Google Patents

Abstract

本发明公开了一种面向密集频率结构振动特征值的随机摄动方法，该方法首先对刚度矩阵和质量矩阵进行谱分解，然后对特征值移位，将密频系统化为重频系统。再对修改后的特征值进行摄动分析，获得当结构参数发生扰动后，关于密频结构振动特征值一阶摄动量的矩阵方程。然后由多项式混沌展开方法，构建关于密频特征值一阶摄动量的代理模型。结合摄动方法和代理模型技术，提出了面向密集频率结构振动特征值的近似计算方法，并基于该近似计算方法，进一步得到了密频结构特征值在参数扰动情况下均值和方差的表达式。本发明解决了传统摄动方法中，密频结构特征值无法被结构参数所显式表达，因而无法直接研究其统计特征的困难。

Description

一种面向密集频率结构振动特征值的随机摄动方法

技术领域

本发明适用于固有频率密集分布的结构系统的振动特征值分析，用以求解固有频率密集分布的结构系统在经受各种参数扰动的情况下，其振动特征值的统计学特性和特征值变化范围，本发明可为密集频率结构系统的特征值分析技术、模型更新、结构分析和设计优化提供指导。

背景技术

在实际工程结构中，密集频率结构或者重复频率结构较为常见。而事实上，密集频率结构比重复频率结构更加常见，因为重频结构可以视为密集频率结构的一个特例，亦即密集分布的频率间隔退化到零。这里应该强调，所谓密集频率指的是在很窄的一段频率间隔内有很多个固有频率分布其中。一般来说，具有严格意义上的重复频率结构系统，当它的结构参数经历诸如外界扰动等微小摄动时，就将转化为密集频率结构。局部对称结构，准对称结构及其相应的振动控制系统，它们的振动固有频率和闭环特征值很有可能是密集分布的。一维结构的振动频率往往有较大的间隔，三维结构尤其是三维大型柔性结构，其振动固有频率一般非常密集。对于大型空间柔性结构，如平流层飞艇、大型建筑膜结构来说，其结构振动的一个主要的特征就是固有频率低且非常密集。

Ghosh和Ghanem利用不变子空间方法研究了密频结构的随机特征值问题。Zhao和Xu等人针对重频阻尼振动系统，提出了基于松弛组合近似方法求解重分析问题的快速算法。通过将特征向量表示成为基向量和系数向量的组合形式，避免了求解大规模方程组，简化了复杂求解运算。关于结构实模态特征值的统计特性，Qiu Z.P和Qiu H.C提出了直接方差分析方法(DVA方法)，无需已知或假定结构参数的相关系数矩阵，通过矩阵摄动理论和概率理论便能直接计算得到实模态结构的随机特征值的方差。

研究密集频率结构振动的矩阵摄动方法，在处理实际工程问题中，更加具有普遍意义。矩阵摄动方法是处理结构模态灵敏度分析、结构更新、结构重分析和设计优化的主要方法之一。然而当工程结构系统的固有振动频率密集分布时，现有的孤立特征值摄动法和重频特征值摄动方法都不能直接应用了。从数学角度来说，密集特征值依然还是孤立的特征值，亦即特征方程的单根，只是特征值之间的间隔很小。但是在实际实施结构重分析和结构灵敏度分析时，传统的孤立特征值的摄动方法不能用于密集频率特征值问题，原因主要有两个：一是摄动展开式中由于特征值间隔很小，其展开式的一阶摄动量和二阶摄动量可能都会失效，即便是在扰动很小的情况下，虽然能确保摄动展开式的收敛性，也会由于过大的截断误差而使得处理结果失去可靠性。二是对于固有频率密集的结构系统，其频率聚集组的每个特征向量都不容易算准，因为根据代数中的特征值理论，实对称矩阵的孤立特征值是良态的，但是它的特征向量却很可能是病态的，其病态程度取决于特征值之间的密集程度。

正是由于摄动方法在处理密集频率结构的振动特征值问题中存在的上述困难，目前国内外学者对于摄动方法处理结构振动特征值问题的研究大多集中在孤立特征值问题上。因此，目前关于摄动方法处理密集频率结构的振动特征值问题的研究工作还很少。

发明内容

本发明要解决的技术问题是：克服现有技术的不足，提供一种面向密集频率结构振动特征值的随机摄动方法。根据此方法来计算密集频率结构系统，当其结构参数发生变化和扰动之后，其振动特征值的近似表达式。进而基于该近似表达式，得到密集频率结构的振动特征值在随机参数扰动情况下的均值和方差。本发明解决了传统的随机摄动方法中，密频结构特征值的一阶摄动量无法被显式表达，因而无法进一步直接研究密频特征值的统计学特征的困难，同时本发明解决了基于泰勒级数展开的摄动方法只能处理小范围变化的设计参数的应用瓶颈。本发明方法只需少数的配点计算就能建立较为精确的特征值近似模型，避免了蒙特卡洛方法大规模高成本的样本点计算问题。因此本发明有力地推动了密集频率结构的动态分析、模型更新和设计优化等工作。

本发明采用的技术方案为：一种面向密集频率结构振动特征值的随机摄动方法，适用于具有密集分布的自由振动频率的结构系统。首先利用特征值移位方法，将全部密集特征值移位至它们的平均值，这样原密频系统的特征值问题转化为重频特征值系统基础上的摄动问题。修改后系统的特征值问题转化为重特征值系统基础上，经过两部分摄动叠加后的特征值问题，其中一部分是系统参数矩阵扰动引起的摄动，另一部分是特征值移位引起的摄动。然后根据重频结构的摄动理论，推导了结构系统在刚度矩阵、质量矩阵等参数发生随机扰动后，关于其自由振动特征值的一阶摄动量的矩阵方程。然后引入多项式混沌展开技术，建立面向密频结构自由振动特征值一阶摄动量的代理模型。进一步结合多项式混沌展开技术和矩阵摄动方法，给出了密集频率结构自由振动特征值的近似计算方法。基于该计算方法，得到了密频结构特征值在参数扰动情况下的均值和方差的表达式，其实现步骤如下：

第一步：根据结构振动的特征值理论，给出刚度矩阵和质量矩阵的谱分解；

第二步：利用移位方法，将原密频特征值系统转化为重频特征值系统，修改后系统的特征值问题化为重特征值系统基础上，经过两部分摄动叠加后的特征值问题，其中一部分是系统参数扰动引起的摄动，另一部分是特征值移位引起的摄动。再根据重频结构摄动理论，针对修改后的系统，推导了该系统自由振动特征值的一阶摄动量的矩阵方程；

第三步：针对第二步构建的矩阵方程，其方程的根亦即特征值一阶摄动量无法被结构参数所显式表达的困难，引入多项式混沌展开技术，将上述密频特征值的一阶摄动量表示成为关于系统中结构参数的混沌多项式展式的形式，建立关于密频特征值一阶摄动量的代理模型；

第四步：结合摄动方法和基于多项式混沌展开的代理模型技术，提出了面向密集频率结构自由振动特征值的近似计算方法，基于该近似计算方法，进一步得到了密集频率结构自由振动特征值在参数扰动情况下的均值和方差的表达式。

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

步骤(11)、给出密频结构系统中质量矩阵M₀在刚度矩阵K₀上的谱分解表示其中K₀和M₀分别是n×n维实对称刚度矩阵和质量矩阵，Λ₀和U₀是它的一组密集特征值对角阵及其特征向量矩阵，Λ₀和U₀的维数分别为m×m和n×m，Λ_A和U_A也是原系统的非密集分布的特征值对角阵及其对应的特征向量矩阵；

步骤(12)、给出密频结构系统中刚度矩阵K₀在质量矩阵M₀上的谱分解表示，

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

步骤(21)、给出移位后的重频特征值表达式其中为原系统各个密集频率值，m为原系统中一组密集固有频率的个数；

步骤(22)、推导了修改后系统的自由振动特征值一阶摄动量的矩阵方程其中α表示相应的系数矩阵，Λ₁表示密频特征值一阶摄动项的对角矩阵，δK₀为特征值移位引起的摄动，K₁和M₁分别表示结构刚度矩阵和结构质量矩阵的一阶摄动项；

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

步骤(31)、针对结构参数ξ₁,ξ₂,…,ξ_n满足高斯分布的情况，使用Hermite多项式展式，将表示为式中表示待定的系数矢量，ξ＝(ξ₁,ξ₂,…,ξ_n)为服从高斯分布的随机变量，表示n次多维Hermite多项式；

步骤(32)、将步骤(31)中的Y(ξ)用有限项截断并用二阶Hermite多项式近似展开，得到上式中n为随机变量的维数，c_0,2,c_i,2,c_ii,2,c_ij,2表示最高阶次为二阶的Hermite多项式展开式的待定系数；

步骤(33)、选择好配点组合，依次将不同的配点代入原系统中，生成相应的系统响应函数，最后建立如下的方程，本发明由于利用最高阶次为2阶的Hermite多项式，故配点应当取为的根，即对于n维随机变量来说，采样点的数目，即随机变量取值的组合数为3ⁿ。

$\left[\begin{array}{cccc}{\mathrm{\Γ}}_0\left({\mathrm{\ξ}}_0\right)& {\mathrm{\Γ}}_1\left({\mathrm{\ξ}}_0\right)& \mathrm{...}& {\mathrm{\Γ}}_{s-1}\left({\mathrm{\ξ}}_0\right)\\ {\mathrm{\Γ}}_0\left({\mathrm{\ξ}}_1\right)& {\mathrm{\Γ}}_1\left({\mathrm{\ξ}}_1\right)& \mathrm{...}& {\mathrm{\Γ}}_{s-1}\left({\mathrm{\ξ}}_1\right)\\ .& .& & .\\ .& .& \mathrm{...}& .\\ .& .& & .\\ {\mathrm{\Γ}}_0\left({\mathrm{\ξ}}_N\right)& {\mathrm{\Γ}}_1\left({\mathrm{\ξ}}_N\right)& \mathrm{...}& {\mathrm{\Γ}}_{s-1}\left({\mathrm{\ξ}}_N\right)\end{array}\right]\left[\begin{array}{c}{c}_0\\ c_1\\ .\\ .\\ .\\ c_{s-1}\end{array}\right]=\left[\begin{array}{c}f\left({\mathrm{\ξ}}_0\right)\\ f\left({\mathrm{\ξ}}_1\right)\\ .\\ .\\ .\\ f\left({\mathrm{\ξ}}_N\right)\end{array}\right]$

ξ₀,ξ₁,…,ξ_N为采样点，N为采样点数目，s为待定系数的数目，利用最小二乘法对上式进行回归分析，就能计算出混沌多项式的展开系数c_0,2,c_i,2,c_ii,2,c_ij,2，从而建立关于密频特征值一阶摄动量的代理模型；

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

步骤(41)、结合矩阵摄动方法和基于多项式混沌展开的代理模型技术，建立改进的密集频率结构自由振动特征值的近似计算方法：

${\mathrm{\λ}}^{\mathrm{\η}}={\mathrm{\λ}}_0^{\mathrm{\η}}+{\mathrm{\ϵ\λ}}_1^{\mathrm{\η}}={\mathrm{\λ}}_0^{\mathrm{\η}}+\mathrm{\ϵ}(c_{0,\mathrm{\η}}+\underset{i=1}{\overset{n}{\Σ}}{c}_{i,\mathrm{\η}}{\mathrm{\ξ}}_i+\underset{i=1}{\overset{n}{\Σ}}{c}_{ii,\mathrm{\η}}({\mathrm{\ξ}}_i^2-1)+\underset{i=1}{\overset{n-1}{\Σ}}\underset{j>i}{\overset{n}{\Σ}}{c}_{ij,\mathrm{\η}}{\mathrm{\ξ}}_i{\mathrm{\ξ}}_j)$

原密频结构特征值序列在受到扰动之后，得到一个新的按大小升序排列的特征值序列，其中上标η表示排序为η的特征值；

步骤(42)、根据步骤(41)给出的密集频率特征值的近似计算方法表达式，推导出密频特征值在结构参数发生扰动和变化后的均值和方差。

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

(1)本发明解决了经典摄动方法中，密频结构特征值的一阶摄动量无法被结构参数显式表达，因而无法直接研究其统计特征的困难；

(2)本发明方法只需少数的配点计算就能建立比较精确的近似模型，避免了蒙特卡洛方法大规模的极其耗时的样本点计算的问题；

(3)本发明避免了基于泰勒级数展式的摄动方法在处理密集频率结构的自由振动特征值问题时，只能处理小范围内变化的结构参数所造成的应用限制，因此本发明的应用范围更加广泛。

附图说明

图1为本发明一种面向密集频率结构振动特征值的随机摄动方法的实现流程图；

图2为本发明一种面向密集频率结构振动特征值的随机摄动方法的实施例示意图，其中图2(a)为某型高空飞艇结构的几何模型，图2(b)为某型高空飞艇结构的有限元模型；

图3为本发明一种面向密集频率结构振动特征值的随机摄动方法的实施例示意图，其中图3(a)为某型高空飞艇结构第4阶模态振型示意图，图3(b)某型高空飞艇结构第5阶模态振型示意图。

具体实施方式

下面结合附图以及具体实施方式进一步说明本发明。

本发明提出了一种面向密集频率结构振动特征值的随机摄动方法，其具体实施步骤是：

第一步：根据结构振动的特征值理论，给出刚度矩阵和质量矩阵的谱分解，下面给出具体的过程：

(1)密频结构系统中质量矩阵M₀在刚度矩阵K₀上的谱分解表示：

考虑结构振动特征值问题：

由上式可以得到：

同时又因为对此式两边取逆，得：

由此式得到：

综合以上推导，得到：

这就是M₀在K₀上的谱分解表示。

(2)密频结构系统中刚度矩阵K₀在质量矩阵M₀上的谱分解表示：

考虑广义的特征值问题：

现在我们构造一个与矩阵K₀十分接近的矩阵使得满足：

λ₀和U₀是上式所示重频特征值问题的m重特征值和相应的特征向量，Λ_A和U_A也是它们的特征值对角阵和特征向量矩阵。

由于和K₀十分接近，因此可以将K₀表达为：

$K_0={\stackrel{\‾}{K}}_0+{\mathrm{\ϵ\δK}}_0$

以下推导矩阵和εδK₀的表达式。由方程 可得：

在上式两侧右乘得到：

上式即称为矩阵K₀在M₀上的谱分解表达式。

第二步：利用移位方法，将原密频特征值系统转化为重频特征值系统，修改后系统的特征值问题化为重特征值系统基础上，经过两部分摄动叠加后的特征值问题，其中一部分是系统参数扰动引起的摄动，另一部分是特征值移位引起的摄动。再根据重频结构摄动理论，针对修改后的系统，推导了该系统自由振动特征值的一阶摄动量的矩阵方程。具体实施步骤如下：

(1)不妨令

${\mathrm{\λ}}_0=\left(\underset{i=1}{\overset{m}{\Σ}}{\mathrm{\λ}}_0^i\right)/m,$

那么Λ₀可以表达为：

Λ₀＝λ₀I+εδΛ₀,

其中：

${\mathrm{\ϵ\δ\Λ}}_0={\mathrm{\Λ}}_0-{\mathrm{\λ}}_{0}I={\mathrm{\Λ}}_0-\left(\underset{i=1}{\overset{m}{\Σ}}{\mathrm{\λ}}_0^i\right)/mI,$

(2)将Λ₀＝λ₀I+εδΛ₀式代入中，得到：

$K_0=M_0\left(U_0\right({\mathrm{\λ}}_{0}I+{\mathrm{\ϵ\δ\Λ}}_0\left){U_0}^T\right)M_0+M_0\left(U_A{\mathrm{\Λ}}_A{U_A}^T\right)M_0={\stackrel{\‾}{K}}_0+{\mathrm{\ϵ\δK}}_0,$

其中：

${\stackrel{\‾}{K}}_0=M_0\left({\mathrm{\λ}}_0{U}_0{U_0}^T\right)M_0+M_0\left(U_A{\mathrm{\Λ}}_A{U_A}^T\right)M_0,$

${\mathrm{\ϵ\δK}}_0=M_0\left(U_0\right({\mathrm{\ϵ\δ\Λ}}_0\left){U_0}^T\right)M_0.$

上式确定的满足下式：

密集特征值组Λ₀越密集，Λ₀就越趋近于λ₀I，那么εδΛ₀就越趋近于0，相应的就越与K₀接近。

当矩阵K₀和M₀分别有扰动εK₁和εM₁之后，相应的密集特征值组及其特征向量的摄动问题可改写为：

(K₀+εK₁)U＝(M₀+εM₁)UΛ,

U(M₀+εM₁)U^T＝I.

将代入(K₀+εK₁)U＝(M₀+εM₁)UΛ中，可得：

$({\stackrel{\‾}{K}}_0+\mathrm{\ϵ}{\stackrel{\‾}{K}}_1)U=(M_0+{\mathrm{\ϵM}}_1)U\mathrm{\Λ},$

U(M₀+εM₁)U^T＝I.

其中：

$\mathrm{\ϵ}{\stackrel{\‾}{K}}_1={\mathrm{\ϵ\δK}}_0+{\mathrm{\ϵK}}_1,$

实际上，经过移位之后，密频结构特征值问题在形式上已经跟重频结构特征值问题一致，因此重频结构特征值的摄动理论可以应用在密频特征值上，亦即：

$V={U_0}^T({\mathrm{\δK}}_0+K_1-{\mathrm{\λ}}_0{M}_1)U_0,$

求解如下的特征值矩阵方程：

$\{\begin{array}{c}V\mathrm{\α}={\mathrm{\α\Λ}}_1\\ {\mathrm{\α}}^T\mathrm{\α}=I\end{array}\⇔\{\begin{array}{c}{U_0}^T({\mathrm{\δK}}_0+K_1-{\mathrm{\λ}}_0{M}_1)U_0\mathrm{\α}={\mathrm{\α\Λ}}_1\\ {\mathrm{\α}}^T\mathrm{\α}=I\end{array},$

即得到密频特征值的一阶摄动项Λ₁，α表示相应特征值的系数矩阵。

其中q表示密频特征值的序号(按特征值大小升序排列)，基于以上推导我们获得密频特征值的摄动表达式：

Λ＝λ₀I+Λ₁.

第三步：针对第二步建立的矩阵方程中，方程的根即特征值一阶摄动量无法被结构参数所显式表示的困难，引入多项式混沌展开技术，将上述密频特征值的一阶摄动量表示成为关于系统中结构参数的混沌多项式展式的形式，建立关于密频特征值一阶摄动量的代理模型。具体实施步骤如下：

(1)针对高斯随机场，我们使用齐次Hermite多项式展开，随机响应Y(ξ)可以表示为：

$Y\left(\mathrm{\ξ}\right)=c_0+\underset{i_1=1}{\overset{n}{\Σ}}{c}_{i_1}{H}_1\left({\mathrm{\ξ}}_{i_1}\right)+\underset{i_1=1}{\overset{n}{\Σ}}\underset{i_2=1}{\overset{i_1}{\Σ}}{c}_{i_1{i}_2}{H}_2({\mathrm{\ξ}}_{i_1},{\mathrm{\ξ}}_{i_2})+\underset{i_1=1}{\overset{n}{\Σ}}\underset{i_2=1}{\overset{i_1}{\Σ}}\underset{i_3=1}{\overset{i_2}{\Σ}}{c}_{i_1{i}_2{i}_3}{H}_3({\mathrm{\ξ}}_{i_1},{\mathrm{\ξ}}_{i_2},{\mathrm{\ξ}}_{i_3})+\mathrm{...}$

式中是待定的系数矢量，ξ＝(ξ₁,ξ₂,…,ξ_n)为服从高斯分布的随机变量，表示n次多维Hermite多项式。上式可以被截断，并且用有限项来近似展开，不妨取s项，则上式可以简写为其中为要求解的系数，为j阶广义Wiener-Askey混沌多项式。这里H_n(ξ₁,ξ₂,…,ξ_n)和Γ_j(ξ)一一对应，和一一对应。

(2)我们用二阶Hermite多项式展开构建代理模型，表示随机响应Y(ξ)，则有：

$Y_2\left(\mathrm{\ξ}\right)=c_{0,2}+\underset{i=1}{\overset{n}{\Σ}}{c}_{i,2}{\mathrm{\ξ}}_i+\underset{i=1}{\overset{n}{\Σ}}{c}_{ii,2}({\mathrm{\ξ}}_i^2-1)+\underset{i=1}{\overset{n-1}{\Σ}}\underset{j>i}{\overset{n}{\Σ}}{c}_{ij,2}{\mathrm{\ξ}}_i{\mathrm{\ξ}}_j$

上式中n为随机变量的维数，c_0,2,c_i,2,c_ii,2,c_ij,2表示最高阶次为二阶的Hermite多项式展开式的待定系数。根据上式可以归纳出二阶Hermite随机响应多项式的展开式中待定系数的数目为其中p表示多项式展开式的最高阶次，当p＝2时，s＝(n+2)(n+1)/2。

(3)进行配点的选取，从而计算得出多项式混沌展式的各待定系数

考虑到本文构造的PCE展开式的最高阶次为p＝2，故配点通常取为(p+1)阶Hermite多项式的根，即我们不妨认为系统中的随机变量ξ_i(i＝1,2,…,n)均满足高斯分布那么每一个ξ_i配点的取值有三个对于n维随机变量来说，配点的组合数目为3ⁿ。并且每一个配点组可以表示为：

ξ_β＝(ξ_β1,ξ_β2,…,ξ_βn),β＝1,2,…,3ⁿ

将每一个选定配点组合依次代入中，并将M₁,K₁,δK₀表示成为：

$\left\{\begin{array}{c}{M}_1=M_{\mathrm{\β}1}=\underset{i=1}{\overset{n}{\Σ}}\left(\frac{\∂M}{\∂{\mathrm{\ξ}}_{\mathrm{\β}i}}{\mathrm{\Δ\ξ}}_{\mathrm{\β}i}\right)=\underset{i=1}{\overset{n}{\Σ}}(\frac{M({\mathrm{\μ}}_1,{\mathrm{\μ}}_2,\mathrm{...},{\mathrm{\μ}}_{i-1},{\mathrm{\ξ}}_{\mathrm{\β}i},{\mathrm{\μ}}_{i+1},\mathrm{...},{\mathrm{\μ}}_n)-M({\mathrm{\μ}}_1,{\mathrm{\μ}}_2,\mathrm{...},{\mathrm{\μ}}_{i-1},{\mathrm{\μ}}_i,{\mathrm{\μ}}_{i+1},\mathrm{...},{\mathrm{\μ}}_n)}{({\mathrm{\ξ}}_{\mathrm{\β}i}-{\mathrm{\μ}}_i)}\×2{\mathrm{\σ}}_i)\\ K_1=K_{\mathrm{\β}1}=\underset{i=1}{\overset{n}{\Σ}}\left(\frac{\∂K}{\∂{\mathrm{\ξ}}_{\mathrm{\β}i}}{\mathrm{\Δ\ξ}}_{\mathrm{\β}i}\right)=\underset{i=1}{\overset{n}{\Σ}}(\frac{K({\mathrm{\μ}}_1,{\mathrm{\μ}}_2,\mathrm{...},{\mathrm{\μ}}_{i-1},{\mathrm{\ξ}}_{\mathrm{\β}i},{\mathrm{\μ}}_{i+1},\mathrm{...},{\mathrm{\μ}}_n)-K({\mathrm{\μ}}_1,{\mathrm{\μ}}_2,\mathrm{...},{\mathrm{\μ}}_{i-1},{\mathrm{\μ}}_i,{\mathrm{\μ}}_{i+1},\mathrm{...},{\mathrm{\μ}}_n)}{({\mathrm{\ξ}}_{\mathrm{\β}i}-{\mathrm{\μ}}_i)}\×2{\mathrm{\σ}}_i)\\ {\mathrm{\δK}}_0={\mathrm{\δK}}_{\mathrm{\β}0}=M_0\left(U_0\right({\mathrm{\Λ}}_0-\left(\underset{i=1}{\overset{m}{\Σ}}{\mathrm{\λ}}_0^i/m\right)I\left){U_0}^T\right)M_0\\ =M({\mathrm{\μ}}_1,{\mathrm{\μ}}_2,\mathrm{...},{\mathrm{\μ}}_{i-1},{\mathrm{\μ}}_i,{\mathrm{\μ}}_{i+1},\mathrm{...},{\mathrm{\μ}}_n)\left(U_0\right({\mathrm{\Λ}}_0-\left(\underset{i=1}{\overset{m}{\Σ}}{\mathrm{\λ}}_0^i/m\right)I\left){U_0}^T\right)M({\mathrm{\μ}}_1,{\mathrm{\μ}}_2,\mathrm{...},{\mathrm{\μ}}_{i-1},{\mathrm{\μ}}_i,{\mathrm{\μ}}_{i+1},\mathrm{...},{\mathrm{\μ}}_n)\end{array}\right.,$

然后求解方程式如下：

$\{\begin{array}{c}{U_0}^T({\mathrm{\δK}}_{\mathrm{\β}0}+K_{\mathrm{\β}1}-{\mathrm{\λ}}_0{M}_{\mathrm{\β}1})U_0\mathrm{\α}={\mathrm{\α\Λ}}_1\\ {\mathrm{\α}}^T\mathrm{\α}=I\end{array},$

把所得的值从小到大排列，得到一组解重复以上步骤，将每一个配点顺次代入方程则得到解向量：

$\left[\begin{array}{ccccc}{\mathrm{\λ}}_{11}^q& {\mathrm{\λ}}_{21}^q& {\mathrm{\λ}}_{31}^q& \mathrm{...}& {\mathrm{\λ}}_{3^{n}1}^q\end{array}\right]=\left[\begin{array}{ccccc}{\mathrm{\λ}}_{11}^1& {\mathrm{\λ}}_{21}^1& {\mathrm{\λ}}_{31}^1& \mathrm{...}& {\mathrm{\λ}}_{3^{n}1}^1\\ {\mathrm{\λ}}_{11}^2& {\mathrm{\λ}}_{21}^2& {\mathrm{\λ}}_{31}^2& \mathrm{...}& {\mathrm{\λ}}_{3^{n}1}^2\\ {\mathrm{\λ}}_{11}^3& {\mathrm{\λ}}_{21}^3& {\mathrm{\λ}}_{31}^3& \mathrm{...}& {\mathrm{\λ}}_{3^{n}1}^3\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ {\mathrm{\λ}}_{11}^{\mathrm{\γ}}& {\mathrm{\λ}}_{21}^{\mathrm{\γ}}& {\mathrm{\λ}}_{31}^{\mathrm{\γ}}& \mathrm{...}& {\mathrm{\λ}}_{3^{n}1}^{\mathrm{\γ}}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ {\mathrm{\λ}}_{11}^m& {\mathrm{\λ}}_{21}^m& {\mathrm{\λ}}_{31}^m& \mathrm{...}& {\mathrm{\λ}}_{3^{n}1}^m\end{array}\right]$

考察上式中任意阶的解可以建立以下方程组：

$\left[\begin{array}{cccc}{\mathrm{\Γ}}_0\left({\mathrm{\ξ}}_0\right)& {\mathrm{\Γ}}_1\left({\mathrm{\ξ}}_0\right)& \mathrm{...}& {\mathrm{\Γ}}_{s-1}\left({\mathrm{\ξ}}_0\right)\\ {\mathrm{\Γ}}_0\left({\mathrm{\ξ}}_1\right)& {\mathrm{\Γ}}_1\left({\mathrm{\ξ}}_1\right)& \mathrm{...}& {\mathrm{\Γ}}_{s-1}\left({\mathrm{\ξ}}_1\right)\\ .& .& & .\\ .& .& \mathrm{...}& .\\ .& .& & .\\ {\mathrm{\Γ}}_0\left({\mathrm{\ξ}}_N\right)& {\mathrm{\Γ}}_1\left({\mathrm{\ξ}}_N\right)& \mathrm{...}& {\mathrm{\Γ}}_{s-1}\left({\mathrm{\ξ}}_N\right)\end{array}\right]\left[\begin{array}{c}{c}_0^{\mathrm{\η}}\\ c_1^{\mathrm{\η}}\\ .\\ .\\ .\\ c_{s-1}^{\mathrm{\η}}\end{array}\right]=\left[\begin{array}{c}{\mathrm{\λ}}_{11}^{\mathrm{\η}}\\ {\mathrm{\λ}}_{21}^{\mathrm{\η}}\\ .\\ .\\ .\\ {\mathrm{\λ}}_{N1}^{\mathrm{\η}}\end{array}\right]$

上式由3ⁿ个方程组成，要求解s个待定系数。对于PCE方法，通常来说3ⁿ>s，所以可以基于最小二乘法进行回归分析，求出一组从而建立关于密频特征值一阶摄动量的代理模型。

第四步：结合摄动方法和基于多项式混沌展开的代理模型技术，提出了密集频率结构自由振动特征值的近似计算方法，基于该近似计算方法，进一步得到了密频结构特征值在参数扰动情况下的均值和方差的表达式。具体实施步骤如下：

(1)建立面向密集频率结构自由振动特征值的近似计算方法

经过以上的求解过程，对于密频结构特征值的一阶摄动，其次序为η的值可以表示为：

${\mathrm{\λ}}_1^{\mathrm{\η}}\left(\mathrm{\ξ}\right)=c_{0,\mathrm{\η}}+\underset{i=1}{\overset{n}{\Σ}}{c}_{i,\mathrm{\η}}{\mathrm{\ξ}}_i+\underset{i=1}{\overset{n}{\Σ}}{c}_{ii,\mathrm{\η}}({\mathrm{\ξ}}_i^2-1)+\underset{i=1}{\overset{n-1}{\Σ}}\underset{j>1}{\overset{n}{\Σ}}{c}_{ij,\mathrm{\η}}{\mathrm{\ξ}}_i{\mathrm{\ξ}}_j,$

其中ξ＝(ξ₁,ξ₂,…,ξ_n)为服从高斯分布的随机变量，经过第二步的求解过程，c_0,η,c_i,η,c_ii,η,c_ij,η成为已知量。

这样我们只考虑到一阶摄动项，可以得到：

${\mathrm{\λ}}^{\mathrm{\η}}={\mathrm{\λ}}_0+{\mathrm{\ϵ\λ}}_1^{\mathrm{\η}}={\mathrm{\λ}}_0+\mathrm{\ϵ}(c_{0,\mathrm{\η}}+\underset{i=1}{\overset{n}{\Σ}}{c}_{i,\mathrm{\η}}{\mathrm{\ξ}}_i+\underset{i=1}{\overset{n}{\Σ}}{c}_{ii,\mathrm{\η}}({\mathrm{\ξ}}_i^2-1)+\underset{i=1}{\overset{n-1}{\Σ}}\underset{j>i}{\overset{n}{\Σ}}{c}_{ij,\mathrm{\η}}{\mathrm{\ξ}}_i{\mathrm{\ξ}}_j),$

为了简便起见，在下文中，将上式的上下脚标η去掉，即：

$\mathrm{\λ}={\mathrm{\λ}}_0+{\mathrm{\ϵ\λ}}_1={\mathrm{\λ}}_0+\mathrm{\ϵ}(c_0+\underset{i=1}{\overset{n}{\Σ}}{c}_i{\mathrm{\ξ}}_i+\underset{i=1}{\overset{n}{\Σ}}{c}_{ii}({\mathrm{\ξ}}_i^2-1)+\underset{i=1}{\overset{n-1}{\Σ}}\underset{j>1}{\overset{n}{\Σ}}{c}_{ij}{\mathrm{\ξ}}_i{\mathrm{\ξ}}_j),$

(2)构建密频特征值在结构参数发生扰动和变化后的均值和方差表达式：

对上式求期望得到：

$\begin{array}{c}E\[\mathrm{\λ}\]=E\[{\mathrm{\λ}}_0+{\mathrm{\ϵ\λ}}_1\]=E\[{\mathrm{\λ}}_0+\mathrm{\ϵ}(c_0+\underset{i=1}{\overset{n}{\Σ}}{c}_i{\mathrm{\ξ}}_i+\underset{i=1}{\overset{n}{\Σ}}{c}_{ii}({\mathrm{\ξ}}_i^2-1)+\underset{i=1}{\overset{n-1}{\Σ}}\underset{j>i}{\overset{n}{\Σ}}{c}_{ij}{\mathrm{\ξ}}_i{\mathrm{\ξ}}_j)\]\\ ={\mathrm{\λ}}_0+\mathrm{\ϵ}E\[c_0+\underset{i=1}{\overset{n}{\Σ}}{c}_i{\mathrm{\ξ}}_i+\underset{i=1}{\overset{n}{\Σ}}{c}_{ii}\left({\mathrm{\ξ}}_i^2-1\right)+\underset{i=1}{\overset{n-1}{\Σ}}\underset{j>i}{\overset{n}{\Σ}}{c}_{ij}{\mathrm{\ξ}}_i{\mathrm{\ξ}}_j\]\\ ={\mathrm{\λ}}_0+{\mathrm{\ϵc}}_0\end{array}$

同样的，使用方差算子，并且考虑Hermite多项式的正交性，我们有：

$\begin{array}{c}Var\&lsqb;\mathrm{\&lambda;}\&rsqb;=Var\&lsqb;{\mathrm{\&lambda;}}_0+{\mathrm{\&epsiv;\&lambda;}}_1\&rsqb;=Var\&lsqb;{\mathrm{\&lambda;}}_0+\mathrm{\&epsiv;}(c_0+\underset{i=1}{\overset{n}{\&Sigma;}}{c}_i{\mathrm{\&xi;}}_i+\underset{i=1}{\overset{n}{\&Sigma;}}{c}_{ii}({\mathrm{\&xi;}}_i^2-1)+\underset{i=1}{\overset{n-1}{\&Sigma;}}\underset{j<i}{\overset{n}{\&Sigma;}}{c}_{ij}{\mathrm{\&xi;}}_i{\mathrm{\&xi;}}_j)\&rsqb;\\ ={\mathrm{\&epsiv;}}^{2}Var\&lsqb;(c_0+\underset{i=1}{\overset{n}{\&Sigma;}}{c}_i{\mathrm{\&xi;}}_i+\underset{i=1}{\overset{n}{\&Sigma;}}{c}_{ii}({\mathrm{\&xi;}}_i^2-1)+\underset{i=1}{\overset{n-1}{\&Sigma;}}\underset{j<i}{\overset{n}{\&Sigma;}}{c}_{ij}{\mathrm{\&xi;}}_i{\mathrm{\&xi;}}_j)\&rsqb;\\ ={\mathrm{\&epsiv;}}^2\underset{i=1}{\overset{n}{\&Sigma;}}{c}_i^2+{\mathrm{\&epsiv;}}^{2}Var\&lsqb;\underset{i=1}{\overset{n}{\&Sigma;}}{c}_{ii}({\mathrm{\&xi;}}_i^2-1)\&rsqb;+{\mathrm{\&epsiv;}}^{2}Var\&lsqb;\underset{i=1}{\overset{n-1}{\&Sigma;}}\underset{j<i}{\overset{n}{\&Sigma;}}{c}_{ij}{\mathrm{\&xi;}}_i{\mathrm{\&xi;}}_j\&rsqb;\\ ={\mathrm{\&epsiv;}}^2\underset{i=1}{\overset{n}{\&Sigma;}}{c}_i^2+{\mathrm{\&epsiv;}}^2(\underset{i=1}{\overset{n}{\&Sigma;}}{c}_{ii}^2<{\mathrm{\&Gamma;}}_{ii}^2>)+{\mathrm{\&epsiv;}}^2(\underset{i=1}{\overset{n-1}{\&Sigma;}}\underset{j<i}{\overset{n}{\&Sigma;}}{c}_{ij}^2<{\mathrm{\&Gamma;}}_{ij}^2>)\end{array}$

综上所述，本发明首先采用特征值移位方法，然后结合摄动方法和基于多项式混沌展开的代理模型技术，提出了面向密集频率结构自由振动特征值的近似计算方法，并基于该近似计算方法，进一步建立了密频结构自由振动特征值在参数扰动情况下的均值和方差的表达式。本发明解决了经典摄动方法中，密频结构特征值的一阶摄动量无法被显式表达，因而无法直接研究其统计学特征的困难；避免了基于泰勒级数展式的摄动方法处理密集频率结构的振动特征值问题时，只能处理小范围内变化的结构参数所造成的应用限制；同时本发明只需少数的配点计算就能建立较精确的特征值近似模型，避免了极其耗时的的样本点计算问题。

实施例：

为了更充分的了解该发明的特点及其对工程实际的适用性，本发明以图2的高空飞艇结构系统为例进行密集频率结构的特征值分析验证。图2中的高空飞艇结构长度为16.5m，长细比3.9。采用2D膜单元，在MSC.Patran/Nastran平台上进行带预应力的模态分析。材料属性上，我们选用各向同性材料，弹性模量、材料密度和泊松比的名义值分别取为E＝1096MPa，ρ＝727kg/m³，μ＝0.3。同时将弹性模量和材料密度都处理为随机变量，并且分别满足高斯分布N(1096,54.8²)，N(727,36.35²)。

首先执行带预应力的非线性模态分析，飞艇的内部压强取为500Pa。使用Patran的SOL106求解器，求解前十阶的固有频率。此时得到一对密集分布的特征值，如图3所示。第4阶固有频率(图3(a))和第5阶固有频率(图3(b))的振动频率值分别为8.6027Hz和8.6035Hz，分布非常密集，其具体的特征值和固有频率如表1所示。

表1

然后我们改写由Patran生成的bdf文件，将改写之后的bdf文件提交Nastran进行计算，从而得到含有总体刚度矩阵和总体质量矩阵的f06文件，并提取出f06文件中的总体刚度阵和总体质量阵。进一步，我们将Nastran和Matlab集成在Isight平台上，并执行本发明所提出的方法。

应用本发明所提出的方法，得出此密集频率结构的第4阶和第5阶自由振动特征值的均值和方差如表2所示。为了验证本发明所提出的方法，同样采用Monte-Carlo方法计算了本实施例中的密频结构特征值的均值和方差。在随机数样本取值为10⁵时，由Monte-Carlo方法计算得到的结果和由本发明方法计算得到的结果对比如表2所示。

表2

针对如本实施例所示的高空飞艇结构，和经典的蒙特卡洛方法相比，本发明方法的计算精度是能够令人满意的。不仅如此，本发明所提出的方法，在计算时只需要在原系统中对各随机变量进行有限的几次配点计算，无需如统计方法那样对所有样本点进行重复计算。因此在计算效率方面也有明显的优势。

由实施例可以看出，本发明所提出的方法对于处理大型复杂空间对称或准对称结构的密频特征值问题具有计算精度高、计算耗时少的显著优势。能够有效处理密集频率结构的动态分析，模型更新，不确定性分析和设计优化等问题。经过二次开发，具有成为成熟的商业软件的潜力。以上实施例验证了本方法面向密集频率结构的振动特征值摄动分析的可行性和优越性。

以上仅是本发明的具体步骤，对本发明的保护范围不构成任何限制。

本发明未详细阐述部分属于本领域技术人员的公知技术。

Claims (5)

1.一种面向密集频率结构振动特征值的随机摄动方法，其特征在于：该方法适用于具有密集分布的自由振动频率的结构系统，并且考虑其结构参数发生变化和扰动的情况，包括以下步骤：

第一步：根据结构振动的特征值理论，给出刚度矩阵和质量矩阵的谱分解；

第二步：利用移位方法，将原密频特征值系统转化为重频特征值系统，修改后系统的特征值问题化为重特征值系统基础上，经过两部分摄动叠加后的特征值问题，其中一部分是系统参数扰动引起的摄动，另一部分是特征值移位引起的摄动；再根据重频结构摄动理论，针对修改后的系统，推导了该系统自由振动特征值的一阶摄动量的矩阵方程；

第三步：针对第二步构建的矩阵方程，其方程的根亦即特征值一阶摄动量无法被结构参数所显式表达的困难，引入多项式混沌展开技术，将上述密频特征值的一阶摄动量表示成为关于系统中结构参数的混沌多项式展式的形式，建立关于密频特征值一阶摄动量的代理模型；

第四步：结合摄动方法和基于多项式混沌展开的代理模型技术，提出了面向密集频率结构自由振动特征值的近似计算方法，基于该近似计算方法，进一步得到了密集频率结构自由振动特征值在参数扰动情况下的均值和方差的表达式。

2.根据权利要求1所述的一种面向密集频率结构振动特征值的随机摄动方法，其特征在于：所述第一步具体实现如下：

步骤(11)、给出密频结构系统中质量矩阵M₀在刚度矩阵K₀上的谱分解表示其中K₀和M₀是n×n维实对称矩阵，Λ₀和U₀是它的一组密集特征值对角阵及其特征向量矩阵，Λ₀和U₀的维数分别为m×m和n×m，Λ_A和U_A也是原系统的特征值对角阵和特征向量矩阵；

步骤(12)、给出密频结构系统中刚度矩阵K₀在质量矩阵M₀上的谱分解表示，

3.根据权利要求1所述的一种面向密集频率结构振动特征值的随机摄动方法，其特征在于：所述第二步具体实现如下：

步骤(21)、给出移位后的重频特征值表达式其中为原系统各个密集频率值，m为原系统中一组密集固有频率的个数；

步骤(22)、推导了修改后系统的自由振动特征值一阶摄动量的矩阵方程其中α表示相应的系数矩阵，Λ₁表示密频特征值一阶摄动项的对角矩阵，δK₀为特征值移位引起的摄动，K₁和M₁分别表示结构刚度矩阵和结构质量矩阵的一阶摄动项。

4.根据权利要求1所述的一种面向密集频率结构振动特征值的随机摄动方法，其特征在于：所述第三步具体实现如下：

步骤(31)、针对结构参数ξ₁,ξ₂,…,ξ_n满足高斯分布的情况，使用Hermite多项式展式，将表示为式中表示待定的系数矢量，ξ＝(ξ₁,ξ₂,…,ξ_n)为服从高斯分布的随机变量，表示n次多维Hermite多项式；

步骤(32)、将步骤(31)中的Y(ξ)用有限项截断并用二阶Hermite多项式近似展开，得到上式中n为随机变量的维数，c_0,2,c_i,2,c_ii,2,c_ij,2表示最高阶次为二阶的Hermite多项式展开式的待定系数；

步骤(33)、选择好配点组合，依次将不同的配点代入原系统中，生成相应的系统响应函数，最后建立如下的方程，由于利用最高阶次为2阶的Hermite多项式，故配点应当取为的根，即对于n维随机变量来说，采样点的数目，即随机变量取值的组合数为3ⁿ；

ξ₀,ξ₁,…,ξ_N为采样点，N为采样点数目，s为待定系数的数目，利用最小二乘法对上式进行回归分析，就能计算出混沌多项式的展开系数c_0,2,c_i,2,c_ii,2,c_ij,2，从而建立关于密频特征值一阶摄动量的代理模型。

5.根据权利要求1所述的一种面向密集频率结构振动特征值的随机摄动方法，其特征在于：所述第四步具体实现如下：

步骤(41)、结合矩阵摄动方法和基于多项式混沌展开的代理模型技术，建立改进的密集频率结构自由振动特征值的近似计算方法：

${\mathrm{\&lambda;}}^{\mathrm{\&eta;}}={\mathrm{\&lambda;}}_0^{\mathrm{\&eta;}}+{\mathrm{\&epsiv;\&lambda;}}_1^{\mathrm{\&eta;}}={\mathrm{\&lambda;}}_0^{\mathrm{\&eta;}}+\mathrm{\&epsiv;}\left(c_{0,\mathrm{\&eta;}}+\underset{i=1}{\overset{n}{\&Sigma;}}{c}_{i,\mathrm{\&eta;}}{\mathrm{\&xi;}}_i+\underset{i=1}{\overset{n}{\&Sigma;}}{c}_{ii,\mathrm{\&eta;}}\left({\mathrm{\&xi;}}_i^2-1\right)+\underset{i=1}{\overset{n-1}{\&Sigma;}}\underset{j>i}{\overset{n}{\&Sigma;}}{c}_{ij,\mathrm{\&eta;}}{\mathrm{\&xi;}}_i{\mathrm{\&xi;}}_j\right)$

原密频结构特征值序列在受到扰动之后，得到一个新的按大小升序排列的特征值序列，其中上标η表示排序为η的特征值；

步骤(42)、根据步骤(41)给出的密集频率特征值的近似计算方法表达式，推导出密频特征值在结构参数发生扰动和变化后的均值和方差。