管路中常溫推進劑的兩相充填特性仿真

任孝文,李平,陳宏玉,周晨初

1.西安航天動力研究所 液體火箭發動機技術重點實驗室,西安 710100

2.航天推進技術研究院,西安 710100

在航天推進系統的起動過程中,液體推進劑供應管路中氣體的存在通常會使得系統動力特性變得難以預測。氣體的壓縮性具有貯存和釋放能量的功能,特別對于液體推進器中的快速充填等瞬態流動而言,這一特點增大了數學建模的難度。根據管路末端結構的不同,充填過程可分為受限式充填和排氣式充填,其中前者是對末端閥門關閉的閉端管路進行充填,而后者是對末端敞口或末端接有節流孔的管路進行充填。

針對受限式充填,國內外已有較多的學者進行了研究：陳宏玉、秦艷平等在處理氣液交界面的移動時使用了一維坐標變換法,Bandyopadhyay和Majumdar采用通用流體系統仿真程序(Generalized Fluid System Simulation Package,GFSSP)進行數值仿真,上述學者使用的充填模型均只能對液體區域進行分布參數計算,預存氣體部分使用了集中參數的方法；在兩相分布參數模型方面,程謀森和張育林使用特征線有限差分方法,對液體與氣體計算區域分別采取不同的空間網格,實現了氣液控制方程在時間上的同步性和穩定性；劉昆和張育林將每個網格單元劃分為未充填、充填以及充滿3種狀態,建立了一維有限元狀態變量模型；本文作者團隊在文獻[8]模型的基礎上基于模塊化建模思想建立了一維有限體積的受限式充填模型；Zhou L等使用FLUENT 軟件的流體體積 (Volume of Fluid,VOF)法進行了二維和三維的充填仿真,并針對夾帶有2個氣柱的液體進行了閉端管路的充填仿真和實驗研究,其建立的一維模型使用特征線方法,在氣液界面附近忽略了液體的慣性和壓力損失。

針對管路末端排氣的快速充填過程,已有學者通過實驗手段研究了其中的瞬態流動現象：該過程包含兩個典型的流動特征,即氣體受液體充填鋒面的壓縮而經末端節流孔快速排放的過程和液體充填鋒面到達節流孔產生的水擊振蕩過程；Bucur等實驗研究了貯箱壓力、節流孔徑及充填長度對充填過程的影響,并通過無量綱分析提出了一個兩方程模型和相應的適用準則以預測不同參數下的最大水擊壓力峰值。在數值仿真方面,由于排氣式充填需要同時考慮液體和氣體的壓縮性以及氣體排盡之后控制方程由氣體向液體的切換,是一個復雜的非穩態兩相流動問題,因此針對該過程的仿真建模研究較少：僅有的仿真工作也多集中在氣體未完全排盡時的流動模擬,早期的有Zhou F等建立的數學模型,該模型只能計算液體充填鋒面到達節流孔之前的流動壓力變化；在Tijsseling等近期建立的模型中,仿真程序設定在氣體體積達到初始體積的5%時停止計算,因此其給出的結果未能覆蓋液體到達節流孔時的水擊過程。Zhou L 等針對夾帶有2個氣柱的液體建立的充填模型同樣無法計算管內氣體排盡之后的流動。目前,Modelica標準模型庫中只有純液相的管路模型,僅能計算高壓液體對低壓液體管路的充填,因此該模型難以真實模擬推進劑充填過程中的氣液兩相和慣性充填過程。對AMEsim 商業庫中管路充填模型的控制方程進行分析,其充填過程的壓力由理想氣體狀態方程計算得到,未給出氣體排盡后流動向液相轉換時壓力的計算方法。有鑒于此,針對管路末端排氣的快速充填過程進行相關的仿真建模工作是十分必要的。

為拓展液體火箭發動機瞬態特性模塊化通用仿真模型庫對兩相充填的仿真能力,同時為了彌補上述模型的缺陷,本文基于Modelica模塊化建模思想,在MWorks平臺上自主開發了一維有限體積的兩相充填管路組件模型。不同于公開文獻中常將液體假設為不可壓縮流體,本文模型同時考慮了氣體和液體的壓縮性,在使用VOF法捕捉氣液界面的基礎上采用等效流容方程計算壓力,避免了特征線法在氣液界面位置復雜的插值格式；此外針對存在多個氣液界面的流動問題,例如夾帶有氣柱的液體充填過程,利用本文模型的模塊化優勢可快速進行系統模型的建立和仿真；與本文模型類似的工作在公開文獻中鮮有報道。利用建立的兩相充填管路模型,在實驗驗證的基礎上,本文研究了氣體壓縮性、節流孔徑對充填特性的影響,并對管路中夾帶有單個氣柱的推進劑充填過程進行了仿真研究。

