?

浮式氫能儲運過程中FLH2通道管外降膜流動的海上適應性強化機理

2024-02-23 10:23孫崇正李玉星許潔韓輝宋光春盧曉
化工進展 2024年1期
關鍵詞:降膜流型海況

孫崇正,李玉星,許潔,韓輝,宋光春,盧曉

(1 山東科技大學儲能技術學院,山東 青島 266590;2 中國石油大學(華東)儲運與建筑工程學院,山東省油氣儲運安全重點實驗室,山東 青島 266580;3 國家管網集團北京管道有限公司,北京 100101)

氫能具有清潔、高效的優點。在“碳達峰、碳中和”背景下,全球主要能源消費國相繼把氫能提升到國家能源戰略層面[1]。隨著海上制氫和碳捕獲、利用與封存(carbon capture, utilization and storage,CCUS)技術的發展,利用海上風電資源或天然氣制備氫氣,并通過儲運技術送到氫能源市場,為解決海上風電并網和消納的難題、促進深海天然氣資源的低碳發展提供了可行的思路[2]。國際海洋工程界提出一種參考浮式液化天然氣生產裝置(FLNG)的浮式氫氣液化裝置(FLH2),浮式氫氣液化裝置通過將氫氣直接在海上液化并裝船運輸,簡化了海上氫能的儲運過程,極大地提高了設備使用的靈活性和安全性[3-6]。但受海上風浪的影響,應用于浮式氫氣液化裝置的制冷劑在降膜流動換熱過程中,會出現冷劑分配不均、薄液膜失穩等問題,進而影響液化工藝的性能指標,因此為強化浮式氫氣液化裝置的海上適應性,亟需開展相關領域的研究工作。

繞管式換熱器具有傳熱系數高、冷劑充灌量少、結構緊湊等優點,被廣泛應用于天然氣和氫氣液化領域。繞管式換熱器主要由纏繞管路和殼體兩部分組成,換熱通過殼側制冷劑在繞管管壁和管間的降膜沸騰來實現,為繞管內的工質提供冷量。Jian 等[7]搭建了整體繞管式換熱器的實驗裝置,通過實驗研究尺寸參數對換熱器換熱性能的影響。Zhang 等[8]通過低溫實驗和數值模擬相結合的方法研究了超低溫氦在繞管式換熱器的換熱特性。Chien 等[9]通過實驗研究了光管、翅片管和沸騰強化管的管外降膜沸騰換熱特性。Hu 等[10]通過實驗方法研究了乙烷/丙烷混合制冷劑在繞管式換熱器殼側流動沸騰的傳熱機理。Geir等[11]對應用于氫氣液化工藝中板翅式換熱器和繞管式換熱器的換熱性能進行理論研究,結果表明繞管式換熱器中正-仲氫轉化的 損失較小。Geni?等[12-13]通過實驗研究優化了繞管式換熱器殼側換熱計算方程。Lu 等[14-15]通過數值模擬方法研究了管間距、管外徑、管層數以及中心筒直徑對換熱器性能的影響。Wang等和Jian等[16-18]通過多目標遺傳算法分析了幾何參數對繞管式換熱器性能的影響。Yu 等[19-21]構建了數值模擬模型,研究了螺旋管內烷烴類工質的冷凝流動與換熱特性。Li 等[22-24]分別研究了常規圓管、波紋管、方波管和螺旋強化管內烷烴類介質的冷凝傳熱特性,研究結果表明,相較于常規圓管,波紋管、方波管和螺旋強化管的管型結構均可提高管內冷凝換熱特性。Ding 等[25-27]研究了烷烴類介質在換熱器殼側純氣相流動和降膜流動過程的換熱與壓降性能。

與陸上的應用條件不同,繞管式換熱器在浮式海上平臺應用時,受多自由度傾斜晃蕩海況的影響,換熱器內部會出現液體偏流的現象,導致管外制冷劑分布不均,管壁局部液膜較薄甚至出現干涸現象,引起換熱惡化。為提高繞管式換熱器的海上適應性,Zheng 等[28]對其分配器性能進行研究,提出一種新型氣液分配器結構。Wang 等[29]建立了FLNG 繞管式換熱器的三維數學模型,研究海上晃蕩條件下換熱器傳熱傳質過程。Duan 等[30]建立了分相多股流動態模型,預測海況下FLNG繞管式換熱器不同相態區域之間的相變換熱規律。Li等[31]研究了晃蕩條件下繞管式換熱器盤管內烷烴類介質冷凝過程中的流動與換熱特性。Ren 等[32-33]通過數值模擬研究發現,相較于純氣相流動換熱,晃蕩對降膜流動換熱過程影響較大。

