何德勝,郭 正,鮑福廷,李廣武

(1.西北工業(yè)大學(xué)航天學(xué)院,西安 710072;2.國防科技大學(xué)航天與材料工程學(xué)院,長沙 410073;3.中國航天科技集團公司四院401所,西安 710025)

0 引言

減壓器是一種能夠根據(jù)上游氣體的壓力自動調(diào)節(jié)自身活門的開度,從而使下游出口的氣體壓力保持穩(wěn)定的機械裝置,被廣泛用于航天飛行器燃料箱壓力保持和地面試車臺氣源供氣系統(tǒng),起著重要的控制和調(diào)節(jié)作用。本文研究了一種用于固體沖壓發(fā)動機地面試車臺的減壓器。其在不工作時處于常開狀態(tài),工作時高壓氣體由進氣口進入高壓腔,然后通過主活塞座上的孔道進入低壓腔,當(dāng)?shù)蛪呵粔毫Ω哂陬A(yù)定值時,主膜片所受氣壓力大于彈簧力,從而向下方運動,此時主活塞在副彈簧的推動下也向下運動,使活門開度減小,低壓腔壓力隨之減小。通過這種自動調(diào)整,低壓腔壓力可維持在設(shè)定值。彈性運動部件中的活塞頂桿和膜片硬芯既可獨立運動,也可貼合在一起運動,并可能發(fā)生碰撞,導(dǎo)致減壓器的內(nèi)部結(jié)構(gòu)非常復(fù)雜。在上游壓力開啟作用的極短時間內(nèi),減壓器內(nèi)部結(jié)構(gòu)受到?jīng)_擊,伴隨彈性敏感元件的運動,內(nèi)部氣體流動必然是非定常的,并且與結(jié)構(gòu)的運動相耦合,有可能表現(xiàn)出強烈的非線性波動,嚴(yán)重時甚至導(dǎo)致膜片破裂,使減壓器失去功能。因此工程上迫切需要采用流固耦合數(shù)值模擬技術(shù)對減壓器系統(tǒng)進行數(shù)值分析。目前國外對于管路閥門系統(tǒng)動態(tài)特性的研究已經(jīng)開始大量采用非定常計算流體力學(xué)(CFD)技術(shù),分析閥門內(nèi)部流動參數(shù)的非定常、非線性特性,如文獻[1]。近年來,為了研究閥門開閉過程中的動態(tài)特性,HOSANGADI與CAVALLO[2-3]采用動網(wǎng)格技術(shù)模擬閥門內(nèi)部結(jié)構(gòu)運動與流體的耦合特性,目前其采用的固體運動模型主要是已知運動規(guī)律的剛體模型。國內(nèi)采用CFD技術(shù)研究閥門流動特性仍以定常問題為主[4-5],對于水錘等瞬態(tài)問題的研究主要采用特征線方法[6-7],該方法難以得出三維空間中流場參數(shù)隨時間變化的規(guī)律,存在一定的局限性。

本文模擬減壓器動態(tài)特性采用流固耦合方法,流場計算采用三維非定常Euler模型,固體結(jié)構(gòu)變形采用一維質(zhì)量彈簧模型。

1 數(shù)值方法

1.1 控制方程

某固體沖壓發(fā)動機地面試車臺的減壓器結(jié)構(gòu)外形及計算機網(wǎng)絡(luò)如圖1所示。其主要由高壓腔、低壓腔及彈性運動部件組成,彈性運動部件包括活塞頂桿和膜片硬芯兩部分。

圖1 減壓器結(jié)構(gòu)外形及計算網(wǎng)格Fig.1 Geometry and computational grids of a pressure regulation value

ALE(Arbitrary Lagrangian-Eulerian)有限體積描述下的三維可壓縮非定常Euler方程可表示為如式(1)的積分形式:

其中,Ψ是控制體;?Ψ是控制體邊界;→n為控制體邊界外法向單位向量;守恒變量Q和對流項為

式中 U為流體相對于網(wǎng)格的速度;at為網(wǎng)格運動的法向速度;nx、ny、nz是→n的3個分量。

式中 xt、yt、zt分別是網(wǎng)格運動速度的3個分量。

1.2 離散方法

在控制體(三維網(wǎng)格單元)內(nèi)積分式(1),有

式中 m表示時間層;Sk是第k個積分面的矢量面積。

為得到積分面上對流項通量的二階近似,首先采用泰勒展開法,由單元中心處向積分面中心對基本變量進行重建,然后用通量矢量分裂方法(Van-Leer或Steger分裂法)解決積分面上的Riemann問題,即

