基于结构模态试验的矩阵型动力学模型修正方法

摘要

本发明公布了一种基于结构模态试验的矩阵型动力学模型修正方法。它是利用多次同时改变结构质量分布和刚度后产生的新结构的模态实验结果，运用矩阵计算和代数方程求解，计算出有限元模型的质量阵和刚度阵的修正量。本发明提高了有限元模型的修正精度，简单易行。

说明书

技术领域

本发明涉及一种基于模态试验的矩阵型动力学模型修正方法，尤其是对结构误差位置或单元确定的动力学模型修正。

背景技术

目前，利用结构动力学模态试验结果进行有限元模型修正有两大类方法。一类是参数型修正法，即通过求解灵敏度来修正结构模型的物理参数。另一类是由Berman，J.C.Chen，Friswell等学者提出的对模型的整个质量阵和刚度阵进行修正的矩阵型修正法。而对一个结构系统进行模态试验所得到的试验数据是有限的，往往不足以很好地识别动力学模型。近年来，CHA等人曾提出通过测量原结构动力学系统模态和添加质量后的结构动力学系统模态来修正模型。但附加质量在实际中往往比较困难，因为质量往往是具有一定体积和惯量的，而在模型修正中却忽略了它们；而且附加质量与结构的连接又不可避免引入的刚度，这些都影响了修正精度。另外当测试模态数目较少时，由于识别方程中未知数数目大于方程数，解方程得到的修正量精度也会很差。

发明内容

本发明目的为了克服附加质量的困难、刚度引入以及参数识别方程的有限性对修正精度的影响，提供一种基于模态试验的矩阵型动力学模型修正方法，该修正方法不仅容易改变结构的质量和刚度，消除引入刚度，而且还能有效增加参数识别方程的个数，提高刚度阵和质量阵的修正精度。

本发明为实现上述目的，采用如下技术方案：

本发明基于结构模态试验的矩阵型动力学模型修正方法：

首先对原结构进行模态测试，得到原结构的振型Φ和频率Λ；然后分多次同时改变结构的质量分布和结构刚度，并对产生的新结构进行模态测试，得到第i次改变结构后新结构的振型Φ_aⁱ和频率Λ_aⁱ，其中i为大于等于1的自然数；

测试模态数为N_e，有限元模型自由度与模态测试自由度均为N，则测试得到的原结构的频率Λ和新结构的频率Λ_aⁱ规模均为N_e×N_e，原结构的振型Φ和新结构的振型Φ_aⁱ规模均为N×N_e；原结构有限元模型修正后的质量矩阵M和刚度矩阵K均为未知量，规模为N×N；第i次结构质量矩阵的改变量M_aⁱ和刚度矩阵的改变量K_aⁱ均为已知量，规模为N×N，其中N为大于1的自然数；

此时原结构模态满足，

MΦΛ＝KΦ        (1)

第i次改变结构质量和刚度后，新结构模态满足，

$(M+M_a^i){\Phi }_a^i{\Lambda }_a^i=(K+K_a^i){\Phi }_a^i---\left(2\right)$

将式(1)转置再后乘Φ_aⁱ，得到，

${\mathrm{\Lambda \Phi }}^{T}M{\Phi }_a^i={\Phi }^{T}K{\Phi }_a^i---\left(3\right)$

用Φ^T前乘式(2)得，

${\Phi }^T(M+M_a^i){\Phi }_a^i{\Lambda }_a^i={\Phi }^T(K+K_a^i){\Phi }_a^i---\left(4\right)$

式(3)、式(4)两式相减得，

${\mathrm{\Lambda \Phi }}^{T}M{\Phi }_a^i-{\Phi }^{T}M{\Phi }_a^i{\Lambda }_a^i={\Phi }^T{M}_a^i{\Phi }_a^i{\Lambda }_a^i-{\Phi }^T{K}_a^i{\Phi }_a^i---\left(5\right)$

由于M为原结构质量阵，未知，令：