1 建模方法

本文所研究的流體均為常溫流體,易相變的低溫流體不在考慮范圍之內,同時不考慮流體與壁面的傳熱。假設充填過程中液體與氣體不相溶且兩相之間具有明確的分界面,液體為弱可壓縮流體,氣體為完全可壓縮流體,液體溫度保持恒定,氣體溫度則隨壓力變化。且在任一瞬間,管流的壓強、密度、溫度和速度在任一截面上是均一的且只隨管長方向變化,管路橫截面積為定值,忽略流體重力的影響。

數值離散采用一維有限體積方法,使用圖1所示的狀態網格單元與速度網格單元交錯的方式將管路進行等分,圖中：為網格索引,p 為流體壓力,為流體密度,u為流體速度,管路兩端半個狀態網格單元內流體的速度取上下游邊界的速度,分別用下標“a”和“b”表示。狀態網格單元計算流體的壓力、密度、溫度以及充填率；充填鋒面位置通過VOF法計算得到。

圖1 狀態單元(實線)與速度單元(虛線)的交錯網格Fig.1 Staggered grids of state cells(solid lines)and velocity cells(dashed lines)

1.1 控制方程

對于單相液體,控制方程為

1)質量守恒方程：

式中：為流體密度；為速度。

2)動量守恒方程：

式中：為壓力；為流 阻。

3)狀態方程：=(,),為流體溫度,根據鏈式求導法則,有

忽略液體的溫度變化,則密度只是壓力的函數,式(3)變為

因此,質量守恒方程轉變為如下壓力方程：

對于考慮氣體壓縮性的兩相充填過程,可將氣液混合物看作單相流體,則仍可使用上述壓力方程及動量方程描述流動過程,不同點及關鍵點在于需要將氣液兩相混合物的密度和聲速與充填率進行關聯。

假設充填鋒面垂直于管路軸向,且將氣液兩相混合流體近似認為弱可壓非均質流體,則有如下充填率方程：

式中：為液體的充填率,其定義為控制體內液體所占的體積分數。

根據等效流容方程,氣液兩相混合物的聲速計算式為

式中：為考慮流固耦合作用的液體聲速,其計算表達式為

其中：為純液體的聲速,對于水,=1 200 m/s；為管壁彈性模量；為管壁厚度；為管路內徑；為考慮管路支撐影響的修正量,具體取值可參見文獻[31]。

在進行兩相充填過程的仿真時,管網中的閥門、節流圈等節流元件同樣需要具有處理兩相流動的能力。通常,節流元件長度相比于管路長度較小,因此可忽略充填時間,認為節流元件在一瞬間即充滿液體。以上游管路末端的壓力作為節流元件的入口壓力,出口壓力由下游組件給定,則根據進出口壓力可以分別計算得到流體均為液體時的流量和均為氣體時的流量。

節流元件中氣液之間的流量切換可根據上游兩相充填管路模型的充填率確定：

式中：()為關于充填率的函數。

簡單地,()可隨充填率線性變化：

本文假設氣液界面嚴格垂直于流動方向,因此使用階躍函數：

1.2 對流項離散方法

本文的數值計算基于Modelica語言的面向對象建模的MWorks仿真平臺,由于該平臺能夠使用時間步長自適應的DASSL算法求解常微分方程組,因此建模過程中只需要對空間項進行數值離散,將分布參數的偏微分方程組轉化為網格單元內的常微分方程組。

由于充填鋒面垂直于管路軸向的假設,初始時刻充填率隨管路沿程具有如圖2所示的分布形式,圖中：橫坐標/表示充填鋒面位置相對管路長度的比值,=1表示充滿液體,=0表示充滿氣體。

圖2 初始時刻充填率分布Fig.2 Distribution of filling rates at initial time

可以看到,初始時刻液體充填率隨管路的沿程分布是一個典型的初始間斷,因此數值計算的關鍵在于準確捕捉間斷的傳播。為此,控制方程中對流項的數值離散采用總變差減小格式(Total Variation Diminishing,TVD),以保證間斷的傳播具有較小的耗散性和較好的迎風性。

在以上數學模型的基礎上,本文基于MWorks平臺,使用Modelica語言開發了常溫液體兩相充填管路的組件模型,將一階迎風格式、二階迎風格式以及7種總變差減小格式全部集成在了該組件模型中,以方便對離散格式進行比較和選擇,圖3所示為組件模型的參數設置界面。