其中,q表示基本變量,q=(ρ,u,v,w,p)T基本變量的重建值由式(6)計算:

式中 φij是限制器。

采用傳統(tǒng)的格林公式方法求解變量梯度,即

以三維網(wǎng)格單元為控制體,近似計算式(7)右端的面積分,可得單元中心處的變量梯度:

上述公式中的變量定義見圖2。為了抑制流場中物理量間斷處可能出現(xiàn)的數(shù)值振蕩,本文改進了Barth和Jespersen[8]的通量限制器。

1.3 對虛擬擋板面的積分

采用動網(wǎng)格技術(shù)真實模擬彈性敏感元件的動態(tài)過程,然而活門關(guān)閉時與活門座接觸,微微開啟時與活門座分離,這一過程中計算域拓撲發(fā)生變化,動網(wǎng)格處理比較復(fù)雜且費時。接觸問題是目前動網(wǎng)格技術(shù)的難點,一些涉及動網(wǎng)格技術(shù)的商業(yè)軟件都未很好解決這一問題。為實現(xiàn)對活門關(guān)閉∕開啟過程的數(shù)值模擬,本文提出了可變邊條零厚度虛擬擋板技術(shù)。即在網(wǎng)格生成時人為地在活門與活門座之間布置一個無厚度的擋板,當(dāng)活門關(guān)閉時給該擋板面雙向賦壁面條件,當(dāng)活門部分開啟時則根據(jù)實際開度計算出開啟部分所占的面積比,計算通量時認為開啟部分是內(nèi)部積分面,作通氣處理。

圖2 變量定義方式示意圖Fig.2 Definition of variables

式(4)是對流場中一般網(wǎng)格單元的離散形式,等號右邊的求和是對四面體單元的4個三角形表面進行的。如果網(wǎng)格單元的某一個表面位于虛擬擋板上,不失一般性,假設(shè)圖2中單元i的p1p2p3表面位于虛擬擋板上,則式(4)對于單元i應(yīng)改寫為

式(9)中等號右邊的兩項分別為3個一般表面的通量求和以及位于虛擬擋板上的表面的通量。Sb、Fb分別是位于虛擬擋板上的表面的面積和其中心處的對流項。

圖3是p1p2p3表面在虛擬擋板上的示意圖,圖3中預(yù)置的虛擬擋板寬度為2層網(wǎng)格,某時刻物體間的實際距離為陰影部分所示,此時三角形p1p2p3上陰影部分作通氣處理,而其他部分作壁面處理。于是有

式中 Sf、Sw分別是通氣部分和壁面部分的面積。

對流項Fb按下面的方法計算:

式中 Ff是作為通氣的內(nèi)部單元面時該表面上的對流項,由式(5)計算得到;Fw是作為壁面時該表面上的對流項。

Fw對表面的通量積分由式(12)計算:

對壁面采用無粘流的無穿透條件,即U·→n=0,則式(12)簡化為

計算通量積分時對流項F應(yīng)取積分面上的平均值,可以證明,對于三角形積分面,其中心處的物理量值是整個面上平均值的二階近似,因此對于二階格式,可采用三角形中心處的對流項。對于具有虛擬擋板表面的網(wǎng)格單元,應(yīng)該針對該表面通氣和壁面2種狀態(tài)分別求解變量梯度,計算對流項Ff和Fw時分別采用相應(yīng)的梯度進行基本變量的重建。

圖3 虛擬擋板示意圖Fig.3 Sketch of virtualbaffle faces

時間離散采用四步Runge-Kutta方法。動網(wǎng)格技術(shù)采用彈簧近似方法實現(xiàn),詳細描述可參見文獻[9]。

2 彈性敏感元件的運動方程

減壓器運動機構(gòu)包括主彈簧、副彈簧、閥芯(頂桿)、膜片,組成一個質(zhì)量彈簧阻尼系統(tǒng)。圖4為系統(tǒng)運動模型及參數(shù)定義示意圖。

圖4 彈性敏感元件運動模型Fig.4 Kinetic model of the elastic sensitive parts

鑒于系統(tǒng)運動的非保守性,運用拉格朗日方程來推導(dǎo)運動方程。根據(jù)減壓器結(jié)構(gòu)及運動特點,建立運動模型時考慮以下問題:

(1)由于膜片厚度很薄,當(dāng)橫向位移相對其較大時,需要考慮大位移效應(yīng)(幾何非線性)。為此將膜片看成非線性彈簧K2。采用有限元法計算主閥膜片位移變形與反力的關(guān)系。

(2)橫向位移較大的情況下,膜片內(nèi)部應(yīng)力很大,可能超出其屈服限。因此,采用彈塑性模型來描述膜片的應(yīng)力應(yīng)變關(guān)系。

(3)安裝有初始變形,在計算時考慮膜片預(yù)變形。(4)假設(shè)漲圈提供摩擦力恒定。

(5)為了完整描述系統(tǒng)動態(tài)運動過程,引入阻尼系數(shù)c來描述系統(tǒng)阻尼力。當(dāng)考慮有摩擦力存在,阻尼力為次要地位,因此阻尼系數(shù)c很小(c=0.1)。

(6)閥芯(頂桿)與膜片調(diào)整塊存在分離面,從分離面處將系統(tǒng)運動分為上下兩部分,單獨建模。

(7)閥芯與向下限位的撞擊用非線性彈簧描述,即:接觸→壓縮→回彈→脫離,當(dāng)脫離接觸,不再有彈簧力。

(8)與膜片相連的調(diào)整塊向上運動,當(dāng)與閥芯頂桿接觸后與其一起運動。當(dāng)向下運動脫離接觸,可繼續(xù)運動,如果達到下限程,同樣用非線性彈簧描述與下限位的碰撞。

取初始位置x2=0,向下運動L5=0.8741E-3(m)即到下限位。根據(jù)拉格朗日方程,可推導(dǎo)出上部分運動方程如式(14):

式中 L3=196/K3=3.94×10-3m為副彈簧預(yù)壓縮量;f為摩擦力;Fq2為作用在閥芯上的氣動力;FX為與下部接觸力,當(dāng)上下部分分離后,FX=0。

取初始位置x1=0向下運動L4=0.7×10-3+0.874 1×10-3=1.574 1×10-3m,即到下限位。根據(jù)拉格朗日方程,可推導(dǎo)出下部分運動方程如下:

式中 L1=1 050/K1=10.38×10-3m,為主彈簧預(yù)壓縮量;Fq1為作用在膜片上的氣動力;FX為與上部接觸力,當(dāng)上下部分分離后,FX=0。

膜片上的氣動力Fq1分為兩部分:

式中 Fq1c為φ＜25 mm范圍的氣動力;Fq1w為25mm＜φ＜35mm范圍的氣動力,折算系數(shù)ξ=2.12。

減壓器彈簧-質(zhì)量-阻尼系統(tǒng)的運動方程式(14)、式(15)是非線性的,采用紐馬克方法求解。

3 計算結(jié)果與分析

計算網(wǎng)格采用非結(jié)構(gòu)網(wǎng)格,幾何建?；窘咏鎸嵲O(shè)計外形。整個減壓器表面網(wǎng)格見圖5。減壓器入口連接高壓進氣管路,高壓氣通過電磁閥進入進氣管路,電磁閥開啟時間從10～200ms不等。認為入口壓力在電磁閥開啟時間段內(nèi)呈線性增大,電磁閥開啟導(dǎo)致的快速增壓過程使得進氣管路中形成激波沖擊。為了避免反射波向入口回傳,干擾入口邊界條件,本文采用一維流動模型入口出口邊界條件。具體做法是在三維計算域的入口和出口處連接足夠長的一維等截面管路數(shù)值模型,交界面上物理量的傳遞通過延拓2層網(wǎng)格單元實現(xiàn)。該方法示意見圖5。

工作氣體為空氣,氣體總溫300 K。計算條件為上游高壓氣壓力分別為21、15、8 MPa,電磁閥開啟時間取20 ms和50 ms 2種情況。

圖5 計算網(wǎng)格及出入口一維流邊界條件Fig.5 Whole 3D model and 1D pipe flow boundary conditions in in let and outlet

