基于氮氣擠代的推進劑供應系統動態特性研究

王 謙,秦 侃,郝常樂,張安靜,羅 凱,黨建軍

(西北工業大學 航海學院,陜西 西安,710072)

0 引言

隨著作戰需求的提升與科技水平的進步,世界各國都把增強魚雷綜合性能作為水下武器裝備的重要發展方向之一[1-2]。能源供應系統作為魚雷動力系統的重要組成部分,對魚雷航速、航程和航深等性能都有較大影響。典型魚雷的能源供應系統原理[3]如圖1 所示。

圖1 某型魚雷能源供應系統原理圖Fig.1 Principle of torpedo energy supply system

擠代式推進劑供應系統由高壓氣艙、電爆活門、管路以及推進劑艙組成。當魚雷入水后,海水電池一方面激活點火器點燃固體藥柱,產生大量燃氣后使互鎖閥打開;另一方面激活電爆活門,使高壓氣艙與推進劑艙連通。此時,高壓氣體快速流出高壓氣艙并向推進劑艙充氣,擠壓其中的燃料填充管路空隙并在燃料泵前建立一定壓力,以確保燃料泵供出足量的燃料進入燃燒室進行穩定燃燒[4]。若在固體藥柱燃盡之前,燃料泵前未建立起指定壓力,會使動力系統偏離設計點工況運行,甚至導致啟動失敗等情況發生[5]。由此可見,擠代式推進劑供應系統[6-7]的動態特性與魚雷動力系統的啟動加速性能密切相關。圖1 中的擠代式供應系統與泵供式供應系統采用級聯的方式,確保燃料能夠快速穩定地向燃燒室供應,但在擠代過程中,管內壓力的劇烈變化會對氣體產生強擾動,在流場中產生激波從而影響系統的動態穩定性。因此,有必要建立一種方法用于研究擠代式推進劑供應系統的動態特性。

國內針對水下能源供應系統的相關研究較少。羅凱等[8]研制了一種水下熱動力能源供應系統,采用比例器對三組元燃料的動態供應進行精確控制,三組元供應比例最大誤差不超1.5%;李代金等[9]建立了由燃燒室、霧化噴嘴以及旋轉燃燒室組成的水下熱動力能源供應系統的數學模型,討論了柱塞燃料泵以及噴嘴與燃燒室的動態匹配關系,并提出了設計和選擇時應當避免的問題。目前尚無針對擠代式推進劑供應系統動態性能的相關公開研究報道。

在管內激波流動研究方面,Roy[10]采用黎曼解算器對管內流體流動進行了理論研究與數值計算,并針對經典Sod 激波管問題,采用Python 代碼對其密度、壓力、馬赫數和熵值進行了近似求解。Chen 等[11]采用高階點隱式求解方法研究了激波管內激波與邊界層的相互作用,得到了二維激波管中激波與邊界層的非定常相互作用過程。Qiu等[12]提出了可壓縮格子玻爾茲曼方法(lattice Boltzmann method,LBM)仿真激波/邊界層在激波管中的相互作用,結果表明加權無波動無自由參數耗散差分格式(weighted non-oscillatory and non-freeparameter dissipation difference,WNND)-Z 方案和WNND-L 方案都能很好地捕捉渦和激波的結構。目前國內外針對激波管內流動的數值計算較為豐富,能夠為擠代式推進劑供應系統內部的激波捕捉以及流動仿真提供理論依據。