圖3 組件參數設置界面Fig.3 GUI of parameter settings of component

2 數值格式比較及驗證

2.1 不同離散格式數值計算結果對比

針對圖4所示典型的貯箱-管路系統,進行常溫液體對管路充填過程的仿真計算。液體工質選用水；管路內徑=10 mm,貯箱壓力保持在2 MPa,貯箱出口閥門流量系數為0.707；管路末端節流孔內徑=8 mm,流量系數為0.65；初始時刻,管內氣體與大氣相通且壓力為0.1 MPa,管路預先充液長度=0 m,待充填長度=5 m；0.1 s時刻閥門開啟,閥門從關閉狀態到完全打開狀態用時5 ms；管路等分為50段。

圖4 典型貯箱-管路系統示意圖Fig.4 Schematic of typical tank-pipe system

由于二階迎風格式存在較大的色散容易導致計算發散,因此在進行不同離散格式的對比時只給出了一階迎風格式與7種TVD 格式的計算結果,如圖5～圖7所示,分別為管路中部第25個網格的充填率隨時間的變化、管路末端壓力及溫度隨時間的變化。

圖5 管路中部第25個網格內的充填率Fig.5 Filling rate of the 25th cell in middle of pipe

圖7 管路末端溫度Fig.7 Temperature at downstream end of pipe

如圖5所示,與TVD 格式相比較,在計算充填鋒面這一間斷的傳播過程中,一階迎風格式具有嚴重的數值耗散,這在圖6中表現為管路末端壓力由大氣壓平滑過渡至穩態壓力,而事實上由于節流孔的存在,在液體充填鋒面到達節流孔位置時會出現水擊現象；在上游貯箱壓力和管路內徑一定的條件下,水擊壓力的大小取決于節流孔的大小,在本節算例節流孔徑給定的情況下,水擊壓力峰值的不同完全是由于數值耗散引起的,根據快速充填過程中充填鋒面垂直于管路軸向的假設,數值耗散越小則計算結果越準確。

圖6 管路末端壓力Fig.6 Pressure at downstream end of pipe

圖7所示為充填過程中管路末端溫度的變化,整個過程溫度出現了2次上升,第1個溫度峰在0.1 s附近,是由于閥門瞬間開啟導致液體開始充填所引起的,第2個較高的溫度峰是由于液體充填鋒面對氣體的壓縮而引起的；同樣可以看到,由于一階迎風格式較大的耗散性,其計算得到的溫度值偏低。

表1給出了一階迎風格式與7種TVD 格式的計算資源消耗對比,其中以一階迎風格式的計算時間為參照,并定義其為1個CPU 時間。綜合比較圖5～圖7及表1可以看到,各種離散格式中,具有TVD 性質的SUPERBEE格式數值耗散最小,QUICK 格式次之；Van Albada、Min-Mod、QUICK 以及UMIST 格式的CPU 消耗時間與一階迎風格式相差不大,Van Leer格式的計算耗時最大,需要35個CPU 時間,這可能是因為Modelica語言在編譯絕對值時會產生大量的狀態事件而占用了較多的計算資源；綜合以上數值耗散和CPU 耗時兩方面的因素,在同時滿足計算準確性和實時性的要求下,推薦使用具有TVD 性質的QUICK 格式進行管路的兩相充填計算。

表1 不同離散格式CPU 耗時對比Table 1 Comparison of CPU time among different schemes

2.2 模型驗證

Bucur等針對圖4所示的貯箱-管路充填過程進行了實驗研究。選取其中一個實驗工況對本文建立的兩相充填模型進行有效性確認,參數設置如下：管路內徑=39 mm,管壁厚度=10 mm,管路材料為石英玻璃,彈性模量=2.5 GPa；貯箱壓力為0.4 MPa,末端節流孔內徑=5 mm；初始時刻,管內氣體與大氣相通,壓力為0.1 MPa,管路預先充液長度=9.56 m,待充填長度=0.55 m；0.3 s時刻閥門開啟,閥門從關閉狀態到完全打開狀態用時20 ms。

圖8給出了圖4中閥門Valve出口位置壓力的計算結果與實驗數據(EXP)的對比,圖中為充填鋒面位置相對充填長度的比值,計算結果能較好地反映管路末端排氣的充填水擊過程；圖中第1個壓力峰值由氣體排盡后充填鋒面第1次到達節流孔位置時產生的水擊引起,第2個壓力峰值對應于圖中先減小后增大的凹坑位置,這表明在第1個水擊壓力后充填鋒面發生了倒流,之后又重新到達節流孔位置時產生了第2 次水擊,液體開始流出節流孔,直至到達穩態。