應用于浮式氫氣液化裝置的新型繞管式換熱器通道管內由于填充了正-仲氫轉化催化劑,導致通道管內流體區域狹小,而海況對有自由表面的管外流體降膜流動過程影響較大,進而影響其傳熱傳質過程,降低FLH2的海上適應性。因此,為強化FLH2裝置的海上適應性,本文基于高速攝像系統和晃蕩平臺系統以及3D 建模打印技術,同時基于耦合的VOF 和Level Set 兩相流模型和動網格技術,搭建了降膜流動浮式可視化實驗裝置和數值模擬模型。采用浮式微觀可視化實驗和數值模擬相結合的方法,研究FLH2中降膜流動過程的海上適應性強化機理,提出降膜流動強化措施。

1 降膜流動可視化實驗裝置

1.1 實驗裝置建立

氫氣液化工藝中制冷循環由預冷循環和深冷循環兩部分組成,深冷循環以氦等工質的膨脹制冷為主,預冷循環中冷劑種類較多,包括氮氣,甲烷、乙烷、丙烷、正丁烷等烷烴類工質,將深冷循環中的氦氣和原料氫氣從常溫25℃冷卻至-193℃左右,預冷循環的穩定運行對整個氫氣液化系統非常重要。預冷循環的非共沸制冷劑在降膜換熱過程中,存在薄液膜界面,海況對含有自由液膜的流動過程影響較大,通過影響降膜流動狀態與液膜分布形態,進而影響預冷循環中繞管式換熱器的換熱效率,最終降低氫氣液化系統的液化率與比功耗,因此本文主要研究海況對烷烴類工質降膜流動與液膜的影響情況。但是常規低溫液相制冷劑溫度較低且組分配比復雜,采用液相混合冷劑作為實驗工質時需要對換熱器進行較為復雜的保冷和吹掃試壓等操作,不利于觀測和研究不同管型結構對降膜流動特性的影響,因此本文選用常溫下為液體的工質進行替代。常溫下水的密度為1007kg/m3、黏度為0.8904×10-3Pa·s、表面張力為72.10mN/m,乙醇密度為787.6kg/m3、黏度為0.9626×10-3Pa·s、表面張力為22.10mN/m,戊烷密度為620.7kg/m3、黏度為0.2204×10-3Pa·s、表面張力為15.39mN/m。與常規低溫液相制冷劑相比,戊烷的物性相差較小。同時戊烷的量綱為1 的伽利略數Ga1/4(559.3)與常規低溫液相制冷劑的Ga1/4(558.9)數值相近,而與乙醇(178.3)和水(497.5)的伽利略數Ga1/4則相差較大。因此優選戊烷作為降膜流動可視化實驗工質,解決了液氫用預冷循環制冷劑中組分較輕不利于可視化實驗,而常規的降膜流動可視化實驗工質(水與乙醇)的物性與制冷劑相差較大的問題。

降膜流動浮式微觀可視化實驗流程與設備實物如圖1 所示,實驗裝置包括:測試流體循環系統、高速顯微攝像與圖像采集處理系統、海上傾斜晃蕩平臺系統。其中測試流體循環系統主要由動力裝置、儲罐、調節閥組和降膜流動測試裝置組成。測試流體由動力系統增壓,經流量計、調節閥組后,流入降膜流動測試裝置,在測試裝置內部降膜流動后進入儲罐,測試流體從儲罐底部流回動力裝置的入口,實現流體循環。

圖1 降膜流動浮式微觀可視化實驗流程(a)與設備實物圖(b)

降膜流動測試裝置安裝在海上傾斜晃蕩平臺上,在平臺操作系統的作用下,隨著平臺的運動模擬海上的傾斜晃蕩工況。測試裝置外殼為透明亞克力材料,由于其內部流體流動速度較快,普通攝像機不能對其流動過程進行精確捕捉,因此采用高速攝像顯微技術,將高速攝像機安裝在海上傾斜晃蕩平臺上,對海況條件下測試裝置內部降膜流動過程進行動態捕捉。不同的管型結構尺寸參數對管外降膜流動過程影響較大,為了更加準確地研究不同管型降膜流動流型演化機理,基于3D建模打印技術,制作了降膜流動可視化實驗測試管,包括光圓管、螺紋波管、方波管、波紋管(圖2)。

