CN105912508A - 一种改进的基于代理模型的重频结构振动特征值的随机摄动方法 - Google Patents

Abstract

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

Description

一种改进的基于代理模型的重频结构振动特征值的随机摄动 方法

技术领域

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

背景技术

在工程实际中常会出现重频系统或者密频系统，如平流层飞艇、飞机、大型火箭等。特别是对于空间柔性结构，其自由振动的主要特点之一就是固有频率重复或者密集出现。重频是实际工程结构中固有频率出现重叠的现象，在工程中，若不考虑重频的影响，将重频模态处理为孤立模态，便会引起模态的遗漏。具有重特征值的系统称为退化系统，退化系统有两个重要的特点：一是结构参数变化后，原来的一组重复特征值可能分离为一组非重复特征值；二是结构参数变化后，特征向量可能产生跳跃现象，即相同编号的新旧两个特征向量之差不是一个小量。产生这种现象的原因，主要是在结构具有重复频率的情况下，原系统的特征向量的选取具有一定的随意性。

Kaminski M和Solecka M将结构特征值和特征向量用混沌多项式(PCE)方法进行展开，研究了线性随机系统的受迫振动响应分析。X.W.YANG和S.H.CHEN将Padeapproximation应用到矩阵摄动理论中，求出了特征向量和特征值变化量的表达式。关于结构实模态特征值的统计特性，Qiu Z.P和Qiu H.C提出了直接方差分析方法(DVA方法)，无需已知或假定结构参数的相关系数矩阵，通过矩阵摄动理论和概率理论便能直接计算得到实模态结构的随机特征值的方差。Zhao和Xu等人针对重频阻尼振动系统，提出了基于松弛组合近似方法求解重分析问题的快速算法。通过将特征向量表示成为基向量和系数向量的组合形式，避免了求解大规模方程组，简化了复杂求解运算。Palej和Krowiak利用解析的方法进行了不同自由度下，重频弹簧质量块系统的模态分析。Vessel和Ram等人给出了两种含不确定参数的阻尼振动系统重复特征值导数的计算方法。Wang和Zhang研究了将两个相近频率处理成重复频率带来的误差。Chen和Jiao构造出一个弹簧质量块系统，并且提出能够通过调整系统的参数，产生任意类型的重复频率振动系统。Xia和Shi等人将最大化结构的第一阶独立或者重复固有频率为优化目标，研究了基于水平集方法的结构拓扑优化问题。Reimherr进行了具有重复频率特征值的方程的回归分析。Fan和Xiang等人通过引入子域方法，将随机摄动方法的处理范围扩展到了参数具有较大范围不确定性的领域以及重复特征值的领域。

矩阵摄动方法作为一种能够进行快速结构重分析和结构快速灵敏度分析的实用工具，具有计算效率高、计算成本低和易于实现等特点，已经同时在理论基础和工程应用中受到了广泛的关注并取得了长足的进步。但目前国内外学者对于摄动方法的研究大多集中在基于泰勒级数展式的摄动方法上，当处理重频特征值时，由于特征值的导数不唯一，其分析过程遇到困难。因此，目前关于使用摄动方法处理重复频率结构特征值问题的研究工作还很少。

发明内容

本发明要解决的技术问题是：克服现有技术的不足，提供一种改进的基于代理模型的重频特征值随机摄动方法，用于计算重复频率结构在参数发生扰动和变化后，其特征值的近似表达式，并基于该近似计算方法，进一步得到了重频结构特征值在参数扰动情况下的均值和方差。本发明解决了传统摄动方法中，重频结构特征值的一阶摄动量无法被显式表达，因而无法进一步研究其统计特征的困难。同时本发明方法只需少数的配点计算就能建立比较精确的近似模型，也避免了蒙特卡洛方法大规模高成本的样本点计算问题。因此本发明极大的方便了重复频率结构的动态分析、模型更新和优化设计等工作。

本发明采用的技术方案为：一种改进的基于代理模型的重频结构振动特征值的随机摄动方法，适用于具有重复自由振动频率的结构系统。首先根据矩阵摄动理论，推导了重复频率结构系统在刚度矩阵、质量矩阵等参数发生扰动和变化后，关于其自由振动特征值的一阶摄动量的矩阵方程。然后引入多项式混沌展开技术，建立针对重频结构自由振动特征值一阶摄动量的代理模型。进一步结合矩阵摄动方法和多项式混沌展开技术，给出了重复频率结构自由振动特征值的近似计算方法。基于该近似计算方法，得到了重频结构特征值在参数扰动情况下的均值和方差的表达式，其实现步骤如下：

