基于多復域的頻響函數靈敏度分析

田 宇， 曹芝腑， 姜 東，

（1.南京林業大學 機械電子工程學院，南京 210037；2.東南大學 空天機械動力學研究所，南京 211189）

參數靈敏度分析在工程結構優化[1]、模型修正[2-4]和參數識別[5-8]中具有廣泛應用。通過參數的靈敏度分析，可以量化工程結構中各設計參數對結構性能的影響程度和影響規律，指導結構設計和分析[9-10]。復雜結構的影響因素比較多，通過靈敏度分析可以將對結構性能影響較大的參數作為變量，從而降低結構分析設計難度，提升分析效率。

直接求導法通過對頻響函數的模態展開式進行求導，得到頻響函數靈敏度微分方程，該方法計算量較大，且存在由模態截斷導致的頻響函數分析誤差。頻響函數殘差法是在實數域上進行參數攝動，在函數偏導數急劇震蕩處或攝動量過小時，計算精度容易受攝動步長的影響。Lin[11]運用加權靈敏度法，在建模存在較大誤差的情況下，依然可以得到收斂性好的頻響函數靈敏度分析結果。Zhu等[12]基于Sherman-Morrison-Woodbury公式，提出了一種頻響函數的再分析方法，利用該方法可以得到多參數同時攝動時的頻響函數靈敏度，提高分析效率。黃修長等[13]基于頻響函數綜合的子結構方法推導了結構動力系統的靈敏度，提高了結構建模與優化設計的效率。周俊賢等[14]基于頻響函數、模態、靜力位移的靈敏度分析，對結構的局部損傷進行識別。Esfandiari等[15]提出了一種基于頻率響應函數主成分變化和靈敏度分析的模型修正方法，用不完全測量的結構響應結合主成分分析得到靈敏度，該靈敏度分析方法比直接使用頻響數據的結果更準確。Wang等[16]提出了一種基于頻響函數靈敏度分析的優化方法，研究了頻響函數在一定頻帶寬度約束下的桁架結構的動態尺寸優化問題。

復變量求導法通過提取復數域的虛部響應得到參數靈敏度，避免差分運算，可以得到高精度的一階導數。Wang等[17]將復變量求導方法應用于彈簧質量系統的固有頻率和振型的一階導數的求解。Gao等[18]通過求解奇異函數和強非線性函數的偏導數，證明了復變量求導法得到的導數解精度遠高于差分法。Garza等[19]運用復變量求導方法計算結構動響應的一階和二階靈敏度，與前向差分法和中心差分法的靈敏度分析結果相比較，說明復變量求導法不易受攝動步長影響。在頻響函數靈敏度分析中，運用復變量求導法計算頻響函數靈敏度更加高效精確。但由于頻響函數本身含有復數，對設計變量再進行虛部攝動時，需要對方程本身的復域與攝動量的復域進行解耦，產生了不同復域下的虛數單位同時運算的問題。

本文提出了一種基于多復域的頻響函數靈敏度分析方法，建立復數在實域中的等價矩陣格式，得到多復域下結構動力學方程與靈敏度方程的解耦形式，實現頻響函數關于結構參數的靈敏度求解，為頻域內結構動響應的靈敏度分析提供更加準確的數值求解方法。

1 理論基礎

1.1 復變量求導法

復變量求導法是一種計算精度高且不受攝動步長影響的數值靈敏度求解方法，復數的定義為

式中:C1為所有單復數的集合；i1為虛數單位，=-1；R為所有實數的集合。在復變量求導法中，通過構造設計參數的虛部攝動量，利用復數域泰勒展開公式，可以得到其對應的一階靈敏度

式中:θ為設計參數；h=Δh×θ為參數攝動量，Δh為攝動系數。

則式（2）的虛部響應可以表示為

式中，Im（f）為函數f的虛部響應，最后得到θ的一階偏導數可表示為

式中，O（h2）為攝動量的二階截斷。由式（4）可以看出，計算函數的一階靈敏度時，對實函數的自變量攝動i h，求出復變函數的值，取其虛部再除以h，即可得到函數的一階靈敏度，避免了極小數相減的問題。該過程可以解決一般解析分析無法解決的強非線性和隱函數偏導數計算問題，且一步即可完成計算，計算量比常規差分法更小。計算系統頻響函數靈敏度時依然可以使用此方法。

1.2 復域的矩陣等價

在采用復數進行分析時，通常采用式（1）的形式構造計算方程。為了將復數的實部和虛部進行分離，可以采用實數域的矩陣形式得到復數域的等價格式。對同一個運算過程，矩陣形式的計算結果重新化為多項式形式后，與使用式（1）的多項式形式計算出的結果完全一致[20]。任何單復數都可以用一個2×2的實數矩陣表示。定義如下