考慮到國內外對水下擠代式推進劑供應系統動態特性的實驗研究有限,且由于系統模型內產生的激波問題,以及管內產生的流動分離和渦流等多動態結構,加大了流體計算軟件計算此類問題的復雜性[13-15]。針對此問題,開發了一種一維瞬態可壓縮求解器用于解決擠代式推進劑供應系統的動態特性仿真問題。建立一維程序分析推進劑供應系統內部的非定常流動狀態,采用計算流體力學(computational fluid dynamics,CFD)仿真和經典文獻進行對比,驗證了程序的可靠性;采用一維程序研究了不同活門直徑、高壓氣艙壓力和推進劑艙容積對系統模型的平衡時間和穩定壓力的影響,分析了三者對系統動態特性產生影響的規律,并給出了相關因素對系統動態特性影響的設計參考。通過對系統動態特性的研究顯示一維程序可以更加便捷高效地得到所需系統動態特性下的模型結構參數,具有較好的工程應用價值。

1 推進劑供應系統數值模型

推進劑供應系統由高壓氣艙、管路、活門和推進劑艙等部分組成,其基本結構如圖2 所示。

圖2 推進劑供應系統結構簡圖Fig.2 Structure of propellant supply system

1.1 高壓氣艙和管路數值模型

1.1.1 計算模型

由于魚雷空間的限制,高壓氣艙不易過大且需要有一定量的氣體儲備,因此設計為細長罐結構。同時,考慮到管內流動以軸向流動為主,高壓氣艙內的氣體流動與管路中的流動類似,均可簡化為一維流動,其一維連續性方程為

式中:ρ為密度;t為時間;ux為軸線方向的速度分量。

一維動量方程為

式中:τij為垂直于i軸的微元表面上受到的沿j方向的切應力;fx為軸線方向上的單位質量力。

一維能量方程為

式中:E為流體微團的總能,即內能、動能和勢能之和;T為溫度;λ為熱傳導系數;h為焓。

采用有限體積法將高壓氣艙和管路離散為儲存流體信息的網格單元,每個網格單元都視為1 個控制體,網格單元之間的交界面視為控制面。圖3為典型的儲存流體信息的控制體(j),以及相鄰控制體的交界面

圖3 儲存流體信息的控制體Fig.3 Control body for storing fluid information

每一個網格單元內的氣體數據都被儲存于數據結構之中,包括邊界條件、時間步長、網格單元中心位置坐標、網格單元交界面位置坐標、網格單元內流體信息(密度、速度、溫度、壓力、內能和音速)變化的平均值以及狀態通量(質量通量、動量通量和能量通量)的時間導數。為保證程序的2 階精度,采用網格中心點數據來代表網格單元內氣體的各種數據進行計算。

網格單元的交界面上儲存著流體經過此壁面的質量通量、動量通量和能量通量,因此流體在高壓氣艙和管路中的流動就可轉化為在多個控制體內流體信息的變化以及控制面上的通量交換。僅考慮理想氣體的一維流動,忽略粘性項且假設管路區域中任何狀態變化都是連續的,則根據方程(1)～(3)計算出的網格單元內的質量、動量和能量,可用來更新網格單元內的流體信息,即

式中:umass,j、umomentum,j和uenergy,j分別為j單元網格上的質量、動量和能量;vj為單元網格上流體速度;ej為網格單元上流體比內能。擠代系統中的氮氣可被認為是熱完全氣體,則其比內能與溫度成正比,表示為

式中,cv為氣體的定容比熱容。繼而根據已求解出的密度與溫度通過理想氣體狀態方程可求出壓力

式中,R為氣體常數。

聲速為

式中:aj為網格單元上的音速;γ為比熱比。通過假設每個網格單元內流體信息的線性變化,將網格單元內流體信息的平均值插值到控制體的交界面上,從而得到每個控制體左右兩端的流體信息。采用AUSMDV(advection upstream splitting methoddifference vector)方法[16]計算交界面的狀態通量。

1.1.2 邊界條件

高壓氣艙與管路的邊界條件統一采取在外部擴建2 層虛擬網格的方法,以保證流體的流動在空間分布和時間推進上都具有2 階精度。設置的2 層網格雖在實際的氣體流動物理網格之外,但卻包含應用于邊界條件的各種數據信息。圖4為增加了虛擬網格后,每個儲存氣體各種數據信息的網格單元結構。

