薛睿淵，俞樹榮，張希恒

(蘭州理工大學 石油化工學院，甘肅 蘭州 730050)

隨著人類對地震認識的加深，對核電站中的系統(tǒng)和設備提出了新的抗震要求，主要包括對現(xiàn)役和在建核電站進行抗震裕度分析[1-2]和對關鍵設備在高頻激勵下破壞形式的研究[3]。對于核電站中被劃分為抗震重要等級的管道而言，一個準確的動力學有限元模型(FEM)可解決其在服役過程中無法進行試驗的問題，對于其健康檢測、損傷診斷也是必不可少的。由于管道在加工制造和安裝過程中的不確定性，其FEM必然存在一些參數(shù)的取值難以準確確定，如支撐為系統(tǒng)提供的約束強度、管道連接件的剛度等。對這些參數(shù)的不準確取值會導致實際結構與FEM之間存在偏差。在工程實踐中，利用基于模態(tài)的FEM修正方法對機械結構模型中的結構參數(shù)進行反演分析是一種普遍方法。郭勤濤等[4]總結了FEM修正方法的發(fā)展、應用和技術難點。費慶國等[5]指出模態(tài)頻率和反共振頻率是FEM修正過程中優(yōu)選的兩種參數(shù)。為減少利用模態(tài)進行反演分析的計算量，近年來代理模型法得到了廣泛運用[6-7]。但對于核電站中完整的管道而言，準確提取足夠數(shù)量的模態(tài)存在很大難度，因此不能直接使用基于模態(tài)的FEM修正方法來反演其模型中的結構參數(shù)。

核電站中相對容易獲得的與管道相關的振動測量數(shù)據(jù)是在其服役或模擬運行試驗過程中監(jiān)測到的地震響應。在橋梁工程領域有大量直接利用從環(huán)境振動試驗中測得的加速度響應進行橋梁壽命和損傷監(jiān)測的實例[8-10]，但本質(zhì)使用的試驗數(shù)據(jù)依然是從加速度響應中提取的模態(tài)，提取模態(tài)的數(shù)量和精度對反演過程有很大限制。文獻[11-12]報道了直接利用時域振動響應進行參數(shù)反演、損傷識別的方法，但這些方法在機械結構中的實際應用還未得到探索。Lin等[13]將結構在基礎激勵下的響應與激勵的比值定義為一種新形式的頻響函數(shù)(FRF)，然后提出基于FRF殘差的可直接使用基礎激勵下響應數(shù)據(jù)的FEM修正方法，推導過程中激勵與響應均為位移，且假設只有固定節(jié)點受力。該方法在航天器FEM修正方面得到廣泛運用[14-15]。而管道的相對位移測量難度大、精度低，且地震過程中系統(tǒng)的每個節(jié)點均受到慣性力的作用。

本文參考文獻[13]，將地震響應與地震激勵的比值定義為一種新形式的FRF，以FRF和管道振動控制方程誤差為目標特征值，推導二者可通用的基于靈敏度的FEM修正模型。然后基于兩種方法利用地震模擬試驗測量數(shù)據(jù)對某管道FEM中的結構參數(shù)進行反演分析，對比二者在實踐運用中的可行性與差異。

1 修正數(shù)學模型推導

基于靈敏度的FEM修正過程的第1步是建立理論和試驗特征值之間的偏差與修正參數(shù)之間的數(shù)學關系，然后通過迭代求解對修正參數(shù)進行識別。這種關系可表示為：

ε=b-Adθ

(1)

式中：ε為FEM的建模誤差；b為理論與試驗特征值之間的偏差；dθ為修正參數(shù)的變化向量；A為修正參數(shù)取值發(fā)生波動時對目標特征值的影響程度，稱之為修正參數(shù)對目標特征值的靈敏度矩陣。

構造dθ關于ε的罰函數(shù)，并對修正參數(shù)進行加權同時引入?yún)?shù)的初始估計信息，可得：

J(dθ)=‖ε‖2=(b-Adθ)T(b-Adθ)+

(θ-θ0)TW(θ-θ0)

(2)

式中：J為ε的罰函數(shù)；θ、θ0和W分別為修正參數(shù)的當前估計向量、初始估計向量和加權矩陣。為向變化較大的參數(shù)分配較小的權重，W可取參數(shù)估計值協(xié)方差矩陣的逆。求解J關于dθ最小值的方法是計算J關于dθ的偏微分并令其等于0，可得：

dθ=[ATA+W]-1(ATb-W(θ-θ0))

(3)

