利用有限元技术,$n$ 自由度的离散保守力学系统的自由响应可以由下面的二阶微分方程来刻画
其中 $M_a,K_a \in {\bf R}^{n \times n}$ 分别是分析质量矩阵与刚度矩阵,$q(t)$ 是 $n\times 1$ 位移向量,$f(t)$ 是 $n\times 1$ 外力向量. 通常,$M_a$ 是实对称正定矩阵,$K_a$ 是实对称半正定矩阵,方程 (1.1) 通常被称为有限元分析模型. 考虑方程 (1.1) 的齐次部分,并假定方程 (1.1) 的位移响应是谐波的,即
那么结构特征问题可表示为
其中 $\lambda_j=\omega_j^2$ 是第 $j$ 个特征值,$x_j$ 是第 $j$ 个特征向量. 众所周知,特征值与特征向量在物理学上可解释为自然频率的平方及模态振型. 设
其中 $\Lambda_1=\mbox{diag}\{\lambda_1,\cdots,\lambda_p\},$ $\Lambda_2=\mbox{diag}\{\lambda_{p+1},\cdots,\lambda_n\},$ $ X_1=[x_1,\cdots,x_p] ,$ $X_2= [x_{p+1},\cdots,x_n].$ 易见 (1.2)式的 $n$ 个特征值-特征向量关系式可以组装为如下的矩阵方程
无阻尼振动模态的一个重要性质是它们关于质量矩阵的正交性,即
在现代结构的动力学分析中,许多努力均致力于推导准确的结构动力学有限元模型. 许多应用,如验证创新结构设计,评估现有结构对地震或强风的影响,结构控制的实施和结构的健康监测等都需要准确的结构模型. 但有限元模型可能有一些误差或不确定性,例如应用有限元技术将连续的结构离散为一个在节点上互相关联的离散结构会产生离散误差. 而其它的不确定性可能是由于工程师对单元的选择,对几何特性的假设及对结构的约束条件,边界条件和连接条件等进行简化处理所导致. 当用动力测试对分析模型进行验证时,通常自然频率和模态,与实测结果会不一致. 结构动力学模型修正即为利用结构现场实测的振动信息修正结构有限元分析模型,使得修正后结构分析的模态参数与试验值趋于一致.
在进行模型修正时,总期望测量的低阶模态数据能够融于修正模型而不改变原模型的高阶模态数据,这样的修正称为无溢出模型修正. 设 $\Lambda_{1} \in {{\bf R}}^{p\times p},\ X_{1}\in {{\bf R}}^{n \times p}$ 是已知的初始模型的前 $p$ 个特征对. 剩余的特征信息 $\Lambda_{2} \in{{\bf R}}^{(n-p) \times (n-p)},\ X_{2}\in {{\bf R}}^{n \times (n-p)}$ 是未知的且需保持不变的. 之所以要关注无溢出现象是由于测量的模态数据一般是不完整的,而原模型的高阶模态数据或是可接受的而不必在修正中引进新的振动,或是工程师对这些高阶模态数据根本就一无所知.
在过去的40多年,结构动力学模型修正问题受到了相当多的关注. 应用测量响应数据修正无阻尼结构系统的质量矩阵与刚度矩阵的许多模型修正技术,已有文献[1-10] 等进行了讨论,这些方法均为直接法且计算效率较高. 然而,这些方法不能保证额外的,虚假的模态不被引入到有意义的频率范围[11],即,这些方法不能保证修正是无溢出的. 最近,在假定质量矩阵是准确的情况下, 文献[12]提出了一个修正无阻尼结构系统的无溢出直接修正法. 文献[13-15] 提出了修正阻尼结构系统的无溢出方法. 文献[16] 发展了保持矩阵正定性的无溢出直接方法. 注意到用迭代方法进行模型修正在这些年被关注较少,本文的目的就是建立一个修正有限元模型的无溢出迭代方法. 在数学上,同时修正质量矩阵与刚度矩阵可以归结为如下的反特征值问题.
问题 IEP 设 $\Upsilon=\mbox{diag}\{\mu_1,\cdots,\mu_p\}\in {\bf R}^{p \times p}$ 和 $Z=[z_1,\cdots,z_p]\in{\bf R}^{n \times p}$ 分别是测量的特征值与特征向量矩阵,其中 $p\ll n.$ 求对称矩阵 $M$ 和 $K$ 使得
通过计算可知方程 (1.5) 与 (1.6) 有 $n(n+1)$ 未知数而只有 $n^2$ 个方程,这意味着由方程 (1.5) 与 (1.6) 求得的矩阵 $M$ 和 $K$ 不唯一. 众所周知 $M_a$ 和 $K_a$ 是 $M$ 和 $K$ 的较好的近似矩阵. 因此,要得到改进的模型,即为求满足方程 (1.5) 与 (1.6) 的,且与 $M_a$ 和 $K_a$ 偏离程度最小的矩阵 $M$ 和 $K$. 这样,我们需进一步求解下面的最佳逼近问题.
问题 Ⅱ 设 ${{\cal S}}_{MK}$ 是问题 IEP 的解集,求 $(\hat{M},\hat{K}) \in {{\cal S}}_{MK}$ 使得
文献[17] 考虑了用迭代方法求解阻尼模型修正问题,他们建立了一个广义牛顿方法并应用于逆二次特征值问题的对偶问题,该方法可保证平方收敛且具有较高的计算效率. 表面看来,本文研究的问题是文献[17] 的特殊情形 ($C=0$),但主要的差别是本文考虑了无溢出现象. 本文组织如下: 在第二节,提出了求解问题 IEP 和问题 Ⅱ 的一个有效迭代方法,通过此方法,在不计舍入误差的情况下,通过选取特殊的初始矩阵对,经有限步迭代就可得到问题 Ⅱ 的极小 Frobenius 范数对称解 $(\hat{M},\hat{K})$; 在第三节,给出了两个数值例子用来检验本文算法的有效性$\!;$ 第四节对本文结论进行了总结.
本文记 ${\bf R}^{m \times n}$ 为 $m \times n$ 矩阵空间,${\bf R}^{n \times n}$ 中的对称矩阵全体记为 ${\bf SR}^{n \times n}$. $A^{\top},$ $\mbox{tr}(A)$ 和 $R(A)$ 分别表示矩阵 $A$ 的转置,迹及列空间. $I_n$ 表示 $n$ 阶单位矩阵. 矩阵空间 $ {\bf R}^{m \times n}$ 的内积定义为
则 ${\bf R}^{m \times n}$ 是一个 Hilbert 空间,此内积诱导的矩阵范数 ${\| \cdot \|}$ 即为 Frobenius 范数. 给定矩阵 $A=[a_{ij}] \in {\bf R}^{m \times n}$,$ B=[b_{ij}] \in {\bf R}^{p \times q}$,$A$ 和 $B$ 的 Kronecker 积定义为 $A \otimes B=[a_{ij}B]\in {\bf R}^{mp \times nq}.$ 此外,对 $m\times n$ 矩阵 $A=[a_1,a_2,\cdots,a_n]$,其中 $a_j,j=1,\cdots,n,$ 是 $A$ 的第 $j$ 个列向量,拉直函数 $\mbox{vec}(A)$ 定义为 $\mbox{vec}(A)=[a_1^\top,a_2^\top,\cdots,a_n^\top]^\top.$ 设 $A,B$ 和 $X$ 是三个适当维数的矩阵,则可得如下的等式[18]
由 (1.3) 和 (1.4)式,可得
由 (2.1) 式及注意到 $X_2$ 是列满秩阵,则方程 (1.5) 可等价地表为
其中 $H=M_a^{-1}K_a M_a^{-1}-X_1\Lambda_1 X_1^\top,\ L=M_a^{-1}-X_1X_1^\top.$ 在工程实践,方程 (1.5) 中的矩阵 $X_2$ 和 $\Lambda_2$ 通常是未知的,注意到在方程 (2.2) 中没有出现矩阵 $X_2$ 和 $\Lambda_2$. 对给定的对称矩阵对 $(M_a,K_a)$,由 (1.6) 和 (2.2) 式可得
注意到
设
那么求解矩阵逼近问题 Ⅱ 等价于求解下列矩阵方程的极小 Frobenius 范数对称解
一旦求得方程 (2.3) 的极小 Frobenius 范数解 $(\tilde{M}^*,\tilde{K}^*)$,则可计算最佳逼近问题 Ⅱ 的唯一解,此时,唯一解可表示为
现在,我们给出求解矩阵方程 (2.3) 的迭代方法.
算法 2.1
S1. 输入矩阵 $M_a,K_a,Z,X_1,\Upsilon,\Lambda_1,$ 并任意选取 $n \times n$ 对称矩阵 $\tilde{M}_1$ 及 $\tilde{K}_1$.
S2. 计算
S3. 若 $R_s=0,$ 或 $R_s \neq 0 $ 但 $ P_s=0,\ Q_s=0,$ 则停止; 否则,$s:=s+1$.
S4. 计算
S5. 转 S3.
由算法 2.1 易知,$\tilde{M}_s,P_s \in {\bf SR}^{n \times n},$ $\tilde{K}_s,Q_s \in {\bf SR}^{n \times n},$ $s=1,2,\cdots.$
定义2.1 设 $Y,Z \in {\bf R}^{m\times n}$,若 ${\rm{tr}}(Y^\top Z)=0$,则称矩阵 $Y,Z$ 是相互正交的.
关于算法 2.1,可得下面的基本性质.
引理2.1 由算法 2.1 产生的数列 $\{R_i\},\{P_i\}$ 和 $\{Q_i\}$ 满足关系式
引理 2.1 的证明见附录 A.
引理2.2 若 $(\tilde{M}^*,\tilde{K}^*)$ 是方程 (2.3) 的任一解,则对任意的对称矩阵对 $(\tilde{M}_1,\tilde{K}_1),$ $\tilde{M}_1,$ $\tilde{K}_1 \in {\bf SR}^{n \times n}$,可得
其中数列 $ \{\tilde{M}_i\},\{\tilde{K}_i\},\{R_i\},\{P_i\}$ 及 $\{Q_i\}$ 由算法 2.1 给出.
引理 2.2 的证明见附录 B.
定理2.1 对任意给定的初始矩阵对 $(\tilde{M}_1,\tilde{K}_1),$ $\tilde{M}_1,\tilde{K}_1 \in {\bf SR}^{n \times n}$,则在不计舍入误差的情况下,经有限步迭代就可得到方程 $(2.3)$ 的一个对称解.
证 若 $R_l\neq 0,\ l=1,2,\cdots,n^2+np$,据引理 2.2 可知 $P_l\neq 0$ 或 $Q_l\neq 0$. 这样由算法 2.1 可以计算 $R_{n^2+np+1}$ 和 $(\tilde{M}_{n^2+np+1},\tilde{K}_{n^2+np+1})$,由引理 2.1 可得
及
因此,$\{R_1,R_2,\cdots,R_{n^2+np}\}$ 构成矩阵空间 $\Phi$ 的一组正交基,其中
这意味着 $R_{n^2+np+1}=0$,即 $(\tilde{M}_{n^2+np+1},\tilde{K}_{n^2+np+1})$ 是方程 $(2.3)$ 的对称解.
引理2.3 方程 $(2.3)$ 有对称解 $(\tilde{M},\tilde{K})$ 当且仅当下列矩阵方程组
是相容的.
证 设 $(\tilde{M}^*,\tilde{K}^*)$ 是方程 (2.3) 的一个对称解,则 $\tilde{M}^* Z \Upsilon -\tilde{K}^* Z=A,\tilde{M}^* H -\tilde{K}^* L=0,$ 且
即 $(\tilde{M}^*,\tilde{K}^*)$ 是方程组 (2.7) 的解. 反之,若方程组 (2.7) 有解,设为 $\tilde{M}=U,\ \tilde{K}=V$. 记
则 $\tilde{M}^* $ 和 $ \tilde{K}^*$ 是对称阵,且
由此可知,$(\tilde{M}^*,\tilde{K}^*)$ 是方程 (2.3) 的一个对称解.
引理2.4 [19] 设相容线性方程组 $Ax=b$ 有解 $\tilde{x} \in R(A^\top)$,那么 $\tilde{x}$ 是此线性方程组的唯一极小 Frobenius 范数解.
应用 Kronecker 积及拉直函数,易知矩阵方程组 (2.7) 等价于
设 $G \in {\bf R}^{n \times p}$,$W \in {\bf R}^{n \times n}$ 为任意矩阵,则
显然,若选取
则由算法 2.1 产生的矩阵 $\tilde{M}_s$ 和 $\tilde{K}_s$ 满足
据引理 2.4 可知,若由式 (2.8) 选取初始矩阵对,其中 $G \in {\bf R}^{n \times p},$ $W \in {\bf R}^{n \times n}$ 为任意矩阵,那么由算法 2.1 所得的矩阵对 $(\tilde{M}^*,\tilde{K}^*)$ 即为方程 (2.3) 的极小 Frobenius 范数对称解. 总结上面的讨论即得下面的结果.
定理2.2 若据式 $(2.8)$ 选取初始矩阵对,其中 $G \in {\bf R}^{n \times p},$ $W \in {\bf R}^{n \times n}$ 为任意矩阵,特别地,可取 $\tilde{M}_1=0$ 及 $\tilde{K}_1=0$,那么在不计舍入误差的情况下,经有限步迭代即可得到方程 $(2.3)$ 的唯一极小 Frobenius 范数对称解,从而可得问题 Ⅱ 的唯一解,其唯一解由式 $(2.4)$ 给出.
在这一节,我们将给出两个数值例子来验证算法的有效性,所有的测试均使用计算软件 MATLAB 6.5.
例3.1 考虑一个$5\,$-$\!\!$自由度系统,其质量矩阵与刚度矩阵分别为
用来产生实验数据的仿真模型为 $M=M_a,$
注意到 $K_a$ 与 $K$ 的区别仅在 $(3,3)$ 和 $(4,4)$ 位置元素. 仿真模型的特征解用来产生实验模态数据,假设测量的特征值与特征向量矩阵分别为 $ \Upsilon $ 和 $Z,$ 其中
分析模型的前两阶模态数据为
据算法 2.1 并选取 $\tilde{M}_1=0,\ \tilde{K}_1=0,$ 经 $77$ 步迭代,可得方程 (2.3) 的极小 Frobenius 范数对称解 $(\tilde{M}^*,\tilde{K}^*)$:
相应的误差为
图 1 给出了误差随迭代步变化的情况. 因此,由式 (2.4),问题 Ⅱ 的最优逼近解为
修正的质量矩阵与刚度矩阵的相对误差为
$(\Upsilon,Z)$ 和 $(\Lambda_2,X_2)$ 的相对误差为
因此,指定的特征值 (矩阵 $\Upsilon$ 的对角元 ) 和特征向量 (矩阵 $Z$ 的列向量) 已嵌入新模型 $ \hat{M}Z\Upsilon=\hat{K}Z.$ 进一步,从 (3.1) 式可看到,修正系统的剩余谱数据与原模型的剩余谱数据一致,这意味着修正是无溢出的,并且矩阵 $ \hat{M} $ 和矩阵 $ \hat{K} $ 仍为对称正定矩阵.
下面的测试例子来自文献 [20-21],这里我们作了些许修改 (原来的分析质量矩阵为 $M_a=4I_n.$)
例3.2 假设分析质量矩阵与刚度矩阵分别为
实验仿真模型为 $M=M_a$ 及 $K=0.85K_a,$ 并取仿真模型的前四阶模态数据作为测量的实验数据.
选取初始迭代矩阵为 $\tilde{M}_1=0$ 及 $\tilde{K}_1=0,$ 并且若误差满足 $\|R_s\|\leq 1.0e-07$ 时,矩阵对 $( \tilde{M}_s,\tilde{K}_s)$ 即被认为是方程 (2.3) 的极小 Frobenius 范数对称解. 数值结果见表 3.1.
其中 "It." 表示迭代步数,
$\hat{M}=M_a+\tilde{M}_s,\ \hat{K}=K_a+ \tilde{K}_s.$ 有趣的是,从表 3.1 看到,在精度要求不太高的情况下,迭代步数并不一定随着矩阵阶数的提高而增加.
本文给出了一个利用测量模态数据同时修正无阻尼结构系统的质量矩阵与刚度矩阵的无溢出迭代方法. 利用此方法,在不计舍入误差的情况下,通过选取特殊的初始矩阵对,经有限步迭代,即可得到问题 Ⅱ 的最佳逼近解 $(\hat{M},\hat{K}).$ 给出的两个数值例子验证了算法的有效性.
A 引理 2.1 的证明
证 因为 ${\rm{tr}}(R_j^\top R_i)={\rm{tr}}(R_i^\top R_j),$ $ {\rm{tr}}(P_j^\top P_i)={\rm{tr}}(P_i^\top P_j)$ 及 ${\rm{tr}}(Q_j^\top Q_i)={\rm{tr}}(Q_i^\top Q_j)$,所以我们只需证明
我们用数学归纳法证明这个结果,证明分两步.
第一步 我们首先证明
$i=1$ 时,由算法 2.1 并注意到 $P_1 \in {{\bf SR}}^{n \times n},$ $Q_1 \in {{\bf SR}}^{n \times n}$,可得
其中 $\delta_1=\frac{\|R_{1}\|^2}{\|P_{1}\|^2+\|Q_{1}\|^2}$.
应用已证明的结果: $ {\rm{tr}}(R_2^\top R_1)=0,$ 可得
假设式 (A.1) 在 $i=t-1$ 时成立,对 $i=t$,我们有
其中 $\delta_t=\frac{\|R_{t}\|^2}{\|P_{t}\|^2+\|Q_{t}\|^2}$.
因此$\!,\ i=t$ 时,式 (A.1) 成立,由数学归纳法可知对所有的 $i,$ 式 (A.1) 均成立. 第二步 假设
我们将证明
据 (A.1),上面的等式在 $i=d$ 时成立. 对 $i=1,$ 我们有
对 $1<i<d,$ 可得
其中 $\delta_{d} =\frac{\|R_{d}\|^2}{\|P_{d}\|^2+ \|Q_{d}\|^2}.$
由上面的结果可得 $ {\rm{tr}}(R_{d+1}^\top R_i)=0.$ 因此,对 $1\leq i<d,$ 我们可得
其中 $\xi = \frac{\|P_{i}\|^2+\|Q_{i}\|^2}{\|R_{i}\|^2}.$ 据数学归纳法,由第一步和第二步可知结论 (2.5) 成立.
B 引理 2.2 的证明
证 我们用数学归纳法证明关系式 (2.6). 当 $i=1$ 时,
现在假设 $1 \leq i \leq t-1$ 时结论 (2.6) 成立,那么我们可得
由数学归纳法可知引理 2.2 成立.