基于Newmark－β法的換熱管振動反問題的求解方法

徐 波，吳遠(yuǎn)峰，蘇文獻(xiàn)，許 伍

（上海理工大學(xué) 能源與動力工程學(xué)院，上海 200093）

換熱器是化工、石油、動力、食品、航天以及其他工業(yè)部門的通用工藝設(shè)備，在國民經(jīng)濟(jì)和工業(yè)生產(chǎn)領(lǐng)域中起著至關(guān)重要的作用［1］.它的機(jī)械結(jié)構(gòu)性能以及傳熱性能的好壞將直接影響產(chǎn)品質(zhì)量、實際生產(chǎn)中能源的綜合利用率，以及整個系統(tǒng)運(yùn)行的經(jīng)濟(jì)性、安全性和可靠性.在各種形式的換熱器中，管殼式換熱器由于結(jié)構(gòu)簡單、用材范圍廣、清洗方便、適應(yīng)性好、生產(chǎn)成本低、運(yùn)行安全可靠、流量大、耐高壓高溫等特點在產(chǎn)量和應(yīng)用范圍等方面處于主導(dǎo)地位.在世界各國工程師的共同努力下，換熱器的諸多領(lǐng)域雖然取得了可喜的進(jìn)步，但仍有一些問題至今未能解決，特別是機(jī)理復(fù)雜的流體誘導(dǎo)振動問題，成了制約換熱器進(jìn)一步發(fā)展的關(guān)鍵因素.

目前諸多學(xué)者在預(yù)防換熱器管束振動正問題方面做了大量的試驗和研究，取得了很多進(jìn)展，提出了不少理論見解和試驗判據(jù).很多方法都可以有效避免管束的共振，但是由于流動的復(fù)雜性以及振動因素的不確定性，仍然不能從根本上解決流體誘導(dǎo)換熱器的管束振動.幾十年來，許多學(xué)者對振動反問題產(chǎn)生了極大的興趣，也提出了一些巧妙的反問題求解方法［2－4］，BUSBY等［5］提出了對受正弦和隨機(jī)激勵的間隙支撐單自由度系統(tǒng)沖擊力的估算方法，在時域范圍內(nèi)成功識別出了激勵的位置；WU等［6］使用時域分析技術(shù)研究了若干個同時作用的激勵的分離識別問題，采用了所謂的反卷積譜分析法，該方法對處理非分散現(xiàn)象是非常有用的，但是不能用于處理離散的彎曲波［7］.

求解這些反問題的主要困難是問題本身的病態(tài)和求解過程的數(shù)值病態(tài)，這也就導(dǎo)致了反問題求解方法對測量信號的噪聲擾動特別敏感.采用正則化技術(shù)——奇異值分解和優(yōu)化技術(shù)，該病態(tài)問題能得到部分的解決.由于求解精度問題，本文將采用WHITON［8］和DOYLE［9］提出的一些通用方法，在時域內(nèi)，采用Newmark－β法直接識別換熱管所受的載荷，完成載荷的識別.

1 Newmark－β法

1959年，NEWMARK提出了一種計算振動動力響應(yīng)方法的假設(shè)，該假設(shè)可以看成是線性加速度方法的推廣.Newmark－β法一般用于求解正問題，即根據(jù)加載的力求解結(jié)構(gòu)的位移、速度和加速度.由于其良好的收斂性，計算結(jié)果精度很高，適用性很廣.Newmark－β法所采用的方程為：

計算步驟可以簡單歸納為兩個步驟.

步驟一 初始計算

①計算結(jié)構(gòu)的剛度矩陣K、質(zhì)量矩陣M和阻尼矩陣C.

②根據(jù)初始條件u0和，計算.

③選取合適的時間步長Δt和參數(shù)α和β，計算積分常數(shù).

④計算有效剛度矩陣K′.

步驟二 對于每一時間步長

① 計算t＋Δt時刻的有效載荷向量P′.

式中：Pt＋Δt為t＋Δt時刻的載荷向量.

②求解t＋Δt時刻的位移.

③計算t＋Δt時刻的加速度和速度.

2 Newmark－β反分析法原理

Newmark－β法是在線性加速度的基礎(chǔ)上發(fā)展起來的一種計算結(jié)構(gòu)動響應(yīng)的數(shù)值算法，該算法具有很好的數(shù)值穩(wěn)定性和計算精度.由于廣泛的適用性，在一定條件下可以保證算法的無條件穩(wěn)定性，這些條件可以反向利用在逆問題上的求解.

基于時間推進(jìn)法，可以知道與時刻t相關(guān)的量是已知的，求t＋Δt時刻的量，因此可以令t時刻的有效載荷向量對式（7）做如下變形：

式中：為t＋Δt時刻的測點位移.

在每一個微小的時間步內(nèi)，對于線彈性系統(tǒng)，根據(jù)疊加原理，可得

在t＋Δt時刻，由式（11）可知，P′t為t時刻的節(jié)點有效載荷向量，為已知量.求解方程式（11）可以得到t＋Δt時刻位移分量

設(shè)在t＋Δt時刻，節(jié)點i上作用的單位節(jié)點載荷為，則有

