嚴(yán) 偉，朱景文，周衛(wèi)雪

（1.藍(lán)箭航天空間科技股份有限公司，北京 100176；2.北京航天長征飛行器研究所，北京 100076）

0 引言

金屬波紋管是一種具有褶合輪廓的柔性管材，在管路系統(tǒng)中通過褶合處的變形補(bǔ)償管道連接處的相互位移，同時吸收振動能量。這使得其在保持局部剛性的同時，整體又具有柔韌性，因此，被廣泛應(yīng)用于液體火箭發(fā)動機(jī)的推進(jìn)劑輸送系統(tǒng)中。

液體火箭發(fā)動機(jī)常采用低溫介質(zhì)液氧作為助燃劑。當(dāng)液氧在波紋管中輸送時，流動狀態(tài)通常是復(fù)雜的非定常湍流。同時，受褶合區(qū)域復(fù)雜幾何的影響，液氧的流動會呈現(xiàn)出一定的周期性，進(jìn)而誘導(dǎo)波紋管結(jié)構(gòu)發(fā)生振動［1-2］。流體誘導(dǎo)振動的強(qiáng)度通常不大，影響可忽略，但如果振動頻率接近波紋管固有頻率，將引起結(jié)構(gòu)共振，此時迅速增大的交變應(yīng)力會在很短時間內(nèi)令波紋管疲勞失效［3-5］。

美國早在20 世紀(jì)60 年代就開展了波紋管流體誘導(dǎo)振動的相關(guān)研究。在土星5 運(yùn)載器第2 次飛行期間，其中1 臺2 級J-2 發(fā)動機(jī)發(fā)生故障導(dǎo)致提前關(guān)機(jī)；同一次飛行，3 級的J-2 發(fā)動機(jī)在地球軌道上起動失敗。這些異常最后查明是由于流體誘導(dǎo)振動導(dǎo)致波紋管疲勞失效引起的。這些失效案例促使美國對流體誘導(dǎo)波紋管振動問題進(jìn)行了廣泛研究，調(diào)查分析了超過50 個不同種類的撓性軟管和管路樣本［6］。

早期的流動可視化試驗，觀察到波紋管褶合處有漩渦不斷脫落，確認(rèn)了流體在波紋管內(nèi)壁的渦脫落是誘導(dǎo)振動的原因?；跍u脫落理論，美國國家航空航天局（NASA）建立了流體誘導(dǎo)波紋管振動的數(shù)學(xué)模型，提出了估算波紋管在給定幾何尺寸和來流速度下的流體誘導(dǎo)振動應(yīng)力的半經(jīng)驗方法。該方法將波紋管每半個褶合等效為集中質(zhì)量，各集中質(zhì)量之間采用彈簧連接，通過計算波紋管固有頻率以及渦脫落頻率，求得激起的波紋管模態(tài)階數(shù)。在此基礎(chǔ)上估算波紋管應(yīng)力，同時編制了計算波紋管流體誘導(dǎo)交變應(yīng)力的計算程序。該方法對高壓撓性軟管以及剛度較大的自由金屬波紋管的預(yù)測較為準(zhǔn)確。應(yīng)用該方法對土星/阿波羅運(yùn)載器上所有使用的波紋管進(jìn)行了計算，以判斷流體誘導(dǎo)振動問題是否存在。此后，NASA 完善了試驗數(shù)據(jù)，進(jìn)一步修正了模型中的中的系數(shù)，提高了精度［7］。

在國內(nèi)，李中來等［8］、馬明等［9］、羅宏瀚［10］分別針對波紋管的疲勞失效問題展開了研究。其中，羅宏瀚基于商用軟件ANSYS 建立了流體誘導(dǎo)波紋管振動的計算模型，比較研究了波紋管參數(shù)、流動等性質(zhì)對振動的影響，定性分析了不同條件下振動的特征，并根據(jù)計算結(jié)果提出了波紋管流固耦合的減振方法。

本文利用二維大渦模擬模型，對波紋管內(nèi)液氧流動進(jìn)行非定常數(shù)值仿真，發(fā)展了波紋管流動的渦脫落幅頻特性預(yù)測方法；研究了液氧流動速度、波紋管幾何對渦脫落頻率的影響，并歸納出了計算渦脫落頻率的經(jīng)驗式，彌補(bǔ)了已有研究在工程適用性上的不足。

1 數(shù)學(xué)模型

液氧在液體火箭發(fā)動機(jī)波紋管內(nèi)的流動本質(zhì)上是一種復(fù)雜的三維、非定常、高雷諾數(shù)湍流。由于波紋管的結(jié)構(gòu)具有軸對稱性（如圖1 所示），液氧內(nèi)部的流動可近似視為軸對稱流動，故此可將三維流動問題簡化為二維。

同時引入流動不可壓假設(shè)，并忽略重力和熱傳導(dǎo)影響，則可用如下二維、黏性、不可壓控制方程來描述波紋管內(nèi)液氧的流動。