$P^i={\Phi }^{T}M{\Phi }_a^i---\left(6a\right)$

$Q^i={\Phi }^T{M}_a^i{\Phi }_a^i{\Lambda }_a^i-{\Phi }^T{K}_a^i{\Phi }_a^i---\left(6b\right)$

其中矩阵Pⁱ、Qⁱ的维数均为N_e×N_e；

式(5)化为，${\mathrm{\Lambda P}}^i-P^i{\Lambda }_a^i=Q^i---\left(7\right)$

由于Λ和A_aⁱ都是对角阵，式(7)可展开为N_e×N_e个代数方程，

$({\lambda }_j-{\lambda }_{\mathrm{ak}}^i)P_{\mathrm{jk}}^i=Q_{\mathrm{jk}}^i---\left(8\right)$

j、k为Pⁱ、Qⁱ矩阵中元素的下标，λ_j为原结构的第j个测试特征值，λ_akⁱ为第i次结构质量和刚度改变后新结构的第k个测量特征值；

原结构有限元模型修正前的质量矩阵和刚度矩阵分别为M₀和K₀，为已知量，规模为N×N，原结构有限元模型质量矩阵和刚度矩阵的修正量分别为ΔM和ΔK，为未知量，规模为N×N；则有，

M＝M₀+ΔM    (9a)

K＝K₀+ΔK    (9b)

原结构有限元模型质量矩阵的修正如下：

将式(9a)代入式(6a)后得到，

${\Phi }^T\mathrm{\Delta M}{\Phi }_a^i=P^i-{\Phi }^T{M}_0{\Phi }_a^i---\left(10a\right)$

设$P^i-{\Phi }^T{M}_0{\Phi }_a^i=P^i---\left(10b\right)$

将式(10a)行拉直后再转置得到Aⁱδm＝rⁱ        (11a)

其中$A^i={\Phi }^T\otimes {\left({\Phi }_a^i\right)}^T---\left(11b\right)$

δm＝[ΔM₁₁…ΔM_1N|ΔM₂₁…ΔM_2N|……|ΔM_N1…ΔM_NN]^T    (11c)

$r^i={[R_{11}^i···R_{1N}^i|R_{21}^i···R_{2N}^i|······|R_{N1}^i···R_{\mathrm{NN}}^i]}^T---\left(11d\right)$

其中Aⁱ矩阵维数为N_e²×N²；

分G次改变原结构的质量和刚度所产生的新结构进行模态测试，然后将得到的G个方程(11a)联立求解得：

$\left(\begin{array}{c}{A}^1\\ A^2\\ ·\\ ·\\ ·\\ A^G\end{array}\right)\mathrm{\delta m}=\left(\begin{array}{c}{t}^1\\ t^2\\ ·\\ ·\\ ·\\ t^G\end{array}\right)---\left(12a\right)$

ΔM＝[unvec_N×N(δm)]^T    (12b)

即$\mathrm{\Delta K}={\left(\begin{array}{cccccc}{\mathrm{\delta m}}_1& {\mathrm{\delta m}}_{N+1}& ·& ·& ·& \delta m_{(N-1)\times N+1}\\ ·& ·& & & ·& ·\\ ·& ·& & & ·& ·\\ ·& ·& & & ·& ·\\ {\mathrm{\delta m}}_N& {\mathrm{\delta m}}_{2N}& ·& ·& ·& {\mathrm{\delta m}}_{N^2}\end{array}\right)}^T---\left(12c\right)$

原结构有限元模型刚度矩阵的修正包括如下两种方法：

(1)ΔK＝(M₀+ΔM)ΦΛΦ^T(ΦΦ^T)^-1-K₀

(2)用Λ^-1前乘式(3)可得，

${\Lambda }^{-1}{\Phi }^{T}K{\Phi }_a^i={\Phi }^{T}M{\Phi }_a^i---\left(14\right)$

用(Λ_aⁱ)^-1后乘式(4)可得，

