陳宏玉，劉紅軍，劉 上

（西安航天動力研究所，陜西 西安710100）

0 引言

液體火箭發(fā)動機關(guān)機時，推進劑供應(yīng)系統(tǒng)由于閥門關(guān)閉或渦輪泵減速等因素影響，不可避免地會在管路中形成水擊和壓力振蕩。這種水擊壓力波動會造成諸如推進劑泄漏、結(jié)構(gòu)破壞和推力輸出精度降低等一系列問題，在發(fā)動機設(shè)計中需要對其進行深入分析。

水擊問題可由一維瞬變流動基本方程描述，它們是一組非線性雙曲型偏微分方程，一般情況下不存在解析解，只能采用數(shù)值計算方法求解。特征線法由于其具有較高的計算精度，數(shù)值方法簡便，運算穩(wěn)定性較好和易于編程等特點廣泛應(yīng)用在水擊問題的分析中[1-4]。但由于推進劑供應(yīng)管路系統(tǒng)是包括液路、氣路和容腔的復(fù)雜管網(wǎng)系統(tǒng)，因此采用特征線法在時間步長協(xié)調(diào)上比較困難，邊界條件的處理也比較復(fù)雜，而且非線性迭代解算收斂較慢[5]。常用求解水擊問題的方法還有有限元法[6-8]、有限體積法[9]和無網(wǎng)格法[10]等。

Quarteroni,Canuto和郭本瑜等[11-12]提出并發(fā)展的譜方法 (Spectral Method)提供了除有限差分和有限元法之外的第三類數(shù)值計算方法,開辟了數(shù)值計算的全新領(lǐng)域。譜方法的最大優(yōu)點是所謂“無窮階收斂性”，即如果原問題的解充分光滑，那么用適當(dāng)?shù)淖V方法所求得的近似解將以基函數(shù)個數(shù)的任意次冪速度收斂于精確解。本文采用譜方法中的Fourier譜方法就推進劑供應(yīng)系統(tǒng)管路中閥門關(guān)閉形成的水擊問題進行數(shù)值計算，給出了數(shù)學(xué)模型，求解原理、求解方法及仿真結(jié)果，并將該結(jié)果與已發(fā)表的采用特性線法和有限元法的求解結(jié)果進行了比較。

1 一維瞬變流動控制方程

假設(shè)管道內(nèi)推進劑流動是一維絕熱有摩擦的液體瞬變流，摩擦損失均按準(zhǔn)穩(wěn)態(tài)處理，則描述管道內(nèi)液體瞬變流的基本方程包括連續(xù)方程和運動方程，是一組擬線性雙曲型偏微分方程，具體形式如下[4]：

式中：p為x斷面的流體平均壓力；u為x斷面的流體平均流速；ρ為流體密度；a為壓力波傳播速度；f為管壁摩擦損失系數(shù)；g為質(zhì)量場加速度；θ為管流軸線方向與加速度方向的夾角。

式 (1)和 (2)即為管內(nèi)非定常流的基本方程，又稱波動方程。在一般流體管道中，有a＞＞u，故可略去式 (1)和 (2)中第二項，對于水平直管還可忽略重力場影響。在此基礎(chǔ)上，令

那么，可得無量綱波動方程：

式中：p0為管道內(nèi)液體初始壓力。

2 Fourier譜方法求解

2.1 Fourier譜方法

譜方法源于經(jīng)典的Ritz-Galerkin方法，以整體無限光滑的函數(shù)系 (例如：Fourier級數(shù)，Chebyshev多項式，Legendre多項式等)作為基函數(shù)，對定義在計算域中的函數(shù)進行近似，然后結(jié)合加權(quán)余量法對偏微分方程進行求解。推進劑供應(yīng)系統(tǒng)管路中流體瞬變特性具有時間周期性，因此可通過選用特殊的周期函數(shù) (Fourier級數(shù))作為基函數(shù)，將未知函數(shù)分解展開成空間變量基與時間函數(shù)系數(shù)的分離變量形式的解式。這里將液體無量綱壓力φ(y,τ)作Fourier正弦級數(shù)展開，液體無量綱流速U(y,τ)作Fourier余弦級數(shù)展開。