第一步：根据矩阵摄动理论，针对具有重复自由振动频率的结构系统，推导了在其刚度矩阵、质量矩阵等结构参数发生扰动和变化后，关于该系统自由振动特征值的一阶摄动量的矩阵方程，这一步是建立重复频率结构自由振动特征值近似算法的基础，后续的算法构建都是基于此方法实现的；

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

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

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

步骤(11)、给出退化系统具有m重的特征值和相应m个两两正交的特征向量wⁱ的特征方程其中K₀和M₀分别表示重复频率结构系统的原刚度矩阵和原质量矩阵；

步骤(12)、根据矩阵摄动理论，建立关于重复频率结构自由振动特征值一阶摄动量的矩阵方程其中wⁱ(i＝1,2,…,m)为原重复频率系统的特征向量，α_j(j＝1,2,…,m)是待定系数，K₁和M₁分别是系统刚度矩阵和质量矩阵的一阶摄动量，和分别为重频系统自由振动的原特征值和扰动后的一阶摄动量；

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

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

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

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

$\left[\begin{array}{cccc}{\mathrm{\Γ}}_0\left({\mathrm{\ξ}}_0\right)& {\mathrm{\Γ}}_1\left({\mathrm{\ξ}}_0\right)& ...& {\mathrm{\Γ}}_{s-1}\left({\mathrm{\ξ}}_0\right)\\ {\mathrm{\Γ}}_0\left({\mathrm{\ξ}}_1\right)& {\mathrm{\Γ}}_1\left({\mathrm{\ξ}}_1\right)& ...& {\mathrm{\Γ}}_{s-1}\left({\mathrm{\ξ}}_1\right)\\ .& .& & .\\ .& .& ...& .\\ .& .& & .\\ {\mathrm{\Γ}}_0\left({\mathrm{\ξ}}_N\right)& {\mathrm{\Γ}}_1\left({\mathrm{\ξ}}_N\right)& ...& {\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，从而建立关于重频特征值一阶摄动量的代理模型；

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

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

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

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

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

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

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

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

(3)本发明避免了基于泰勒级数展式的摄动方法处理重复频率结构的自由振动特征值问题时，由于特征值的导数非唯一而造成的分析困难，因此本发明的应用范围更加广泛。

附图说明

图1为本发明一种改进的基于代理模型的重频结构振动特征值的随机摄动方法的实现流程图；

图2为本发明一种改进的基于代理模型的重频结构振动特征值的随机摄动方法的实施例示意图，其中图2(a)为穹顶结构侧视图，图2(b)为穹顶结构俯视图。

具体实施方式

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

本发明提出了一种改进的基于代理模型的重频结构振动特征值的随机摄动方法，其具体实施步骤是：

第一步：根据矩阵摄动理论，针对具有重复自由振动频率的结构系统，推导了在其刚度矩阵、质量矩阵等结构参数发生扰动和变化后，关于该系统自由振动特征值的一阶摄动量的矩阵方程，这一步是建立重复频率结构自由振动特征值近似算法的基础，后续的算法构建都是基于此方法实现的，下面给出具体的过程：

(1)给出重频结构自由振动的特征值方程

设重频结构(退化系统)的某个特征值为m重特征值，相应有m个两两正交的特征向量wⁱ(i＝1,2,…,m)，并且满足方程：

$K_0{w}^i={\mathrm{\λ}}_0^i{M}_0{w}^i,$

不仅如此，wⁱ的线性组合也是该重特征值的特征向量，即：

$\begin{array}{c}{u}_0^i={\mathrm{\α}}_1{w}^1+{\mathrm{\α}}_2{w}^2+...+{\mathrm{\α}}_m{w}^m\\ =\[w^1,w^2,...,w^m\]\·{\[{\mathrm{\α}}_1,{\mathrm{\α}}_2,...,{\mathrm{\α}}_m\]}^T=W\·\mathrm{\α}\end{array}$

其中α是待定的向量。

(2)根据矩阵摄动理论，建立关于重复频率结构自由振动特征值一阶摄动量的矩阵方程

受到使用环境的变化和材料批次的差异等外界扰动时，退化系统的结构参数发生变化，相应的有质量和刚度的变化发生，可以表示为：

$\left\{\begin{array}{c}K=K_0+{\mathrm{\ϵK}}_1\\ M=M_0+{\mathrm{\ϵM}}_1\end{array}\right.$

扰动后的系统不再出现重特征值，即应该有m个不同的特征值其对应的特征向量为uⁱ(i＝1,2,…,m)。类似于孤立的特征值情况，λⁱ和uⁱ仍可以表示成为ε的级数，

${\mathrm{\λ}}^i={\mathrm{\λ}}_0^i+{\mathrm{\ϵ\λ}}_1^i+{\mathrm{\ϵ}}^2{\mathrm{\λ}}_2^i+...$

$u^i=u_0^i+{\mathrm{\ϵu}}_1^i+{\mathrm{\ϵ}}^2{u}_2^i+...$

根据摄动方法，将以上四式同时代入Ku＝λMu，展开之后略去O(ε³)项之后，比较ε的同次幂系数可以得到：

$K_0{u}_1^i+K_1{u}_0^i={\mathrm{\λ}}_0^i{M}_0{u}_1^i+{\mathrm{\λ}}_0^i{M}_1{u}_0^i+{\mathrm{\λ}}_1^i{M}_0{u}_0^i$

将代入上式，得到：

$K_0{u}_1^i+K_1\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_j{w}^j={\mathrm{\λ}}_0^i{M}_0{u}_1^i+{\mathrm{\λ}}_0^i{M}_1\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_j{w}^j+{\mathrm{\λ}}_1^i{M}_0\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_j{w}^j$

用[w^k]^T左乘上式两边，则有：

${\[w^k\]}^T{K}_0{u}_1^i+\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\[w^k\]}^T{K}_1{w}^j{\mathrm{\α}}_j={\mathrm{\λ}}_0^i{\[w^k\]}^T{M}_0{u}_1^i+{\mathrm{\λ}}_0^i\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\[w^k\]}^T{M}_1{w}^j{\mathrm{\α}}_j+{\mathrm{\λ}}_1^i\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\[w^k\]}^T{M}_0{w}^j{\mathrm{\α}}_j$