圖2 降膜流動可視化實驗測試3D打印管

1.2 可視化實驗研究結果

通過可視化實驗發現流體的雷諾數(Re)影響管外降膜流動流型,隨著雷諾數的增加,降膜流動流型依次為滴狀流流型、柱狀流流型、扇狀流流型及其過渡流型。相鄰流型之間的Re與Ga有關,表達式如式(1)所示。滴狀流流型下冷劑流量不足,管壁表面會出現干涸區域,換熱效率較低,扇狀流流型下制冷劑流量較大,冷劑大量浪費,導致液化系統比功耗較低,而柱狀流流型下制冷劑的經濟性較好?;诟咚亠@微攝像與圖像處理系統,通過微觀可視化實驗得到管外徑8mm、管間距4mm 光圓管、螺紋波管、波紋管和方波管的降膜流動流型演化規律,如圖3所示,Re與工質流量直接相關,本文所研究的Re范圍為0~900??梢钥闯龀R巿A管和波紋管比螺紋波管和方波管的柱狀流范圍較大,流體更易形成柱狀流流型。流型劃分圖的提出為工程中合理選擇制冷劑流量范圍提供理論依據。

圖3 不同結構的管外降膜流動流型劃分

式中,a是降膜流型劃分因子;Ca是毛細數,Γ是換熱管單側單位長度降膜流動液體質量流量;μ是液體的動力黏度;σ是液體的表面張力;ρ是液體的密度;g是重力加速度。

水平靜止條件下光圓管外降膜流動狀態如圖4所示,可以看出,在水平靜止條件下降膜流動流型較為穩定,液膜均勻覆蓋在管壁表面,管間液柱在液體入口中心,液柱形態穩定,液膜均勻鋪展。在海況3°傾斜條件下光圓管外降膜流動狀態如圖5所示,實驗結果表明,海況3°條件對降膜流動流型影響較小,管間液柱和液膜鋪展都較為穩定,因此在海況3°條件下,常規光圓管的降膜流動流型有較好的海上適應性。

圖4 水平靜止條件下光圓管外降膜流動狀態

圖5 海況3°條件下光圓管外降膜流動狀態

基于高速攝像系統和晃蕩平臺系統,研究了海況6°條件下光圓管外降膜流動的動態過程。如圖6(a)所示,傾斜一側流型以扇狀流流型為主,流體流量過剩,而另一側沒有液膜覆蓋;隨著時間的推移,如圖6(b)~(e),另一側流體流量不足,降膜流動流型以滴柱狀流為主,降膜流動和液膜鋪展的穩定性均較差;最后如圖6(f)所示,管間降膜流動流型又恢復到圖6(a)的形態。海況9°條件下光圓管外降膜流動狀態如圖7 所示,可以看出,管間降膜流動狀態與海況6°條件下流動狀態相似,管間液體偏移更為明顯,液滴快速傾斜降落,在傾斜側聚并形成液扇。這是因為海況主要是對重力起主導作用的降膜流動過程影響較大,傾斜海況下,重力的軸向分力會產生軸向加速度,從而影響液膜沿軸向的鋪展過程,導致管間降膜流型多以局部扇狀流和滴狀流為主,因此在海況6°和9°條件下常規光圓管的海上適應性較差。

圖6 海況6°條件下光圓管外降膜流動狀態

圖7 海況9°條件下光圓管外降膜流動狀態