将式(14)、式(15)两式相减得，

${\Lambda }^{-1}{\Phi }^{T}K{\Phi }_a^i-{\Phi }^{T}K{\Phi }_a^i{\left({\Lambda }_a^i\right)}^{-1}={\Phi }^T{K}_a^i{\Phi }_a^i{\left({\Lambda }_a^i\right)}^{-1}-{\Phi }^T{M}_a^i{\Phi }_a^i---\left(16\right)$

设$S^i={\Phi }^T{K}_a^i{\Phi }_a^i{\left({\Lambda }_a^i\right)}^{-1}-{\Phi }^T{M}_a^i{\Phi }_a^i---\left(17a\right)$

$U^i={\Phi }^{T}K{\Phi }_a^i---\left(17b\right)$

Λ^-1、(Λ_aⁱ)^-1为对角矩阵，将式(17a)和式(17b)代入式(16)并展开后，可得

$(\frac{1}{{\lambda }_j}-\frac{1}{{\lambda }_{\mathrm{ak}}^i})U_{\mathrm{jk}}^i=S_{\mathrm{jk}}^i---\left(18\right)$

将式(9b)代入式(17b)，可以得到

令$U^i-{\Phi }^T{K}_0{\Phi }_a^i=T^i---\left(19b\right)$

等式(19a)右边已知，将式(19a)行拉直后再转置得到

Aⁱδk＝tⁱ            (20a)

其中$A^i={\Phi }^T\otimes {\left({\Phi }_a^i\right)}^T---\left(20b\right)$

δk＝[ΔK₁₁…ΔK_1N |ΔK₂₁…ΔK_2N|……|ΔK_N1…ΔK_NN]^T    (20c)

$t^i={[T_{11}^i···T_{1N}^i|T_{21}^i···T_{2N}^i|······|T_{N1}^i···T_{\mathrm{NN}}^i]}^T---\left(20d\right)$

将分G次改变原结构的质量和刚度所产生的新结构进行模态测试，然后将得到的G个方程(20a)联立求解得

$\left(\begin{array}{c}{A}^1\\ A^2\\ ·\\ ·\\ ·\\ A^G\end{array}\right)\mathrm{\delta k}=\left(\begin{array}{c}{t}^1\\ t^2\\ ·\\ ·\\ ·\\ t^G\end{array}\right)---\left(21a\right)$

其中i＝1、2、3……G；

ΔK＝[unvec_N×N(δk)]^T    (21b)

即$\mathrm{\Delta K}={\left(\begin{array}{cccccc}{\mathrm{\delta k}}_1& {\mathrm{\delta k}}_{N+1}& ·& ·& ·& \delta k_{(N-1)\times N+1}\\ ·& ·& & & ·& ·\\ ·& ·& & & ·& ·\\ ·& ·& & & ·& ·\\ {\mathrm{\delta k}}_N& {\mathrm{\delta k}}_{2N}& ·& ·& ·& {\mathrm{\delta k}}_{N^2}\end{array}\right)}^T---\left(21c\right)$

将修正矩阵ΔM和ΔK代入式(9a)、(9b)即可得原结构有限元模型修正后的质量矩阵M和刚度矩阵K。

本发明同时改变结构的质量和刚度实施起来简单易行，并且运用多次模态试验结果可以有效增加参数识别方程个数，提高有限元模型的修正精度；此外本发明用矩阵计算直接修正有限元模型的质量阵和刚度阵，对使用者工作要求较低，修改过程程序化较强。

附图说明

图1：本发明方法流程图。

具体实施方式

下面结合附图对发明的技术方案进行详细说明：

实际工程中所建立的有限元模型自由度数远大于模态试验测试自由度数，但可以采用Kammer与Serep提出的模态缩聚方法使有限元模型自由度与模态试验测试自由度一致。故本发明假设有限元模型自由度与模态测试自由度位置、数量一致。

