航空动力学报
JOURNAL OF AEROSPACE POWER
1998年  第13卷  第1期科技期刊

重特征值模态灵敏度分析的非线性摄动法

西北工业大学 崔博文**　沈允文
【摘要】  利用非线性摄动等式，推导了实对称矩阵的重特征值模态灵敏度计
算公式。该方法的优点是保持了原特征系统矩阵的对称性和带宽。数值算例表
明该方法精度高，易于计算和实施。
主题词：  特征值  模态分析  摄动法
分类号：  O241.6

    最近，Ojalvo［1］，Mills Curran［2］和Dailey［3］等推广了Nelson［4］的方法。
但他们在求系数矩阵时忽略了二阶摄动项，这势必影响其结果准确性。陈［5］文
提出了一个改进方法，但并未论及系数矩阵的求法。Kim［6］虽则考虑了二阶摄动
项，但在求解过程中却采用分块矩阵法。但是正如作者所述，当重特征值数较大时
这种方法可能失效。本文在前述作者方法基础上，考虑二阶摄动项，利用非线性摄
动等式进行求解。

1  理论基础
    实对称矩阵的特征问题为：

(K-λM)x=0　　　　　　　　　　　　(1)

xTMx=1　　　　　　　　　　　　　　(2)

式中：K，M为结构刚度矩阵、质量矩阵；λ，x为特征值及对应的特征向量。

    设特征值问题有m重特征值λi，即λ1=λ2…=λm=λ，相应的特征向量为xi(i=1,
2,…m)。令X=［x1 x2 … xm］，那末相应于重特征值的特征值问题可表示为：

KX=MXΛ　　　　　　　　　　　　　　(3)

XTMX=I　　　　　　　　　　　　　　　(4)

式中Λ=λI∈Rmxm。由于X不唯一，可能使得X不连续且不可导，因而其特征向
量的导数可能不存在。为获取特征向量导数，必须获得满足式(3)，(4)且连续
可导的特征向量(表示为Z)。由于Z和X都满足式(3)，(4)，因而必定存在一正交
变换矩阵A∈Rmxm使下式成立

Z=XZA X∈Rnxm　　　　　　　　　　　　　　(5)

当结构某参数p发生小变化时，结构质量阵M和刚度阵K均发生变化。设M和K的改
变量分别为△M和△K，这时的特征值问题为：

(K+△K)(Z+△Z)=(M+△M)(Z+△Z)(Λ+△Λ)　　　　　　　　(6)

式中：△Λ和△Z分别表示对应于参数p变化后的特征值和特征向量的变化量。

    将式(6)展开，取掉高阶摄动项后有：

(K-λM)△Z=(λ△M-△K)Z+MZ△Λ　　　　　　　　　　　　(7)

式(7)前乘XT，并将Z=XA代入，注意到XT(K-λM)=0，于是可得到下式：

XT(λ△M-△K)XA+XTMXA△Λ=0　　　　　　　　　　　　　　(8)

令D=XT(△K-λ△M)X，则式(8)可表示为：

DA=A△Λ　　　　　　　　　　　　　　　　　　　　　　　(9)

式(9)是一个简单的特征值问题。由于特征值的重数m远小于总自由度数n，因而
从式(9)可很容易获得△Λ和A，进而可获得相应于重特征值λi的特征向量Z=XA。

    利用推广了的Nelson方法［1］，假定对应于参数p变化的特征向量的变化量
△Z具有：

△Z=V+ZC　　　　　　　　　　　　　　　　　　　　　　　(10)

式中：C为系数矩阵∈Rmxm。为求解矩阵V，将式(10)代入式(7)，并定义G=K-λm，
F=(λ△M-△K)Z+MZ△Λ，这时有：

GV=F　　　　　　　　　　　　　　　　　　　　　　　　　(11)

　　因矩阵G为一奇异阵，其秩为(n-m)。因而方程(11)为一奇异线性方程组。
该方法的求解可参见文献［1］。

为便于说明，式(10)可变化为：　　　　　(12)

同样对应的特征问题(6)可表示为：

［(K+△K)-(λi+△λi)(M+△M)］(zi+△zi)=0　　　　　　(13)

展开式(13)，考虑到(K-λiM)zi=0，可得下式：

(K-λiM)△zi+（△K-λi△M）(zi+△zi)-△λi(M+△M)(zi+△zi)=0　　　　(14)

把式(12)代入式(14)，并前乘zTK，并略去三阶摄动项可得到下式：

Cki=(zTk(△K-λi△M)（zi+vi）-△λizTkMvi-△λizTk△Mzi)/(△λi-△λk) (k≠i)　　(15)

当k=i时，考虑到质量正交性有：

(zi+△zi)T(M+△M)(zi+△zi)=1　　　　　　　　　　　　　(16)

展开式(16)，略去高阶摄动项后有：

△zTiMzi+zTiM△zi+zTi△Mzi=0　　　　　　　　　　　　　(17)

把式(12)代入式(17)后，并化简可得到：

Cii=-(1/2)(vTiMzi+zTiMvi+zTi△Mzi)　　　　　　　　　　(18)

