南京航空航天大学学报  2017, Vol. 49 Issue (5): 606-611   PDF    
基于测量试验数据修正有限元模型质量矩阵
戴华, 魏伟     
南京航空航天大学理学院,南京,210016
摘要: 对无阻尼结构系统有限元模型质量矩阵修正问题,以该矩阵修正量的F-范数为目标函数,并以待修正质量矩阵应具有的性质,如满足正交关系,对称性,半正定性和稀疏性作为约束条件,数学上形成带约束的矩阵最佳逼近问题。给出了问题有解的条件,基于循环投影方法,提出了求解矩阵最佳逼近问题的数值方法。数值结果说明了所给方法的有效性。
关键词: 无阻尼结构系统     模型修正     矩阵最佳逼近问题     投影方法    
Mass Matrix Correction of Finite Element Model Using Measured Test Data
DAI Hua, WEI Wei     
College of Science, Nanjing University of Aeronautics & Astronautics, Nanjing, 210016, China
Abstract: For problem of correcting mass matrix of finite element model of undamped structural systems, the desired mass matrix properties, including orthogonality relation, symmetry, positive semi-definiteness and sparsity, are imposed as side constraints to form mathematically the optimal matrix approximation problems. The solvability conditions for the problems are presented. Based on the cyclic projection method, numerical methods are proposed for solving the matrix nearness problems. Numerical results show that the proposed methods are efficient.
Key words: undamped structural system     model updating     matrix nearness problem     projection method    

振动对于航空、航天、船舶、机械、电子、大型桥梁等许多工程结构常常是造成恶性破坏的直接原因。因此,在大型结构设计中,振动设计与控制是至关重要的。对工程结构进行定量、准确的动力学分析,解决工程结构中普遍存在的振动问题,首先必须建立结构的动力学模型。有限元方法具有分析速度快、设计周期短、与结构试验相比成本低等优点,已成为结构动力系统建模最常用、最有效的方法之一。由有限元方法可建立无阻尼结构系统的运动微分方程

$ {\mathit{\boldsymbol{M}}_a}\mathit{\boldsymbol{\ddot x}}\left( t \right) + {\mathit{\boldsymbol{M}}_a}\mathit{\boldsymbol{x}}\left( t \right) = \mathit{\boldsymbol{f}}\left( t \right) $

式中:x(t)∈Rn为位移向量;n为有限元模型的自由度;Ma, KaRn×n分别为系统有限元模型的质量和刚度矩阵; f(t)∈Rn为载荷向量。通常,有限元模型的质量和刚度矩阵都是对称半正定矩阵,并且具有稀疏结构。结构系统的固有频率ω或特征值λ(λ=ω2)以及相应的振型或特征向量x由如下特征方程确定

$ {\mathit{\boldsymbol{K}}_a}\mathit{\boldsymbol{x}} = \lambda {\mathit{\boldsymbol{M}}_a}\mathit{\boldsymbol{x}} $ (1)

应用特征值问题的数值方法[1],可求得广义特征值问题(1)的特征值λi(a)和相应的特征向量xi(a)

随着试验模态分析技术、信号处理和动态测试技术的发展,由振动试验可获得实际结构系统的一些低阶频率ωi(e)(i=1, 2, …, mn),即特征值λi(e)=(ωi(e))2以及相应的振型xi(e)。由于实际结构的复杂性,有限元方法必须对结构的边界条件、连接条件和约束条件进行简化处理,还必须对结构的几何特性做一些假设,因而有限元方法所建立的结构有限元模型往往不能准确地反映实际结构的动态特性,有限元模型的计算结果λi(a)xi(a)与实测结果λi(e)xi(e)之间存在着误差。另一方面,受试验条件和成本的限制,实测信息的不完整导致试验建模无法获得完备的数学模型。有限元模型修正就是将有限元建模和试验建模结合起来,先用有限元方法建立结构的有限元模型,然后依据系统动态特性的试验测量数据修正有限元模型,使得修正后结构模型的低阶动态特性参数与试验值趋于一致,并且使有限元模型在保持拓扑结构的同时改变量尽可能小。数学上,有限元模型修正问题归结为求矩阵,使其满足低阶动态特性、保持物理和结构性质的同时,并使矩阵的改变量达到最小。修正的有限元动力学模型可用于结构响应分析和控制。

四十多年来,结构有限元模型修正的理论和方法取得了许多成果[2-3],并且在工程实际中得到了重要应用[4-5]。因为有限元模型质量和刚度矩阵是实际结构的较好近似,所以最早发展并被广泛使用的有限元模型修正方法是以系统矩阵或矩阵元素为修正对象的矩阵型和元素型修正方法。Berman等[6-7]利用测量的部分振型修正质量矩阵M,要求修正的质量矩阵M满足对称性MT=M和正交性条件XeTMXe=Im,并使质量矩阵的改变量最小,应用Lagrange乘子法,导出了修正质量矩阵M的表达式。Zhang和Lee等[8-9]分别利用矩阵变换和Moore-Penrose广义逆给出了修正质量矩阵M的显式表示。但这些方法修正的质量矩阵M不仅改变了解析质量矩阵Ma的带状或稀疏结构,而且未必是半正定的。为保持修正质量矩阵的稀疏结构,Wei和Zhang[10], Cha[11]以质量矩阵的待修正元素为对象, 发展了修正质量矩阵的元素型修正方法。Yuan和Dai[12]考虑了解析质量矩阵的一个主子矩阵准确、修正其余元素的问题,利用广义奇异值分解给出了修正质量矩阵的表达式。这类方法具有灵活、方便等特点,但难以保证修正质量矩阵M的半正定性。