圖4 根據實際網格流體信息復制的虛擬網格圖Fig.4 Virtual grid copied from actual grid fluid information

若邊界條件為壁面時,網格中的參數采用以單元壁面為對稱軸的方式復制對稱面網格單元中的流體信息,但要注意的是,此時速度的方向應與對稱面單元的速度相反;若邊界條件為入口時,單元中的質量流量和溫度從上一級流入的網格單元中獲取;若邊界條件為出口時,單元中的溫度與壓力將傳遞到下一級的網格單元中。

1.2 活門計算模型

在高壓氣體向推進劑艙充氣擠代推進劑的過程中,活門某一截面處氣體流速會達到音速從而發生壅塞現象,因此,通過活門的質量流量和出口溫度按照壅塞和非壅塞2 種狀態來處理。

將高壓氣艙出口處的壓力和溫度賦值為活門入口處的壓力和溫度,分別記作PL,valve和TL,valve,將管路入口處的壓力賦值為活門出口處的壓力PR,valve,則活門出口處的溫度為

式中:A為活門橫截面積;Cd為流量系數,其受活門直徑、壓力和溫度的影響,根據經驗Cd的取值范圍為0.71～0.93。

1.3 推進劑艙計算模型

考慮到推進劑艙為直徑與高近似相等的圓柱狀結構,可采用集總參數法忽略推進劑艙的實際體積尺寸,將其視為一點,用該點上流體信息隨時間變化的一元函數關系來表征推進劑艙內部整體的變化情況,如圖5 所示。

圖5 推進劑艙控制體Fig.5 Propellant camber control body

選取推進劑艙為控制體(control volume,CV),管路出口處的接觸面為控制面(control surface,CS),顯然推進劑艙內恒有質量流入而無質量流出,根據式(1)和式(3)可得

式中:V為體積;e為單位質量流體所具有的能量;Acs為管路出口處的橫截面積。

又因

由此,可以獲得時間推進當中t+?t時刻下推進劑艙內部氣體流動的質量與能量。根據式(20)～(22)即可得到推進劑艙內部能量、溫度和壓力隨時間的分布變化情況

1.4 時間推進

對于每一個網格單元,其內部各狀態參量從時間n推進到時間n+1,采用歐拉預測校正格式迭代計算出每個網格內的質量、動量和能量,再由單元體內已知信息推導出內部各狀態參量,即

式中:x為離散控制體單元的空間步長;包含控制體j單元內的質量、動量和能量;包含控制體前后j±1/2單元壁面上的質量通量、動量通量和能量通量;是由上一次迭代計算出的控制體內部的狀態參數得到的控制體內儲存的質量、動量和能量。

為了保持穩定,時間推進步長限制為

式中: Δtallowed是對所有網格單元設定的最小值;α是收斂條件判斷數,在仿真中通常限定α<0.5。由于流動過程為可壓縮流動,在Δtsignal的計算上考慮音速a的影響,對于每個網格單元取值為

2 模型驗證

2.1 管路流動模型驗證

針對經典一維Sod 激波管問題,分別采用一維程序與L1d 程序[17]進行求解,進行交互對比來驗證一維程序在流體流動方面的可靠性。

將激波管管徑設為D=0.01 m,長度為1 m,其內填充氮氣,分為高壓和低壓2 個部分。壁面邊界條件應用于管路兩端,交換邊界條件應用于管路中間,忽略黏性效應,管路內部初始狀態的設定如表1 所示。

表1 計算域模型初始狀態Table 1 Initial state of computing domain model

取流動時間t=0.6 ms 時2 種方式下的溫度、壓力、密度和速度在軸向分布的對比情況,如圖6所示。結果顯示,一維程序在仿真流體流動方面與Sod 激波管的程序仿真相近,誤差較小,通過對比分別采用一維程序、5 階目標本質無振蕩(targeted essentially non-oscillatory,TENO)格式、5 階加權本質無振蕩(weighted essentially non-oscillatory,WENO)格式、5 階加權緊致非線性方案(weighted compact nonlinear schemes,WCNS)格式、5 階Steager-Warming 格式和加入人工黏性項,計算經典Sod 激波管問題后得到的數值計算情況,證明了文中數值方法在仿真管內流動方面的精確性,因此所建立的一維程序可以用于管內流動狀態的仿真計算。

