基于隨機攝動方法的廣義隨機復特征值問題研究

仇翯辰, 樊維超

(1. 中國商飛復合材料中心，典型結構部,上海 201324)(2. 中國商飛北京民用飛機技術研究中心，民用飛機結構與復合材料北京市重點實驗室，北京 102211)

0 引言

20世紀90年代初期，邱志平教授將非概率不確定理論引入結構振動問題中，研究了不確定參數對結構固有頻率的影響[1]。Ghosh將結構特征值和特征向量用混沌多項式(PCE)方法進行展開，開展了線性隨機系統的不確定性分析[2]。

在進行實際工程結構的分析和設計優化時，當設計參數發生變化和擾動，結構力學行為(如固有振動特性)也會發生相應的變化，那么為了獲得最優設計方案，迭代計算或者說反復的結構重分析是無法避免的。然而，對于大型復雜結構來說，諸如Monte-Carlo模擬和分層抽樣等迭代計算方法的時間消耗是巨大的，計算成本過高。因此，能夠進行快速靈敏度分析和結構重分析的矩陣攝動方法，便受到學者們的極大關注，近年來得到蓬勃發展。較高的計算效率和較少的時間消耗是矩陣攝動方法的顯著特征，因此學者們紛紛投身于該領域的研究之中。陳塑寰系統地闡述了矩陣攝動理論在結構振動分析中的應用[3]，在后來的著作中，陳塑寰又論述了結構動力學的矩陣攝動理論，開展了孤立特征值的攝動分析和線性振動虧損系統的矩陣攝動分析等研究[4]?？紤]具有實特征值或者復特征值的隨機結構，其特征值的統計量(如均值和方差)的計算對于結構的不確定性傳播分析、動力學重分析和設計優化具有重要指導意義和參考價值。通過對不確定性傳播進行準確的定量化，攝動之后實/復特征值的變化范圍能夠被準確地預測，那么攝動之后的特征值(或固有頻率)和原設計使用工況之間的矛盾則能夠被避免，這樣結構系統的性能和可靠性將得到充分保障。對于傳統的矩陣攝動理論，求結構隨機特征值的統計量的常用方法是將剛度矩陣和質量矩陣在隨機結構參數的均值附近進行Taylor展開，而隨機特征值的方差由特征值的敏感度矩陣、隨機結構參數的標準差矩陣和結構參數的相關系數矩陣求得。然而，在工程實際中，結構參數的相關系數矩陣是很難獲得的，一般只能通過假設或者大量的樣本統計試驗獲得，其成本極高，工程可實現性較差。因此使用傳統的矩陣攝動方法進行不確定性傳播分析受到了嚴重的制約，亟待建立一種改進的工程適用性強的矩陣攝動方法。

針對具有復特征值的隨機結構系統，提出一種改進的不確定傳播分析方法。通過該方法，在考慮結構參數不確定性的基礎上，隨機結構的復特征值的方差可以被方便地求出。

對于隨機實模態結構，在之前的研究工作[6]中建立了一種改進的隨機攝動方法來克服瓶頸。該方法將結構振動特征值和特征向量的一階攝動展式與概率理論相結合，那么結構隨機特征值的方差就可以由剛度矩陣的攝動項、質量矩陣的攝動項和原始特征向量直接推導出來。上述方法在對實模態結構隨機特征值變化范圍的計算中，不涉及結構參數的相關系數矩陣，在工程實際中是易于實現的。具體操作上，只需要知道隨機結構參數本身的方差(或標準差)，而不是結構參數之間的相關系數矩陣，這使得計算隨機結構特征值的統計學性質變得更為便利。

實模態結構隨機攝動方法被進一步改進并推廣到復模態特征值領域。將改進的復模態結構隨機攝動方法與概率理論相結合，推導出了隨機復特征值方差關于隨機結構參數的顯式表達式。根據本文提出的方法，隨機復特征值的方差能夠由隨機結構參數的統計量直接求出，那么相應的復模態結構的動力學重分析和不確定性分析也能夠很方便地完成。值得一提的是，非對稱系統也可能是虧損系統，換句話說，虧損系統中沒有完備的特征向量系。在這種情況下，應用于非虧損系統的矩陣攝動方法將不適用于虧損系統，因此本文研究只針對非虧損系統。