为了保证修正的质量矩阵既满足对称半正定性,又具有稀疏结构,本文以质量矩阵改变量的F-范数为目标函数,并以待修正质量矩阵应具有的性质,如满足正交关系或特征方程,对称性,半正定性和稀疏性作为约束条件,数学上将问题归结为带约束的矩阵最佳逼近问题。应用矩阵的QR分解,给出了该问题有解的条件。基于循环投影方法[13],提出了求解这类问题的数值方法,提供了两个数值例子以说明所给方法的有效性。

为了表述方便,对使用的记号作如下说明:rank(A)和A+分别表示矩阵A的秩和Moore-Penrose广义逆。对A, BRm×n, (A, B)=tr(BTA)和A×B分别表示矩阵AB的内积和Hadamard积,基于此内积,Rm×n构成Hilbert空间,并且由此内积导出的长度即为矩阵的F-范数‖*‖FSRn×n表示n阶实对称矩阵的全体,SR0n×n表示n阶对称半正定矩阵的集合,ORn×n表示n阶正交矩阵的全体。对ASRn×nA≥0表示A为实对称半正定矩阵。sparse(M)=sparse(Ma)表示矩阵MMa具有相同的稀疏结构。In表示n阶单位矩阵。记Λe=diag(λ1(e), …, λm(e)),Xe=[x1(e), …, xm(e)],并设Λe>0,rank(Xe)=mMa=(maij), Ka=(kaij)∈SR0n×n,并且Ma正定。

1 利用正交关系修正质量矩阵

为了使修正质量矩阵M既满足如下正交性条件

$ \mathit{\boldsymbol{X}}_e^T\mathit{\boldsymbol{M}}{\mathit{\boldsymbol{X}}_e} = {\mathit{\boldsymbol{I}}_m} $ (2)

又是对称半正定的,并保持稀疏结构,本文将修正质量矩阵M的问题归结为如下带约束的矩阵逼近问题

$ \begin{array}{*{20}{c}} {\min \frac{1}{2}\left\| {\mathit{\boldsymbol{M}} - {\mathit{\boldsymbol{M}}_a}} \right\|_{\rm{F}}^2}\\ {{\rm{s}}{\rm{.}}\;{\rm{t}}{\rm{.}}\;\;\mathit{\boldsymbol{X}}_e^T\mathit{\boldsymbol{M}}{\mathit{\boldsymbol{X}}_e} = {\mathit{\boldsymbol{I}}_m},{\mathit{\boldsymbol{M}}^{\rm{T}}} = \mathit{\boldsymbol{M}} \ge 0}\\ {{\rm{sparse}}\left( \mathit{\boldsymbol{M}} \right) = {\rm{sparse}}\left( {{\mathit{\boldsymbol{M}}_a}} \right)} \end{array} $ (3)

设测量的振型矩阵XeQR分解为

$ {\mathit{\boldsymbol{X}}_e} = \mathit{\boldsymbol{Q}}\left( \begin{array}{l} \mathit{\boldsymbol{R}}\\ {\bf{0}} \end{array} \right) $ (4)

式中:Q=[Q1, Q2]∈ORn×n, Q1Rn×mRm阶可逆上三角矩阵。将式(4)代入式(2),可导出满足正交条件(2)的对称矩阵为

$ \mathit{\boldsymbol{M}} = \mathit{\boldsymbol{Q}}\left( {\begin{array}{*{20}{c}} {{{\left( {\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{R}}^{\rm{T}}}} \right)}^{ - 1}}}&{{\mathit{\boldsymbol{M}}_{12}}}\\ {\mathit{\boldsymbol{M}}_{12}^{\rm{T}}}&{{\mathit{\boldsymbol{M}}_{22}}} \end{array}} \right){\mathit{\boldsymbol{Q}}^{\rm{T}}} $ (5)

式中:M12Rm×(nm)M22SR(nm)×(nm)是任意的。

为了求得矩阵方程(2)的对称半正定解,需要如下引理。

引理1[14] 设ESRn×n, FRn×k, GSRk×k$ \boldsymbol{H }= \left( {\begin{array}{*{20}{c}} \boldsymbol{E}&\boldsymbol{F}\\ {{\boldsymbol{F}^{\rm{T}}}}&\boldsymbol{G} \end{array}} \right)$,则H≥0当且仅当E≥0, GFTE+F≥0,且rank[E, F]=rank(E)。

由引理1,即得矩阵方程(2)的对称半正定解M

$ \mathit{\boldsymbol{M}} = \mathit{\boldsymbol{Q}}\left( {\begin{array}{*{20}{c}} {{{\left( {\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{R}}^{\rm{T}}}} \right)}^{ - 1}}}&{{\mathit{\boldsymbol{M}}_{12}}}\\ {\mathit{\boldsymbol{M}}_{12}^{\rm{T}}}&{\mathit{\boldsymbol{M}}_{12}^{\rm{T}}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{R}}^{\rm{T}}}{\mathit{\boldsymbol{M}}_{12}} + \mathit{\boldsymbol{G}}} \end{array}} \right){\mathit{\boldsymbol{Q}}^{\rm{T}}} $ (6)

式中:M12Rm×(nm)GSR0(nm)×(nm)是任意的。

为了给出问题(3)解的存在性和唯一性,需要如下引理。

引理2(最佳逼近定理)[15] 设H是一个Hilbert空间,WH中的一个非空闭凸集,则对gH,存在唯一的向量g*W,使得

$ \left\| {\mathit{\boldsymbol{g}} - {\mathit{\boldsymbol{g}}^ * }} \right\| = \mathop {\min }\limits_{\mathit{\boldsymbol{x}} \in \mathit{\boldsymbol{W}}} \left\| {\mathit{\boldsymbol{g}} - \mathit{\boldsymbol{x}}} \right\| $ (7)

式中:x表示Hilbert空间中向量x的长度或由长度导出的范数。

问题(7)的解g*通常称为gW中的最佳逼近或称g*gW中的投影,记为g*=PW(g)。