圖6 一維程序與L1D 代碼的驗證對比Fig.6 Validation comparison of one-dimensional program and L1D code

2.2 系統模型整體驗證

文中采用Fluent 軟件仿真計算的結果與一維程序計算結果進行交叉驗證。

建立如圖7 所示的推進劑供應系統模型,模型共分為高壓氣艙段、活門段、管路段和推進劑艙段4 個部分。給出如表2 所示的典型推進劑供應系統模型各部分尺寸及初始狀態設定。

表2 系統模型尺寸及初始狀態設定Table 2 System model size and initial state setting

圖7 典型推進劑供應系統模型幾何尺寸Fig.7 Geometry of typical propellant supply system model

為驗證一維程序對模型在某一時刻數值仿真的準確性,如圖8 是一維程序求解得到的0.5、1.0和1.5 ms 流動時刻下,模型軸線處的流體參數與CFD 仿真所得結果的對比情況。

圖8 不同時間下系統壓力軸線變化圖Fig.8 System pressure axis changes under different time

從圖中可以看出,各個時刻內高壓氣艙段壓力的變化曲線與CFD 仿真結果一致。而推進劑艙中CFD 仿真的管路內出現了壓力的震蕩,其主要是由于管徑變化[18-20]使活門后產生了如圖9 所示的膨脹波導致管路內部壓力的震蕩,而一維程序雖然無法很好地捕捉到其震蕩特性,但由圖8(d)～(f)可以看出,其對于震蕩減弱后的數值捕捉良好,并且其震蕩數值均值與一維程序仿真結果的相對誤差小于10%。

圖9 速度云圖Fig.9 Velocity cloud picture

總體來說,一維程序基本可以反映系統模型在軸線方向上流體信息的變化狀態。

為進一步驗證一維程序對模型在某一截面處流體參數隨時間變化的準確性,圖10 是采用數值方法求解得到的高壓氣艙前10 mm 位置和推進劑艙中部處的流體參數與一維程序仿真所得數據結果的對比情況。

圖10 高壓氣艙與推進劑艙流體參數對比Fig.10 Comparison of fluid parameters between high-pressure chamber and propellant chamber

從圖10(a)～(c)中可以看出,高壓氣艙段一維程序仿真得到的結果與CFD 仿真結果基本一致,能夠反映高壓氣艙的泄氣過程。而從圖10(e)和圖10(f)中可以看出,在推進劑艙段的仿真曲線中兩者仿真的壓力和質量流量效果也接近一致,而圖10(d)中一維仿真程序仿真的溫度與CFD 仿真結果有些許偏差,其主要原因是CFD 仿真了流體在管路中流動與壁面因黏性摩擦而產生的熱量[21],激波較強時溫度的真實值較之以熱完全氣體計算出的相應數值相差較大。因此,CFD 仿真相比于一維程序有更高的溫度,兩者總體趨勢接近,相對誤差不超過7%。總體來說,一維程序可以仿真系統模型各個截面流體信息隨時間的變化情況。

綜上所述,一維程序在仿真系統模型流動方面有較高的準確性,可以節省CFD 數值計算因網格數過大所需的過長時間,因此可以采用一維程序來求解系統模型內流體參數的變化。

3 推進劑供應系統動態特性研究

采用一維程序選取表2 中不同活門直徑(D=4、6、8 mm)下的系統模型進行數值計算,計算過程中通過改變高壓氣艙壓力以及推進劑艙容積,獲得如圖11 所示的系統模型內流動達到穩定時,系統的平衡時間和穩定壓力云圖。