式中，X1為單復數x0+x1i對應的矩陣表達形式，其每個元素與C1中的每個元素一一對應。

雙復數的定義在單復數的基礎上，含有兩個不同方向的虛數單位。在所有的多復數中，雙復數的性質是已知和研究最多的，并已經應用于多個應用領域，如分形和量子理論[21-22]。雙復數可以用2×2的復數矩陣或4×4的實數矩陣表示，定義如下

式中，X2中的每個元素與C2中的每個元素一一對應，任何雙復數都可以用一個2×2的復數矩陣或4×4的實數矩陣表示。

每增加一個虛數單位，就要對矩陣進行一次擴維。含有n個虛數單位的參數，需要2n×2n維實數矩陣進行等價。

表示多復數的實數矩陣在文獻[23]中稱為柯西-黎曼矩陣，這種形式的柯西-黎曼矩陣與多復數之間存在一一對應關系。因此，多復數的算術運算（+，-，×）等價于矩陣表示上的算術運算。

1.3 多復域頻響函數靈敏度分析

為了求解結構頻響函數關于設計參數的靈敏度，考慮頻域內結構的動力學方程中已有的i復域，引入僅考慮設計參數攝動的j復域，從而實現多復域下的頻響函數靈敏度分析。

頻域內結構的動力學方程為

式中:M，C，K分別為系統的質量矩陣、阻尼矩陣和剛度矩陣；i和j為虛數單位；ω為圓頻率；θ為設計參數；x（ω）和F（ω）為位移與激勵的傅里葉變換。

首先，對i復域內的結構矩陣進行等價轉換。利用復數的等價矩陣表達，式（7）中各項的等價實數矩陣為

式中:xRe，xIm分別為x（ω）的實部和虛部；FRe，FIm分別為F（ω）的實部和虛部。 將式（8） ～式（12）代入式（7）可得

式（13）為四個等式構成的方程組，其中只有兩個有效方程，舍去x矩陣與F矩陣的第二列，式（13）可簡化為

然后，考慮設計參數θ所在的j復域，再次利用復數的矩陣表達對j復域內的參數進行等價轉換。將式（14）展開為j復域下的表達式

式中，ΔM，ΔC和ΔK分別為M，C，K的攝動量。利用雙復數式（6），將式（15）轉化為實數矩陣的表達

最后，通過在i復域和j復域上的二次轉換，得到擴維后的控制方程式（22），即為頻響函數對設計參數θ的靈敏度計算公式。

若對某設計參數的攝動只導致質量矩陣、阻尼矩陣和剛度矩陣中的某一個或兩個矩陣發生變化，則只要將不發生變化的攝動矩陣ΔM，ΔC或ΔK變為零矩陣代入式（22）即可。

1.4 實現流程

基于多復域的結構頻響函數靈敏度分析方法（見圖1），按照如下步驟實現:

圖1 多復域法頻響函數靈敏度分析流程圖Fig.1 Flow chart of sensitivity analysis of frequency response function in multicomplex domain

步驟1對結構進行動力學建模，得到結構的質量矩陣、阻尼矩陣、剛度矩陣和外部力向量，從而得到結構的頻響函數方程；

步驟2構造設計參數的虛數攝動量，該虛數方向與頻響函數自身虛數的方向不同，使設計參數從實數變為復數，結構質量矩陣、阻尼矩陣、剛度矩陣也變成對應復數量，頻響函數函數方程變為多復數方程；

步驟3將攝動后的頻響函數方程運用多復數理論進行擴維，得到等價的實數矩陣，并利用數值分析方法對該實數矩陣求解；

步驟4提取實數矩陣的運算結果，得到結構的實部響應和虛部響應，并可計算得到頻響函數對設計參數的的實部靈敏度和虛部靈敏度。

2 算例研究

2.1 多自由度系統

為驗證多復域方法靈敏度分析方法的有效性，用五自由度彈簧質量阻尼系統開展研究，結構如圖2所示。

圖2 五自由度彈簧質量阻尼系統Fig.2 Five degrees of freedom spring mass damping system

該系統兩端固定，共有5個質量參數，9個剛度參數和6個阻尼參數，分析模型選取的參數值為m1=5 kg，m2=2 kg，m3=1 kg，m4=2 kg，m5=5 kg，c1=c2=c3=c4=c5=c6=0.5 N·s/m，k1=k2=k5=k6=100 000 N/m，k3=k4=50 000 N/m，k13=k24=k35=40 000 N/m。計算結構的幅頻曲線，如圖3所示。