如果引理2中Hilbert空间H的非空闭凸集W= $\bigcap\limits_{i = 1}^k {{\boldsymbol{W}_i}} $,其中Wi(i=1, 2, …, k)都是Hilbert空间H中的闭凸集,并且对gH,仅能计算g在每个Wi上的投影PWi(g)。Boyle和Dykstra[13]给出了计算gW中投影的一个迭代方法,称为循环投影方法,并证明了循环投影方法的收敛性。

由引理2可得如下结论。

定理1  设测量的振型矩阵XeQR分解为式(4),如果存在M12Rm×(nm)GSR0(nm)×(nm),使得式(6)中的质量矩阵M满足sparse(M)=sparse(Ma),则问题(3)的可行域是非空闭凸集,从而问题(3)存在唯一解。

$ {S_{{M_1}}} = \left\{ {\mathit{\boldsymbol{M}} \in {\mathit{\boldsymbol{R}}^{n \times n}}\left| {\mathit{\boldsymbol{X}}_e^{\rm{T}}{\mathit{\boldsymbol{M}}_a}{\mathit{\boldsymbol{X}}_e} = {\mathit{\boldsymbol{I}}_m}} \right.,{\mathit{\boldsymbol{M}}^{\rm{T}}} = \mathit{\boldsymbol{M}}} \right\} $
$ {S_{{M_2}}} = \left\{ {\mathit{\boldsymbol{M}} \in {\mathit{\boldsymbol{R}}^{n \times n}}\left| {{\rm{sparse}}\left( \mathit{\boldsymbol{M}} \right) = {\rm{sparse}}\left( {{\mathit{\boldsymbol{M}}_a}} \right)} \right.} \right\} $
$ {S_{{M_3}}} = \left\{ {\mathit{\boldsymbol{M}} \in {\mathit{\boldsymbol{R}}^{n \times n}}\left| {{\mathit{\boldsymbol{M}}^{\rm{T}}} = \mathit{\boldsymbol{M}} \ge 0} \right.} \right\} $

SM3=SR0n×n,问题(3)的可行域为SM1SM2SM3,并且是Rn×n中的闭凸集。因此问题(3)等价于求MaSM1SM2SM3中的最佳逼近。下面依次考虑MaSM1, SM2, SM3上的投影。

对解析质量矩阵MaSRn×n,记

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}^{\rm{T}}}{\mathit{\boldsymbol{M}}_a}\mathit{\boldsymbol{Q}} = \left( {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{M}}_{a11}}}&{{\mathit{\boldsymbol{M}}_{a12}}}\\ {\mathit{\boldsymbol{M}}_{a12}^{\rm{T}}}&{{\mathit{\boldsymbol{M}}_{a22}}} \end{array}} \right)}\\ {{\mathit{\boldsymbol{M}}_{aij}} = \mathit{\boldsymbol{Q}}_i^{\rm{T}}{\mathit{\boldsymbol{M}}_a}{\mathit{\boldsymbol{Q}}_j}\;\;\;\;i,j = 1,2} \end{array} $ (8)

MSM1,由式(5,8)和F-范数的正交不变性,可得

$ \begin{array}{*{20}{c}} {\left\| {\mathit{\boldsymbol{M}} - {\mathit{\boldsymbol{M}}_a}} \right\|_{\rm{F}}^2 = \left\| {\left( {\begin{array}{*{20}{c}} {{{\left( {\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{R}}^{\rm{T}}}} \right)}^{ - 1}}}&{{\mathit{\boldsymbol{M}}_{12}}}\\ {\mathit{\boldsymbol{M}}_{12}^{\rm{T}}}&{{\mathit{\boldsymbol{M}}_{22}}} \end{array}} \right) - {\mathit{\boldsymbol{Q}}^{\rm{T}}}{\mathit{\boldsymbol{M}}_a}\mathit{\boldsymbol{Q}}} \right\|_{\rm{F}}^2 = }\\ {\left\| {{{\left( {\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{R}}^{\rm{T}}}} \right)}^{ - 1}} - {\mathit{\boldsymbol{M}}_{a11}}} \right\|_{\rm{F}}^2 + 2\left\| {{\mathit{\boldsymbol{M}}_{12}} - {\mathit{\boldsymbol{M}}_{a12}}} \right\|_{\rm{F}}^2 + }\\ {\left\| {{\mathit{\boldsymbol{M}}_{22}} - {\mathit{\boldsymbol{M}}_{a22}}} \right\|_{\rm{F}}^2} \end{array} $

当且仅当M12=Ma12, M22=Ma22$ \mathop {\min }\limits_{\boldsymbol{M} \in \boldsymbol{S}{\boldsymbol{R}^{n \times n}}} \boldsymbol{M} - \boldsymbol{M}_{a{\rm{F}}}^2$取得最小值。因此, 解析质量矩阵MaSM1中的投影为

$ {P_{{S_{{M_1}}}}}\left( {{\mathit{\boldsymbol{M}}_a}} \right) = \mathit{\boldsymbol{Q}}\left( {\begin{array}{*{20}{c}} {{{\left( {\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{R}}^{\rm{T}}}} \right)}^{ - 1}}}&{{\mathit{\boldsymbol{M}}_{a12}}}\\ {\mathit{\boldsymbol{M}}_{a12}^{\rm{T}}}&{{\mathit{\boldsymbol{M}}_{a22}}} \end{array}} \right){\mathit{\boldsymbol{Q}}^{\rm{T}}} $

由式(4)可得Xe+=R-1Q1T。利用Moore-Penrose广义逆的性质和式(8),上式可化为

