高壓氫氣減壓器內部流場仿真分析

(1.重慶大學 機械傳動國家重點實驗室，重慶 400044； 2.重慶大學 機械與運載工程學院，重慶 400044；3.重慶凱瑞動力科技有限公司，重慶 401120)

引言

氫燃料電池汽車作為新能源汽車之一，由于其使用氫氣作為燃料并通過化學反應轉化成電能供給汽車使用，反應產物主要為水，不產生任何污染性氣體，因此具有零排放的優點。目前，氫燃料電池汽車多采用氣態儲氫方式，為了提高汽車的續航里程，一般采用以下兩種方法增加載氣量：一種是增加氣瓶的容積，另一種是增加氣瓶的壓力。在氣瓶的安裝空間受限制的條件下，增加氣瓶的壓力是提高氫氣儲量的主要方式。

為保證一次加氫后行使距離達200 km，通常儲氫壓力需35 MPa以上，而要保證一次加氫后行駛距離500 km，儲氫壓力則要求高達70 MPa[1]。在新能源汽車未來的發展中，傳統的單級減壓閥無法在高壓、大壓力比的工況下實現壓力穩定調節，會出現振動大、能耗高等一系列問題。因此近年來，許多學者對適用于高壓、大壓力比等環境或能實現多級降壓的調壓器做了研究。儲景瑞等[2]設計了一種高壓氣動減壓閥，并研究了輸出壓力控制精度及響應速度的影響因素。徐志鵬等[3]研制了一種高壓氣動比例減壓閥，并研究了結冰特性對減壓閥性能的影響。張春[4]研究了先導式高壓氣動減壓閥主要參數對其能量和靜動態特性的影響。劉佳等[5]研究了高壓差迷宮式調節閥，采用數值模擬方法著重分析了不同流道形式對壓降控速性能的影響。唐騰飛等[6]提出了一種適用于迷宮式調節閥節流碟片流道優化的設計流程，可以有效改善流道的流動特性。彭健等[7]對多級降壓調節閥的閥芯結構進行了改進，改善了閥門的空化現象。

隨著氫燃料電池汽車研究的開展，不少學者研究并設計了一些新型式的氫氣減壓閥，訚耀保等[1,8-9]研究了氫能源汽車兩級超高壓氣動減壓閥，分析其工作機理、特性，并對其內部流場分別在35 MPa，5 MPa壓力下進行了二維流體仿真，流場中氫氣為音速或超音速流動狀態，速度較大，主要出現在節流位置下游，減壓過程主要在節流位置處實現。CHEN Fuqiang等[10-11]研究了一種由1個多級減壓器和1個多級消聲器兩部分組成的氫燃料電池電動汽車(FCEV)2步高壓減壓系統，對35 MPa迷宮式多級減壓器進行流體仿真，并分析了多孔孔板結構參數對其內部流動特性的影響，研究發現壓力降低和速度增大主要反映在閥門開口的節流部件上，孔板層數增加或者孔徑減小都會降低湍流強度。QIAN Jinyuan等[12]對氫燃料電池氫氣減壓過程中通過多級特斯拉閥的反向氫流動進行了數值研究。李純杰[13]研究了減壓閥的減壓機理，設計了一種高壓氫氣減壓閥，對其進行了二維數值模擬計算。JIN Zhijiang等[14]為實現加氫站氫穩定減壓，設計了一種用于加氫站多級減壓閥，并通過分析過熱蒸汽和氫氣2種流體流動特性，驗證其在加氫站中應用的可行性。

值得注意的是，盡管減壓閥的工作氣體為氫氣，但為了簡化計算，在以往的仿真研究中多采用理想氣體做為工作介質。盡管真實氣體與理想氣體仿真結果會存在誤差，但仿真結果也具有一定參考價值。陳晨[15]設計了一種用于燃料電池的氫氧電氣比例減壓閥，并針對燃料電池雙路壓差控制問題，采用理想氣體仿真，對減壓閥進行了控制特性研究。朱旺[16]分析對比了理想氣體與真實氣體狀態方程的仿真結果，研究發現理想氣體與真實氣體的壓力與速度分布基本相同，溫度會受一些影響。