圖3 五自由度系統幅頻曲線Fig.3 Amplitude-frequency curve of a five degree of freedom system

2.1.1 靈敏度分析

利用多復域法計算彈簧質量阻尼系統頻響函數對結構參數的靈敏度。為了便于比較兩種方法，選取合適的分析頻段，先計算出系統的幅頻曲線見圖3。

由系統幅頻曲線可知，該系統的一階模態頻率在110 Hz附近，為方便討論，取頻率范圍108～112.5 Hz進行研究，分別給出系統頻響函數幅值和頻響函數對特定質量參數、阻尼參數、剛度參數的靈敏度。

記結構頻響函數為Hlp，表示在第p個自由度上施加單位簡諧力時第l個自由度的穩態位移頻率響應。使用有限差分法和多復域法計算結構頻響函數H11的幅值及實部、虛部響應對不同參數的靈敏度，攝動系數設為10-6。以頻響函數的實部靈敏度為橫坐標，虛部靈敏度為縱坐標，畫出靈敏度的奈奎斯特圖。結構頻響函數幅值及頻響函數的實部和虛部分別對質量參數、阻尼參數、剛度參數的靈敏度如圖4所示。

圖4中Am表示頻響函數的幅值，圖4（a）、圖4（c）、圖4（e）分別為結構頻響函數幅值對質量參數m4、阻尼參數c1和剛度參數k1的靈敏度曲線，圖4（b）、圖4（d）、圖4（f）為位移頻率響應H11分別對質量參數m4、阻尼參數c1和剛度參數k1的實部和虛部靈敏度的奈奎斯特圖。奈奎斯特圖上每一點都是對應一個特定頻率點下的頻響函數靈敏度，橫坐標和縱坐標分別為頻率響應的實部靈敏度和虛部靈敏度。頻響函數靈敏度更大的點計算誤差也越大，更需要被關注，而頻率-靈敏度圖中絕大部分都為數值較小的點，誤差大的部分難以被重點關注。奈奎斯特圖中數值較小的點都集中在原點附近，數值越大的點距原點越遠，變化較大的靈敏度在圖中得到了很好的反映，精確度和誤差也因此更容易比較。通過選取特定點計算響應對質量、阻尼和剛度的靈敏度，發現多復域方法是一種有效的方法，可以得出較為精確的靈敏度曲線，但精確度需要進一步探究。

圖4 兩種不同方法計算得到的頻響函數靈敏度曲線Fig.4 Sensitivities of the FRF calculated by two different methods

2.1.2 攝動系數影響分析

為給出多復域法計算結果的精確誤差，計算有限差分法和多復域法不同方法在攝動系數從10-1～10-7時與精確值之間的誤差，以解析法的計算結果為精確值，給出頻響函數靈敏度精確值的計算表達式（23）。

定義每次攝動后的靈敏度分析結果誤差e為

表1 兩種分析方法對于精確值的誤差Tab.1 The error of the two analytical methods for the exact value

可以看出，多復域方法的誤差明顯小于有限差分法，具有更高的精度。畫出兩種方法在不同攝動系數時的靈敏度奈奎斯特圖，如圖5所示，通過圖片更加直觀地反應兩種方法的誤差。

圖5 攝動系數對頻響函數靈敏度的影響Fig.5 Influence of perturbation of two different methods on the sensitivities of the FRF

在攝動系數變小的過程中，有限差分法得出的剛度靈敏度曲線逐漸歸一，誤差逐漸減小。但在攝動系數過小之后，誤差又開始逐漸增大。運用多復域方法計算剛度靈敏度在攝動系數到達10-2后，曲線迅速平穩，可以看出多復域方法收斂快，精度高。

2.2 桁架結構

進一步說明多復域靈敏度分析方法在復雜結構中的適用性，本文以GARTEUR AG-11[24]桁架結構為例開展頻響函數靈敏度分析。

圖6所示GARTEUR AG-11桁架結構的幾何尺寸為寬3 m，長3 m×5 m，采用桿單元建立該結構的有限元模型，共劃分為36個單元和32個節點，每個節點有兩個自由度Zix和Ziy，節點編號與自由度編號如圖所示。其中31和32號節點為邊界節點，約束其全部自由度。結構的材料分別參數為:彈性模量E=7.5 GPa，密度ρ=2 800 kg/m3，橫截面積A=0.004 m2。結構阻尼矩陣采用比例阻尼，設C=αM+βK，利用加權最小二乘法計算得出合理的α和β，取α=β=0.000 1。

圖6 GARTEUR AG-11桁架結構Fig.6 GARTEUR AG-11 truss structure

計算桁架結構在10號節點x向自由度Z10x的原點頻響函數，記為H，其幅值曲線如圖7所示。