1 復模態結構系統振動的基本方程

N維線性系統的復模態受迫振動方程可以表示為：

(1)

復模態結構系統的自由振動方程可以表示為：

(2)

不妨假設X=xeλt，將其代入式(2)，則獲得相應的右特征值問題如下：

(Mλ2+Cλ+K)x=0

(3)

相應的伴隨特征值問題為(Mλ2+Cλ+K)Ty=0，對伴隨特征值方程進行轉置，得到左特征值方程如下：

yT(Mλ2+Cλ+K)=0

(4)

其中向量x和y分別表示右特征值和左特征值，λ表示振動方程的特征值。

這樣，式(3)和式(4)可以被分別改寫為：

(Aλ+B)u=0

(5)

vT(Aλ+B)=0

(6)

其中

(7)

經過狀態變換之后，式(2)可以改寫為

(Aλ+B)u=0

(8)

以及

vT(Aλ+B)=0

(9)

式(8)和式(9)具有相同的特征值，均可以由以下方程求出：

det(Aλ+B)=0

(10)

式(10)是一個在復數域中具有2N個特征值的代數方程，這些特征值可以被表示為λi(i=1,2,…,2N)?？紤]其中第i個特征值λi，其左狀態特征向量vi和右狀態特征向量ui均滿足如下方程：

(Aλi+B)ui=0

(11)

viT(Aλi+B)=0

(12)

同時由特征向量的正交歸一化條件有

vjTAui=δij

(13)

vjTBui=-λiδij

(14)

其中δij為Kronecker delta函數。

將u=Tx和v=Ty代入式(13)和式(14)，得到下式：

同時注意到式(7)，進一步將上式改寫為：

yjT[M(λi+λj)+C]xi=δij

(15)

yjT[-Mλiλj+K]xi=-λiδij

(16)

式(15)和式(16)也被稱為原二階系統的模態正交化條件。

2 復模態結構振動的矩陣攝動方法

實際工程中，結構參數的變化是具體體現在質量矩陣、阻尼矩陣和剛度矩陣的變化之中的?？紤]參數在其均值周圍的小擾動，可以將質量矩陣、阻尼矩陣和剛度矩陣表示為：

(17)

其中ε是一個小參數，并且對于原始系統有ε=0；M0,C0,K0分別表示質量矩陣、阻尼矩陣和剛度矩陣的均值；εM1,εC1,εK1分別表示M0,C0,K0的一階攝動項。進一步，式(13)和式(14)中的A和B可以表示為：

(18)

需要補充的是，本節討論的是原始特征值為特征方程(10)中單根的情況。根據矩陣攝動方法，復模態結構振動系統的特征值和特征向量可以展成攝動級數的形式：

(19)

將式(19)和式(18)代入式(11)得到：

(20)

相似的，將式(19)和式(18)代入式(12)有：

(21)

展開式(20)并且忽略O(ε3)的項，比較等式左右兩邊ε同次冪項的系數得到：

(22)

對于式(21)進行同樣的操作有：

(23)

(24)

將式(24)代入式(22)得：

(25)

(26)

考慮式(13)和式(14)所表示的模態正交歸一化條件，式(26)能夠被簡化為:

(27)

(28)

式(28)可以進一步改寫為：

(29)

3 隨機復特征值的不確定性傳播分析

由于原材料性能的差異，制造加工的誤差和載荷環境的變化，實際工程結構中的質量矩陣、阻尼矩陣和剛度矩陣一般被認為是隨機結構參數。這樣一來，便導致了隨機結構中特征值和特征向量的隨機性。

本文中，矩陣攝動方法被引入到具有隨機參數的復模態結構的特征值分析之中，同時也為非對稱系統的不確定性傳播分析打下了基礎?？紤]結構參數的隨機性質，剛度矩陣K，質量矩陣M，阻尼矩陣C，復特征值λi，左復特征向量y和右復特征向量x可以被表示為如下的形式：

(30)

根據概率理論，對式中的各隨機參數求期望，可以得到下式：

(31)

(32)

對式兩邊取期望算子，有

(33)

另外，復特征值λi的方差可以表示成如下形式：

(34)