上述研究中，關于適用于70 MPa高壓的減壓器研究較少，且研究的減壓器多為迷宮式，加工制造難度大，以錐閥為主減壓部件的減壓器數值模擬采用簡化的二維模型，無法觀測到整體流場的分布狀態。本研究對一種以錐閥、節流孔為節流元件的兩級高壓氫氣減壓器開展了三維數值模擬，分析不同閥芯開度對其流場分布特點及流場物理參數的影響。本研究中溫度對殼體產生的熱應力遠小于壓力作用產生的機械應力，可忽略不計，因此采用理想氫氣來研究其在流場中的流動特性，為氫燃料電池汽車供氫系統在高壓、大壓力比等復雜條件下，實現壓力穩定調節提供技術支持。

1 基本原理

1.1 高壓氫氣減壓器兩級減壓原理

高壓氫氣減壓器的減壓原理為節流降壓，即當氣體經過突縮管時，有效流動橫截面積急劇減小，氣體呈加速流動狀態，根據能量守恒定律，壓力能轉化為動能，使得動能增大，壓力能降低，因此，高壓氣體流經突縮管后可以實現節流降壓。圖1為高壓氫氣減壓器兩級減壓示意圖，高壓氫氣從高壓腔室進入，經過閥門實現第一次減壓后進入中壓腔室，再經過節流孔完成第二次減壓，最終低壓氫氣進入低壓腔室。通過兩級減壓，實現較大減壓比的高壓氣體減壓[17]。其中，閥門開口大小由閥芯的位置確定。

1.2 流體流動控制方程

流體流動也需要滿足質量守恒定律、動量守恒定律、能量守恒定律。守恒定律的基本控制方程[18]列出如下。

圖1 兩級減壓示意圖Fig.1 Schematic diagram of two-stage decompression

1) 質量守恒方程

質量守恒方程即為連續性方程：

(1)

式中，ρ—— 流體的密度

u—— 流體的速度

t—— 時間

上述連續性方程對于可壓縮或不可壓縮流動，黏性或無黏性流動，定常或非定常流動均可適用。當流體為定?？蓧嚎s流動時，式(1)可變為：

▽·(ρu)=0

(2)

需要提及的是，由于空氣與氫氣的性質、參數不同，在保證其他參數相同的情況下，空氣的質量流量會大于氫氣的質量流量，但依然滿足質量守恒關系，因此滿足本研究關于流量的趨勢性分析。

2) 動量守恒方程

按照動量守恒定律，可以推導出動量守恒方程如下：

ρu·▽u=-▽p+▽·(τ)+ρg+F

(3)

式中，p—— 靜壓

g—— 重力加速度

F—— 除重力之外的外部體積力。計算中不考慮重力和重力之外的體積力的影響，因此ρg=F=0。

τ—— 黏性剪切應力張量，可表示為：

(4)

式中，μ—— 動力黏度

u′ —— 速度u的轉置矩陣

I—— 單位張量，其余符號意義同上

3) 能量守恒方程

流體流動過程中能量守恒，表達式如下：

▽·[u(ρE+p)]=▽·[keff▽T+(τeff·u)]

(5)

式中，E—— 能量

T—— 溫度

keff—— 有效傳熱系數

τeff—— 有效黏性剪切應力張量

方程左側表示單位時間凈流入能量，右側分別表示傳熱和黏度耗散的能量轉換。

4) 湍流模型

計算流體力學軟件Fluent中包含多種湍流模型，其中k-ε雙方程模型應用較為廣泛。根據湍流黏度計算方法以及控制湍流擴散的參數k，ε的不同可分為Standardk-ε，RNGk-ε和Realizablek-ε3種模型。Standardk-ε模型對工程實際中的可壓縮流動問題具有良好的求解精度，且相對節省計算量、容易收斂，該模型是一個由大量實驗現象總結得到的半經驗公式，湍流動能方程是通過精確的方程推導出來的，而湍流耗散率方程是由物理推理、數學上模擬相似的原形方程得到。