式(3)中當修正參數(shù)性質(zhì)不同時，有必要對參數(shù)進行歸一化處理以改善修正過程的數(shù)值條件。此外，當參與修正的頻點計算得到的[ATA]接近奇異時，式(3)難以求得符合實際的解，為提高修正過程的成功率，此處利用Marquardt方法求解式(3)，即將一個正定的對角矩陣與[ATA]相加。規(guī)定κi+1=dθi/θi，i為當前迭代步數(shù)，數(shù)值條件改善后，用于反演分析的迭代數(shù)學模型如下：

κi+1=[(Aθi)T(Aθi)+W+λI]-1·

((Aθi)Tb-W(κi-κ0))

(4)

式中：I為單位矩陣；λ為Marquardt方法中需要設置的常數(shù)。修正過程λ的取值需根據(jù)每一步的迭代結果進行調(diào)整，在保證修正過程順利完成的前提下盡可能地提高達到收斂的速度。相對于以往修正方法中使用最小二乘法對迭代公式進行求解，這樣雖增加了計算量，但大幅度提高了順利完成修正過程的可能性。

將地震響應與地震激勵的比值定義為一種新形式的FRF，修正過程以這種FRF的殘差為目標特征值時，b和A的數(shù)學模型如下：

b=-B(Ht-Ha)

(5)

(6)

式中：B、M和e為修正模型的動剛度矩陣、質(zhì)量矩陣和單位向量；Ht和Ha分別為實測與理論的FRF向量，可按下式進行計算：

(7)

Ha=-B-1Me

(8)

當以振動控制方程誤差為目標特征值時，b和A具有以下形式：

(9)

(10)

對于上述兩種方法而言，當p為加速度響應時：

B=M-ω-2K

(11)

當p為位移響應時：

B=K-ω2M

(12)

式中：K為修正模型的剛度矩陣；ω為用于反演分析的頻率向量。

綜上利用地震響應對核電管道模型中的結構參數(shù)進行反演分析的流程如圖1所示。

可以觀察到上述兩種FEM修正方法在推導過程中均不涉及任何模態(tài)信息，僅需測量實際結構的地震響應即可完成FEM中結構參數(shù)的反演分析。此外，上述修正過程未考慮阻尼，這是因為目前流行的兩階段FEM修正方法[16]已證明阻尼可在建立起準確的動剛度矩陣后通過不同方式進行單獨識別，因此本文只對比這兩種方法在反演和動剛度矩陣相關參數(shù)過程中的差。

2 試驗數(shù)據(jù)采集與管道模型建立

2.1 地震模擬試驗

為驗證上述兩種FEM修正方法在實踐運用中的可行性和差異，搭建如圖2所示的管道結構并對其進行地震模擬試驗。圖2所示管道結構總長3.5 m，由兩段重3.78 kg/m、φ48 mm×3.5 mm規(guī)格的鍍鋅焊接管通過管接頭連接而組成，管接頭附近焊接一重為5 kg的法蘭盤作為負重。管道與支架通過螺母連接，支架底板通過8個M20螺栓固定在振動臺臺面。采用電液伺服地震模擬振動臺為系統(tǒng)提供水平y(tǒng)方向的人工地震波激勵，12個加速度計(A1～A12)吸附在管道上測量試驗過程中的加速度響應，同時在地震臺臺面上吸附1個加速度計(A13)監(jiān)測臺面振動加速度時程。A13測得的人工地震波如圖3所示，測點A1～A12測量數(shù)據(jù)減去A13測量數(shù)據(jù)即可得到修正過程所需的相對加速度響應。試驗管道在人工地震波激勵下僅表現(xiàn)出第1階固有頻率，為17.33 Hz。

圖1 結構參數(shù)反演分析流程圖Fig.1 Flow chart of structure parameter inversion analysis

圖2 試驗結構及測點布置Fig.2 Test structure and arrangement of test point

圖3 人工地震波Fig.3 Artificial seismic wave

2.2 建立管道FEM