首先对原结构进行模态测试，得到原结构的振型Φ和频率Λ；然后分若干次同时改变结构的质量分布(如在不同位置固定安装已知轻重的质量块)和结构刚度(如增加固定约束等)，并对产生的新结构进行模态测试，得到第i次改变结构后新结构的振型Φ_aⁱ和频率Λ_aⁱ。

假设测试模态数为N_e，有限元模型自由度与模态测试自由度均为N，则测试得到的原结构的频率Λ和新结构的频率Λ_aⁱ规模均为N_e×N_e，原结构的振型Φ和新结构的振型Φ_aⁱ规模均为N×N_e。原结构的质量阵M和刚度阵K均为未知量，规模为N×N。第i次结构质量阵的改变量M_aⁱ和刚度阵的改变量K_aⁱ均为已知量，规模为N×N。

此时原结构模态满足，

MΦΛ＝KΦ                            (1)

第i次改变结构质量和刚度后，新结构模态满足，

$(M+M_a^i){\Phi }_a^i{\Lambda }_a^i=(K+K_a^i){\Phi }_a^i---\left(2\right)$

将式(1)转置再后乘Φ_aⁱ，得到，

$\Lambda {\Phi }^{T}M{\Phi }_a^i={\Phi }^{T}K{\Phi }_a^i---\left(3\right)$

用Φ^T前乘式(2)得，

${\Phi }^T(M+M_a^i){\Phi }_a^i{\Lambda }_a^i={\Phi }^T(K+K_a^i){\Phi }_a^i---\left(4\right)$

式(3)、式(4)两式相减得，

$\Lambda {\Phi }^{T}M{\Phi }_a^i-{\Phi }^{T}M{\Phi }_a^i{\Lambda }_a^i={\Phi }^T{M}_a^i{\Phi }_a^i{\Lambda }_a^i-{\Phi }^T{K}_a^i{\Phi }_a^i---\left(5\right)$

由于M为原结构质量阵，未知，设：

$P^i={\Phi }^{T}M{\Phi }_a^i---\left(6a\right)$

$Q^i={\Phi }^T{M}_a^i{\Phi }_a^i{\Lambda }_a^i-{\Phi }^T{K}_a^i{\Phi }_a^i---\left(6b\right)$

其中矩阵Pⁱ、Qⁱ的维数均为N_e×N_e；

式(5)化为，${\mathrm{\Lambda P}}^i-P^i{\Lambda }_a^i=Q^i---\left(7\right)$

由于Λ和Λ_aⁱ都是对角阵，式(7)可展开为N_e×N_e个代数方程，

$({\lambda }_j-{\lambda }_{\mathrm{ak}}^i)P_{\mathrm{jk}}^i=Q_{\mathrm{jk}}^i---\left(8\right)$

j、k为Pⁱ、Qⁱ矩阵中元素的下标，λ_j为原结构的第j个测试特征值，λ_akⁱ为第i次结构质量和刚度改变后新结构的第k个测量特征值。由于总能通过调整结构的附加质量或刚度实现λ_j不等于λ_akⁱ，因此由式(6b)、式(8)可以求出Pⁱ中的元素P_jkⁱ；

原结构有限元模型的质量矩阵和刚度矩阵分别为M₀和K₀，为已知量，规模为N×N。由于建模时对实际结构的简化和建模人员技术水平的限制，有限元模型的质量阵和刚度阵都存在误差。与实际结构的质量阵和刚度阵相比，原结构有限元模型质量阵和刚度阵的修正量分别为ΔM和ΔK，为未知量，规模为N×N。则有，

M＝M₀+ΔM                (9a)

K＝K₀+ΔK            (9b)

将式(9a)代入式(6a)后整理得，

${\Phi }^T\mathrm{\Delta M}{\Phi }_a^i=P^i-{\Phi }^T{M}_0{\Phi }_a^i---\left(10a\right)$

设$P^i-{\Phi }^T{M}_0{\Phi }_a^i=R^i---\left(10b\right)$