連續(xù)方程為

圖1 波紋管幾何模型Fig.1 Schematic diagram of bellows geometric model

x方向動量方程為

y方向動量方程為

式中：u、v分別為流動沿x、y方向的速度分量；p為壓強(qiáng)；ρ、μ分別為液氧的密度和動力黏性系數(shù)；t為時間。

本文采用大渦模擬模型來描述波紋管內(nèi)的湍流運(yùn)動。它的基本思想為：大尺度脈動（或大尺度湍流渦旋）的結(jié)構(gòu)依賴于個別的流動問題，不存在通用模型，因此，直接通過求解Navier-Stokes 方程獲得；而小尺度脈動近似各向同性，其對大尺度運(yùn)動的作用可以使用模型假設(shè)［11］。然而，大渦模擬是針對三維問題提出的，為了能在二維問題中應(yīng)用該模型，本文參考文獻(xiàn)［12］中建立的二維大渦模擬格式進(jìn)行數(shù)值計算。

2 計算方法

2.1 網(wǎng)格劃分

采用有限容積法對計算區(qū)域進(jìn)行離散，生成如圖2 所示的結(jié)構(gòu)化網(wǎng)格，并且對近壁面處網(wǎng)格進(jìn)行加密。在波紋管褶合區(qū)域內(nèi)，使用C 型網(wǎng)格。

圖2 網(wǎng)格劃分示意圖Fig.2 Schematic diagram of meshing

2.2 定解條件

在進(jìn)口截面給定速度邊界條件

在出口截面給定壓力邊界條件

在壁面處給定無滑移固壁邊界

在對稱邊界給定

式中：下標(biāo)inn 為進(jìn)口截面；out 為出口截面；wall 為固壁；axis 為對稱軸。初始條件給定波紋管褶合內(nèi)速度為零，其余部分同進(jìn)口邊界條件。

3 仿真結(jié)果及分析

3.1 仿真結(jié)果

進(jìn)口速度uinn=30 m/s 時，流動充分發(fā)展后在通道褶合附近8 個不同時刻的速度云圖如圖3 所示，時間差Δt=1.2×10-4s。由圖可知：主流區(qū)的速度波動較小；當(dāng)液氧流經(jīng)褶合的上坡側(cè)時，流道擴(kuò)張導(dǎo)致速度下降，并較大的旋渦產(chǎn)生、脫落；當(dāng)液氧流經(jīng)褶合的下坡側(cè)時，流道收縮，速度增大，并且受回流的影響，之前在上坡側(cè)出現(xiàn)的脫落渦被卷吸進(jìn)入褶合內(nèi)部，與褶合內(nèi)原有的渦相互作用。該過程不斷重復(fù)，誘導(dǎo)波紋管結(jié)構(gòu)發(fā)生振動，振動的頻率等于褶合上坡側(cè)的渦脫落頻率。

圖3 不同時刻波紋管內(nèi)的速度分布（uinn=30 m/s）Fig.3 Contours of velocity magnitude（uinn=30 m/s）at different time in the bellows

為了定量分析渦脫落的頻率，在C點(diǎn)（見圖1 標(biāo)注）設(shè)置數(shù)字探針，監(jiān)控其速度的脈動。當(dāng)流動充分發(fā)展后，截取0.1 s 之內(nèi)的速度信號如圖4 所示。橫坐標(biāo)表示時間，縱坐標(biāo)表示C點(diǎn)的速度大小。由圖4 可見，C點(diǎn)速度隨時間高頻振蕩，平均速度為21.57 m/s。

對速度-時間信號進(jìn)行快速傅里葉分析（FFT），得到其幅頻特性如圖5 所示。橫坐標(biāo)表示信號的頻率成分，縱坐標(biāo)表示相應(yīng)頻率下的振動幅值。

圖4 C 點(diǎn)速度信號（uinn=30 m/s）Fig.4 Velocity signal at point C（uinn=30 m/s）

圖5 C 點(diǎn)速度信號的幅頻特性（uinn=30 m/s）Fig.5 Amplitude-frequency characteristic of the velocity signal at point C（uinn=30 m/s）

由圖可知：圖4 所示速度-時間信號的基頻為2 030 Hz，對應(yīng)速度分量波動的幅值為0.684 2 m/s，約占平均速度的3%；二階頻率對應(yīng)分量的波動幅值為0.225 7 m/s，約占平均速度的1%；三階及其以上頻率分量強(qiáng)度很小，可忽略。

為了分析流動速度對渦脫落頻率的影響，本文設(shè)計了多組對比算例，對不同進(jìn)口速度uinn下的流場進(jìn)行了仿真，并分別對C點(diǎn)的速度-時間信號進(jìn)行FFT 處理，不同進(jìn)口速度下的基頻見表1?？梢?，渦脫落運(yùn)動的基頻與液氧流動速度近似成正比，液氧流動速度越大，渦脫落頻率越高。

