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

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

0 引言

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

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

1 數(shù)值方法

1.1 控制方程

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

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

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

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

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

式中 xt、yt、zt分別是網(wǎng)格運(yùn)動(dòng)速度的3個(gè)分量。

1.2 離散方法

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

式中 m表示時(shí)間層;Sk是第k個(gè)積分面的矢量面積。

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

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

式中 φij是限制器。

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

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

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

1.3 對(duì)虛擬擋板面的積分

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

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

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

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

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

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

對(duì)流項(xiàng)Fb按下面的方法計(jì)算:

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

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

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

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

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

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

2 彈性敏感元件的運(yùn)動(dòng)方程

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

圖4 彈性敏感元件運(yùn)動(dòng)模型Fig.4 Kinetic model of the elastic sensitive parts

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

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

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

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

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

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

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

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

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

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

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

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

膜片上的氣動(dòng)力Fq1分為兩部分:

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)束語

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

本文提出的虛擬擋板技術(shù)成功模擬了活門由閉到開的瞬態(tài)過程,表明該技術(shù)模擬物體間由零距離開始分離的實(shí)際過程是成功的,且對(duì)最小網(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] 徐克鵬,蔡虎,崔永強(qiáng),等.大型汽輪機(jī)主汽調(diào)節(jié)閥的實(shí)驗(yàn)與數(shù)值分析[J].動(dòng)力工程,2003,23(6):2785-2789.

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

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

[7] 聶萬勝,陳新華,戴德海,等.姿控推進(jìn)系統(tǒng)發(fā)動(dòng)機(jī)關(guān)機(jī)的管路瞬變特性[J].推進(jìn)技術(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] 郭正.包含運(yùn)動(dòng)邊界的多體非定常流場(chǎng)數(shù)值模擬方法研究[D].長(zhǎng)沙:國(guó)防科技大學(xué),2002.