由于Φ^T和Φ_aⁱ都是矩形矩阵，不能从式(10a)中直接求解ΔM，故将式(10a)行拉直后再转置得到Aⁱδm＝rⁱ                        (11a)

其中$A^i={\Phi }^T\otimes {\left({\Phi }_a^i\right)}^T---\left(11b\right)$

δm＝[ΔM₁₁…ΔM_1N|ΔM₂₁…ΔM_2N|……|ΔM_N1…ΔM_NN]^T        (11c)

$r^i={[R_{11}^i···R_{1N}^i|R_{21}^i···R_{2N}^i|······|R_{N1}^i···R_{\mathrm{NN}}^i]}^T---\left(11d\right)$

其中Aⁱ矩阵维数为N_e²×N²。

将分G次改变原结构的质量和刚度所产生的新结构进行模态测试，然后将得到的G个方程(11a)联立求解得

$\left(\begin{array}{c}{A}^1\\ A^2\\ ·\\ ·\\ ·\\ A^G\end{array}\right)\mathrm{\delta m}=\left(\begin{array}{c}{r}^1\\ r^2\\ ·\\ ·\\ ·\\ r^G\end{array}\right)---\left(12a\right)$

ΔM＝[unvec_N×N(δm)]^T                        (12b)

即$\mathrm{\Delta K}={\left(\begin{array}{cccccc}{\mathrm{\delta m}}_1& {\mathrm{\delta m}}_{N+1}& ·& ·& ·& \delta m_{(N-1)\times N+1}\\ ·& ·& & & ·& ·\\ ·& ·& & & ·& ·\\ ·& ·& & & ·& ·\\ {\mathrm{\delta m}}_N& {\mathrm{\delta m}}_{2N}& ·& ·& ·& {\mathrm{\delta m}}_{N^2}\end{array}\right)}^T---\left(12c\right)$

这样测试模态总数远大于测试模态自由度数时，从式(12b)中解出的有限元质量阵的修正矩阵ΔM使修正后的质量阵M能较好地逼近真实值。

一旦修正后的质量阵M由式(9a)得到后，有限元刚度阵的修正矩阵ΔK亦可以由式(1)和(9b)推导得到ΔK＝(M₀+ΔM)ΦΛΦ^T(ΦΦ^T)^-1-K₀            (13)

有限元刚度阵的修正矩阵ΔK也可以用与质量阵修正方法相近似的方法进行直接修正计算。

用Λ^-1前乘式(3)可得，

${\Lambda }^{-1}{\Phi }^{T}K{\Phi }_a^i={\Phi }^{T}M{\Phi }_a^i---\left(14\right)$

用(Λ_aⁱ)^-1后乘式(4)可得，

将式(14)、式(15)两式相减得，

${\Lambda }^{-1}{\Phi }^{T}K{\Phi }_a^i-{\Phi }^{T}K{\Phi }_a^i{\left({\Lambda }_a^i\right)}^{-1}={\Phi }^T{K}_a^i{\Phi }_a^i{\left({\Lambda }_a^i\right)}^{-1}-{\Phi }^T{M}_a^i{\Phi }_a^i---\left(16\right)$

设$S^i={\Phi }^T{K}_a^i{\Phi }_a^i{\left({\Lambda }_a^i\right)}^{-1}-{\Phi }^T{M}_a^i{\Phi }_a^i---\left(17a\right)$

$U^i={\Phi }^{T}K{\Phi }_a^i---\left(17b\right)$

由于Λ^-1、(Λ_aⁱ)^-1为对角矩阵，将式(17a)和式(17b)代入式(16)并展开后，可得

$(\frac{1}{{\lambda }_j}-\frac{1}{{\lambda }_{\mathrm{ak}}^i})U_{\mathrm{jk}}^i=S_{\mathrm{jk}}^i---\left(18\right)$

当λ_j不等于λ_akⁱ，未知数Uⁱ中的未知数U_jkⁱ可以被解出。