湍流動能k的輸運方程為：

Gb-ρε-YM+Sk

(6)

湍流耗散率ε的輸運方程為：

(7)

式中，xi，xj—— 流體流動方向坐標(i,j=1,2,3)

μ—— 動力黏度

k,ε—— 湍流動能和湍流耗散率

Gk—— 由于平均速度梯度產生的湍流動能產生項

Gb—— 由于浮力而產生的湍流動能產生項

YM—— 可壓縮湍流中的波動膨脹對整體耗散率的貢獻

C1ε,C2ε,C3ε—— 常數

σk,σε——k和ε的普朗特(Prandtl)數

Sk,Sε—— 用戶自定義源項

μt—— 湍流(或渦流)黏度，可通過k和ε聯合計算出：

(8)

式中，Cμ—— 常數

上述方程中的常數默認值，C1ε=1.44，C2ε=1.92，Cμ=0.09，σk= 1.0，σε=1.3。

1.3 焦耳-湯姆遜效應

焦耳-湯姆遜效應(Joule-Thomson-effect)[19]描述的是當氣體流經截面突然減小的節流元件時產生壓降，溫度隨之產生變化的現象。該效應表明真實氣體在節流膨脹后溫度會降低或升高，導致冷效應(正效應)或熱效應(負效應)。

對于理想氣體，經絕熱節流過程后，其溫度保持不變[19]，這是由于理想氣體的溫度僅與焓值有關，而在氣體節流過程前后焓值應相等，故理想氣體的溫度在節流前后不變。

2 CFD數值模擬

2.1 物理模型假設

氣體在減壓器中的流動過程是一個復雜的變質量的過程。為了抓住問題的本質，作出如下假設：

(1) 工作介質為理想氣體；

(2) 減壓器中的減壓過程為絕熱過程[19]；

(3) 氣體流動為定常流動；

(4) 數值模擬考慮黏性耗散；

(5) 考慮壓縮性效應；

(6) 忽略氣體重力對流動的影響。

2.2 仿真模型建立

為了能更好的運用高壓氫氣減壓器，采用仿真軟件Fluent對其內部流場進行數值模擬。圖2所示為減壓器三維模型，根據三維模型的內腔結構特征，利用ANSYS Workbench中SCDM模塊通過體積抽取命令快速抽取流體域，由于生成的流體域具有對稱性，同時為了節省仿真時間，避免仿真過程中出口回流現象的影響，將流體域模型進行如下處理：① 以對稱面為中心平面，取其1/2模型來進行仿真研究； ② 去除流體域模型的圓角、倒角以及反饋腔等形狀復雜且對計算結果影響不大的區域；③ 將實際出口邊界位置延長至30 mm處，以解決出口回流現象，提高仿真結果的準確性；流體域及其簡化后的模型如圖3所示。

圖2 高壓氫氣減壓器的三維模型

圖3 1/2流體域模型Fig.3 1/2 fluid domain model

2.3 網格劃分

使用ANSYS Workbench中自帶的Mesh模塊進行流體域網格劃分，圖4所示為流體域的網格劃分模型，考慮到內部流場結構比較復雜，因此不規則區域采用Automatic Method自動生成四面體網格。對于入口端、出口端的規則區域，為了避免網格數量過多，減小仿真模型計算負荷，均采用Sweep方法生成網格。由于閥口處壓力梯度較大，該處網格質量較低會導致計算結果發散，因此對閥口進行局部網格加密處理，設置局部網格尺寸為0.05 mm，以保證網格質量。本研究中，仿真模型的網格質量以偏斜度(Skewness)作為評價指標，文中涉及的所有流體計算域模型網格偏斜度均小于0.75。

網格數量對仿真計算有著重要的影響，網格數量越多，計算結果并非越精確，且需要更多的計算資源。為確保仿真計算結果的準確性，減少計算資源的消耗，對流體域網格進行網格無關性驗證。本研究以60%閥芯開度為例，通過分析不同網格數量對質量流量變化的影響來進行驗證，當質量流量不隨網格數量增加而變化時，即可以忽略網格對仿真計算結果精度的影響。由于70 MPa高壓數值模擬需要較高的時間成本，因此選擇較低的壓力進行網格無關性分析，即入口壓力5 MPa，出口壓力1 MPa。