基于海上傾斜晃蕩平臺、高速攝像和3D 打印技術,分別研究了螺紋波管、方波管和波紋管的降膜流動海上適應性。水平靜止條件下方波管外降膜流動狀態如圖8 所示,可以看出,液柱多從方波管的大管徑處降落,液柱的流動過程較為穩定,液膜鋪展性也較好。但當處于傾斜9°海況時,如圖9(a)~(f)所示,由于大管徑和小管徑處突變的影響,液體無法在軸向均勻鋪展,出現了大量的干涸區域,因此方波管外降膜流動的海上適應性較差。水平靜止0°條件下螺紋波管外降膜流動狀態如圖10 所示。從圖10(a)、(b)、(f)可以看出,流體多從螺紋波管的螺紋縫隙降落,會出現少量的干涸區域。當螺紋波管處于傾斜6°海況時,如圖11 所示,與常規光圓管相似,在傾斜一側以扇狀流流型為主,而另一側以滴狀流流型為主,受螺紋縫隙阻隔的影響,滴狀流流型下管壁表面出現了大量的干涸區域,海上適應性較差[圖11(c)~(f)]。因此通過可視化實驗發現,在海況條件下,由于管壁表面突變導致了較強的阻隔效應,螺紋波管和方波管外降膜流動和液膜鋪展的均勻性、穩定性均較差。

圖8 水平靜止條件下方波管外降膜流動狀態

圖9 海況9°條件下方波管外降膜流動狀態

圖10 水平靜止條件下螺紋波管外降膜流動狀態

圖11 海況6°條件下螺紋波管外降膜流動狀態

圖12 給出了水平靜止0°條件下波紋管外降膜流動過程,可以看出,液柱多從波紋管的突出部分降落,液柱較為穩定,液膜鋪展性也較好。當處于傾斜6°海況時,如圖13 所示,管壁表面仍未出現干涸區域,管間流型多以柱狀流流型為主,雖然液柱傾斜,但管間流型依然穩定,液膜的均勻穩定性較好。當海上傾斜角度增加到9°時,如圖14所示,與6°海況類似,波紋管管間液柱雖有少量傾斜,但整體流型較為穩定,且液膜鋪展的均勻性和穩定性均較好。因此在海況條件下,從管間流型穩定性、液膜分布均勻穩定性、管壁干涸情況、液膜縱向偏移和柱狀流流型區間等多方面考慮,相較于常規光圓管、螺紋波管和方波管,推薦采用波紋管作為FLH2中降膜流動換熱段的管型結構,以強化浮式氫氣液化裝置的海上適應性。

圖12 水平靜止條件下波紋管外降膜流動狀態

圖13 海況6°條件下波紋管外降膜流動狀態

圖14 海況9°條件下波紋管外降膜流動狀態

2 波紋管降膜流動數值模擬模型建立

2.1 控制方程與計算方法

為了進一步定量研究波紋管的管外降膜流動特性,得到薄液膜厚度分布等微觀尺度參數,本文建立了波紋管降膜流動數值模型。模型的連續性方程如式(2)。

動量方程如式(3)~式(6)所示。

F為動量源項,如式(7)所示,包括重力源項(FG)、表面張力源項(Fσ)、氣液拖曳力源項(FLG),如式(8)、式(9)。

本文采用耦合的VOF[volume of fluid,如式(10)]和Level Set兩相流模型進行降膜流動數值模擬。數值模擬模型采用壓力速度耦合方法選擇PISO,動量方程空間離散采用精度較高的QUICK,壓力離散選用PRESTO!,體積分數選用Geo-Reconstruct,Level Set采用QUICK。

2.2 網格劃分及模型驗證

波紋管管外降膜流動數值模型的網格劃分如圖15 所示,圓周角定義為α。在管壁周圍區域進行網格加密,貼近管壁處網格尺寸為0.015mm,以捕捉管壁附近的薄液膜等微觀尺寸參數。模型為質量較高的六面體結構化網格,時間步長設置為1×10-5s。入口邊界條件設置為質量流量,左右兩側邊界為對稱邊界條件(symmetry)。為了驗證波紋管管外降膜流動數值模擬結果的可靠性,模擬得到的戊烷工質(密度620.7kg/m3、黏度0.2204×10-3Pa·s、表面張力15.39mN/m)降膜流動流型與可視化實驗結果進行了對比,如圖16所示,通過對比發現,構建的降膜流動數值模擬模型與微觀可視化實驗結果相似性較好,可以準確模擬管外降膜流動流型。

圖15 網格劃分示意圖

圖16 降膜流動流型模擬結果與實驗結果對比

3 數值模擬結果

