采用流固耦合方法數(shù)值分析減壓器動態(tài)響應①

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

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

0 引言

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

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

1 數(shù)值方法

1.1 控制方程

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 計算結果與分析

計算網(wǎng)格采用非結構網(wǎng)格,幾何建模基本接近真實設計外形。整個減壓器表面網(wǎng)格見圖5。減壓器入口連接高壓進氣管路,高壓氣通過電磁閥進入進氣管路,電磁閥開啟時間從10～200ms不等。認為入口壓力在電磁閥開啟時間段內(nèi)呈線性增大,電磁閥開啟導致的快速增壓過程使得進氣管路中形成激波沖擊。為了避免反射波向入口回傳,干擾入口邊界條件,本文采用一維流動模型入口出口邊界條件。具體做法是在三維計算域的入口和出口處連接足夠長的一維等截面管路數(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ù)時間都是結合在一起運動,只有開啟壓力大于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 結束語

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

本文提出的虛擬擋板技術成功模擬了活門由閉到開的瞬態(tài)過程,表明該技術模擬物體間由零距離開始分離的實際過程是成功的,且對最小網(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ù)值計算理論研究——模型求解與應用[J].南華大學學報(自然科學版),2005,19(2):6-9.

[7] 聶萬勝,陳新華,戴德海,等.姿控推進系統(tǒng)發(fā)動機關機的管路瞬變特性[J].推進技術,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].長沙:國防科技大學,2002.