圖8 模型驗證Fig.8 Validation of present model

由于氣液界面嚴格垂直于管路軸向的假設,在氣體完全排盡后,流體的聲速完全切換為了純液相的聲速,由此造成了計算得到的水擊振蕩頻率大于實驗數據的結果；實際在充填過程中,氣體不會立刻排盡,會有部分氣體以氣泡形式殘存于液體中導致流體聲速減小,這是造成實驗數據中水擊振蕩頻率相對較小的原因。

3 仿真及分析

3.1 與集中參數充填模型的對比

采用本文開發的常溫液體兩相充填管路的組件模型與液體火箭發動機瞬態特性模塊化通用仿真模型庫中已有的集中參數模型分別對圖4所示的貯箱-管路系統進行充填過程的仿真計算。根據2.1 節的結論,對流項離散選用TVD_QUICK 格式；集中參數模型忽略初始時刻管路中原有氣體的影響,并假設充填過程中氣液界面氣體一側的壓力始終與大氣環境一致,為0.1 MPa；系統模型及參數設置同2.1節。

圖9所示為不同節流孔徑與管路內徑之比/(簡稱節流孔徑比)下,2種管路充填模型的計算結果對比。隨著節流孔徑比的減小,2種模型計算結果的差距逐漸增大。由于集中參數模型不考慮管路內預存氣體的影響,因此計算得到的水擊壓力峰值偏大,這一特點在節流孔徑比小于0.5時尤其明顯。造成這一差別的原因在于實際充填過程中,管路內原有的氣體會對液體的順利充填造成一定的阻礙作用；特別對于管路末端節流孔徑較小的情況,在向外界或下游排放受阻的同時,氣體受到液體充填鋒面的壓縮,導致氣體壓力增大,進一步減弱了液體的充填速度。根據茹科夫斯基最大水擊壓強公式=,充填速度減小后,當液體與末端節流孔撞擊時,產生的水擊壓強也將減小。綜上所述,由于考慮了管路內預存氣體壓縮性的影響,相比于集中參數模型,本文建立的兩相充填管路模型的仿真結果更符合實際過程。

圖9 不同節流孔徑比下的管路末端壓力Fig.9 Pressure at downstream end of pipe under different d/D

3.2 節流孔徑對充填水擊的影響

為進一步探索管路末端節流孔徑對充填水擊的影響,以圖4 所示的貯箱-管路系統為研究對象,采用本文建立的常溫液體兩相充填管路的組件模型對不同節流孔徑比的管路結構進行充填過程的仿真計算。管路預先充液長度=5 m,待充填長度=5 m。保持待充填長度與預先充液長度的比值/=1不變,末端節流孔內徑在0～39 mm 范圍內變化,對應于節流孔徑比在0～1.0之間。0.1 s時刻閥門開啟,閥門從關閉狀態到完全打開狀態用時5 ms。其余參數設置與2.2節相同。

從共計進行的35次仿真計算結果來看,在不同的節流孔徑比下,根據充填過程中表現出的壓力峰值大小以及壓力峰值出現時刻與液體充填率的對應關系,管路內的壓力振蕩形式可以劃分為3種模式,分別是水擊效應忽略模式、水擊效應微弱模式以及水擊效應主導模式。下面將分別展開討論。

3.2.1 水擊效應忽略模式

對于管路末端沒有氣體排出或節流孔徑比足夠小的情況,液體充填過程中管路內的壓力振蕩周期較長,同時沒有表現出明顯的水擊壓力,如圖10所示。此時管路內的壓力振蕩主要由預存氣體的壓縮性造成,在圖中具體表現為：管路內的壓力振蕩集中于充填率小于1的時間段內。

特別地,對于/=0 的極限情況,即對預存氣體閉端管路的充填,如圖10(a)所示,液體始終無法充滿管路,液體充填鋒面隨壓力波動往復振蕩,在上游貯箱壓力和管壁摩擦的作用下,液體充填率最終穩定在0.6左右。隨著節流孔徑比的增大,充填過程中管路內原有的氣體經節流孔排出,但是在節流孔徑比足夠小的情況下,管路內氣體量短時間內無法排放完全,壓力振蕩同樣主要由氣體的緩沖作用主導。然而,隨著液體充填率的升高,管路內氣體量減小,充填過程的壓力峰值減小；液體慣性的影響變大,壓力振蕩頻率升高,如圖10(b)及圖10(c)所示。

圖10 水擊效應忽略模式仿真結果Fig.10 Simulation results of type of negligible water hammer effect