$ {P_{{S_{{M_1}}}}}\left( {{\mathit{\boldsymbol{M}}_a}} \right) = {\mathit{\boldsymbol{M}}_a} + {\left( {\mathit{\boldsymbol{X}}_e^ + } \right)^{\rm{T}}}\mathit{\boldsymbol{X}}_e^ + - {\mathit{\boldsymbol{X}}_e}\mathit{\boldsymbol{X}}_e^ + {\mathit{\boldsymbol{M}}_a}{\mathit{\boldsymbol{X}}_e}\mathit{\boldsymbol{X}}_e^ + $ (9)

对给定的解析质量矩阵MaSRn×n,使用如下指标集SMa刻画矩阵Ma的稀疏结构

$ \begin{array}{*{20}{c}} {{S_{{M_a}}} = \left\{ {\left( {i,j} \right)\left| {{m_{aij}} = 0,} \right.} \right.}\\ {\left. {{\mathit{\boldsymbol{M}}_a} = \left( {{m_{aij}}} \right) \in \mathit{\boldsymbol{S}}{\mathit{\boldsymbol{R}}^{n \times n}}} \right\}} \end{array} $ (10)

由指标集SMa的定义,容易证明如下结论。

引理3  设T=(tij)∈SRn×n,则问题

$ \begin{array}{*{20}{c}} {\min \frac{1}{2}\left\| {\mathit{\boldsymbol{M}} - \mathit{\boldsymbol{T}}} \right\|_{\rm{F}}^2}\\ {{\rm{s}}{\rm{.}}\;{\rm{t}}{\rm{.}}\;{\rm{sparse}}\left( \mathit{\boldsymbol{M}} \right) = {\rm{sparse}}\left( {{\mathit{\boldsymbol{M}}_a}} \right)} \end{array} $ (11)

有唯一解M=(mij)∈SRn×n,且mij= $\left\{ \begin{array}{l} 0,\left( {i,j} \right) \in {S_{{M_a}}}\\ {t_{ij}},\left( {i,j} \right) \notin {S_{{M_a}}}。\end{array} \right. $

由引理3可知,问题(11)的M就是TSM2中的投影,即

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{M}} = {P_{{S_{{M_2}}}}}\left( \mathit{\boldsymbol{T}} \right) = \left( {{m_{ij}}} \right)}\\ {{m_{ij}} = \left\{ \begin{array}{l} 0\;\;\;\;\left( {i,j} \right) \in {S_{{M_a}}}\\ {t_{ij}}\;\;\;\left( {i,j} \right) \notin {S_{{M_a}}} \end{array} \right.} \end{array} $ (12)

为了计算一个矩阵在SM3中的投影,需要如下引理。

引理4[16]  设ARn×n$ \boldsymbol{\hat A} = \frac{1}{2}$(A+AT)的谱分解为$\boldsymbol{\hat A} $=Udiag(θ1, …, θn)UT,其中θ1, …, θn为矩阵$\boldsymbol{\hat A} $的特征值, UORn×n,则存在唯一的对称半正定矩阵A+,使得

$ {\left\| {{\mathit{\boldsymbol{A}}_ + } - \mathit{\boldsymbol{A}}} \right\|_{\rm{F}}} = \mathop {\min }\limits_{\mathit{\boldsymbol{X}} \in \mathit{\boldsymbol{SR}}_0^{n \times n}} {\left\| {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{A}}} \right\|_{\rm{F}}} $ (13)

并且A+=Udiag(β1, …, βn)UT,其中βi=max{θi, 0}(i=1, 2, …, n)。

对矩阵ARn×n,问题(13)的解称为矩阵A的正逼近,记为A+。实际上,A+ARn×nSM3中的投影,即A+=PSM3(A)=PSR0n×n(A)。

于是,可给出求解问题(3)的数值算法如下:

算法1  修正质量矩阵的循环投影方法

输入:MaSRn×nXeRn×m和收敛准则ε,并置M1, 0=MaI0, j=0(j=1, 2, 3)。

输出:修正的质量矩阵M

i=1, 2, …,执行以下运算直到收敛

j=1, 2, 3,计算

$ {\mathit{\boldsymbol{M}}_{i,j}} = {P_{{S_{{M_j}}}}}\left( {{\mathit{\boldsymbol{M}}_{i,j - 1}},{\mathit{\boldsymbol{I}}_{i - 1,j}}} \right) $
$ {\mathit{\boldsymbol{I}}_{i,j}} = {\mathit{\boldsymbol{M}}_{i,j}} - \left( {{\mathit{\boldsymbol{M}}_{i,j - 1}} - {\mathit{\boldsymbol{I}}_{i - 1,j}}} \right) $

Mi+1, 0=Mi, 3, 如果‖Mi+1, 0Mi, 0Fε,则M=Mi+1, 0,停止;否则继续迭代。

定理2[13]  如果问题(3)的可行域非空,则对任意j(j=1, 2, 3),算法1所产生的矩阵序列{Mn, j}n=1收敛于问题(3)的唯一解。

由算法1修正的质量矩阵不仅满足正交性条件(2),而且是对称半正定的,与解析质量矩阵Ma具有相同的稀疏结构。

2 利用特征方程修正质量矩阵

假定解析刚度矩阵Ka准确,并设测量的特征值Λe和相应的振型Xe满足

$ \mathit{\boldsymbol{X}}_e^{\rm{T}}{\mathit{\boldsymbol{K}}_a}{\mathit{\boldsymbol{X}}_e} = {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_e} $ (14)

基于待修正质量矩阵M应满足如下特征方程

$ \mathit{\boldsymbol{M}}{\mathit{\boldsymbol{X}}_e}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_e} = {\mathit{\boldsymbol{K}}_a}{\mathit{\boldsymbol{X}}_e} $ (15)

和对称半正定性,戴华[17]将修正质量矩阵的问题归结为如下带约束的矩阵逼近问题