在式(34)中，使用了復數/復向量方差的定義和期望的定義來完成推導。

將式(30)代入式(11)，采用與第1節相同的方法，得到：

(35)

(36)

由式(36)可得

(37)

簡便起見，下文的推導中將移除符號的上腳標i。同時，注意到只有特征值λd、左特征向量yd和右特征向量xd為復數形式，同時各個結構參數矩陣的一階攝動項Mr,Cr,Kr均為實數矩陣。所以假設：

(38)

其中λdr,ydr,xdr和λdy,ydy,xdy分別表示λd,yd,xd的實部和虛部。

注意到E[Mr]=E[Cr]=E[Kr]=0，所以隨機結構參數Mr,Cr,Kr兩兩乘積的數學期望均為零。因此，將式(35)和式(38)代入式(37)，得到

Var(λr)=E[(λr)·(λr)H]

(39)

為表達式的簡明起見，令：

(40)

將式(40)代入式(39)可得：

(41)

并且

(42)

值得一提的是，如果特征值λ∈R，則特征值和相應的左/右特征向量的虛部均為零，亦即，λdy=0,ydy=xdy=0,λd=λdr。因此，可以得到

式(42)將退化為：

(43)

顯而易見，式(43)與隨機結構的實特征值的方差表達式相一致。

4 數值算例

4.1 算例一：兩自由度振動系統

一個二自由度的振動系統如圖 1所示，該系統滿足c=1,k=9,m=1，并且有c1=c2=c3=c。由達朗貝爾原理，得到該系統的運動微分方程

圖1 兩自由度振動系統

根據上式，右狀態特征向量u，矩陣A和矩陣B可以表示為：

由式(11)，得到：

(44)

式(44)的特征方程等價于：

(mλ2+3cλ+3k)(mλ2+cλ+k)=0

(45)

進一步，式(45)的特征值可以推導成如下的形式：

(46)

其中ξ=c/2mω，ω2=k/m?？紤]到該二自由度振動系統的特征值是兩對共軛復根，將c=1,k=9,m=1代入式(46)得到：

λ1,2=-0.5±2.958i,λ3,4=-1.5±4.975i

(47)

考慮式(39)，由本節方法計算出的復特征值方差見表1。

表1 所提方法計算出的復特征值的方差

為了驗證本文中所提出的方法，在本數值算例中，Monte-Carlo Simulation (MCS)被使用來計算復特征值的方差。此外，為了確定MCS方法中使用的樣本容量，首先測試不同樣本大小下MCS方法計算出的方差結果，見表2。

表2 確定MCS方法的樣本容量

如表 2所示，隨著樣本容量的增大，MCS方法計算的結果趨于穩定。特別的，當樣本容量大于1×106，由MCS方法計算出的復特征值的方差收斂。因此，在本算例中取MCS方法的樣本容量為1×106是可信和合理的。進一步，由本文方法計算出的復特征值方差與由MCS方法計算出的復特征值方差對比見表3。

表3 由MCS方法和本文方法計算得到的復特征值方差的比較

注意到由本文方法計算得到的方差比MCS方法計算得到的要小。分析之后，主要有兩個原因：第一，本文中所提出的方法是基于一階矩陣攝動方法建立起來的，所以與精確解肯定會有一定的偏差，在采用二階或者更高階的矩陣攝動方法之后，這一偏差將會顯著減小；第二，本文中的隨機剛度矩陣Kr，隨機阻尼矩陣Cr，隨機質量矩陣Mr被假設為獨立變量，亦即

在本方法中不考慮不同隨機參數矩陣之間的相關性。也就是說，在本文方法的計算中只計及每一個隨機參數矩陣自身的二次項，而不計及每一個隨機參數矩陣的交叉項。這一求解策略同樣導致由本方法計算得到的復特征值的方差相比MCS方法得到的結果要小。

多項式混沌展式(PCE)方法同樣在本算例中被使用，來驗證本節中所提出方法的可行性和準確性。根據本算例中不確定結構參數(m,c,k)的分布類型，選定Hermite多項式作為正交多項式基函數[7-9]。本算例中的隨機變量的維數是三，Hermite多項式的最高階次被分別設為1和2，以增強對比效果。由于使用Hermite多項式作為PCE方法中的正交基函數，如果Hermite多項式的最高階次為p，則(p+1)次Hermite多項式的根一般被取做相應的配點。對于p=1，二階Hermite多項式的根為ξ=±1，相應的配點分布如圖