表1 不同uinn下C 點(diǎn)速度波動的基頻Tab.1 Base frequency of C-point velocity signal at different uinn

如圖3 所示，渦脫落運(yùn)動的動態(tài)過程，渦脫落頻率除了與液氧流動速度有關(guān)，還可能受波紋管褶合處幾何參數(shù)的影響。

本文基于控制變量法，在波紋管進(jìn)口速度不變的條件下，分別研究了r1、r2（如圖1 所示）對渦脫落頻率的影響。

某一較小r1、進(jìn)口速度為30 m/s 時，波紋管內(nèi)液氧流動充分發(fā)展后，在通道褶合附近8 個不同時刻的渦量云圖如圖6 所示，時間差Δt=1.0×10-4s。從渦量云圖上也可以觀察到周期性的渦脫落運(yùn)動，同時，由于褶合上坡側(cè)曲率半徑r1較?。ㄏ鄬τ趫D3所示），加快了渦脫落的過程。

圖6 某較小r1條件下不同時刻波紋管內(nèi)的渦量分布Fig.6 Contours of vorticity under a small r1 at different time in the bellows

對不同r1、r2下C點(diǎn)的速度-時間信號分別進(jìn)行FFT 處理，得到速度波動的基頻與波紋管r1、r2的關(guān)系，結(jié)果見表2。其中，r0為特征半徑，利用它對r1、r2進(jìn)行無量綱處理。由表2 可見，r1、r2越大，速度波動的基頻越小，意味著渦脫落頻率隨r1、r2增大而減小。

表2 不同r1、r2下C 點(diǎn)速度信號的基頻Tab.2 Base frequency of C-point velocity signal at different r1，r2

3.2 結(jié)果分析

綜上所述，波紋管內(nèi)渦脫落運(yùn)動的基頻fst與進(jìn)口速度uinn正相關(guān)，與r1、r2負(fù)相關(guān)，據(jù)此將fst與uinn、r1、r2之間的函數(shù)關(guān)系歸納為如下形式：

基于式（4），對表1 和表2 中數(shù)據(jù)進(jìn)行最小二乘分析，確定比例系數(shù)ξ=0.26。

分析經(jīng)驗式（4）可知，當(dāng)r1、r2一定時，fst與uinn成正比。經(jīng)驗式計算值與表1 中仿真數(shù)據(jù)對比如圖7所示，數(shù)據(jù)點(diǎn)與擬合直線貼合較好，表明經(jīng)驗式正確反映了波紋管內(nèi)渦脫落運(yùn)動的基頻與流動速度之間的關(guān)系。

圖7 經(jīng)驗式計算值與表1 數(shù)據(jù)對比Fig.7 Comparison of empirical formula calculations with Tab.1

當(dāng)uinn一定時，fst與r1、r2之和成反比。經(jīng)驗式計算值與表2 中仿真數(shù)據(jù)對比如圖8 所示。圖中曲面表示uinn=30 m/s 時fst與r1、r2之間的函數(shù)。從圖中可以看出，表2 中數(shù)據(jù)基本上均勻分布在曲面兩側(cè)，且與函數(shù)曲面貼合較好，這表明經(jīng)驗式合理反映了波紋管內(nèi)渦脫落運(yùn)動的基頻與r1、r2之間的關(guān)系。

圖8 經(jīng)驗式計算值與表2 數(shù)據(jù)對比Fig.8 Comparison of empirical formula calculations with Tab.2

在流體力學(xué)中，通常用斯特羅哈數(shù)St來表征渦脫落頻率。用St代替經(jīng)驗式（4）中的比例系數(shù)ξ，可得

已有試驗研究證明，流體在波紋管內(nèi)產(chǎn)生的渦脫落運(yùn)動，與圓柱繞流時的渦脫落運(yùn)動具有相近的斯特羅哈數(shù)。對于圓柱繞流，一般其St大約在0.2附近。這與本文通過仿真歸納出的系數(shù)0.26 較接近。

4 結(jié)束語

本文針對液氧在波紋管內(nèi)的非定常流動，采用大渦模擬方法進(jìn)行了數(shù)值仿真，研究了渦脫落幅頻特性與流動速度、波紋管幾何特征之間的關(guān)系，并歸納了波紋管流動渦脫落頻率的經(jīng)驗式。主要結(jié)論如下：

1）液氧在波紋管內(nèi)流動，會在褶合上坡側(cè)產(chǎn)生脫落渦；

2）渦脫落運(yùn)動的基頻與來流速度成正比，與褶合特征長度成反比，使用經(jīng)驗式（4）可以估算渦脫落頻率。

后續(xù)計劃開展相關(guān)試驗，驗證或修正經(jīng)驗式（4），為液體火箭發(fā)動機(jī)波紋管的防疲勞設(shè)計提供技術(shù)支撐。