用無質(zhì)量的管單元模擬管件，支撐以線性彈簧模擬，建立試驗結構的集中質(zhì)量FEM如圖4所示，由15個單元、16個節(jié)點構成，模型中節(jié)點與試驗結構測點一一對應。由于試驗管道質(zhì)量較小，分配于各節(jié)點的質(zhì)量只能在小范圍內(nèi)變化，而這種幅度的變化對管道的動態(tài)特性幾乎沒有影響，因此可認為FEM中質(zhì)量矩陣是近似確定的。管道左右兩端支架如圖5a、b所示，由30 mm厚的鋼板焊接而成，具有足夠的剛度將振動臺提供的激勵傳遞至管道，但考慮到支架與管道連接部位在制造和安裝過程的不確定性，支架對系統(tǒng)提供的約束強度是隨機的，即FEM中線性彈簧剛度的取值是未知的。此外系統(tǒng)左右兩端管道外伸出支架的長度是不一致的。單元8對應的管接頭位置如圖5c所示，管接頭制造安裝過程中的不確定性、管接頭與管件之間螺紋連接松動及外徑不均勻等因素造成單元8的彈性模量和慣性矩取值也是未知的。綜上，圖4中取值未知的結構參數(shù)包括左右兩端線性彈簧的平動剛度Kt1、Kt2，轉動剛度Kr1、Kr2及單元8的抗彎剛度W。W為單元8彈性模量與慣性矩的乘積。

為反演分析過程的順利完成，需從上述參數(shù)中進一步挑選修正參數(shù)，且還需從測得的地震響應中挑選用于反演分析的頻點，因此需對上述參數(shù)進行靈敏度分析。各參數(shù)對FRF殘差的靈敏度如圖6所示，對振動控制方程誤差的靈敏度如圖7所示。

圖4 試驗結構的集中質(zhì)量有限元模型Fig.4 Lumped mass FEM of test structure

a——左端支架；b——右端支架；c——管接頭圖5 試驗結構中參數(shù)不確定的部件Fig.5 Component with uncertain parameter in test structure

a——Kt1；b——Kt2；c——Kr1；d——Kr2；e——W圖6 結構參數(shù)對FRF殘差的靈敏度Fig.6 Sensitivity of structure parameter to FRF residual

a——Kt1；b——Kt2；c——Kr1；d——Kr2；e——W圖7 結構參數(shù)對振動控制方程誤差的靈敏度Fig.7 Sensitivity of structure parameter to vibration control equation error

由圖6、7可知，上述各參數(shù)對兩種特征值的靈敏度曲線有類似的趨勢，各參數(shù)的最大靈敏度均出現(xiàn)在共振頻率附近，靈敏度在5 Hz以下也具有較高取值，意味著5 Hz以下測量數(shù)據(jù)受噪聲影響較大，靈敏度取值較低的頻率范圍為5～13 Hz。不同的是各參數(shù)對FRF殘差的靈敏度是對振動控制方程誤差靈敏度的100倍。值得注意的是，彈簧的平動剛度對兩種特征值的靈敏度均遠小于其他參數(shù)，參考文獻[17]的處理方式，反演過程中應將Kt1、Kt2的取值設置為無限大。綜上選擇模型中的修正參數(shù)為Kr1、Kr2和W。

3 參數(shù)的反演分析

3.1 反演條件準備

評價一個FEM修正方法性能的兩個重要指標是方法的魯棒性和對修正參數(shù)初始估計的依賴性。Kwon等[18]的研究成果表明用于反演分析的頻點應在靈敏度較低的頻率區(qū)域內(nèi)選擇。參與修正的頻點避開共振頻率的另一重要意義在于可最小化忽略阻尼對反演結果的影響，因為阻尼對共振區(qū)域的響應幅值最為敏感。本次試驗測得的響應數(shù)據(jù)在頻域受噪聲影響隨頻率的升高而降低，因此選取表1所列4組依次降低的頻率組合用于參數(shù)的反演分析，以對比兩種方法的魯棒性。

表1 參與反演過程的頻率組Table 1 Frequency group used for inversion process

反演分析前還需對修正參數(shù)的允許變化范圍及初始取值進行估計。由于沒有任何工程經(jīng)驗可指導修正參數(shù)的取值，應允許各參數(shù)在較大范圍內(nèi)變化。修正參數(shù)的初始估計列于表2。初始估計Ⅰ是考慮到管道支架由鋼板焊接組成，管接頭的存在使管件壁厚增大，因此判斷彈簧的轉動剛度和管接頭部位的抗彎剛度應大于管單元的轉動剛度(5.9×105N·m/rad)和抗彎剛度(2.86×104N·m6)，其對應的理論基頻為19.6 Hz。初始估計Ⅱ則是以理論和實測固有頻率吻合為目標，經(jīng)多次調(diào)整而確定的較為準確的初始估計，對應理論基頻為17.7 Hz。

表2 修正參數(shù)的初始估計Table 2 Initial estimation of updating parameter

由于試驗過程中管道轉動自由度方向沒有外界激勵，故可利用靜力凝聚法消去FEM中的轉動自由度以解決試驗過程各測點轉動自由度響應未被測量的問題。按照圖 1所示流程編制程序實現(xiàn)修正參數(shù)的迭代求解。