2所示，對應的PCE展式可以寫成：

(48)

(49)

圖2 p=1時PCE方法對應的配點分布

確定完變量空間的配點之后，利用基于最小二乘法的回歸分析來求解上述Hermite多項式的待定系數。結合本算例中隨機結構參數m,c,k的分布類型，得出了p=1和p=2時的復特征值λ1、λ2、λ3、λ4的多項式混沌展式見表4。進一步，基于多項式混沌展式，復特征值λ1、λ2、λ3、λ4的方差也能被直接求出，見表5。

圖3 p=2時PCE方法對應的配點分布

表4 算例一中的多項式混沌展式(變異系數=0.05)

續表4

表5 本文方法和PCE方法計算復特征值方差的結果比較

在上述算例研究中，將隨機結構參數的變異系數設定為一個固定的值(CV=0.05)。為了不失一般性，在以下的研究中，將考慮隨機結構參數取不同大小的變異系數，來進一步研究本節方法的可行性和適用性。

對于本算例中所有的結構參數，不妨認為它們的變異系數同步從0.01變化到0.1。同時，注意到本數值算例中的特征值是成對出現的共軛復根。簡明起見，以下只討論Var(λ1)和Var(λ3)，對于Var(λ2)和Var(λ4)來說，由于共軛復根的關系，很顯然其結果與Var(λ1)和Var(λ3)是相似的。考慮隨機結構參數取不同的變異系數，同時比較三種不同的不確定性分析方法(本節所提方法、MCS方法和PCE方法)所得的Var(λ1)和Var(λ3)，其結果如圖 4和圖 5所示。

圖4 不同變異系數下三種方法計算出的Var(λ1)結果比較

圖5 不同變異系數下三種方法計算出的Var(λ3)結果比較

盡管變異系數在增加(從0.01到0.1)，隨機結構參數的不確定性在增大，但由本文提出方法計算出來的結果與MCS方法和PCE方法計算出的結果在總體趨勢上保持一致。即使當變異系數增大到0.1時，本文方法計算出的結果與MCS方法和PCE方法之間的偏差仍然是可以接受的，這進一步驗證了所提方法的可行性和適用性。另外，發現由本節方法計算得到的Var(λ1)和Var(λ3)相比MCS方法和PCE方法的結果偏小，這與表 3和表 5的結果對比一致，相關原因已經在上文中討論過。

比較所提方法、PCE方法和MCS方法所需的計算時間，本文所提方法在計算效率方面具有明顯的優勢，如圖6所示。

圖6 三種方法的計算耗時比較

隨著樣本大小的增大，MCS方法的計算耗時顯著增大。因此當處理某些大型復雜結構時，MCS方法的計算耗時將是無法承受的。對于PCE方法，盡管計算耗時相對于MCS方法已經減少很多，但是由于該方法需要在各個配點處進行計算，所以在計算耗時上相比于本文所提方法仍然大了幾乎一個數量級。反觀本文所提方法，其計算成本相對較低，這不僅提高了結構不確定性分析和結構重分析的計算效率，并且能夠縮短具有復模態結構的設計優化周期。

為了更加直觀地比較三種方法的計算耗時，形成圖 6，這樣本文方法在計算效率方面的優勢被更直接和充分地展現。需要指出的是，本文所提出的方法、PCE方法和MCS方法都是在同一臺Intel Core i7-2600@3.40GHz電腦上被執行的。

4.2 算例二：五自由度彈簧質量塊系統

如圖7所示，一個五自由度彈簧質量塊力學系統被作為說明算例，來進一步驗證本文所提出的方法。同時，假設振動只在豎直平面內發生。

圖7 五自由度彈簧質量塊系統

該系統中質量矩陣M的元素可以表示為：

(50)

其中Ji(i=4,5)表示慣性矩，L為正方形面板的邊長，θi(i=4,5)表示轉動角。

該系統中剛度矩陣K的各個元素可以寫成：

(51)

而阻尼矩陣C的各元素為：

