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

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

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

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

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

1 Newmark－β法

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

計(jì)算步驟可以簡(jiǎn)單歸納為兩個(gè)步驟.

步驟一 初始計(jì)算

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

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

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

④計(jì)算有效剛度矩陣K′.

步驟二 對(duì)于每一時(shí)間步長(zhǎng)

① 計(jì)算t＋Δt時(shí)刻的有效載荷向量P′.

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

②求解t＋Δt時(shí)刻的位移.

③計(jì)算t＋Δt時(shí)刻的加速度和速度.

2 Newmark－β反分析法原理

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

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

式中：為t＋Δt時(shí)刻的測(cè)點(diǎn)位移.

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

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

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

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

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

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

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

令

則，方程組可以改寫成

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

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

2.1 算例分析

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

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

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

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

圖2 加載點(diǎn)的識(shí)別載荷隨時(shí)間的變化曲線Fig.2 Identification load curve of loading node ith time

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

3 結(jié)語

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

［1］鄭津洋，董其伍，桑芝富.過程設(shè)備設(shè)計(jì)［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.