$ \begin{array}{*{20}{c}} {\min \frac{1}{2}\left\| {\mathit{\boldsymbol{M}} - {\mathit{\boldsymbol{M}}_a}} \right\|_{\rm{F}}^2}\\ {{\rm{s}}{\rm{.}}\;{\rm{t}}{\rm{.}}\;\mathit{\boldsymbol{M}}{\mathit{\boldsymbol{X}}_e}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_e} = {\mathit{\boldsymbol{K}}_a}{\mathit{\boldsymbol{X}}_e},{\mathit{\boldsymbol{M}}^{\rm{T}}} = \mathit{\boldsymbol{M}} \ge 0} \end{array} $ (16)

$ {D_M} = \left\{ {\mathit{\boldsymbol{M}} \in {\mathit{\boldsymbol{R}}^{n \times n}}\left| {\mathit{\boldsymbol{M}}{\mathit{\boldsymbol{X}}_e}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_e} = {\mathit{\boldsymbol{K}}_a}{\mathit{\boldsymbol{X}}_e},{\mathit{\boldsymbol{M}}^{\rm{T}}} = \mathit{\boldsymbol{M}} \ge 0} \right.} \right\} $

DM是问题(16)的可行域。

为了给出DM中元素的表达式,需要如下引理。

引理5[18]  设X, BRn×m,如果rank(X)=m,并且XQR分解X=Q$\left( {\begin{array}{*{20}{c}} \boldsymbol{R}\\ \boldsymbol{0} \end{array}} \right) $,其中Q=[Q1, Q2]∈ORn×n, Q1Rn×mRm阶可逆上三角矩阵,则AX=B有解ASR0n×n当且仅当

$ {\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{B}} = {\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{X}} \ge 0,{\rm{rank}}\left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{B}}} \right) = {\rm{rank}}\left( \mathit{\boldsymbol{B}} \right)。$

在有解的情况下,AX=B的对称半正定解可表示为A=BX++[(BX+)T+(InXX+)B(XTB)+BT](InXX+)+Q2GQ2T,其中GSR0(nm)×(nm)是任意的。

定理3  设测量的特征值Λe和相应的振型Xe满足式(14),且XeQR分解为式(4),则问题(16)的可行域非空,并且问题(16)存在唯一解

$ \mathit{\boldsymbol{M}} = {\mathit{\boldsymbol{M}}_0} + {\mathit{\boldsymbol{Q}}_2}{\left[ {\mathit{\boldsymbol{Q}}_2^{\rm{T}}\left( {{\mathit{\boldsymbol{M}}_a} - {\mathit{\boldsymbol{M}}_0}} \right){\mathit{\boldsymbol{Q}}_2}} \right]_ + }\mathit{\boldsymbol{Q}}_2^{\rm{T}} $ (17)

其中M0=KaXeΛe-2XeTKa

证明  由引理5知,矩阵方程MXeΛe=KaXe存在对称半正定解,并且其通解可表示为

$ \mathit{\boldsymbol{M}} = {\mathit{\boldsymbol{K}}_a}{\mathit{\boldsymbol{X}}_e}\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_e^{ - 2}\mathit{\boldsymbol{X}}_e^{\rm{T}}{\mathit{\boldsymbol{K}}_a} + {\mathit{\boldsymbol{Q}}_2}\mathit{\boldsymbol{GQ}}_2^{\rm{T}} $ (18)

式中:GSR0(nm)×(nm)是任意的,从而问题(16)的可行域是非空闭凸集。由引理2可知,问题(16)存在唯一解。

对任意MDMM可表示为

$ \mathit{\boldsymbol{M}} = {\mathit{\boldsymbol{M}}_0} + {\mathit{\boldsymbol{Q}}_2}\mathit{\boldsymbol{GQ}}_2^{\rm{T}} = {\mathit{\boldsymbol{M}}_0} + \mathit{\boldsymbol{Q}}\left( {\begin{array}{*{20}{c}} {\bf{0}}&{\bf{0}}\\ {\bf{0}}&\mathit{\boldsymbol{G}} \end{array}} \right){\mathit{\boldsymbol{Q}}^{\rm{T}}} $

由F-范数的正交不变性,可得

$ \begin{array}{*{20}{c}} {\left\| {\mathit{\boldsymbol{M}} - {\mathit{\boldsymbol{M}}_a}} \right\|_{\rm{F}}^2 = \left\| {\mathit{\boldsymbol{Q}}\left( {\begin{array}{*{20}{c}} {\bf{0}}&{\bf{0}}\\ {\bf{0}}&\mathit{\boldsymbol{G}} \end{array}} \right){\mathit{\boldsymbol{Q}}^{\rm{T}}} - \left( {{\mathit{\boldsymbol{M}}_a} - {\mathit{\boldsymbol{M}}_0}} \right)} \right\|_{\rm{F}}^2}=\\ {\left\| {\left( {\begin{array}{*{20}{c}} {\bf{0}}&{\bf{0}}\\ {\bf{0}}&\mathit{\boldsymbol{G}} \end{array}} \right) - {\mathit{\boldsymbol{Q}}^{\rm{T}}}\left( {{\mathit{\boldsymbol{M}}_a} - {\mathit{\boldsymbol{M}}_0}} \right)\mathit{\boldsymbol{Q}}} \right\|_{\rm{F}}^2 = }\\ {\left\| {\mathit{\boldsymbol{Q}}_1^{\rm{T}}\left( {{\mathit{\boldsymbol{M}}_a} - {\mathit{\boldsymbol{M}}_0}} \right){\mathit{\boldsymbol{Q}}_1}} \right\|_{\rm{F}}^2 + 2\left\| {\mathit{\boldsymbol{Q}}_1^{\rm{T}}\left( {{\mathit{\boldsymbol{M}}_a} - {\mathit{\boldsymbol{M}}_0}} \right){\mathit{\boldsymbol{Q}}_2}} \right\|_{\rm{F}}^2 + }\\ {\left\| {\mathit{\boldsymbol{G}} - \mathit{\boldsymbol{Q}}_2^{\rm{T}}\left( {{\mathit{\boldsymbol{M}}_a} - {\mathit{\boldsymbol{M}}_0}} \right){\mathit{\boldsymbol{Q}}_2}} \right\|_{\rm{F}}^2} \end{array} $