(52)

假設本算例中所有的結構不確定參數均滿足高斯隨機分布，同時有k3=k4=k5=k6=kG，c3=c4=c5=c6=cG。并且，令k1,k2和kG分別滿足高斯分布N[15 000,(0.05×15 000)2],N[30 000,(0.05×30 000)2]和N[1 000,(0.05×1 000)2]；令c1,c2和cG分別滿足高斯分布N[6,(0.05×6)2],N(9,(0.05×9)2]和N[40,(0.05×40)2]；令m1,m2,m3和L分別滿足高斯分布N[300,(0.05×300)2],N[750,(0.05×750)2],N[1 500,(0.05×1 500)2]和N[1,(0.05×1)2]。

首先，將各個結構參數的均值取為它們的名義值，則該振動系統的名義參數矩陣可以表示如下：

(53)

進一步得到振動系統的名義復特征值如下(共5對共軛復根，其中有兩對在數值上相同)：

(54)

考慮各個結構參數的不確定性，根據本節所提出的計算方法，得到了該振動系統復特征值的方差(表6)。簡明起見，每一對不同的共軛特征值中，僅僅取一個特征值進行計算。

表6 由MCS方法和本文方法計算得到的復特征值方差的比較 (變異系數=0.05)

如表6所示，這一部分的計算是將算例中所有10個不確定參數的變異系數均取為0.05。作為對比算法，經計算，本算例中蒙特卡洛方法的樣本容量取為106。通過比較兩種方法的計算結果及其相對誤差發現，本節所提出方法的計算結果完全能夠滿足工程應用的要求。

在前面的討論中，研究了在不確定結構參數取固定變異系數的情況下，本節所提出方法的計算精度。為了進一步研究所提出方法的計算能力和計算精度，在下面的討論中，考慮不確定結構參數取變化的變異系數，同時假設所有的不確定結構參數的變異系數同步變化。簡明起見，有針對性地取Var(λ1)進行研究，對于Var(λ3)、Var(λ5)和Var(λ7)，其討論過程是相似的。

圖8比較了由所提出方法和MCS方法計算得到的復特征值的方差[i.e.Var(λ1)]?？梢钥闯?，隨著變異系數的增大，由所提出方法計算得到的復特征值的方差也在逐漸增大；由MCS方法計算得到的結果與之具有相同的趨勢，但在計算結果的數值上偏大一些。當不確定結構參數的變異系數取較小值時，兩種方法所得結果非常接近，誤差很小；當不確定結構參數的變異系數取值接近0.1時，兩種方法計算結果的偏差逐漸增大。值得一提的是，當變異系數取0.1時，所提方法計算結果與MCS方法計算結果之間的相對誤差僅為7.5%，在工程上依然是比較令人滿意的。這也驗證了本文所提方法的計算精度。

圖8 由兩種方法計算得到的Var(λ1)隨變異系數的變化比較

通過對以上兩個數值算例的分析和研究，對于具有隨機結構參數的復模態結構系統，其不確定性傳播分析的重要意義和必要性被充分地認識，結構參數的不確定性及其產生的影響也被深入地研究。不僅如此，數值算例表明，本文所提方法能夠滿足工程實際對于算法可行性和算法適用性的要求。

5 結論

本文將實模態結構隨機攝動方法改進并推廣到復模態特征值的不確定分析領域，改進的復模態隨機攝動方法與概率理論相結合，建立了隨機復特征值的方差關于隨機結構參數的高效算法，開展了隨機結構復特征值的不確定傳播分析。

通過對多個數值算例的研究，本文所提出方法的可行性、適用性和計算精度得到驗證。傳統的矩陣攝動方法，為了進行結構特征值的不確定傳播分析，需要預先假設或者通過大量的樣本統計試驗獲得結構參數的相關系數矩陣，在工程上往往難以實現，所提出的改進的隨機攝動方法不涉及結構參數的相關系數矩陣計算，因此克服了結構參數相關系數矩陣對于采用攝動方法進行隨機特征值不確定性分析的制約。所提出的方法在工程上易于實現，并且使得計算隨機結構特征值的統計學性質變得更為便利，為相應的動力學重分析和結構設計優化打下良好的基礎。