基于動網格技術,編寫了UDF(user defined function)程序導入數值模擬模型中,研究不同角度的海況條件對降膜流動過程的影響規律。在傾斜6°海況條件下,常規光圓管外降膜流動隨時間的變化情況如圖17 所示。在48.581ms 時,第一根管壁表面開始覆蓋液膜,液膜逐漸鋪展,到65.581ms時,第一根管壁表面完全覆蓋,液柱開始逐漸形成,隨后在100.301ms時,液柱逐漸向下拉伸,與第二根管壁表面接觸,但在傾斜海況和重力的作用下,傾斜一側液柱先降落,重力的軸向加速度導致傾斜側管壁表面液膜較厚,另一側管壁表面液膜鋪展性較差。在174.791ms 時,出現了大量干涸區域,干涸區域面積隨著時間的增加雖有所減小,但仍然存在(t=210.110~252.681ms)。

圖17 海況6°條件下光圓管外降膜流動隨時間變化過程

水平靜止0°條件下波紋管外降膜流動隨時間變化過程如圖18 所示。在25.086ms 時,液體開始降落,在35.086ms 時,降膜流動液體逐漸在軸向匯聚形成波峰,隨后59.876ms 時,第一根管壁表面已經完全鋪展液膜,并在突出部分匯聚液體。在103.016ms 時,第二根管壁表面液膜覆蓋約一半,到123.016ms,液膜幾乎完全覆蓋波紋管,第二根管只有底部少量管壁表面出現干涸區域。隨著時間的推移,到149.341ms時,液膜完全鋪展覆蓋在波紋管壁表面。

圖18 水平靜止條件下波紋管外降膜流動隨時間變化過程

圖19 給出了海況6°條件下波紋管外降膜流動隨時間變化過程,在31.303ms時,液膜開始覆蓋管壁,逐漸沿著軸向和周向鋪展,到51.303ms時,液膜完全覆蓋第一根波紋管壁表面。在71.303ms時,液體以液柱的形式開始降落,隨后到91.303ms,液柱逐漸與第二根管壁表面相接觸,在重力軸向加速度的作用下,傾斜一側液體較多,覆蓋的區域也較寬。在117.880ms時,傾斜一側管壁幾乎完全覆蓋,而另一側管壁還出現少量干涸區域,但在183.223ms時,波紋管壁表面被液膜完全覆蓋。

圖19 海況6°條件下波紋管外降膜流動隨時間變化過程

海況9°條件下波紋管外降膜流動隨時間變化過程如圖20 所示。在20.701ms 時,液體降落到第一根管壁表面,在65.712ms時,液體逐漸向軸向鋪展,完全覆蓋第一根波紋管并在突出部位匯聚。隨后108.652ms時,與海況6°時相似,管間液體在海況的影響下逐漸偏移。在138.652ms時,液膜幾乎覆蓋了第二根波紋管的表面,但在傾斜另一側仍舊有少量的干涸區域,到152.852ms,干涸區域逐漸縮小。隨著時間的推移,在251.412ms時,液膜完全覆蓋第二根波紋管管壁表面。通過圖19 和圖20研究發現,在海況條件下,波紋管管間液柱雖有少量偏移,但降膜流動流型總體穩定,液膜也較為均勻地覆蓋在管壁表面,因此波紋管有較好的海上適應性。

圖20 海況9°條件下波紋管外降膜流動隨時間變化過程

為了與水平靜止0°時的降膜液膜厚度進行定量對比,采用齊次坐標技術來描述空間的各點坐標及其變換,將通過動網格技術傾斜的降膜流動數值模型恢復到水平靜止0°狀態,轉換方法如式(11)~式(14)所示。

通過式(11)~式(14)的計算分析,得到了水平靜止0°、海況6°和9°傾斜條件下,波紋管外降膜流動液膜厚度沿縱向管長分布。如圖21 所示,在圓周角為45°時,水平靜止、海況6°和9°的主波峰厚度分別為0.482mm、0.519mm、0.546mm,海況6°和9°的主波峰縱向位置分別偏移了1.237mm 和1.425mm。圓周角增加到80°時,水平靜止0°、海況6°和9°的主波峰厚度分別為0.532mm、0.557mm、0.535mm,主波峰的縱向位置分別偏移了1.328mm和1.895mm。進一步增加圓周角至110°時,傾斜海況6°的波峰位置偏移了1.565mm,傾斜海況9°的波峰位置偏移了1.853mm。圓周角為145°時,水平靜止0°、海況6°和9°的主波峰厚度分別為0.688mm、0.701mm、0.766mm,傾斜6°和9°海況的主波峰位置分別偏移了1.140mm 和1.378mm。傾斜海況導致了波峰沿縱向管長的偏移,傾斜角度越大,偏移量越大;但隨著圓周角增加到145°時,波峰縱向位置偏移量逐漸減小,可以看出,隨著流體沿波紋管降膜鋪展,波紋管的管型結構對傾斜流體的修正作用不斷增強,偏移程度逐漸降低,進一步驗證了波紋管較好的海上適應性。