3.2 利用加速度響應進行反演

采用初始估計Ⅰ，基于FRF殘差和振動控制方程誤差各參數(shù)反演分析結果列于表3、4。將修正后FEM對應的基頻與試驗結果的相對偏差σ小于1%作為反演結果可靠的一個必要條件。由表3、4可知，當對修正參數(shù)的初始估計不準確時，只有使用測量精度較高的地震響應基于振動控制方程誤差的方法才能獲得可靠的反演結果，即工況4-1(表4中的工況1)。這證明基于振動控制方程誤差的方法對修正參數(shù)的初始估計依賴性低。

當采用初始估計Ⅱ時，基于FRF殘差和振動控制方程誤差各參數(shù)反演分析結果列于表5、6。

表4 采用初始估計Ⅰ基于振動控制方程誤差的反演結果Table 4 Inversion result based on vibration control equation error under initial estimation Ⅰ

表5 采用初始估計Ⅱ基于FRF殘差的反演結果Table 5 Inversion result based on FRF residual under initial estimation Ⅱ

表6 采用初始估計Ⅱ基于振動控制方程誤差的反演結果Table 6 Inversion result based on vibration control equation error under initial estimation Ⅱ

由表5、6可見，反演結果對應的基頻與實測基頻之間的相對偏差隨測量數(shù)據(jù)受噪聲影響程度的增加而增加。當采用相對準確的初始估計Ⅱ時，基于FRF殘差的方法在利用表1所列的前3組頻率組合對應的地震響應可獲得較為可靠的反演結果，即工況5-1、5-2和5-3。基于振動控制方程誤差的方法利用表1所列的前兩組頻率組合時能獲得可靠的反演結果，即工況6-1和6-2。這說明基于FRF殘差的方法具有更強的魯棒性。而工況4-1與6-1的反演結果完全相同進一步證明了當實測地震響應具有足夠的精度時，修正參數(shù)的初始估計對基于振動控制方程誤差方法反演結果的影響較小。值得注意的是，在工況4-1、5-1和6-2中各參數(shù)反演結果存在明顯差異，但對應理論基頻與試驗結果的相對偏差均在1%以內(nèi)，因此有必要通過驗證修正后FEM計算所得理論響應的功率譜密度(PSD)曲線與實測響應PSD曲線之間的吻合程度來判斷反演結果的可靠性。

采用Rayleigh阻尼模擬該管道的耗能機制。取管道阻尼比為2%，上下截止頻率為50 Hz和5 Hz，依據(jù)文獻[19]中的公式計算該管道對應的理論質(zhì)量和剛度阻尼系數(shù)分別為1.097和1.14×10-4。一般地，理論阻尼系數(shù)對應的響應與實測響應之間總是存在一定偏差，由于已建立了該管道的質(zhì)量和剛度矩陣，參考文獻[20-21]，采用手動調(diào)整的方法，以各測點理論和實測響應峰值相吻合為目標對理論阻尼系數(shù)進行修正。最終確定圖4所示FEM中各單元修正后的質(zhì)量阻尼系數(shù)為1.018，與管接頭位置對應的單元8的修正后的剛度阻尼系數(shù)為1.428×10-4，其他各單元修正后的剛度阻尼系數(shù)為8.08×10-5。由阻尼修正結果可知管接頭的存在使該管道阻尼分布表現(xiàn)出一定的非比例特點。相對2%阻尼比下的理論剛度阻尼系數(shù)，管接頭部位的剛度阻尼系數(shù)較大而管道阻尼系數(shù)較小。

測點7、8、10的理論與實測PSD曲線的對比示于圖8。由圖8可見，與各工況對應的理論PSD曲線與實測結果在波形上基本吻合，但理論曲線峰值與實測結果之間依然存在差異。這是由大地脈動和振動臺驅動油泵的振動引起的試驗結構初始狀態(tài)不為0導致。工況4-1對應的理論與實測結果之間吻合程度最高，說明工況4-1所示的反演結果是最可靠的。對應理論基頻十分接近的參數(shù)組合計算所得的PSD曲線在幅值上存在明顯差異，說明固有頻率不能作為反演分析結果可靠的唯一判據(jù)，特別是在測得固有頻率階數(shù)較少的情況下。

a——測點7；b——測點8；c——測點10圖8 理論與實測PSD曲線對比Fig.8 Comparison of theoretical and measured PSD curves