从而║MMaF2=min当且仅当║GQ2T(MaM0)Q2F2=min。由引理4可知,当G=[Q2T(MaM0)Q2]+时,║MMaF2取最小值。将G=[Q2T(MaM0)Q2]+代入式(18)即得式(17)。

定理3给出的修正质量矩阵M不仅满足特征方程(15),而且是对称半正定矩阵。袁永新和戴华[19]利用奇异值分解给出了问题(16)解的表达式。但这些表达式不能保持解析质量矩阵的稀疏结构。

为了使修正质量矩阵M既满足特征方程(15),又是对称半正定的,并保持稀疏结构,本文将修正质量矩阵M的问题归结为如下带约束的矩阵逼近问题

$ \begin{array}{*{20}{c}} {\min \frac{1}{2}\left\| {\mathit{\boldsymbol{M}} - {\mathit{\boldsymbol{M}}_a}} \right\|_{\rm{F}}^2}\\ {{\rm{s}}{\rm{.}}\;{\rm{t}}{\rm{.}}\;\mathit{\boldsymbol{M}}{\mathit{\boldsymbol{X}}_e}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_e} = {\mathit{\boldsymbol{K}}_a}{\mathit{\boldsymbol{X}}_e},{\mathit{\boldsymbol{M}}^{\rm{T}}} = \mathit{\boldsymbol{M}} \ge 0}\\ {{\rm{sparse}}\left( \mathit{\boldsymbol{M}} \right) = {\rm{sparse}}\left( {{\mathit{\boldsymbol{M}}_a}} \right)} \end{array} $ (19)

由定理3和引理2,可得如下结论。

定理4  设测量的特征值Λe和振型Xe满足式(14),且XeQR分解为式(4)。如果存在矩阵GSR0(nm)×(nm),使得M=KaXeΛe-2XeTKa+Q2GQ2T满足sparse(M)=sparse(Ma),则问题(19)的可行域非空,从而问题(19)存在唯一解。

$ {S_1} = \left\{ {\mathit{\boldsymbol{M}} \in {\mathit{\boldsymbol{R}}^{n \times n}}\left| {\mathit{\boldsymbol{M}}{\mathit{\boldsymbol{X}}_e}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_e} = {\mathit{\boldsymbol{K}}_a}{\mathit{\boldsymbol{X}}_e},{\mathit{\boldsymbol{M}}^{\rm{T}}} = \mathit{\boldsymbol{M}} \ge 0} \right.} \right\} $
$ {S_2} = \left\{ {\mathit{\boldsymbol{M}} \in {\mathit{\boldsymbol{R}}^{n \times n}}\left| {{\rm{sparse}}\left( \mathit{\boldsymbol{M}} \right) = {\rm{sparse}}\left( {{\mathit{\boldsymbol{M}}_a}} \right)} \right.} \right\} $

则问题(19)的可行域DM=S1S2

对给定的矩阵MaSRn×n,由定理3可得MaS1上的投影为

$ {P_{{S_1}}}\left( {{\mathit{\boldsymbol{M}}_a}} \right) = {\mathit{\boldsymbol{M}}_0} + {\mathit{\boldsymbol{Q}}_2}{\left[ {\mathit{\boldsymbol{Q}}_2^{\rm{T}}\left( {{\mathit{\boldsymbol{M}}_a} - {\mathit{\boldsymbol{M}}_0}} \right){\mathit{\boldsymbol{Q}}_2}} \right]_ + }\mathit{\boldsymbol{Q}}_2^{\rm{T}} $ (20)

因为S2=SM2,对给定的矩阵T=(tij)∈SRn×n,由式(12)容易计算TS2上的投影PS2(T)=PSM2(T)。

于是,给出求解问题(19)的数值算法如下:

算法2  修正质量矩阵的循环投影方法

输入:Ka, MaSRn×n, XeRn×m, ΛeRm×m和收敛准则ε,置M1, 0=MaI0, j=0(j=1, 2)。

输出:修正的质量矩阵M

i=1, 2, …,执行以下运算直到收敛。

j=1, 2,计算

$ {\mathit{\boldsymbol{M}}_{i,j}} = {P_{{S_j}}}\left( {{\mathit{\boldsymbol{M}}_{i,j - 1}} - {\mathit{\boldsymbol{I}}_{i - 1,j}}} \right) $
$ {\mathit{\boldsymbol{I}}_{i,j}} = {\mathit{\boldsymbol{M}}_{i,j}} - \left( {{\mathit{\boldsymbol{M}}_{i,j - 1}} - {\mathit{\boldsymbol{I}}_{i - 1,j}}} \right) $

Mi+1, 0=Mi, 2, 如果‖Mi+1, 0Mi, 0Fε,则M=Mi+1, 0,停止;否则继续迭代。

由文献[13]中定理2,可得算法2的收敛性结果。

定理5  如果问题(19)的可行域非空,则对任意j(j=1, 2),算法2所产生的矩阵序列{Mn, j}n=1收敛于问题(19)的唯一解。

3 数值结果

本节给出两个数值例子,以说明算法1和算法2的有效性。算法1和算法2用MATLAB R2015b编程, 程序在主频为2.90 GHz的计算机上实现。收敛准则ε=10-7。用M$\boldsymbol{\tilde M} $分别表示精确的质量矩阵和计算的修正质量矩阵,Iter和CPU分别表示执行算法的迭代步数和计算时间(单位: s),记