圖11 不同活門直徑下的系統動態特性分析Fig.11 Analysis of system dynamic characteristics under different valve diameters

通過不同活門直徑條件下,高壓氣艙壓力對系統平衡時間的影響可以看出: 由圖11(a)所示,在D=4 mm 時,平衡時間隨著高壓氣艙壓力的升高而緩慢增大,高壓氣艙的壓力每升高2 MPa,最終平衡的時間增加0.006～0.01 s,而當活門直徑增大到D=8 mm 時如圖11(c)所示,平衡時間基本不會再隨著高壓氣艙壓力的升高而變化。這可能是由于在活門直徑較小時,系統內部不平衡的壓差導致通過活門產生的入射激波與反射激波相互干擾而導致系統達到平衡的時間延長,而增大活門直徑后,活門處產生的入射激波大幅減小[22],對平衡時間的干擾也逐漸減弱。

類似的,不同活門直徑條件下,推進劑艙容積變化對系統平衡時間也有相應的影響。從圖11(a)～(c)中可以看出,在活門直徑為D=4 mm 時,推進劑艙的容積每增加2 L,最終平衡時間增加0.06 s,隨后隨著活門直徑的增大,推進劑艙容積的改變量對最終平衡時間的影響逐漸減弱,在活門直徑變為原來的2 倍后,推進劑艙容積的改變量對最終平衡時間的影響減弱為原來的1/5。由式(12)～(14)可以看出,通過活門的質量流量也會隨著活門直徑的增大而增加,從而導致平衡時間可以在一定的容積范圍內更快地達到平衡。

由圖11(d)～(f)中不同的活門直徑條件下,高壓氣艙壓力及推進劑艙容積變化對系統穩定壓力的影響可以看出: 系統穩定壓力隨著高壓氣艙壓力的升高而增加,但增加的幅度又隨著推進劑艙容積的增大而減小,這也符合理想氣體狀態方程中的描述。隨著推進劑艙容積的增大,最終平衡的壓力從一開始的高壓氣艙壓力每增加2 MPa,平衡壓力增加0.12～0.16 MPa,減小到高壓氣艙壓力每增加2 MPa,平衡壓力增加0.1 MPa;隨著高壓氣艙壓力的升高,最終平衡的壓力從一開始的推進劑艙容積每增加2 L,平衡壓力減小0.04 MPa,增大到推進劑艙容積每增加2 L,平衡壓力減小0.08MPa。

不同活門直徑、高壓氣艙壓力以及推進劑艙容積條件下系統的平衡時間與穩定壓力各有不同,說明不同因素的變化情況對不同條件下系統的動態特性均有影響。下文通過一維程序計算得到的不同位置處的流體信息,進一步分析三者對系統動態特性的影響。

3.1 高壓氣艙壓力對系統動態特性的影響

改變表2 中高壓氣艙壓力,分別取P=10、15、20 MPa 后進行一維程序仿真,得到高壓氣艙末端與推進劑艙處的流體壓力特性曲線如圖12 所示。

圖12 不同高壓氣艙壓力下系統動態特性對比Fig.12 Comparison of system dynamic characteristics under different pressures of high-pressure chamber

由圖12(a)可以看出,高壓氣艙壓力分別為P=10、15、20 MPa 時,系統平衡時間分別為0.85、0.86、0.86 s,可見高壓氣艙內壓力的變化對高壓氣艙與推進劑艙達到平衡時間的影響較弱。根據小孔出流規律,更高的壓差會增大活門出口處的質量流量,加速系統內整體壓力的平衡,因此高壓氣艙壓力的變化對系統平衡時間不敏感。由圖12(b)可以看出,隨著高壓氣艙內壓力的上升,由于系統為密閉容器沒有與外界聯通,因此推進劑艙內最終穩定的壓力也隨之上升,在高壓氣艙壓力分別為10、15、20 MPa 的條件下,系統最終穩定壓力分別為0.64、0.95、1.28 MPa。