由工況4-1反演結果可知管道FEM中左右兩端線性彈簧轉動剛度不同，這意味著管道在兩端支架的外伸長度不一致，且固定螺母通過手工加工、手動擰緊導致支架為系統(tǒng)提供了不同強度的轉動約束。單元8抗彎剛度也小于其他管單元，說明兩段管件連接部位相對于其他部位相當于存在一定程度的損傷，管接頭的存在并未增加該部位的剛度。

3.3 利用位移響應進行反演

Lin等[13]提出在基于FRF殘差的修正方法中，當測量數(shù)據(jù)為加速度時應將加速度轉換為位移完成反演過程。為研究上述兩種方法采用不同類型的響應時對反演結果的影響，將表1所列4組頻率對應的加速度響應在頻域轉換為位移對管道FEM中的修正參數(shù)進行反演分析。采用初始估計Ⅰ利用兩種方法對各參數(shù)反演結果列于表7、8，采用初始估計Ⅱ利用兩種方法對各參數(shù)的反演結果列于表9、10。此時動剛度矩陣采用式(12)的形式。

由表7～10可知，采用位移響應進行反演分析時，只有對參數(shù)的初始估計足夠精確才能獲得可靠的反演結果。工況9-1、9-2和9-3可獲得基本相同的反演結果，說明利用位移響應可提高基于FRF殘差修正方法反演結果的穩(wěn)定性。工況9-4、10-3和10-4依然無法獲得可靠的反演結果，說明利用位移響應對提高修正方法的魯棒性沒有貢獻。此外值得注意的是，在表9、10中低頻數(shù)據(jù)對應反演結果的可靠性高于高頻，說明利用位移響應進行反演分析降低了測量精度較高的高頻響應數(shù)據(jù)對應反演結果的準確性，從而增加了兩種修正方法對修正參數(shù)初始估計的依賴性。

表7 利用位移響應情況下采用初始估計Ⅰ基于FRF殘差的反演結果Table 7 Inversion result based on FRF residual under initial estimation Ⅰ by using displacement response

表8 利用位移響應情況下采用初始估計Ⅰ基于振動控制方程誤差的反演結果Table 8 Inversion result based on vibration control equation errorunder initial estimation Ⅰ by using displacement response

表9 利用位移響應情況下采用初始估計Ⅱ基于FRF殘差的反演結果Table 9 Inversion result based on FRF residual under initial estimation Ⅱ by using displacement response

表10 利用位移響應情況下采用初始估計Ⅱ基于振動控制方程誤差的反演結果Table 10 Inversion result based on vibration control equation errorunder initial estimation Ⅱ by using displacement response

4 結論

本文將地震響應與地震激勵的比值定義為一種新的FRF，分別以FRF殘差和管道振動控制方程誤差為目標特征值，利用模擬地震響應對某一管道系統(tǒng)FEM中的結構參數(shù)進行反演分析，對比二者修正過程得到以下結論。

1) 基于FRF殘差的FEM修正方法具有更強的魯棒性，而基于振動控制方程誤差的FEM修正方法對修正參數(shù)的初始估計依賴性較低，在合適條件下，二者識別結果具有類似的精度。

2) 采用位移響應可提高基于FRF殘差修正方法反演結果的穩(wěn)定性，但導致兩種方法在利用受噪聲影響較小的高頻響應進行修正時反演結果準確性降低，反而增加了修正過程對參數(shù)初始估計的依賴性。

3) 利用地震響應對核電管道FEM中的結構參數(shù)進行反演是一種簡便可行且更加科學的方法，一方面避免了附加的模態(tài)提取試驗和傳統(tǒng)的FRF測量試驗，另一方面解決了不同參數(shù)組合對應相同模態(tài)時反演結果可靠性難以評價的困難。

4) 本文方法并不限于地震響應，也可直接利用其他類型基礎激勵下的響應來完成核電管道中未知結構參數(shù)的反演分析。對于已安裝狀態(tài)的管道系統(tǒng)而言，可實時監(jiān)測其服役過程中安裝樓板位置的振動和管道在樓板激勵下的響應來反演其結構參數(shù)。對于典型支撐和連接設備應建立其修正后模型的數(shù)據(jù)庫，在類似系統(tǒng)的建模中直接使用。

核電管道在服役過程中的任何振動測量數(shù)據(jù)都有潛力作為反演其結構參數(shù)的試驗數(shù)據(jù)。管道的自激振動是比較常見的，如文獻[22]中提到的主蒸汽隔離閥氣流誘發(fā)的管道振動，只要能構造出管道的振動控制方程，亦可按照本文方法進行修正，難點在于如何確定激勵的相關信息，這將是下一步的研究方向。