$ {\mathit{\boldsymbol{E}}_1} = {\left\| {\mathit{\boldsymbol{X}}_e^{\rm{T}}\mathit{\boldsymbol{\tilde M}}{\mathit{\boldsymbol{X}}_e} - \mathit{\boldsymbol{I}}} \right\|_{\rm{F}}} $
$ {\mathit{\boldsymbol{E}}_2} = {\left\| {\mathit{\boldsymbol{\tilde M}}{\mathit{\boldsymbol{X}}_e}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_e} - {\mathit{\boldsymbol{K}}_a}{\mathit{\boldsymbol{X}}_e}} \right\|_{\rm{F}}},{\rm{Err}} = \frac{{{{\left\| {\mathit{\boldsymbol{\tilde M}} - \mathit{\boldsymbol{M}}} \right\|}_{\rm{F}}}}}{{{\mathit{\boldsymbol{M}}_{\rm{F}}}}} $

例1  考虑图 1所示的无阻尼弹簧-质点系统,其中mi=1 kg(i=1, 2, …, 5),ki=0.5 N/m(i=1, 2, …, 6),则系统精确的质量矩阵与刚度矩阵分别为

图 1 无阻尼弹簧-质点系统 Figure 1 Undamped spring-mass system

$ \mathit{\boldsymbol{M}} = {\rm{diag}}\left( {1,1,1,1,1} \right) $
$ \mathit{\boldsymbol{K}} = {\mathit{\boldsymbol{K}}_a} = \left( {\begin{array}{*{20}{c}} 1&{ - 0.5}&{}&{}&{}\\ { - 0.5}&1&{ - 0.5}&{}&{}\\ {}&{ - 0.5}&1&{ - 0.5}&{}\\ {}&{}&{ - 0.5}&1&{ - 0.5}\\ {}&{}&{}&{ - 0.5}&1 \end{array}} \right) $

取广义特征值问题Kx=λMx的二个最小特征值及其相应的特征向量作为测量特征对, 分别记为Λe, Xe, 即

$ {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_e} = {\rm{diag}}\left( {0.134\;0,0.500\;0} \right) $
$ {\mathit{\boldsymbol{X}}_e} = \left( {\begin{array}{*{20}{c}} {0.288\;7}&{ - 0.500\;0}\\ {0.500\;0}&{ - 0.500\;0}\\ {0.577\;4}&0\\ {0.500\;0}&{0.500\;0}\\ {0.288\;7}&{0.500\;0} \end{array}} \right) $

容易验证:ΛeXe满足式(14)。M的估计为Ma=M+μEr×M,其中Er为随机生成的五阶对称矩阵, 其元素随机分布在区间[-1.0, 1.0], 参数μ表示对M相对扰动的大小。利用算法1和算法2计算的修正质量矩阵$\boldsymbol{\tilde M} $与精确的质量矩阵M完全相同。数值结果见表 1

表 1 例1的数值结果 Table 1 Numerical results of Example 1

例2  考虑一端固定、另一端自由的离散变截面杆[20],基于线性形函数的解析质量矩阵与刚度矩阵为

$ \mathit{\boldsymbol{M}} = \left( {\begin{array}{*{20}{c}} {0.4}&{0.1}&{}&{}&{}&{}&{}&{}&{}&{}\\ {0.1}&{0.4}&{0.1}&{}&{}&{}&{}&{}&{}&{}\\ {}&{0.1}&{0.4}&{0.1}&{}&{}&{}&{}&{}&{}\\ {}&{}&{0.1}&{0.4}&{0.1}&{}&{}&{}&{}&{}\\ {}&{}&{}&{0.1}&{0.6}&{0.2}&{}&{}&{}&{}\\ {}&{}&{}&{}&{0.2}&{0.8}&{0.2}&{}&{}&{}\\ {}&{}&{}&{}&{}&{0.2}&{0.8}&{0.2}&{}&{}\\ {}&{}&{}&{}&{}&{}&{0.2}&{0.8}&{0.2}&{}\\ {}&{}&{}&{}&{}&{}&{}&{0.2}&{0.8}&{0.2}\\ {}&{}&{}&{}&{}&{}&{}&{}&{0.2}&{0.4} \end{array}} \right) $
$ \mathit{\boldsymbol{K}} = {\mathit{\boldsymbol{K}}_a} = \left( {\begin{array}{*{20}{c}} 2&{ - 1}&{}&{}&{}&{}&{}&{}&{}&{}\\ { - 1}&2&{ - 1}&{}&{}&{}&{}&{}&{}&{}\\ {}&{ - 1}&2&{ - 1}&{}&{}&{}&{}&{}&{}\\ {}&{}&{ - 1}&2&{ - 1}&{}&{}&{}&{}&{}\\ {}&{}&{}&{ - 1}&4&{ - 3}&{}&{}&{}&{}\\ {}&{}&{}&{}&{ - 3}&6&{ - 3}&{}&{}&{}\\ {}&{}&{}&{}&{}&{ - 3}&6&{ - 3}&{}&{}\\ {}&{}&{}&{}&{}&{}&{ - 3}&6&{ - 3}&{}\\ {}&{}&{}&{}&{}&{}&{}&{ - 3}&6&{ - 3}\\ {}&{}&{}&{}&{}&{}&{}&{}&{ - 3}&3 \end{array}} \right) $

这里刚度矩阵和质量矩阵的元素缩放了一定比例。取广义特征值问题Kx=λMx的4个最小特征值及其相应的特征向量作为测量特征对, 分别记为Λe, XeM的估计为Ma=M+μEr×M,其中Er为随机生成的10阶对称矩阵, 其元素随机分布在区间[-1.0, 1.0]。利用算法1和算法2计算的修正质量矩阵$\boldsymbol{\tilde M} $与精确的质量矩阵M完全相同。数值结果见表 2

