航天器液氧輸送管路閥控瞬變的數(shù)值模擬

航天器液氧輸送管路閥控瞬變的數(shù)值模擬

以氫氧為燃料的推進器供應(yīng)系統(tǒng)在打開、關(guān)閥等閥控操作過程中不可避免地會涉及到變換流的問題。變變流的控制直接關(guān)系到系統(tǒng)的安全和可靠性。由于液氧密度遠高于液氫，在閥控過程中產(chǎn)生的水壓過大，對閥控的要求更高。同時，液氧屬于低溫液體，其管道的組成明顯不同于室溫液體。因此，萃取能耗系統(tǒng)的瞬變流是系統(tǒng)設(shè)計和控制程序制定中應(yīng)考慮的主要問題之一。對于在設(shè)計階段的液氧傳輸系統(tǒng)，由于試驗方法無法獲得液體傳輸系統(tǒng)的瞬變性能，因此通過計算數(shù)值計算預(yù)測瞬變過程是非常重要的。閥控是低溫泵送管路控制的主要方面。由于送送管的在線控制是流量控制的熱點問題，文獻介紹了計算不帶波紋管的液氫傳輸系統(tǒng)水壓的問題。在文獻中，我們討論了聯(lián)合波動方程和閥方程對管道水壓降致的優(yōu)化控制，以及管道參數(shù)的設(shè)計和動態(tài)流量控制對復(fù)雜管道系統(tǒng)壓力的優(yōu)化控制。目前，關(guān)于低溫液體輸送管路數(shù)值的研究很少。由于高速液體輸送系統(tǒng)中有大量的波形管道、波幅補償和管道由許多短管組成，因此很難正確計算軍閥控制引起的瞬變。本文通過變化液氧輸送管路中最常用兩種特性閥門的動作指數(shù)和動作時間條件,綜合考慮了非穩(wěn)態(tài)摩阻和波紋軟管對水力瞬變的影響,數(shù)值預(yù)測了一條液氧輸送管路在幾種閥控條件下的瞬變過程,并通過對比分析提出優(yōu)化方案,為液氧輸送系統(tǒng)的瞬變控制提供參考.1系統(tǒng)的物理模型和計算方法1.1多層預(yù)熱系統(tǒng)密度計算液氧輸送系統(tǒng)包括貯槽、真空多層絕熱硬管、真空多層軟管、控制閥、、流量計等,管路通徑從75~250mm不等.真空多層絕熱硬管由短鋼管串聯(lián)組合而成,鋼管一端固定,一端自由伸縮,鋼管之間采用波紋補償連接.真空多層絕熱軟管長度一般分為1.5和3m兩種,每隔一定距離或在彎頭處或設(shè)置一根波紋軟管.文中以通徑為150mm的液氧輸送管路為例進行計算,包括液氧貯槽、真空多層絕熱硬管、真空多層絕熱軟管和控制閥4種決定系統(tǒng)瞬變的主要部件,管路總長420m,其中波紋軟管總長31.5m,如圖1所示.G1、G2分別為上、下游液氧貯罐,V為末端控制閥,p1、p2、f分別為末端閥前壓力、閥后壓力和管路末端流量.1.2物理模型和計算方法1.2.1維管流數(shù)值計算真空多層絕熱硬管中液氧瞬變流問題屬于有壓管道的非恒定流問題.閥控動作時,壓力和流速等流動參數(shù)變化很快,而傳熱過程很慢,快慢過程同時進行計算在仿真上屬于剛性問題,因此不考慮傳熱的影響.液氧在管路內(nèi)的流速一般不超過5m/s,而壓力波在硬管中的波速為幾百m/s,忽略一維管流瞬變流動控制方程中的對流項,以壓力p和質(zhì)量流率6)m表示的硬管內(nèi)瞬變流動的控制方程為其中:ρ為液體密度;g為重力加速度;x為沿管長的流向坐標(biāo);A為管道橫截面積,A=A(x,t);α為管道軸線即x的正方向與水平面的夾角,高度沿x方向增加時α為正;p為管截面平均壓力,p=p(x,t);6)m為管截面質(zhì)量流率,6)m=ρAv,v為管截面平均流速,v=v(x,t);a為壓力波波速;f為摩擦阻力因數(shù).一維管道瞬變流動數(shù)值計算結(jié)果的精確性與摩擦阻力因數(shù)f密切相關(guān).摩擦阻力因數(shù)f決定于摩擦阻力模型.當(dāng)前管道的摩阻模型有半穩(wěn)摩阻模型和非穩(wěn)摩阻模型兩種.半穩(wěn)摩阻模型認為瞬時黏性切應(yīng)力與穩(wěn)態(tài)黏性切應(yīng)力相等.非穩(wěn)摩阻模型發(fā)展起來以前,通常以半穩(wěn)摩阻模型計算瞬變過程,這種替代對黏度較大的流體或快速瞬變過程,計算與實驗差距較大.對于非穩(wěn)摩阻模型,當(dāng)前主要包括兩類.一類是基于瞬時加速度的Brunone模型,Brunone模型考慮時變加速度和位變加速度對摩阻的影響,將摩阻與流體特性聯(lián)系起來,對非穩(wěn)湍流有很好的適用性,在高Re數(shù)流態(tài)模擬上得到廣泛應(yīng)用;另一類是以Zielke模型為代表的加權(quán)函數(shù)模型.Zielke模型認為非穩(wěn)層流摩阻是流動加速度和歷史加速度變化的一個加權(quán)函數(shù).Zielke模型在Re<10000的范圍內(nèi)具有很好的模擬效果,但因沒有考慮對流項的影響,對強對流條件下的瞬變模擬容易產(chǎn)生誤差.雖然液氧的黏度較低,但由于液氧管路較長,存在關(guān)閥等快瞬過程,且液氧流動中的Re>10000,因此文中選擇基于瞬時加速度的摩阻模型進行計算.以質(zhì)量流率表示的瞬時加速度模型(MIAB)的非穩(wěn)摩阻因數(shù)的表達式為式中:fq為Darchy-Weisbach摩擦因數(shù),液氧為牛頓流體,對Re≤2320的層流,fq=64/Re;對Re>2320的紊流,按Colebrook-White公式計算.當(dāng)前,數(shù)值求解有壓管道的非恒定流問題最廣泛的方法是特征線法.特征線法是利用特征線將偏微分方程轉(zhuǎn)化為常微分方程進行求解.以下求解一維管流非恒定流動的方法是根據(jù)文獻、文獻推導(dǎo)而得.將式(2)代入方程(1),并令A(yù)'=A/(1+0.5k),a'=a/(1+0.5k)1/2,k'=0.5k/(1+0.5k)1/2,整理得式(4)、(5)的求解與符號Φ有關(guān),取網(wǎng)格比θ=Δt/Δx=1/a,Φ=±1對應(yīng)的特征線差分網(wǎng)格如圖2(a)、(b)所示.非線性摩擦項采用具有二階精度的線性隱式近似,在特征線上積分即可獲得對應(yīng)的差分方程,其中S、R點通過插值法計算,即可求得P點的流動參數(shù).1.2.2波紋管彈性模數(shù)計算假定:1)在瞬變流動中,波紋軟管長度可在彈性范圍內(nèi)自由伸縮,而波紋補償則忽略其軸向伸長,只考慮其徑向變形;2)波紋管及波紋補償?shù)淖冃尉趶椥苑秶畠?nèi);3)忽略節(jié)點耦合,考慮泊松耦合的影響,根據(jù)質(zhì)量守恒和動量守恒,建立波紋管內(nèi)瞬變流動的控制方程:式中:Ai為以波紋管內(nèi)直徑Di計算的橫截面積;δx為波距;Ω為波紋管單波容積;ζ為單波阻力系數(shù);K'為波紋管有效體積彈性模數(shù),表示流體和波紋體積變化的彈性影響.對于波紋管的有效體積彈性模數(shù),鑒于波紋管波高遠小于內(nèi)徑,在徑向應(yīng)變中波高的應(yīng)變大大小于半徑的應(yīng)變,因此認為波紋管產(chǎn)生軸向變形時,波峰和波谷同時且以相同的增量徑向擴張,而波紋之間的壓差導(dǎo)致軸向伸長,則可通過波峰處的周向應(yīng)變計算半徑變化,通過軸向剛度計算單波長度變化,得到波紋管有效體積彈性模數(shù)的計算式如式(7):式中:Kl為液氧的體積彈性模量;Ks為波紋管彈性模量,Ks=-Ω/W.其中式中:Ro為波紋管外半徑;Ri為波紋管內(nèi)半徑;q為波距;Kr為周向應(yīng)力系數(shù);Acu為單個波紋的金屬橫截面面積,Acu=nδ(0.571q+2h);E為管材在介質(zhì)溫度下的楊氏彈性模量;Dm為波紋管平均直徑,Dm=Db+h+nδ;Db為波紋管直邊段內(nèi)徑,Db=Do;h為波高;δ為波紋管壁厚;n為厚度為δ的波紋管層數(shù);fiu為波紋管單波軸向彈性剛度,且fiu=(1.7DmEδ3n)/(h3Cf),Cf為波紋管單波軸向彈性剛度的修正系數(shù)通過流場數(shù)值模擬建立波紋補償及波紋管單波阻力系數(shù)隨流速的變化關(guān)系.流場計算時,利用單波周期性邊界條件和RNGk-ε雙方程模型,動量方程、湍動能方程、耗散率方程均采用二階迎風(fēng)差分格式;選擇simple算法,采用強化壁面函數(shù)處理和yplus自適應(yīng)網(wǎng)格加密技術(shù),控制殘差ε≤10-6,并要求數(shù)值計算結(jié)果與網(wǎng)格無關(guān).管路中波紋補償長度很短,而且直線安裝,兩端固定,發(fā)生瞬變流時,波紋補償沿軸線方向幾乎沒有伸縮變形,只有徑向膨脹.計算中,將波紋補償作為局部阻力元件處理.而波紋軟管的數(shù)學(xué)模型通過等效轉(zhuǎn)化后,可以利用特征線法求解.式(8)與以內(nèi)徑為Di的等徑管道的瞬變控制方程(1)形式相同,其求解方法與利用特征線法求解一維管道瞬變流動控制方程相同,不同的是此時fe為一確定數(shù)值,而不是表達式.取網(wǎng)格比θ=Δt/Δx=1/a,數(shù)值求解該方程組恒穩(wěn)定.1.2.3貯罐出口壓力平衡方程假定液流方向從上游流向下游為正,貯罐出、入口P點處壓力為式中:“-”號適用貯罐位于管道上游,“+”號適用貯罐位于管路下游,H(t)、p0(t)分別為t時刻貯罐液面高度和氣枕壓力,ζ為貯罐出口處的局部阻力系數(shù),A為貯罐出入口管道橫截面積,ρ為液體密度,g為重力加速度.式(9)求解過程以上游貯罐為例,由貯罐出口P點壓力平衡方程(9)和P點下游管段的負特征線相容方程得到任意第n時間步的一元二次型方程式(10)、(11)中,CM為負特征壓力線,由與P點下游相鄰節(jié)點S的管道特性及點S在(n-1)Δt時刻的流動參數(shù)按式計算,R=(fqaΔt)/(2ρDA2)表示節(jié)點所在管段的阻力系數(shù);B=a/A表示節(jié)點所在管段的特征阻抗.ρgH(t)、p0(t)分別由nΔt時刻貯罐液位高度和液面氣枕體力計算.式(11)可以用線性理論法求解.方法是使阻力損失項線性化:形成迭代方程:C=[CM-(ρgH(nΔt)+p0(nΔt))],k為迭代次數(shù).1.2.4局部阻力及等配比特性閥式中:pi為閥上游入口壓力;po為閥下游出口壓力;(pi-po)是流體通過閥門時壓力的瞬時降落;m6)為閥門開度為τ時從入口到出口的質(zhì)量流率;AD為以閥門公稱通徑計算的截面積;ζ(τ)為門開度為τ時的阻力系數(shù).閥門局部阻力系數(shù)ζ決定于開度τ下流道的幾何形狀、流道壁面的相對粗糙度和雷諾數(shù).因受局部障礙的強烈擾動,閥內(nèi)流體在較小雷諾數(shù)時就已進入湍流粗糙區(qū),雷諾數(shù)的變化對湍動程度的實際影響很小,可以認為閥門局部阻力系數(shù)ζ只取決于開度,與Re無關(guān).根據(jù)閥門流量系數(shù)的定義,可以得到任意開度τ下閥門的局部阻力系數(shù):其中:Δpe=100kPa;ρe=1000kg/m3;Kv(1)表示開度為100%時的閥門流量系數(shù),單位為t/h;f(R,τ)為閥門的特性.等百分比特性閥:f(R,τ)=R(τ-1).其中:R為閥門可調(diào)比.R=qvmax/qvmin.閥桿行程l=0時,qv=qvmin;閥桿行程最大,即l=L時,qv=qvmax.τ=l/L.閥門在開閉過程中,開度τ是時間的函數(shù),關(guān)閥時,τ=(1-t/tc)m,開閥時,τ=(t/to)m.tc和to分別是關(guān)閥和開閥時間,m是閥動作過程中開度與時間的關(guān)系指數(shù)即閥動作指數(shù).求解控制閥瞬變流的方法是將式(12)與閥門進出口連接管道的特征線相容方程聯(lián)合,整理得到一元二次方程:式中:Cp為正特征壓力線,由與P點上游相鄰節(jié)點A的管道特性及點A在(n-1)Δt時刻的流動參數(shù)計算.式(13)同樣可通過線性理論法求解.當(dāng)液氧管路組成、布局、組件尺寸和流動特性、邊界壓力及液氧物性等參數(shù)確定后,即可對管路壓力和質(zhì)量流率進行求解.2瞬變預(yù)測模型的驗證按照上述建模及求解方法,對具有試驗數(shù)據(jù)的某液氧輸送管路進行關(guān)閥瞬變計算.穩(wěn)態(tài)下液氧流速為2.3m/s,在0時刻關(guān)閉管路始端控制閥,閥門特性為快開特性,閥動作指數(shù)m=2.0,關(guān)閥時間tc為4s.管路中部監(jiān)測瞬變壓力.計算結(jié)果與試驗數(shù)據(jù)對比如圖3所示.由圖3可知,計算數(shù)據(jù)與測試數(shù)據(jù)基本吻合,實測與計算所得波峰波谷的相位基本一致,說明按照前文所述方法進行瞬變預(yù)測具有一定的精確度和可信度.3初始液氧分離計算條件:上游貯罐液面到參考面高差6m,液面壓力0.6MPa,下游貯罐液面到參考面高差34m,液面壓力0.11MPa.管內(nèi)液氧溫度80K,對應(yīng)飽和蒸汽壓30.120kPa,初始恒定流狀態(tài)下,管內(nèi)單相液氧流動.瞬變過程中,當(dāng)管內(nèi)壓力低于飽和蒸汽壓即發(fā)生液柱分離.管路末端控制閥公稱通徑DN100,可調(diào)比R為50,額定流量系數(shù)為160.變化閥門特性、閥動作指數(shù)和動作時間,討論對閥控動瞬變的影響.計算過程中取時間步長Δt=0.001s,網(wǎng)格比θ=Δt/Δx=1/a,控制相對誤差ε≤5×10-4.3.1截止閥對壓力的影響對控制閥進行關(guān)閥和開閥瞬變計算,tc/to=4s,m=1.0;閥門分別為等百分比調(diào)節(jié)閥(簡稱調(diào)節(jié)閥)和快開截止閥(簡稱截止閥).閥門從0時刻開始動作.計算得到的主要瞬變過程如圖4和圖5所示.圖4、圖5表明,調(diào)節(jié)閥進行閥控動作產(chǎn)生的波動明顯比截止閥小,調(diào)節(jié)閥能夠有效控制水擊壓力.調(diào)節(jié)閥關(guān)閉過程中,閥前水擊壓力最高為0.8MPa.而截止閥關(guān)閉過程中,閥前產(chǎn)生了1.2MPa的水擊壓力.調(diào)節(jié)閥關(guān)閥過程中,通過閥門的流量均勻減小.截止閥關(guān)閥過程中,前3s流量只減少約10%,在關(guān)閥的最后1s,流量迅速從1140L/min降低到0.關(guān)閥過程中流量均勻減少是調(diào)節(jié)閥能夠有效抑制關(guān)閥水擊的原因.通過積分計算,關(guān)閥的4s時間內(nèi)通過調(diào)節(jié)閥的總流量是38.6L,通過截止閥的流量是78.1L.很明顯采用調(diào)節(jié)閥比采用截止閥更容易控制液氧輸送總量精度.開閥過程中,通過調(diào)節(jié)閥的流量呈均勻上升狀態(tài),過渡過程平穩(wěn)緩和.由于閥門本身的阻力相對整條管路的阻力很小,開啟閥門不能立即提升管路流量,在開閥過程提升管路流量的速度上,截止閥并不比調(diào)節(jié)閥具有明顯優(yōu)勢.3.2開閥和關(guān)閉過程中閥前壓力波動圖6、7是管路末端調(diào)節(jié)閥動作時,閥指數(shù)m分別為0.5、1.0、1.5、2.0,tc/to=4s時計算到得瞬變過程.m=0.5,表示關(guān)閥動作先慢后快,開閥動作先快后慢;m=1.0表示開、關(guān)過程中閥桿勻速動作;m=1.5和2.0表示關(guān)閥動作先快后慢,開閥時則相反.計算表明,m≥1.0時,閥動作指數(shù)越大,閥前壓力波動幅度越大.當(dāng)m=0.5時,其瞬變曲線表現(xiàn)出與其余三者截然不同的行為,閥門關(guān)閉產(chǎn)生的閥前壓力比m=2.0時小,但比m=1.0時大,與m=1.5時相當(dāng).m=1.0時,閥前壓力波動最平緩.閥門關(guān)閉期間通過閥門的流量:m=0.5時為52L;m=1.0時為33L;m=1.5時為31L;m=2.0時為26L.m=0.5的流量要遠遠大于其余三者.關(guān)閥指數(shù)越大,通過閥門的流量越少,理論上越有利于制氧輸送總量精度的控制.但指數(shù)增大帶來的總流量精度控制并不顯著,當(dāng)總量精度控制在100L以下時,優(yōu)化調(diào)節(jié)閥關(guān)閥指數(shù)對總流量精度控制的貢獻可以忽略.綜合考慮壓力波動和輸送總量精度控制因素,在0.5、1.0、1.5、2.0幾種調(diào)節(jié)閥指數(shù)中,采用m=1.0的方案對瞬變過程控制最有利.3.3末端截止閥關(guān)閉情況圖8~10為m=1.0條件下,閥門動作時間tc變化時,管路末端控制閥開、關(guān)閥動作對瞬變過程的影響.預(yù)測表明,閥動作時間越長,引起的壓力波動越小.對調(diào)節(jié)閥,變化閥動作時間對水擊壓力波幅值影響很大,延長關(guān)閥時間,可以有效降低閥前水擊壓力,關(guān)閥時間從3s延長到6s,最高水擊壓力從0.92MPa降低到0.71MPa,降低水0.21MPa.延長關(guān)閥時間在大幅降低水擊壓力的同時,也延長了第一個水擊波到達的時間.末端截止閥關(guān)閉時,關(guān)閥時間從4s延長到7s,最高水擊壓力從1.19MPa降低到1.14MPa,只降低了0.05MPa.對截止閥,延長關(guān)閥時間更多地表現(xiàn)為延長第一個水擊波到達的時間,而不是水擊壓力的降低.對末端調(diào)節(jié)閥,隨著關(guān)閥時間的延長,關(guān)閥期間通過閥門的總流量緩慢增大.在不同關(guān)閥時間下從閥門開始動作到完全關(guān)閉,通過閥門的流量依次為:tc=3s時為33L;4s時為39L;5s時為47L;6s時為56L.如果輸送總量控制精度為100L,則均在總流量控制精度范圍內(nèi).對末端截止閥,關(guān)閥期間通過閥門的流量隨著關(guān)閥時間的延長急劇增加,從閥門開始動作到完全關(guān)閉,通過閥門的流量依次為:tc=4s時為78L;5s時為97L;6s時為117L;7s時為136L.當(dāng)關(guān)閥時間為6s時,關(guān)閥期間通過閥門的總量就超過了100L的控制精度要求.因此,末端截止閥延長關(guān)閉時間對輸送