考察上式，同时根据以及正交关系式[w^k]^TM₀w^j＝δ_kj，可得：

$\underset{j=1}{\overset{m}{\mathrm{\Σ}}}({\[w^k\]}^T{K}_1{w}^j-{\mathrm{\λ}}_0^i{\[w^k\]}^T{M}_1{w}^j){\mathrm{\α}}_j={\mathrm{\λ}}_1^i{\mathrm{\α}}_j,(k=1,2,...,m)$

不妨令a_kj＝[w^k]^TM₁w^j，b_kj＝[w^k]^TK₁w^j，则上式又可以改写成：

$\underset{j=1}{\overset{m}{\mathrm{\Σ}}}(b_{jk}-{\mathrm{\λ}}_0^i{a}_{kj}){\mathrm{\α}}_j={\mathrm{\λ}}_1^i{\mathrm{\α}}_k,k=1,2,...,m$

上式的矩阵形式为其中矩阵D的元素为(k,j＝1,2,…,m)。

上式的矩阵方程是一个m阶的标准特征值方程，求解此特征值问题可以得到相应于重特征值的一阶摄动以及待定系数α^j。仍然按照以前的规定从小到大升序排列，若上式矩阵方程无重特征值，则所求的和α^j是唯一的，简便起见假定上式矩阵方程无重复特征值。求解之，等价于求解：

$\det \[D-{\mathrm{\λ}}_1^j{I}^{m\×m}\]=0$

可以得到将其带回可以得到α_j(j＝1,2,…,m)，求得和α_j(j＝1,2,…,m)之后，根据可以得到与重特征值相应的特征向量

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

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

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

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

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