由于氣體壓縮緩沖作用占主導,充填過程受阻,液體充填至節流孔位置時速度較小,因此在充滿的瞬間并未出現水擊現象。

綜上,在圖10代表的壓力振蕩形式中,氣體的緩沖作用起主導作用,水擊壓力可以忽略,因此稱為水擊效應忽略模式。經計算,在本文討論參數范圍內,節流孔徑比在0～0.110之間時,管路中的兩相充填過程處于水擊效應忽略模式。

3.2.2 水擊效應微弱模式

在充填過程中,當氣體緩沖作用與水擊效應在管路內均有體現時,管內的壓力振蕩形式可稱為水擊效應微弱模式,如圖11所示。從圖中可以看到,以充填率=1為分界線,管內壓力變化可以分為2個階段：充填率＜1的階段內,在充填過程的早期,管內氣體量較多,管內壓力變化呈現出周期較長的振蕩形式,與水擊效應忽略模式類似,同時可以看到隨著節流孔徑比的增大,該階段的振蕩周期逐漸縮短,壓力峰值也明顯減小。之后,當氣體排放完全,管內液體充填率=1的瞬間,液柱與節流孔固體部分撞擊產生水擊壓力。與氣體導致的壓力振蕩周期相比,液體與節流孔撞擊產生的水擊壓力振蕩周期明顯縮短。

圖11 水擊效應微弱模式仿真結果Fig.11 Simulation results of type of weak water hammer effect

經計算,節流孔徑比在0.110～0.174 之間時,管路中的兩相充填過程處于圖11所示的水擊效應微弱模式。

一個值得注意的事實是,2.2節中圖8算例對應的節流孔徑比同樣為0.128,理應與圖11(a)表現出相似的充填性質,但兩者結果卻相差較大。其原因在于兩者的物理模型存在差別：/是不同的,前者為0.0575,后者為1。在節流孔徑比一定的條件下,/直接反映了液體慣性對充填過程的影響,該比值越小表明液體慣性的影響越大而氣體緩沖作用相對越小。圖11算例中,氣體緩沖作用遠大于圖8算例,充填過程氣阻明顯,當液體充填至管路末端位置時產生的水擊壓力峰值較小(圖11中第2個壓力波峰),之后流動未出現明顯的倒流(圖11中第3個壓力波峰很小),因此充填率曲線沒有如圖8那樣出現凹坑。由此可見,/也是影響兩相充填過程的重要因素,是后續工作中值得細致研究的內容。

3.2.3 水擊效應主導模式

在本文討論參數范圍內,當節流孔徑比超過0.174時,管路中的兩相充填過程處于水擊效應主導模式,如圖12所示。在該孔徑比范圍下,節流孔對氣體的流阻較小,充填開始后管內預存的氣體在液體鋒面的推動作用下經節流孔快速排出,氣體不再表現出如水擊效應忽略模式和水擊效應微弱模式中的氣體緩沖作用,因此當液體充填至節流孔位置時具有較大的速度,導致與節流孔撞擊而產生較大的水擊壓力。

對比圖12(a)與圖12(b),當/=0.179時,在充填率＜1的階段內,氣體對液體的充填過程有些許阻力,管路末端壓力上升較為平滑；而在節流孔徑比進一步增大至0.231后,氣體排放非常快,幾乎對液體的充填過程沒有阻力,導致液體快速撞擊至節流孔。從充填時間也可以看到,隨著節流孔徑比的增大,液體充填率達到1所用的時間逐漸縮短。

圖12 水擊效應主導模式仿真結果Fig.12 Simulation results of type of dominant water hammer effect

3.2.4 壓力峰值變化

在工程實踐中,通常對充填過程中出現的水擊壓力峰值較為關心,因此將不同節流孔徑比下的壓力峰值()仿真結果總結對比,可以得到如圖13所示的曲線。從圖中可以看到,在貯箱壓力一定的條件下,隨著節流孔徑比的增大,充填水擊壓力峰值呈現出較為復雜的變化規律。

為便于分析,將圖13中曲線按照3種壓力振蕩模式進行劃分。節流孔徑比在0～0.110范圍內時,管路兩相充填過程中的壓力振蕩屬于水擊效應忽略模式,可以看到,在該模式下,管路內的壓力峰值隨節流孔徑比的增大而逐漸減小。如3.2.1節所述,在水擊效應忽略模式下,管內的壓力變化主要由氣體的壓縮性造成,在液體充填鋒面的推動下,隨著節流孔徑比的增大,氣體向下游或外界排放受到的阻力減小,因此氣體受壓程度也隨之減小,從而導致管內壓力峰值減小。