圖6是高壓沖擊某時刻的馬赫數(shù)云圖。圖7為開啟持續(xù)時間20 ms不同開啟壓力下高壓腔和低壓腔壓力變化歷程,圖7中顯示,開啟壓力為15MPa和8MPa下在穩(wěn)定工作時高壓腔和低壓腔壓力均沒有振動而保持恒定值,開啟壓力為21 MPa時高壓腔和低壓腔壓力均呈現(xiàn)出有限幅值的波動。計算得到的低壓腔壓力在1.7MPa左右。圖8為開啟壓力21 MPa不同持續(xù)時間高壓腔和低壓腔壓力變化歷程。由圖8可見,開啟持續(xù)時間20 ms時高低壓腔壓力出現(xiàn)波動現(xiàn)象,而開啟時間50ms時則沒有波動,壓力呈光滑變化,這說明延長電磁閥開啟持續(xù)時間可減弱彈性元件的振動疲勞。圖9為開啟持續(xù)時間20 ms不同開啟壓力下的閥芯與膜片位移曲線,可見膜片硬芯與閥芯頂桿絕大多數(shù)時間都是結(jié)合在一起運動,只有開啟壓力大于15 MPa時在工作初期有短暫分離。開啟壓力為15 MPa和8 MPa且在穩(wěn)定工作時閥芯開度穩(wěn)定在各自的恒定值,開啟壓力為21 MPa時,則在較長時間內(nèi)閥芯開度處于小幅振動狀態(tài)。圖10為開啟壓力21MPa、持續(xù)時間50ms的位移曲線,圖10中顯示,同樣的開啟壓力,將持續(xù)時間延長到50 ms就能避免閥芯開度的振動,使開度穩(wěn)定在恒定值。

圖6 高壓沖擊某時刻的馬赫數(shù)云圖Fig.6 Transient Mach contours during the impact of high p ressure gas

圖7 不同開啟壓力下高低壓腔壓力變化歷程對比Fig.7 History of p ressure in high and low p ressure room at different upstream p ressures

圖8 入口21MPa不同持續(xù)時間高低壓腔壓力變化歷程Fig.8 History of pressure in high and low pressure room at different upstream p ressurizing duration under 21MPa

圖9 不同開啟壓力下持續(xù)時間20m s閥芯與膜片位移Fig.9 Displacements ofmoving assembly at different upstream pressures and 20ms p ressurizing duration

圖10 開啟壓力21MPa持續(xù)時間50ms閥芯與膜片位移Fig.10 Displacements of moving assembly at 21M Pa upstream pressure and 50ms pressurizing duration

4 結(jié)束語

采用流固耦合數(shù)值模擬技術(shù)研究了減壓器工作過程的動態(tài)特性,得到了減壓器在不同上游壓力和不同開啟時間條件下彈性敏感元件的運動規(guī)律。結(jié)果表明,開啟壓力為21 MPa時,如果開啟持續(xù)時間為20 ms,則高壓腔和低壓腔壓力以及彈性元件的運動均呈現(xiàn)出有限幅值的波動;而開啟時間50 ms時則沒有波動,壓力呈光滑變化趨勢,閥芯開度穩(wěn)定在恒定值。這說明延長電磁閥開啟持續(xù)時間可減弱彈性元件的振動疲勞。

本文提出的虛擬擋板技術(shù)成功模擬了活門由閉到開的瞬態(tài)過程,表明該技術(shù)模擬物體間由零距離開始分離的實際過程是成功的,且對最小網(wǎng)格尺度沒有限制。

[1] Fejtek I,Waller G,Wong R.Computational study of the flowfield of an aircraft outflow valve[R].AIAA 2005-4843.

[2] Ahuja V,Hosangadi A,Shipman J,etal.Simulations of instabilities in com plex valve and feed systems[R].AIAA 2006-4758.

[3] Cavallo P A,Hosangadi A,Ahuja V.Transient simulations of valvemotion in cryogenic systems[R].AIAA 2005-5152.

[4] 徐克鵬,蔡虎,崔永強,等.大型汽輪機主汽調(diào)節(jié)閥的實驗與數(shù)值分析[J].動力工程,2003,23(6):2785-2789.

[5] 王冬梅,陶正良,賈青,等.高壓蒸汽閥門內(nèi)流場的二維數(shù)值模擬及流動特性分析[J].動力工程,2004,24(5):690-697.

[6] 樊希葆,謝水波,陳澤昂,等.水泵站停泵水錘數(shù)值計算理論研究——模型求解與應(yīng)用[J].南華大學(xué)學(xué)報(自然科學(xué)版),2005,19(2):6-9.

[7] 聶萬勝,陳新華,戴德海,等.姿控推進系統(tǒng)發(fā)動機關(guān)機的管路瞬變特性[J].推進技術(shù),2003,24(1):6-8.

[8] Barth T J,Jespersen DC.The design and application of upwind schemes on unstructured meshes[R].AIAA 89-0366.

[9] 郭正.包含運動邊界的多體非定常流場數(shù)值模擬方法研究[D].長沙:國防科技大學(xué),2002.