圖7 桁架結構頻響函數幅值Fig.7 Amplitude of FRF of truss structure

2.2.1 靈敏度分析

桁架橫截面積A是桁架結構的重要參數，決定桁架的剛度和質量。分析桁架頻響函數對橫截面積的靈敏度，可以為結構的優化設計和模型修正提供參考。桁架結構的橫截面積是較為復雜的參數，橫截面積的變化會引起系統質量矩陣、阻尼矩陣、剛度矩陣的變化。

在用公式表達彈簧質量阻尼系統的靈敏度分析的基礎上，給出適用于更一般連續系統的方法。考慮連續結構如圖8所示的簡單桿桁架結構系統的質量、阻尼和剛度攝動。

圖8 桁架單元Fig.8 Truss element

圖8為桁架結構的一個單位元素，單元兩側的節點用i和j表示，每個節點兩個自由度，分別為x向和y向的平動自由度。則集中質量矩陣表示為

式中:上標t為單元編號；ρ為密度；A為橫截面積；l為長度。單元剛度矩陣表示為

式中:E為彈性模量，每個子陣[kij]t維數為2×2。給定傾斜角度θ，則轉換矩陣為

當總質量擾動為Δmi時，在節點i處的質量單元矩陣的單元局部坐標的變化量為

在全局坐標系中，質量單元矩陣的相應變化為

上面的形式可以寫為

式中，上標T為矩陣轉置，其中

由于質量矩陣組裝性質，描述總質量變化的矩陣ΔM為

式中，Yij中的元素按節點編號排列，只有第（2i-1）個和第（2i-2）個對角元素非零。

如果將多個集中質量添加到多個節點，則由于集中質量只修改質量矩陣對角線上的元素，頻響函數的靈敏度可以由此計算得出。類似式（32）要求在矩陣ΔM的對應位置添加變量值。

假設長度l保持不變，橫截面積A作為變量，通過式（25）和式（26）影響質量矩陣和剛度矩陣。我們取單元t為例考慮橫截面積變化情況:A^=A+ΔA。

單元t的橫截面積變化ΔA時，其質量矩陣和剛度矩陣均發生變化。為了簡化，設ΔA＞0，表示質量單元和剛度單元變化的矩陣為

則同理可知，將ΔA代入有限元分析計算，得出的質量矩陣M、剛度矩陣K即為所需的ΔM，ΔK。

在得到質量、阻尼、剛度矩陣的變化量ΔM，ΔC，ΔK后，代入式（22），采用多復域方法，得到Z10x的原點頻響函數H關于橫截面積A的實部和虛部響應靈敏度奈奎斯特圖，如圖9所示。從圖中可以看出，本文方法得到的結果與有限差分法的結果一致，對于復雜桁架結構，橫截面積A的變化同時引起質量矩陣、阻尼矩陣和剛度矩陣的變化，該情況下多復域方法仍然有效。

圖9 攝動系數為10-5時多復域方法的靈敏度分析結果Fig.9 Sensitivity analysis results of multicomplex domain method with perturbation of 10-5

2.2.2 靈敏度影響因素

為探究多復域方法對復雜結構進行靈敏度分析時的準確性和有效性，考慮攝動量對該方法計算精度的影響。分別使用有限差分法和多復域法，將攝動系數從10-1變化到10-8，計算頻響函數H對設計參數橫截面積A的靈敏度，并用式（23）和式（24）計算出精確值和誤差e，結果如圖10所示。

圖10 攝動系數對有限差分法與多復域法計算頻響函數靈敏度的影響Fig.10 The influence of perturbation coefficient on the sensitivity of FRFs by finite difference method and multicomplex domain method

多復域法和有限差分法在攝動量變小的過程中，誤差都逐漸收斂。但有限差分法計算結果在攝動系數到達10-6后開始產生波動，在10-8處產生了劇烈變化，精度明顯下降，而與之相對應的多復域法仍然穩定。

3 結 論

本文提出了一種多復域的結構頻響函數靈敏度分析方法，通過實數域的矩陣格式對復數的實部和虛部進行分離，實現了復變量求導法在結構頻域響應靈敏度分析中的應用，并通過多自由度系統和GARTEUR結構進行參數研究，證明了該方法的有效性。

相較于傳統的實數域有限差分法，本文給出的一種多復域靈敏度分析方法保留了復數域靈敏度分析方法的優點，計算精度高，不受攝動步長影響，可以為頻響函數靈敏度分析提供更加準確的結果。本文僅針對方程中存在兩個復域的情況進行了研究，若存在更多復域時，可對方程進行多次擴維，用以一次攝動多個參數或求更高階靈敏度，對此該方法依然有效。