圖13 充填過程管路內壓力峰值隨d/D 的變化Fig.13 Pressure peak against d/D of filling process in pipe

當節流孔徑比在0.110～0.174范圍內時,管內壓力振蕩處于水擊效應微弱模式,此時水擊壓力峰值隨節流孔徑比的增大呈現出先減小后升高的趨勢。該模式下,壓力峰值主要由液體與節流孔撞擊產生的水擊造成,但是氣體壓縮性的影響不可忽略,兩者孰更占優直接影響著管路內壓力峰值的變化規律。將節流孔徑比范圍進行細分,在0.110～0.141范圍內,充填過程的壓力峰值隨節流孔徑比的增大而呈減小趨勢,該階段氣體壓縮性影響較大,如圖11(a)所示,因此壓力峰值變化規律與水擊效應忽略模式類似；在0.141～0.174之間,液體與節流孔撞擊產生的水擊影響較大,隨著節流孔徑比的增大,液體充填速度加快,導致水擊壓力逐漸上升。

節流孔徑比進一步增大至0.174 以上時,圖13中管內壓力峰值隨節流孔徑比的增大而先升高后降低,在/=0.231時壓力峰值達到最大值1.3 MPa。在/=0.231以前,充填過程仍然會受到氣阻的影響,隨著節流孔徑比的增大,氣阻減小,充填過程加速,導致充填水擊壓力升高；在/=0.231以后,氣體排放非常快,幾乎對液體的充填過程沒有阻力,導致液體快速撞擊至節流孔,但是隨著節流孔徑比的增大,節流孔固體區域減小,節流孔對液體的滯止作用逐漸減小,從而導致水擊壓力隨節流孔徑比的增大而呈減小趨勢。

當節流孔徑比大于0.641后,節流孔徑與管路內徑接近,節流孔流阻大幅降低,管路內最大壓力開始小于箱壓,這表明在該對應結構下的充填過程沒有水擊發生。

3.3 管路中夾帶單個氣柱的推進劑充填過程仿真

由于液體火箭發動機各部件的結構布局較為緊湊,連接其中的管路不可避免地會出現彎折。在發動機起動時,若管路彎折位置積存有一定量的氣體,則閥門打開后就會出現夾氣液體對下游管路的充填過程。在該過程中,液柱與氣柱的相互作用會導致管路內產生較大的充填水擊與壓力振蕩,進而對推進劑的穩定供應造成一定的負面影響。

為研究夾氣流體對充填過程的影響,使用本文開發的氣液兩相充填管路模型對如圖14所示的典型貯箱-管路系統進行仿真,圖中管路內藍色部分為液體,白色部分為氣柱,圖中示出了多個液體段和氣柱；鑒于工程中彎管數量有限,本節選取單個氣柱的情況進行仿真研究,且暫不考慮管路彎曲對流動的影響。

圖14 夾氣流體管路充填示意圖Fig.14 Schematic of filling process of fluid with entrapped air pockets

圖15所示為管路中存在單個氣柱的充填仿真模型,其中Pipe 0為純液相管路模型,Pipe 1～Pipe 3為本文建立的兩相充填管路模型,對流項離散格式使用TVD_QUICK 格式。Valve 1前與貯箱連接的液體管路Pipe 0 及兩相充填管路Pipe 2管路長度均為2 m,Pipe 1長度為0.5 m,待充填管路Pipe 3 長0.55 m；管路內徑均為39 mm。根據網格獨立性測試,Pipe 0、Pipe 2網格數均為20,Pipe 1及Pipe 3網格數為5。閥門Valve 1流量系數設置為0.95,流通面積與管路相同；閥門Valve 2流道直徑為5 mm,Pipe 3管路末端節流圈流道直徑為8.5 mm。

圖15 管路中存在單個氣柱的充填仿真模型Fig.15 Simulation model for filling process of pipe with one entrapped air pocket

計算中貯箱為恒壓源,且壓力穩定在1 MPa,Pipe 3管路末端與大氣相通,壓力穩定在0.1 MPa。初始時刻,Valve 1閥前管路Pipe 0內流體壓力為1 MPa；Pipe 1中充滿標準大氣壓的空氣,即將其每個網格的充填率設置為0；Pipe 2為液體管路,每個網格的充填率設置為1,液體初始壓力為0.1 MPa；Pipe 3為待充填管路,其末端與大氣相通,充填率設置為0。可以看到Valve 1閥后管路為一段夾帶有氣體的液柱,計算起始時刻,Valve 1閥后管路內氣體與液體處于平衡狀態。Valve 2為常開狀態,0.5 s時刻Valve 1打開。