式 (5)是通過φy(y,τ)的Fourier余弦級數(shù)展開，然后逐項積分所得。式 (5)和 (6)中，與時間相關(guān)的系數(shù) φk(τ)和 Uk(τ)為：

式中，k∈［1,2,… ］N 。這樣，非線性無量綱控制方程 (3)和 (4)求解就轉(zhuǎn)換為確定系數(shù)函數(shù)φk(τ)和Uk(τ)的形式。Fourier譜展開的項數(shù)N選取取決于計算精度。

將液體無量綱壓力φ (y,τ)、流速U(y,τ)的Fourier譜展開式 (5)和 (6)代入方程 (3)和 (4)中，得：

值得注意的是，非線性摩擦積分項出現(xiàn)在方程(10)中，令

對摩擦項積分，W的Fourier余弦系數(shù)可表達成離散的Fourier變換形式

式中：yp＝p/K。

式 (9)和 (10)未考慮邊界問題，在具體應(yīng)用時，帶入邊界條件和初值，可得展開系數(shù)φk(t)，Uk(t)(k=0,1,2,…N)。再將其帶入式 (5)和 (6)中，即可求得管路內(nèi)液體瞬變流的近似解。為了適應(yīng)譜方法的高精度并保證時程積分的穩(wěn)定性，采用具有四階精度的Runge-Kutta積分法。

2.2 非物理性振蕩問題

應(yīng)用譜方法求解管道內(nèi)流體瞬變流動時，常會遇到帶有大梯度甚至間斷的問題 (如水擊問題)。在它們附近數(shù)值解會出現(xiàn)大幅度的Gibbs類型的振蕩。為了克服這種困難，學(xué)者們提出很多改進方法，其中包括超粘性方法、濾波法和譜粘性消去法 (SVV)等[13-14]。其中濾波法主要思想是用一個允許低頻波通過的濾波因子使得導(dǎo)致數(shù)值不穩(wěn)定的高頻波大大減弱，那么譜方法在處理水擊和其它間斷處的振蕩基本上被抑制。

文中采用濾波法消除因間斷導(dǎo)致的非物理振蕩，對無量綱壓力和無量綱流量的展開系數(shù)φk(τ)和 Uk(τ)乘以濾波因子 σk。σk為

式中：系數(shù)α和η的取值參考文獻 [14]。

3 算例及結(jié)果分析

為了驗證Fourier譜方法計算水擊問題的有效性，對文獻 [15]中的單管閥門關(guān)斷問題 (見圖1)進行了仿真計算。問題描述：流體介質(zhì)為水，原始參數(shù)：pT=15.0 MPa，pE=14.8 MPa，水平圓管長度L=4.8 m，內(nèi)徑D=0.01 m，ρ=1 000 kg/m3，α=1 200 m/s，CVAV=0.7 A，A為管道流通面積，摩擦系數(shù)f=0.018。

考慮閥門兩種關(guān)閉方式：

(a)管道閥門突然關(guān)閉；

圖1 簡單管路系統(tǒng)Fig.1 A simple pipeline system with valve and tank

按Fourier譜方法 (SFM)的計算結(jié)果與按特征線方法 (MOC)[1]和分段集中參數(shù)有限單元法(FEM)[8]的計算結(jié)果比較見圖2和圖3。計算時，F(xiàn)ourier展開項數(shù)取8，為保證特征線方法和有限單元的精度，特征線方法分段數(shù)取為40，分段集中參數(shù)有限單元法取分段數(shù)為50。由圖可見：采用Fourier譜方法計算單管水擊問題無論是最大水擊壓力還是水擊壓力波的衰減速率都能與特征線方法和分段集中參數(shù)有限單元法保持一致。

濾波法的作用是避免因大梯度參數(shù)變化或間斷問題引起計算時的非物理振蕩。從圖4可以看出，如果不采用濾波技術(shù)，在管內(nèi)流體實際振蕩曲線上疊加一個高頻振蕩，而采用濾波技術(shù)能較好的抑制這一高頻振蕩。

圖4譜濾波后計算結(jié)果對比Fig.4 Contrast of calculated results after filtering by spectral-Fourier method