表 2 例2的数值结果 Table 2 Numerical results of Example 2

例1, 2的数值结果说明由算法1和算法2修正的质量矩阵不仅满足正交性条件或特征方程, 而且是对称半正定矩阵, 还具有要求的稀疏结构。因此,算法1, 2是有效的,并且由表 1, 2可见,如果刚度矩阵已知,由特征方程修正质量矩阵的计算效率更高。

4 结束语

本文研究了有限元模型质量矩阵的修正问题。建立了以待修正质量改变量的F-范数为目标函数并以质量矩阵应具有的性质如满足正交关系或特征方程,对称性,半正定性和稀疏性为约束条件的矩阵最佳逼近问题。给出了问题有解的条件。提出了求解带约束矩阵最佳逼近问题的循环投影方法。所获得的修正质量矩阵不仅满足正交性条件或特征方程,而且是对称半正定的,并与解析质量矩阵具有相同的稀疏结构。本文所发展的修正质量矩阵的理论和方法可应用于修正有限元模型的刚度矩阵或同时修正有限元模型的质量和刚度矩阵。

参考文献
[1] BATHE K J, WILSON E L. Numerical methods in finite element analysis[M]. Englewood Cliffs: Prentice-Hall, 1976.
[2] MOTTERSHEAD J E, FRISWEEL M I. Model updating in structural dynamics: A survey[J]. Journal of Sound and Vibration, 1993, 167: 347–375. DOI:10.1006/jsvi.1993.1340
[3] 朱安文, 曲广吉, 高耀南, 等. 结构动力模型修正技术的发展[J]. 力学进展, 2002, 32(3): 337–348. DOI:10.6052/1000-0992-2002-3-J2001-003
ZHU Anwen, QU Guangji, GAO Yaonan, et al. A survey of the modifying techniques of structure dynamic models[J]. Advances in Mechanics, 2002, 32(3): 337–348. DOI:10.6052/1000-0992-2002-3-J2001-003
[4] FRISWEEL M I, MOTTERSHEAD J E. Finite element model updating in structural dynamics[M]. Dordrecht: Kluwer Academic Publishers, 1995.
[5] 张德文, 魏阜旋. 模型修正与破损诊断[M]. 北京: 科学出版社, 1999.
ZHANG Dewen, WEI Fuxuan. Model updating and damage detection[M]. Beijing: Science Press, 1999.
[6] BERMAN A. Mass matrix correction using an incomplete set of measured modes[J]. AIAA Journal, 1979, 17: 1147–1148. DOI:10.2514/3.61290
[7] BERMAN A, NAGY E J. Improvement of a large analytical model using test data[J]. AIAA Journal, 1983, 21: 1168–1173. DOI:10.2514/3.60140
[8] ZHANG D, ZHANG L. Matrix transformation method for updating dynamic model[J]. AIAA Journal, 1992, 30: 1440–1443. DOI:10.2514/3.11083
[9] LEE E T, RAHMATALLA S, EUN H C. Estimation of parameter matrices based on measured data[J]. Applied Mathematical Modelling, 2011, 35: 4816–4823. DOI:10.1016/j.apm.2011.03.048
[10] WEI F S, ZHANG D W. Mass matrix modification using element correction method[J]. AIAA Journal, 1989, 27: 119–121. DOI:10.2514/3.10106
[11] CHA P D. Correcting system matrices using the orthogonality conditions of distinct measured modes[J]. AIAA Journal, 2000, 38: 730–732. DOI:10.2514/2.1023
[12] YUAN Y X, DAI H. Inverse problems for symmetric matrices with a submatrix constraint[J]. Applied Numerical Mathematics, 2007, 57: 646–656. DOI:10.1016/j.apnum.2006.07.030
[13] BOYLE J P, DYKSTRA R L. A method for finding projections onto the intersection of convex sets in Hilbert spaces[C]//Advances in Order Restricted Statistical Interference. Berlin: Springer, 1986: 28-47.
[14] ALBERT A. Conditions for positive and nonnegative definiteness in terms of pseudoinverses[J]. SIAM Journal on Applied Mathematics, 1969, 17: 430–440. DOI:10.1137/0117040
[15] AUBIN J P. Applied functional analysis[M]. New York: John Wiley & Sons, 1979.
[16] HIGHAM N J. Computing a nearest symmetric positive semidefinite matrix[J]. Linear Algebra and Its Applications, 1988, 103: 103–118. DOI:10.1016/0024-3795(88)90223-6
[17] 戴华. 用振动试验最优校正刚度、柔度和质量矩阵[J]. 振动工程学报, 1988, 1(2): 18–27.
DAI Hua. Optimal correction of stiffness, flexibility and mass matrices using vibration tests[J]. Journal of Vibration Engineering, 1988, 1(2): 18–27.
[18] 张磊. 对称非负定矩阵反问题解存在的条件[J]. 计算数学, 1989, 11(4): 337–343.
ZHANG Lei. The solvability conditions for the inverse problem of symmetric nonnegative definite matrices[J]. Mathematica Numerica Sinica, 1989, 11(4): 337–343.
[19] 袁永新, 戴华. 用振动测量数据最优修正振型矩阵与质量矩阵[J]. 工程数学学报, 2007, 24(4): 631–638.
YUAN Yongxin, DAI Hua. Optimal correction of modal matrix and mass matrix using test data[J]. Chinese Journal of Engineering Mathematics, 2007, 24(4): 631–638.
[20] 田霞, 戴华. 杆有限元模型的特征值反问题[J]. 高等学校计算数学学报, 2007, 29(2): 146–156.
TIAN Xia, DAI Hua. Inverse eigenvalue problems for the finite element model of the rod[J]. Numerical Mathematics—A Journal of Chinese Universities, 2007, 29(2): 146–156.