圖5為不同網格數量對應的質量流量圖，圖中網格數量從17萬增加到162萬。由圖可以看出，隨著網格數量的不斷增加，質量流量qm變化逐漸趨于平緩，當網格數量達到112萬后，質量流量波動僅為0.07%，基本保持不變，此時可以忽略網格數量n對仿真計算結果精度的影響，故選擇112萬的網格數量進行數值模擬仿真。

圖4 網格劃分模型Fig.4 Meshed model

圖5 網格數量-質量流量Fig.5 Grid quantity-mass flow rate

2.4 求解模型和邊界條件

1)求解模型設置

采用基于壓力(Pressure-based)求解器，進行穩態計算。由于氣體在流場中的流動速度較高，流動狀態比較復雜，故選擇湍流模型來進行求解，湍流模型選擇標準k-ε模型，該模型是工程中常用的湍流求解模型，對氣體仿真有較好效果，可靠性較高。設置流體材料屬性，減壓器內的氫氣為等熵可壓縮流體，因此將氫氣設置為理想氣體，采用理想氣體狀態方程，此時Energy選項默認打開。操作壓力(Operating Conditions)設為0。最后選用壓力-速度耦合算法(Coupled)求解，耦合算法求解對運算內存要求較高，但收斂速度更快。

2) 邊界條件設置

由于氣體在流場中為可壓縮流動狀態，入口邊界條件和出口邊界條件分別設置為壓力入口和壓力出口，對稱平面設置為對稱邊界，其他邊界設置為壁面。在低壓壓縮空氣仿真中操作壓力通常保持默認值，而本模型中入口壓力為70 MPa高壓，出口壓力為1 MPa，操作壓力為0 MPa。

在Fluent仿真過程中，當入口壓力較高，入口與出口之間的壓力差較大時，各物理參數在模擬的過程中會出現劇烈波動，導致結果發散。通過大量的模擬計算，本研究提出一種高壓入口邊界的加壓方法，該方法適用于高壓氣體仿真過程，并且能有效的改善仿真收斂效果，提高仿真結果的準確性。方法實施如下：設入口最終需達到的壓力70 MPa為目標壓力pf，每一次入口設置的壓力為pi，則加壓公式可表示為：

pi+1=(pf-pi)×0.1+pi(i=1,2,3…)

(9)

式中，i—— 迭代加壓次數

pi，pi+1—— 入口壓力迭代i和i+1次的壓力值

pf—— 目標壓力

在本研究中，入口初始壓力p1設置為7 MPa，之后按照式(9)進行逐次迭代加壓，直到入口壓力達到目標壓力70 MPa。

值得注意的是，高壓氫氣經減壓器降壓后并不是直接進入燃料電堆，還要經過減壓器出口連接的引射器，通過引射器將氫氣噴入燃料電堆，最終進入電堆的氫氣壓力約為0.2 MPa。仿真的出口平均壓力基本穩定在1 MPa左右，而在整個對稱平面的壓力分布中，最低壓力并不是1 MPa，但這一操作不影響壓力結果及整個壓力場的分布。

3 仿真結果分析

本研究對高壓氫氣減壓器進行了穩態數值模擬求解，求解過程中氣體流動狀態不隨時間變化，為了更好的觀察氣體在流場內部的流動情況，研究壓力、速度等物理參數的分布，數值仿真以閥芯開口量大小(閥芯開度)為主要變量進行研究，根據閥芯開度的不同分別進行了5組算例，各組算例對應的參數見表1。

表1 不同算例參數Tab.1 Parameters of different calculation examples

3.1 壓力場結果分析