基于數值模擬結果,通過修正Hou 等[34-35]提出的液膜厚度計算關聯式(15),得到波紋管外降膜流動平均液膜厚度的計算關聯式。

式中,α為圓周角;ρG和ρL為氣體和液體的密度;μL為液體的黏度;c和n通過研究結果進行擬合得到;S為管間距;管徑D=( )Dmin+Dmax2。本文根據雷諾數Re和管徑D范圍劃分計算關聯式的適用范圍:低(Re/D)<70mm-1,中(Re/D)為70mm-1≤(Re/D) <170mm-1, 高(Re/D) ≥170mm-1。 在 低(Re/D)的條件下,c和n值與Hou模型的數值相同。在中(Re/D)條件下,圓周角<90°時,c為1.15,n為-0.07;圓周角>90°時,c和n值與Hou模型的數值一致。在高(Re/D)條件下,圓周角<90°時,c和n值分別為1.96和-0.07;圓周角>90°時,c和n值分別為1.60和0.07。優化后波紋管平均液膜厚度關聯式計算值與模擬值的誤差分析如圖22所示,可以看出誤差都在±25%以內,大部分誤差在±20%以內。圖23給出了海況下平均液膜厚度計算關聯式的精度,計算關聯式也可以較好地預測海況下波紋管管外降膜流動的平均液膜厚度。平均液膜厚度計算關聯式的提出為工程中合理選擇波紋管型結構和流體雷諾數提供理論依據。

圖22 優化后理論計算關聯式的波紋管液膜厚度與模擬值對比

圖23 海況下理論計算關聯式的波紋管液膜厚度與模擬值對比

4 結論

基于高速攝像系統、晃蕩平臺系統和3D 建模打印技術,搭建了研究常規光圓管、螺紋波管、方波管和波紋管的降膜流動浮式可視化實驗裝置,進行了海況條件下降膜流動的微觀可視化實驗研究;基于耦合的VOF和Level Set兩相流模型以及動網格技術,進行FLH2中降膜流動的數值模擬研究,揭示降膜流動海上適應性的強化機理,主要結論如下。

(1)構建的降膜流動數值模擬模型與微觀可視化實驗結果吻合性較好,可以準確模擬管外降膜流動流型;繪制了常規光圓管、螺紋波管、波紋管和方波管管外降膜流動的流型劃分圖。

(2)通過降膜流動浮式微觀可視化實驗和數值模擬研究結果可知,在海況條件下,從管間流型穩定性、液膜分布均勻穩定性、管壁干涸情況、液膜縱向偏移和柱狀流流型區間等多方面考慮,相較于常規光圓管、螺紋波管和方波管,推薦采用波紋管作為FLH2中降膜流動換熱段的管型結構,以強化浮式氫氣液化裝置的海上適應性。

(3)通過實驗和數值模擬相結合的方法,得到了波紋管降膜流動平均液膜厚度計算關聯式,根據雷諾數和管徑比值劃分計算模型的適用范圍?;贖ou理論模型,擬合了不同雷諾數Re和管徑D下的波紋管外降膜流動的平均液膜厚度計算模型,誤差都在±25%以內,海況條件下液膜厚度計算關聯式也有較高的計算精度。

猜你喜歡
降膜流型海況
水平井油水兩相流型實驗研究
燒堿裝置降膜蒸發器的腐蝕檢查
典型海況下艦載發射箱結構強度仿真分析
惡劣海況下海洋石油116內轉塔式FPSO裝配載優化
板式降膜蒸發器激光焊蜂窩狀波紋傳熱板
極端海況下軟質海崖大規模蝕退計算模型研究
共流型轉子的有限元分析
極限海況下單點系泊系統纜索動張力研究
基于Taitel-Dukler方法的氣液兩相流型邊界計算軟件開發
水平管外R404A降膜蒸發傳熱的實驗研究
91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合