因此，求解第i節(jié)點上施加的載荷的問題就轉(zhuǎn)化為求解載荷系數(shù)的問題.式（10）可以改寫為

求解上式，可以分別得到單獨在節(jié)點i處施加單位載荷時體系的位移向量，再根據(jù)式（14），則方程式（10）的解必定滿足以下關(guān)系：

設(shè)有m個測點，在已知測出t＋Δt時刻的位移（或速度和加速度）值（k＝1，2，…，m），根據(jù)Newmark－β法，有

式（17）實際上是具有n個未知數(shù)（i＝1，2，…，n），m個測點位移（k＝1，2，…，m），可以形成m個方程組成的線性方程組

令

則，方程組可以改寫成

當(dāng)m＝n時，解式（17）就可以求得（i＝1，2，…，n），代入式（14）可以得到，即t＋Δt時刻的振動載荷.若求出每一時刻的振動載荷，則可以得到振動載荷的時程.

當(dāng)m＜n時，測點數(shù)小于要求解的未知數(shù)，故不能識別出振動載荷；當(dāng)m＞n時，測點數(shù)大于要求解的未知數(shù)，可以采用最小二乘法得到高精度的結(jié)果，以消除誤差.

2.1 算例分析

利用Newmark－β反分析法原理，在 MATLAB編程，驗證理論推導(dǎo)過程的正確性.

已知：梁的長度L＝10 000mm，彈性模量E＝197GPa，截面面積S＝1 000mm3，截面慣性矩I＝108mm4，在梁的某位置處施加簡諧載荷F＝1 000sin（t）.?。簳r間步長Δt＝0.1s，參數(shù)α＝1／4，β＝1／2，作用時間t＝10.0s，空間步長ΔL＝500 mm，采用基于Newmark－β法的 MATLAB程序，不考慮阻尼影響，可以得到各個節(jié)點的位移、速度和加速度時間歷程曲線.由于加速度的精度更高，因此顯示加速度時間歷程曲線，如圖1所示.

圖1 各節(jié)點的加速度隨時間的變化曲線Fig.1 Acceleration curve of each node with time

根據(jù)其中的任意一個加速度歷程曲線，按照Newmark－β反分析法，可以得到載荷的識別曲線，如圖2所示.

圖2 加載點的識別載荷隨時間的變化曲線Fig.2 Identification load curve of loading node ith time

從圖2可以看出，在不考慮阻尼和環(huán)境噪聲的影響下，識別載荷和加載載荷基本上是一致的.以上只對一個加載載荷進(jìn)行了計算，同樣地對于在多個節(jié)點上施加不同的載荷時程，也可以按照此方法完成振動載荷的識別.需要指出的是，在實際工程中，位移、速度和加速度的時間歷程曲線是需要進(jìn)行測量獲得的，并作為已知條件進(jìn)行使用，同時選擇合的步長和參數(shù)，采用最小二乘法排除可能的誤差，這樣得到的識別載荷才有一定的精度和可信度.

3 結(jié)語

本文根據(jù)換熱管在激勵載荷作用下的響應(yīng)，提出了應(yīng)用Newmark逆分析法求解在已知該振動響應(yīng)下系統(tǒng)所受的激勵.在激勵載荷識別時，提出了基于時域范圍內(nèi)的識別方法，初步探索了一種較為簡單的管殼式換熱器管束振動反演問題的計算方法，即激勵載荷識別技術(shù)，為今后全面分析整體結(jié)構(gòu)的振動正反問題做鋪墊.

［1］鄭津洋，董其伍，桑芝富.過程設(shè)備設(shè)計［M］.北京：化學(xué)工業(yè)出版社，2010.ZHENG Jingyang，DONG Qiwu，SANG Zhifu.Process equipment design［M］.Beijing：Chemical Industry Press，2010.

［2］PRESS W H，TEUKOLSKY A A，VETTERLING W T，et al.Numerical recipes：theartof scientific computing［M］.Cambridge：Cambridge University Press，1992.

［3］GROETCH C W.Inverse problems in the mathematical sciences［M］.Wiesbaden：Vieweg and Sohn，1993.

［4］HANSEN P C.Regularization tools［J］.Numer Algorithms，1994（6）：1－35.

［5］BUSBY H R，TRUJILLO D M.Solution of an inverse dynamics problem using an eigenvalue reduction technique［J］.Computers ＆ Structures，1987，25：109－117.

［6］WU E，YEH J C.Identification of impact forces at multiple locations on laminated plates［J］.Journal of American Institute of Aeronautics and Astronautics，1994，32：2433－2439.

［7］KIM J T，LYON R H.Cepstral analysis as a toll for robust processing，deverberation and detection of transients［J］.Mechanical Systems and Signal Process，1992（6）：1－15.

［8］JORDAN R W，WHISTON G S.Remote impact analysis by use of propagated acceleration signals：comparison between theory and experiments［J］.Journal of Sound and Vibration，1984，97：53－65.

［9］DOYLE J.Wave propagation in structures：an FFT－based spectral analysis methodology［M］.New York：Springer Verlag，1989.