通過對5組算例分別進行求解，得到不同閥芯開度下對稱平面的壓力分布云圖如圖6所示。由圖6可知，高壓氫氣減壓器在工作過程中，壓力下降主要發生在閥口和節流孔處，在入口腔和出口腔壓力分布都較為均勻。隨著閥芯開度的增大，閥口處的壓力梯度較小，這是因為閥芯開度越大，氫氣流過閥口處的流動阻力越小，因此，該處的壓力梯度隨著閥芯開度增大不斷減小。節流孔處則與閥口處相反，隨著閥芯開度的增加，該處壓力梯度變大。

圖6 不同閥芯開度下對稱平面的壓力分布Fig.6 Pressure distributions in symmetry plane under different valve openings

為準確了解閥芯開度變化對壓力大小的影響，對不同閥芯開度下Z方向上的壓力進行研究(Z方向為入口、出口的軸線方向)，做出不同閥芯開度的壓力圖如圖7所示。由圖7可以看出，在Z=25～35 mm處和Z=40～50 mm處出現壓力降低，該位置即為閥芯與節流孔的位置。對應20%，40%，60%，80%，100%不同閥芯開度，氫氣流經閥口后壓力從70 MPa分別降到9.73, 18.72, 27.61, 35.55, 41.71 MPa；氫氣流經節流孔后壓力分別從9.73, 18.72, 27.61, 35.55, 41.71 MPa降到1.00 MPa。隨著閥芯開度增大，閥口壓降梯度減小，節流孔壓降梯度增大，即第一級減壓比減小，第二級減壓比增大。為了確定兩級減壓比之間的關系，給出了不同閥芯開度下第一、二級的減壓比，見表2。

圖7 不同閥芯開度下Z方向上的壓力Fig.7 Pressure in Z direction under different valve openings

表2 不同閥芯開度下的減壓比R1，R2Tab.2 Decompression ratio R1 and R2 under different valve openings

根據表2中第一、二級減壓比R1，R2的數值，可以得出高壓氫氣減壓器的兩級減壓比關系滿足：

R=R1×R2

(10)

式中，R—— 總減壓比

R1—— 第一級減壓比

R2—— 第二級減壓比

由于5組算例的入口、出口邊界設置均相同，即R保持不變。則根據式(10)可知，當R為定值時，R1與R2成反比，說明第一級閥口減壓作用減弱，會導致第二級節流孔減壓作用增強。

此外，根據圖7中的放大區域可知，隨著閥芯開度的增大，閥芯上游壓力波動越大，對閥芯平衡影響越強，是閥芯受力不穩定的因素之一。

3.2 速度場結果分析

流場中速度的分布可以直接反映氫氣在流場中的流動情況。因此，本研究以高壓氫氣減壓器不同閥芯開度對速度場的影響展開研究，速度分布云圖如圖8所示。

圖8 不同閥芯開度下對稱平面的速度分布Fig.8 Velocity distributions in symmetry plane under different valve openings

圖8中在閥口與節流孔出口兩處，速度有顯著增大的區域，這是因為閥口、節流孔的通流面積變小，氫氣在流經閥口和節流孔時經歷絕熱壓縮過程，壓力降低，導致速度增加。節流孔出口的速度突然增大，這是由于氫氣在高壓差條件下通過節流孔形成了高速噴射流，且由于壓力降低，氫氣迅速膨脹，氫氣流出了節流孔后會繼續加速一段距離，由于出口空間較大，速度降低極為迅速，這與文獻[8]、文獻[13]關于速度場分析的結論基本一致。氫氣高速噴射流長時間沖擊減壓器殼體壁面，會對殼體造成破壞，這為研究殼體優化，減小噴射流提供基礎。

從圖8可以看出，不同閥芯開度下速度分布不同，且隨著閥芯開度的增大，閥口處的速度及其變化區域逐漸減小，而節流孔出口的速度及其變化區域逐漸增大，這說明速度的變化受截面突變程度的影響，截面突變程度越大，速度變化越明顯，且最大速度是出現在截面突變程度最大的位置，與文獻[10]中的結論相符。分析最大速度的變化，圖8中不同的閥芯開度，最大速度出現的位置不同，這是由于閥芯開度增大會導致閥口過流面積增加，從而減小了閥口處截面變化梯度。當閥芯開度為20%時，閥口處的截面變化梯度最大，最大速度出現在閥口處，當閥芯開度為40%，60%，80%，100%時，截面變化梯度最大位置均在節流孔處，因此最大速度均出現在節流孔出口處。