圖5 給出了譜展開項數(shù)不同，利用無濾波的Fourier譜方法計算單管水擊問題時，閥門上游壓力變化。可以看出，當(dāng)Fourier展開項數(shù)由8項增加到16和32項時，由數(shù)值計算引起的高頻振蕩未能得到有效抑制，因此采用Fourier譜方法計算大梯度或間斷問題時，增加展開項數(shù)并不能有效抑制非物理振蕩。

圖5 Fourier展開項數(shù)對數(shù)值振蕩的影響Fig.5 Influence of different Fourier terms on numerical oscillation

4 結(jié)論

推進劑供應(yīng)系統(tǒng)管路水擊問題的Fourier譜方法求解結(jié)果表明：

1）將Fourier譜方法應(yīng)用于推進劑供應(yīng)系統(tǒng)管路水擊這一強非線性問題計算時，求得結(jié)果與特征線法和分段集中參數(shù)有限單元法結(jié)果一致，表明方法和仿真結(jié)果可信。

2）采用濾波技術(shù)能夠較好地抑制數(shù)值解在間斷處產(chǎn)生的非物理振蕩。增加Fourier展開項數(shù)并不能有效抑制非物理振蕩。

3）采用Fourier譜方法求解一維瞬變流動模型時，邊界處理簡單，可以方便地連接集中參數(shù)元件，適用于復(fù)雜流體網(wǎng)絡(luò)的建模與仿真，工程實際應(yīng)用有廣闊前景。

[1]WYLIE E B,STREETER V L.Fluid transient in systems[M].Endlewood Cliffs,USA:Prentice Hall,1993.

[2]RUTH E K,AIM H,BAKER R L,et al.Advanced liquid rocket engine transient model,AIAA1990-2299[R].USA:AIAA,1990.

[3]SAGNICK H D,KRFILL G.Numerical simulation of transients in feed systems for cryogenic rocket engines,AIAA 95-2967[R].USA:AIAA,1995.

[4]林景松,王平陽,高紅,等.液體火箭發(fā)動機關(guān)機水擊的數(shù)值模擬[J].上海航天,2008,25(3):53-57.

[5]張育林,劉昆,程謀森.液體火箭發(fā)動機動力學(xué)理論與應(yīng)用[M].北京:科學(xué)出版社,2005.

[6]KOLCIO K,HELMICKI A J,JAWEED S.Propusion system modeling for condition monitoring and control:partⅠtheoretical foundation,AIAA 94-3227[R].USA:AIAA,1994.

[7]KOLCIO K,HELMICKI A J,JAWEED S.Propusion system modeling for condition monitoring and control:partⅡapplication to the SSME AIAA 94-3228[R].USA:AIAA,1994.

[8]劉昆,張育林.一維可壓縮流的有限元狀態(tài)空間模型[J].推進技術(shù),1999,19(5):62-66.

[9]MAJUMDAR A k.Numerical modeling of fluid transient by a finite volume procedure for rocket propulsion systems[C]//Proceedings of4th ASME/JSME JointFluids Engineering Conference.Honolulu,Hawaii,USA:ASME/JSME,2003:2967-2974.

[10]TIJSSELING A S,BERGANT A.Meshless computation of water hammer[C]//Proceedings of 2nd IAHR International Meeting of the Workgroup on Cavitation and Dynamic Problems in Hydraulic Machinery and Systems.Timisoara,Romania:IAHR,2007:1-34.

[11]CANUTO C,QUARTERONI A.Approximation results for orthogonal polynomials in Sobolev spaces[J].Math Compu,1982,38(1):67-86.

[12]GUO B Y.Spectral methods and their applications[M].Singapore:World Scientific,1998.

[13]MA H P.Chebyschev-Legendre super spectral viscosity method for nonlinear conservation laws[J].SIAM J.Numer.Anal.,1998,35(3):893-908.

[14]SHEN Jie,TANG Tao.Spectral and high-order methods with applications[M].Beijing:Science Press,2006.

[15]劉昆,張育林.推進劑供應(yīng)管道的集中參數(shù)近似模型研究-模態(tài)近似模型[J].推進技術(shù),1998,19(4):41-45.