将式(9b)代入式(17b)，可以得到

设$U^i-{\Phi }^T{K}_0{\Phi }_a^i=T^i---\left(19a\right)$

等式(19a)右边已知。将式(19a)行拉直后再转置得到

Aⁱδk＝tⁱ            (20a)

其中$A^i={\Phi }^T\otimes {\left({\Phi }_a^i\right)}^T---\left(20b\right)$

δk＝[ΔK₁₁…ΔK_1N|ΔK₂₁…ΔK_2N|……|ΔK_N1…ΔK_NN]^T    (20c)

$t^i={[T_{11}^i···T_{1N}^i|T_{21}^i···T_{2N}^i|······|T_{N1}^i···T_{\mathrm{NN}}^i]}^T---\left(20d\right)$

将分G次改变原结构的质量和刚度所产生的新结构进行模态测试，然后将得到的G个方程(20a)联立求解得

$\left(\begin{array}{c}{A}^1\\ A^2\\ ·\\ ·\\ ·\\ A^G\end{array}\right)\mathrm{\delta k}=\left(\begin{array}{c}{t}^1\\ t^2\\ ·\\ ·\\ ·\\ t^G\end{array}\right)---\left(21a\right)$

ΔK＝[unvec_N×N(δk)]^T              (21b)

即$\mathrm{\Delta K}={\left(\begin{array}{cccccc}{\mathrm{\delta k}}_1& {\mathrm{\delta k}}_{N+1}& ·& ·& ·& \delta k_{(N-1)\times N+1}\\ ·& ·& & & ·& ·\\ ·& ·& & & ·& ·\\ ·& ·& & & ·& ·\\ {\mathrm{\delta k}}_N& {\mathrm{\delta k}}_{2N}& ·& ·& ·& {\mathrm{\delta k}}_{N^2}\end{array}\right)}^T---\left(21c\right)$

这样当测试模态总数远大于测试模态自由度数时，从式(21b)中解出的有限元刚度阵的修正矩阵ΔK使修正后的刚度阵K能较好地逼近真实值。

在实际工程中，修正前的有限元模型能够在一定程度上反映结构的真实情况。借助于工程经验，或者根据试验与仿真计算的结果对比，对大多数有限元模型的误差部位或误差元素的位置可以做出初步判断，所以并不需要对有限元模型的质量阵和刚度阵中所有元素进行修正。因此可以分别将式(11a)和式(20a)等式两边做相同初等变换，将需要修正的元素移至方程上部，再通过分块计算可以得到新方程(22)和(23)。

$\left(\begin{array}{cc}{A}_1^i& A_2^i\end{array}\right)\left(\begin{array}{c}\mathrm{\delta M}\\ 0\end{array}\right)=r^i---\left(22\right)$

其中向量δM为向量δm中需要修正的元素组成的向量，矩阵A₁ⁱ为矩阵Aⁱ中与向量v中各元素对应的部分组成的矩阵，矩阵A₂ⁱ为矩阵Aⁱ中与零向量中各元素对应的部分组成的矩阵。

$\left(\begin{array}{cc}{B}_1^i& B_2^i\end{array}\right)\left(\begin{array}{c}\mathrm{\delta K}\\ 0\end{array}\right)=t^i---\left(23\right)$

其中向量δK为向量δk中需要修正的元素组成的向量，矩阵B₁ⁱ为矩阵Bⁱ中与向量δK中各元素对应的部分组成的矩阵，矩阵B₂ⁱ为矩阵Bⁱ中与零向量中各元素对应的部分组成的矩阵。

推出

$A_1^i\mathrm{\delta M}=r^i---\left(24\right)$

$B_1^i\mathrm{\delta K}=t^i---\left(25\right)$

将分G次改变原结构的质量和刚度所产生的新结构进行模态测试，然后将得到的G个方程(24)、方程(25)分别联立求解，即可以得到更精确的有限元模型质量阵和刚度阵中需要修正元素的修正量δM和δK。