此外，研究了不同閥芯開度下Z方向速度v分布的變化，如圖9所示。圖中5種不同開度下閥口處的最大速度分別為2705, 2531, 2388, 2227, 2040 m/s。由此可見，隨著閥芯開度的增大，閥芯處的最大速度在相應減小，這是由于閥芯開度的增大導致閥口處截面突變程度降低導致的。

圖9 不同閥芯開度下Z方向上的速度Fig.9 Velocities in Z direction under different valve openings

3.3 溫度場結果分析

當流體流過高壓氫氣減壓閥中的閥口與節流孔時，屬于絕熱等熵膨脹過程，會導致溫度的急劇下降。以20%閥芯開度為例對溫度場進行分析，圖10中入口、出口溫度基本保持不變，這與理想氣體的焦耳-湯姆遜效應相合符。結合圖10與圖8可知，最低溫度與最大速度位置相對應，最低溫度約為93 K。氫氣的結霜點為14 K[13]，由此可以判斷低溫區的存在不會導致氫氣在閥芯和節流孔處結霜，因此，對高壓氫氣減壓器性能不產生影響。其他開度下的最低溫度與20%開度時差別不大，此處不再贅述。

圖10 20%閥芯開度下的溫度分布Fig.10 Temperature distribution under 20% valve opening

3.4 流量及能耗分析

研究質量流量和閥芯開度的關系，分別對5組不同閥芯開口量s的流場進行仿真，分別為 0.12, 0.24, 0.36, 0.48, 0.60 mm，記錄進出口的質量流量值如圖11所示。由圖11可知，隨著閥芯開口量的增加，質量流量不斷增大，質量流量的變化趨勢與文獻[20]的結果基本一致。

湍流耗散率ε是表征湍流所產生能量消耗的物理量，湍流耗散率ε越大，意味著湍流越大，能量消耗越大。圖12為不同閥芯開度下Z方向上的湍流耗散率，從圖12可以看出，湍流耗散率出現了2個峰值，分別對應閥口處和節流孔處2個位置，其中第1個峰值由于閥芯阻斷了流體而被分為了2個部分。湍流耗散率的峰值位置與渦旋產生的位置基本一致，由此可以從能量上驗證閥口和節流孔處渦旋區域的存在。圖12中，隨著閥芯開度的增加，閥口處的湍流耗散率降低，與之相反的是隨著閥芯開度的增加，節流孔處的湍流耗散率增大。因此，較大的閥芯打開度意味著流量損失較小，與文獻[14]中湍流耗散率的結論相吻合。

圖11 不同閥芯開口量下的質量流量Fig.11 Mass flow rates under different valve openings

圖12 不同閥芯開度下Z方向湍流耗散率Fig.12 Turbulence dissipation rates in Z direction under different valve openings

4 結論

研究了一種可用于氫燃料電池汽車供氣系統的兩級高壓氫氣減壓器，通過數值仿真研究其內部流場物理參數的分布，并以閥芯開度為主要研究變量得出了流場內部壓力、速度、能耗等參數的變化規律：

(1) 在減壓器內部流場中，壓力降低和速度增加主要反映在閥口與節流孔2個位置，隨著閥芯開度增大，閥口的減壓作用減弱，節流孔的減壓作用增強；閥芯開度越大，閥口處最大速度越小，閥口處的速度及其變化區域越小，節流孔出口處的速度及其變化區域越大；

(2) 氫氣流經閥口與節流孔后，入口溫度與出口溫度基本保持不變，符合理想氣體的焦耳-湯姆遜效應，且流場內部最低溫度與最大速度位置相對應，最低溫度高于氫氣的結霜點；

(3) 閥芯開度對質量流量與能耗有較大影響，質量流量隨閥芯開度增大而增加；隨著閥芯開度增大，閥口處的湍流耗散率減小，節流孔處的湍流耗散率增加。