求得矩阵V及系数矩阵C后，可利用式(10)求得重特征值的特征向量的变化量△Z。

    设△K和△M表示由于结构参数p变化△p而产生的刚度矩阵和质量矩阵的增量。
△λi和｛△zi｝为特征值和特征向量的增量，系统振动的特征值和特征向量的导数
可近似表示为：

λ′i=△λi／△p ， z′i=△zi／△p　　　　　　　　　　(19)

2  数值算例

    以两个8自由度梁单元组成的方形悬臂梁(图1)作为本文方法的应用。取梁横截
面y方向尺寸h为设计变量p，位移向量取为：

q=［v1，w1，θy1，θz1，v2，w2，θy2，θz2］T            (20)

悬臂梁的总刚度矩阵K9和总质量矩阵M9分别为：
(21)

(22)

式中：E为杨氏模量，L为梁的长度，ρ为质量密度，h和b分别为横截面尺寸(见图1)。

图1 方形悬臂梁
    运用矩阵理论［7］，可获得Kg和Mg对尺寸h的偏导矩阵K′g和M′g。当结
构尺寸h发生小变化△h时的刚度矩阵和质量矩阵的变化阵可分别表示为：

△Kg=K′g·△h=（Kg/h）△h                 (23)

△Mg=M′g·△h=（Mg/h）△h                 (24)

    为便于计算，令b=h=10，ρ=4.2，L=1，E=1000。由于结构关于y-z对称，求
解特征问题KgX=λMgX可得到四对重特征值，2730.38；61405.83；772445.82；
5997738.62。当△h=0.001时，分别用文献所述方法和本文方法计算了最低阶重
特征值的导数及其对应的特征向量导数，并用差分法检验本文方法计算精度。摄
动矩阵为Kp=Kg+0.001K′g，Mp=Mg+0.001M′g，求解KpXp=MpXpΛp特征问题，满足
约束条件XTpMpXp=I。并令Z′=(Xp-Z)／0.001，Λ′=(λp-λI)／0.001。计算结果见表1。

表1 文献方法、本文方法及差分法所得结果比较
重特征值λ1=λ2=2730.38
特征值
导 数Ojalvo方法Mills方法Dailey方法本文方法差分方法
8949.549371.278949.549371.278949.549371.278949.549371.278949.969829.87
特 征
向 量
导 数-0.00349
0.00014
0.00010
0.01303
-0.00327
0.00005
0.00007
-0.001820
-0.0050
-0.0245
0
0
0.00542
0.00715
0-0.00350
0.00014
0.00010
0.01310
-0.00327
0.00049
0.00007
-0.001810
-0.0050
-0.0245
0
0
0.00542
0.00715
0-0.00349
0.00014
0.00010
0.01303
-0.00327
0.00005
0.00007
-0.001820
-0.0050
-0.0245
0
0
0.00542
0.00715
0-0.00349
0
0
0.01303
-0.00327
0
0
-0.001820
-0.0050
-0.0245
0
0
0.00542
0.00715
0-0.00350
0
0
0.01304
-0.0316
0
0
-0.001780
-0.00497
-0.02449
0
0
0.00566
0.00724
0

    由表1结果比较可知，文献所述方法由于在计算特征向量导数时忽略二阶摄动项，
精度较低，而本文方法则具有良好的计算精度。

3  结  论

    本文提出的计算重特征值模态灵敏度的新方法修正了Ojalvo［1］，Kim［6］等人方
法存在的缺陷。利用非线性摄动等式并考虑二阶摄动项进行模态灵敏度计算，保证了计
算精度。除需要求解一低阶特征值问题外，本文方法仅需要进行矩阵乘法运算及解一高
维线性方程组。由于本文方法保持了矩阵M和K的对称、带宽及稀疏性，因而线性方程组
的求解效率很高。本文方法特别适合于结构刚度矩阵和质量矩阵不能表示为某结构参数p
的显式函数时的情形。通过模态灵敏度可方便地进行结构动力修改，以使结构具有良好
的动态特性。当重特征值的变化(或导数)仍重复时，本文方法会失效。有关这方面的问
题将在以后再进行深入研究。

男  西安市西北工业大学324信箱  710072
参考文献
1Ojalvo I U.Efficient Computation of Modal Sensitivities for Systems with 
Repeated Frequencies.AIAA Journal,1988,26(3):361-366
2Mills-Curran W C.Calculation of Eigenvector Derivatives for Structure 
with Repeated Eigenvalues.AIAA Journal,1988，26(7):867-871
3Dailey R L.Eigenvector Derivatives with Repeated Eigenvalues.AIAA 
Journal,1989,27(4):486-491
4Nelson R B.Simplified Calculation of Eigenvector Derivatives.AIAA 
Journal,1976,14(Sept.):1201-1205
5陈塑寰，宋大同，韩万芝．重特征值的特征向量导数计算的新方法．机械强度，
1995，17(2)：43-47
6Kim K O,Wallerstein D V.Modal Design Sensitivities for Multiple 
Eigenvalues.Comput.Struct.1988，29：755-762
7程云鹏．矩阵论．西北工业大学出版社，1989

1996年12月收稿；1997年3月收到修改稿。

（责任编辑  杨再荣）