3.2 推進劑艙容積對系統動態特性的影響

改變表2 中推進劑艙容積,分別取V=22、24、26 L 進行一維程序仿真,得到高壓氣艙末端與推進劑艙處的流體壓力特性曲線如圖13 所示。

圖13 不同推進劑艙容積下系統動態特性對比Fig.13 Comparison of system dynamic characteristics under different propellant volumes

由圖13(a)可以看出,在推進劑艙容積分別為V=22、24、26 L 的條件下,系統最終平衡時間分別為0.79、0.82、0.86 s,隨著推進劑艙容積的增大,系統最終平衡的壓力也會隨之降低,而通過活門的質量流量依舊不變,高壓氣艙需要更長的時間達到更低的壓力,因此高壓氣艙與推進劑艙達到穩定的時間緩慢增大。由圖13(b) 可以看出,隨著推進劑艙容積的增大,根據理想氣體狀態方程可知,最終整個系統達到穩定時刻的穩定壓力會隨著推進劑艙容積的增大而減小。在推進劑艙容積分別為V=22、24、26 L 的條件下,系統最終穩定的壓力分別為1.43、1.28、1.20 MPa。

3.3 活門直徑對系統動態特性的影響

改變表2 中活門直徑參數,分別取D=4、6、8 mm 后進行一維程序仿真得到高壓氣艙末端與推進劑艙處的流體壓力特性曲線如圖14 所示。

圖14 不同活門直徑下系統動態特性對比Fig.14 Comparison of system dynamic characteristics under different valve diameters

由圖14(a)可以看出,在D=4、6、8 mm 的條件下計算得到的高壓氣艙最終平衡時間分別為2.01、0.82、0.45 s,可見隨著活門直徑的擴大,相同條件下活門出口流出的質量流量增大,系統達到穩定的時間前移。圖14(b)中推進劑艙內穩定的壓力會有不同,這是由于一維程序在流體仿真上無法反映流體在通過突擴管和突縮管內時流體的真實情況,如程序忽略了將活門看做整體進行集總參數的計算時會產生不同的壓力損失和熱損失導致最終的計算偏離理想氣體的條件,以及流體在管內流動的局部損失等,從而會導致最終穩定的壓力有些許偏差。

4 結束語

文中根據典型魚雷能供系統建立了描述推進劑系統特性的一維數值仿真程序,通過與經典文獻以及CFD 仿真的對比,證實了一維程序的可行性。利用一維程序對系統模型不同管徑、高壓氣艙壓力和推進劑艙容積的條件下展開仿真,得出以下結論:

1) 推進劑供應系統擠代過程中,活門管徑的適當增大可以明顯縮短系統達到平衡的時間。在活門管徑由4 mm 變為8 mm 后,活門處質量流量增加,最終使得系統的平衡時間提前至原來的22.4%。

2) 高壓氣艙的壓力對系統達到平衡時間的影響隨著活門管徑的增大而減弱。活門直徑D=4 mm 時,高壓氣艙壓力每升高2 MPa,平衡時間增加0.06 s,穩定壓力增加0.128 MPa;在活門直徑變為原來的2 倍后,高壓氣艙壓力變化對平衡時間基本沒有影響,其主要與不同管徑下活門處產生的膨脹波和激波有關。

3) 推進劑艙容積的增大會使封閉系統內最終平衡的壓力下降,導致高壓氣艙需要更長的泄氣時間以達到更低的壓力。在此條件下,將活門直徑由4 mm 變為原來的2 倍,則會在出口處產生更大的質量流量,使得推進劑艙變化對系統的影響減弱為原來的1/5。

文中對于一維程序仿真擠代式供應系統精度的討論還較為粗略,后續將結合仿真計算與試驗分析,對擠代式系統開展深入的研究,以期獲得更加適用于工程實踐的擠代式供應系統仿真程序。