圖16～圖22給出了圖15仿真模型的計算結果。當0.5 s閥門Valve 1打開后,由于閥門前后的壓差作用,管路內液體和氣體開始向下游流動,在閥門Valve 2、Pipe 3末端節流圈Orifice的節流作用下,以及Pipe 1 中氣體壓縮性的影響下,如圖16中黑色曲線所示,管路內出現了較大的充填水擊壓力和明顯的壓力振蕩。為了更直觀地表現Pipe 1中氣柱的壓縮性,將Pipe 1設置為純液管路進行計算對比,即將其每個網格的初始充填率設置為1,得到的計算結果見圖16紅色曲線。當Pipe 1為純液管路,閥門Valve 1打開時,本質是上游高壓液體對Pipe 1內低壓液體的推動過程,由于閥門Valve 2的節流作用以及Pipe 1內液體的慣性作用,管路內產生了較小的水擊壓力和較高的振蕩頻率。由此可知,夾帶氣柱的液體在充填過程中產生的較大的水擊壓力由氣柱受到壓縮造成,同時由于夾氣會導致聲速降低,因此該情況下振蕩頻率減小。

圖16 管路Pipe 0出口位置壓力曲線Fig.16 Pressure history on outlet of Pipe 0

圖17所示為液體對Pipe 1管路的充填率以及Pipe 1管路出口位置壓力隨時間的變化。由于初始時刻Pipe 1 管路內充滿了標準大氣壓的空氣,0.5 s閥門Valve 1打開后,液體開始充填Pipe 1,同時由于閥門Valve 2流道通徑很小,有很大的節流作用,Pipe 2管路內的液體在Valve 2位置受到較大的流阻,導致Pipe 1管路內的空氣既受到上游液體的推動又受到了下游液體的阻力,氣體受到壓縮導致壓力急劇上升,此時Pipe 1管路出口壓力對應于圖17中的第1個壓力峰值,該壓力峰值達4.87 MPa,幾乎是貯箱壓力的5倍；結合充填率曲線分析也可知該壓力峰值為氣體壓縮導致,因為此時管路內液體的充填率未達到1,即液體并未充滿。隨著充填過程的進行,Pipe 1管路氣體逐漸進入Pipe 2,直至Pipe 1管內充滿液體,如圖17中充填率曲線所示,在1.0 s左右,Pipe 1管內充填率達到1。同時結合圖19,1.0 s 之前,Pipe 1 管路內的壓力波動主要由Pipe 1內殘存的氣柱和進入Pipe 2 管路內氣柱的壓縮性造成,1.0 s之后的壓力波動則由全部進入下游Pipe 2管路的氣柱造成。

圖17 Pipe 1管路充填率及出口位置壓力曲線Fig.17 Filling rate and pressure history on outlet of Pipe 1

圖19 Pipe 2管路充填率及出口位置壓力曲線Fig.19 Filling rate and pressure history on outlet of Pipe 2

圖18所示為Pipe 1管路入口的液體流量隨時間的變化曲線,可以看到Valve 1閥門開啟的瞬間,在前后巨大壓差的作用下,液體以瞬時較高的流量開始向下游流動,而由于Pipe 1管路內氣柱的存在,其表現為一個貯存和釋放能量的彈簧,導致管路內出現流量和壓力振蕩,其頻率大約在10 Hz左右。

圖18 Pipe 1管路入口流量Fig.18 Mass flow rate of Pipe 1 inlet

圖19所示為Pipe 2管路出口的壓力變化以及液體充填率曲線,0.5 s時刻閥門Valve 1開啟的瞬間,Pipe 1管路內的氣體即開始進入Pipe 2中,表現在圖中Pipe 2管路的充填率從1開始下降,如前所述,此時管路內的氣柱受到上下游液體的擠壓而產生高壓,這一高壓迅速傳播至Pipe 2管路出口,并與液體撞擊至Valve 2位置的水擊壓力疊加,導致產生了5.6 MPa的壓力峰值,對應于圖19中第1個壓力峰值,可以看到該壓力峰值要高于圖17中Pipe 1出口位置的壓力峰值,高出的壓力值即為Pipe 2 管路內液體撞擊至Valve 2位置產生的水擊壓力。

1.0 s左右,Pipe 1管內的氣柱全部進入Pipe 2管內,表現為圖19中充填率曲線在1.0 s后不再有較大的波動以及圖20中氣體流量在1.0 s左右降為0,隨著流動向下游的傳播,氣柱逐漸向下游移動,至2.7 s左右氣柱排出Pipe 2管路,此時壓力曲線出現一個較小的凹坑,該凹坑在圖16和圖17中同樣可以觀察到。從圖18和圖20可以看到,該凹坑由液體流量的突然上升而造成,即壓力勢能部分轉化為了動能。