$Y_2\left(\mathrm{\ξ}\right)=c_{0,2}+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{c}_{i,2}{\mathrm{\ξ}}_i+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{c}_{ii,2}({\mathrm{\ξ}}_i^2-1)+\underset{i=1}{\overset{n-1}{\mathrm{\Σ}}}\underset{j>i}{\overset{n}{\mathrm{\Σ}}}{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＝1,2,…,n，那么每一个ξ_i配点的取值有三个对于n维随机变量来说，配点的组合数目为3ⁿ。并且每一个配点组可以表示为：

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

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

$\left\{\begin{array}{c}{M}_{\mathrm{\β}1}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left(\frac{\∂M}{\∂{\mathrm{\ξ}}_{\mathrm{\β}i}}{\mathrm{\Δ\ξ}}_{\mathrm{\β}i}\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(\frac{M(\sqrt{3}{\mathrm{\σ}}_i+{\mathrm{\ξ}}_{\mathrm{\β}i},{\mathrm{\μ}}_j)-M({\mathrm{\ξ}}_{\mathrm{\β}i},{\mathrm{\μ}}_j)}{\sqrt{3}{\mathrm{\σ}}_i}\×2{\mathrm{\σ}}_i)\\ K_{\mathrm{\β}1}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left(\frac{\∂K}{\∂{\mathrm{\ξ}}_{\mathrm{\β}i}}{\mathrm{\Δ\ξ}}_{\mathrm{\β}i}\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(\frac{K(\sqrt{3}{\mathrm{\σ}}_i+{\mathrm{\ξ}}_{\mathrm{\β}i},{\mathrm{\μ}}_j)-K({\mathrm{\ξ}}_{\mathrm{\β}i},{\mathrm{\μ}}_j)}{\sqrt{3}{\mathrm{\σ}}_i}\×2{\mathrm{\σ}}_i)\end{array}\right.,\begin{array}{c}j=1,2,...,n.j\≠i.\\ \mathrm{\β}=1,2,...,3^n\end{array}$

然后求解方程式如下：

$\underset{q=1}{\overset{m}{\mathrm{\Σ}}}({\[w^k\]}^T{K}_{\mathrm{\β}1}{w}^q-{\mathrm{\λ}}_0^q{\[w^k\]}^T{M}_{\mathrm{\β}1}{w}^q){\mathrm{\α}}_q={\mathrm{\λ}}_{\mathrm{\β}1}^q{\mathrm{\α}}_q,(k=1,2,...,m)$

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

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

考察上式中任意阶的解γ＝1,2,…,m，可以建立以下方程组：

$\left[\begin{array}{cccc}{\mathrm{\Γ}}_0\left({\mathrm{\ξ}}_0\right)& {\mathrm{\Γ}}_1\left({\mathrm{\ξ}}_0\right)& ...& {\mathrm{\Γ}}_{s-1}\left({\mathrm{\ξ}}_0\right)\\ {\mathrm{\Γ}}_0\left({\mathrm{\ξ}}_1\right)& {\mathrm{\Γ}}_1\left({\mathrm{\ξ}}_1\right)& ...& {\mathrm{\Γ}}_{s-1}\left({\mathrm{\ξ}}_1\right)\\ .& .& & .\\ .& .& ...& .\\ .& .& & .\\ {\mathrm{\Γ}}_0\left({\mathrm{\ξ}}_{3^n}\right)& {\mathrm{\Γ}}_1\left({\mathrm{\ξ}}_{3^n}\right)& ...& {\mathrm{\Γ}}_{s-1}\left({\mathrm{\ξ}}_{3^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{\λ}}_{3^{n}1}^{\mathrm{\γ}}\end{array}\right]$

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

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

(1)建立改进的重复频率结构自由振动特征值的近似计算方法

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

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

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

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

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

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

$\mathrm{\λ}={\mathrm{\λ}}_0+{\mathrm{\ϵ\λ}}_1={\mathrm{\λ}}_0+\mathrm{\ϵ}(c_0+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{c}_i{\mathrm{\ξ}}_i+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{c}_{ii}({\mathrm{\ξ}}_i^2-1)+\underset{i=1}{\overset{n-1}{\mathrm{\Σ}}}\underset{j>i}{\overset{n}{\mathrm{\Σ}}}{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}{\mathrm{\Σ}}}{c}_i{\mathrm{\ξ}}_i+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{c}_{ii}({\mathrm{\ξ}}_i^2-1)+\underset{i=1}{\overset{n-1}{\mathrm{\Σ}}}\underset{j>i}{\overset{n}{\mathrm{\Σ}}}{c}_{ij}{\mathrm{\ξ}}_i{\mathrm{\ξ}}_j)\]\\ ={\mathrm{\λ}}_0+\mathrm{\ϵ}E\[c_0+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{c}_i{\mathrm{\ξ}}_i+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{c}_{ii}({\mathrm{\ξ}}_i^2-1)+\underset{i=1}{\overset{n-1}{\mathrm{\Σ}}}\underset{j>i}{\overset{n}{\mathrm{\Σ}}}{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}{\mathrm{\&Sigma;}}}{c}_i{\mathrm{\&xi;}}_i+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{c}_{ii}({\mathrm{\&xi;}}_i^2-1)+\underset{i=1}{\overset{n-1}{\mathrm{\&Sigma;}}}\underset{j>i}{\overset{n}{\mathrm{\&Sigma;}}}{c}_{ij}{\mathrm{\&xi;}}_i{\mathrm{\&xi;}}_j)\&rsqb;\\ ={\mathrm{\&epsiv;}}^{2}Var\&lsqb;(c_0+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{c}_i{\mathrm{\&xi;}}_i+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{c}_{ii}({\mathrm{\&xi;}}_i^2-1)+\underset{i=1}{\overset{n-1}{\mathrm{\&Sigma;}}}\underset{j>i}{\overset{n}{\mathrm{\&Sigma;}}}{c}_{ij}{\mathrm{\&xi;}}_i{\mathrm{\&xi;}}_j)\&rsqb;\\ ={\mathrm{\&epsiv;}}^2\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{c}_i^2+{\mathrm{\&epsiv;}}^{2}Var\&lsqb;\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{c}_{ii}({\mathrm{\&xi;}}_i^2-1)\&rsqb;+{\mathrm{\&epsiv;}}^{2}Var\&lsqb;\underset{i=1}{\overset{n-1}{\mathrm{\&Sigma;}}}\underset{j>i}{\overset{n}{\mathrm{\&Sigma;}}}{c}_{ij}{\mathrm{\&xi;}}_i{\mathrm{\&xi;}}_j\&rsqb;\\ ={\mathrm{\&epsiv;}}^2\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{c}_i^2+{\mathrm{\&epsiv;}}^2(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{c}_{ii}^2<{\mathrm{\&Gamma;}}_{ii}^2>)+{\mathrm{\&epsiv;}}^2(\underset{i=1}{\overset{n-1}{\mathrm{\&Sigma;}}}\underset{j>i}{\overset{n}{\mathrm{\&Sigma;}}}{c}_{ij}^2<{\mathrm{\&Gamma;}}_{ij}^2>)\end{array}$

综上，本发明结合摄动方法和基于多项式混沌展开的代理模型技术，提出了改进的重复频率结构自由振动特征值的近似计算方法，并基于该近似计算方法，进一步建立了重频结构自由振动特征值在参数扰动情况下的均值和方差的表达式。本发明解决了经典摄动方法中，重频结构特征值的一阶摄动量无法被显式表达，因而无法直接研究其统计特征的困难；避免了基于泰勒级数展式的摄动方法处理重复频率结构的振动特征值问题时，由于特征值的导数非唯一而造成的分析困难；同时本发明只需少数的配点计算就能建立较精确的特征值近似模型，避免了大规模高成本的样本点计算问题。

实施例：

为了更充分的了解该发明的特点及其对工程实际的适用性，本发明以图2的穹顶结构系统为例进行重复频率结构的特征值分析验证。图2中的穹顶结构由52根杆和21个节点组成，每一根杆件的截面积均为10cm²。穹顶结构相对于Y轴和X轴对称，其杆件的弹性模量、材料密度和泊松比的名义值分别为E＝210GPa,ρ＝7850kg/m³,μ＝0.3。重要的具有代表性的节点的坐标如表1所示。

表1

在本实施例中，将弹性模量和材料密度处理为随机变量，并且分别满足分布N(210,10.5²)和N(7850,392.5²)。应用本发明提出的方法，得出此重复频率结构的自由振动特征值的均值和方差如表2所示。为了验证本发明所提出的方法，同样采用Monte-Carlo方法计算了本算例中的重频结构特征值的均值和方差。在随机数取值为10⁵时，由Monte-Carlo方法计算得到的结果和本发明得到的结果对比如表2所示。

表2

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

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

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

Claims (4)

1.一种改进的基于代理模型的重频结构振动特征值的随机摄动方法，其特征在于：该方法适用于具有重复的自由振动频率的结构系统，并且考虑其结构参数发生变化和扰动的情况，包括以下步骤：

第一步：根据矩阵摄动理论，针对具有重复自由振动频率的结构系统，推导了在其刚度矩阵、质量矩阵结构参数发生扰动和变化后，关于该系统自由振动特征值的一阶摄动量的矩阵方程；

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

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

2.根据权利要求1所述的一种改进的基于代理模型的重频结构振动特征值的随机摄动方法，其特征在于：所述第一步具体实现如下：

步骤(11)、给出退化系统具有m重的特征值和相应m个两两正交的特征向量wⁱ的特征方程其中K₀和M₀分别表示重复频率结构系统的原刚度矩阵和原质量矩阵；

步骤(12)、根据矩阵摄动理论，建立关于重复频率结构自由振动特征值一阶摄动量的矩阵方程其中wⁱ(i＝1,2,…,m)为原重复频率系统的特征向量，α_j(j＝1,2,…,m)是待定系数，K₁和M₁分别是系统刚度矩阵和质量矩阵的一阶摄动量，和分别为重频系统自由振动的原特征值和扰动后的一阶摄动量。

3.根据权利要求1所述的一种改进的基于代理模型的重频结构振动特征值的随机摄动方法，其特征在于：所述第二步具体实现如下：

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

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

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

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

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

4.根据权利要求1所述的一种改进的基于代理模型的重频结构振动特征值的随机摄动方法，其特征在于：所述第三步具体实现如下：

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

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

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

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