圖20 Pipe 2入口液體和氣體流量Fig.20 Mass flow rate of liquid and gas on inlet of Pipe 2

圖21所示為Pipe 3管路充填率及出口位置壓力曲線,圖22給出了Pipe 3管路入口的液體和氣體流量隨時間的變化曲線。初始時刻,Pipe 3為待充填管路,其通過末端節流圈Orifice與外界大氣相通,可以看到閥門Valve 1在0.5 s開啟后,即有管路Pipe 2中的液體進入Pipe 3中,對應于圖21中充填率從0開始上升以及圖22中液體流量的瞬時增大。由于Pipe 3管路末端節流圈的流阻較小,管路內預存氣體向外界的排放較為順利,因此當Pipe 2 中液體充填Pipe 3 管路時,并未對其中預存的氣體產生較大的擠壓作用,反映在圖21中,Pipe 3管路出口位置的壓力在液體充滿前始終保持在初始大氣壓0.1 MPa狀態。圖21中,1.4 s時刻,Pipe 3管路的充填率達到1,表明該時刻Pipe 3管路第一次充滿液體,同時反映在壓力曲線上,該時刻Pipe 3管路出口位置的壓力開始上升,由于末端節流圈流阻較小,液體到達該位置時有小幅壓力波動,并未出現明顯的水擊現象。

圖21 Pipe 3管路充填率及出口位置壓力曲線Fig.21 Filling rate and pressure history on outlet of Pipe 3

2.1 s左右,Pipe 2管路內的液體全部進入了Pipe 3管路內,且Pipe 1管路內的氣柱開始進入Pipe 3管路,圖22中Pipe 3管路入口液體流量開始下降至0,氣體流量開始逐漸增長。由于Pipe 1管路內的氣柱在Valve 1閥門打開后受到上下游液體的擠壓而具有較高的壓力和較小的體積,隨著流動的進行由于沿程流阻、局部流阻以及流動的波動,盡管壓力逐漸降低,但當其移動至Pipe 3管路內時仍具有一定的壓力,表現為圖21中管路出口壓力曲線在充填率開始降低時壓力突然升高,隨著該氣柱從管路末端的逐漸排出,Pipe 3管路內重新充滿液體,至此貯箱后的整個管路內全部充滿液體,流動達到穩態。

圖22 Pipe 3入口液體和氣體流量Fig.22 Mass flow rate of liquid and gas on inlet of Pipe 3

由以上仿真結果及分析可知,夾帶有單個氣柱的液體在充填過程中表現出的壓力振蕩由2個因素耦合造成：①氣柱受到上下游液柱的壓縮而產生的壓力波動；②氣柱下游液柱在節流元件位置產生的水擊壓力振蕩。在氣柱與液柱的耦合振蕩下,在本文參數范圍內,管路內水擊壓力峰值達到上游供應壓力的5倍左右。

本文建立的氣液兩相充填管路模型,能夠捕捉多個氣液交界面的位置移動,如本節單個氣柱與上下游液柱的氣液交界面、充填過程的氣液交界面等共計3個氣液交界面；同時利用本文模型模塊化的優勢,之后還可根據需要對夾帶有多個氣柱的液體充填過程進行仿真計算。

4 結 論

1)綜合數值耗散和CPU 耗時兩方面的因素,在同時滿足計算準確性和實時性的要求下,推薦使用TVD_QUICK 格式進行對流項的離散。

2)與本文模型相比,集中參數模型忽略初始時刻管路中原有氣體的影響,導致計算得到的充填水擊壓力明顯偏高。

3)在本文參數范圍內,節流孔徑比/在0～0.110之間,液體對預存氣體管路的充填過程處于水擊效應忽略模式,/在0.110～0.174之間對應水擊效應微弱模式,/在0.174～0.641之間對應水擊效應主導模式；最大水擊壓力峰值一般發生在水擊效應主導模式下。

4)夾帶單個氣柱的常溫推進劑在充填過程中表現出的壓力振蕩由2個因素耦合造成：①氣柱受到上下游液柱的壓縮而產生的壓力波動；②氣柱下游液柱在節流元件位置產生的水擊壓力振蕩。在氣柱與液柱的耦合振蕩下,管路內水擊壓力峰值達到上游供應壓力的5倍。

5)待充填長度與預先充液長度的比值/直接反映了液體慣性對充填過程的影響,是后續工作中值得進一步研究的內容。