LNG冷箱中降溫過程的動態(tài)模擬并行計算

王 琳,李玉星,朱建魯,王婭婷,盛歡歡,王武昌

(中國石油大學儲運與建筑工程學院,山東青島 266580)

LNG冷箱中降溫過程的動態(tài)模擬并行計算

王 琳,李玉星,朱建魯,王婭婷,盛歡歡,王武昌

(中國石油大學儲運與建筑工程學院,山東青島 266580)

由于單核計算無法承擔LNG冷箱中降溫過程動態(tài)模擬的高負荷,采用共享內存模式下的對稱多核并行計算方法實現(xiàn)帶預冷的氮膨脹液化流程中冷箱降溫過程的動態(tài)模擬。動態(tài)模型中,板翅式換熱器采用一維模型,壓縮機、膨脹機和節(jié)流閥按穩(wěn)態(tài)元件處理;并行方法中,根據(jù)并行機特征和液化流程特點進行計算單元的劃分及邊界耦合,利用路障對每個積分步進行同步控制,通過顯式的函數(shù)調用進行單元間通信。將模擬結果與試驗數(shù)據(jù)進行對比,驗證動態(tài)模型和并行方法的合理性。模擬結果表明:工質的降溫速率主要取決于膨脹機多變效率;該并行方法使模擬進程加速23倍,計算結點效率是單核計算的3.83倍。

天然氣液化;冷箱降溫;動態(tài)模擬;多指令流多數(shù)據(jù)流;多核;并行計算

天然氣液化工藝中,天然氣從常溫冷卻到液化溫度,需要通過多次的熱交換,這些換熱過程主要通過冷箱內不同溫度級別的換熱器來完成[1]。帶預冷的氮膨脹液化流程中,冷箱中工質降溫過程的動態(tài)模擬須進行預冷制冷劑、氮氣制冷劑和原料氣3種工質的物性計算,以及單相換熱、沸騰換熱和冷凝換熱3種方式的換熱計算[2],且物性參數(shù)、換熱特性隨時間變化;多股流板翅式換熱器通道排列,傳熱過程復雜,計算繁復[3];而多級串聯(lián)換熱裝置的結構單元多,流程復雜,進一步增加了動態(tài)數(shù)值模擬的復雜性和計算量。由于受計算機CPU頻率和存儲器訪問速度的限制[4],單核計算機系統(tǒng)性能達到極限,無法滿足動態(tài)數(shù)值模擬的需要。并行計算是解決大數(shù)據(jù)量數(shù)值計算耗時過長的重要途徑,其基礎研究在計算數(shù)學和計算機科學領域開展較早[5]。計算流體力學領域的并行計算技術應用研究取得了很多成果,科研人員在計算機集群網(wǎng)絡環(huán)境下,結合并行機性能特征、應用程序執(zhí)行特點和實際物理模型特征,將一批串行計算程序改造為并行計算程序,解決了工程設計中的一些關鍵氣動力問題[6-9]。但是,計算機集群分布式存儲的特點對通信網(wǎng)絡和并行機性能要求很高,且存在通信開銷大的缺點。因此,筆者根據(jù)帶預冷的氮膨脹液化流程中冷箱的結構特點和計算機性能特征采用一種共享內存模式下的對稱多核并行計算方法進行冷箱降溫動態(tài)模擬的多核并行計算,對動態(tài)模擬結果與試驗數(shù)據(jù)進行對比分析,并對并行計算性能進行測試。

1 數(shù)學模型及數(shù)值方法

冷箱降溫動態(tài)模擬過程中采用Peng-Robinson狀態(tài)方程和van der Waals單流體混合規(guī)則進行工質的相平衡及物性計算。

1.1 板翅式換熱器

采用一維模型,并假設:流動穩(wěn)定,換熱器無漏熱,無流動方向的導熱。動量方程中忽略慣性項與對流項,忽略重力沖量,簡化為摩擦力與壓力沖量的關系式。熱流控制方程[2]經(jīng)簡化后如下:

質量方程為

式中,ρh為熱流體密度,kg·m-3;t為時間,s;uh為速度,m·s-1;z為長度,m。

動量方程為

式中,ph為熱流體壓力,Pa;Dh為水力直徑,m;λh為摩阻系數(shù)。

能量方程為

式中,Ah為橫截面積,m2;hh為熱流比焓,J·kg-1;Qh為熱流體換熱量,J·s-1。

同理可得相同形式的冷流動態(tài)模型。

隔板能量方程為

式中,ρw為隔板密度,kg·m-3;Aw為隔板橫截面積, m2;Tw為隔板溫度,K;cp,w為隔板質量定壓熱容,J· kg-1·K-1;Qc為冷流體換熱量,J·s-1。

式(3)、(4)和冷流的能量方程組成了板式換熱器的動態(tài)換熱控制方程,對其進行求解可以得到板翅式換熱器的動態(tài)參數(shù),具體的離散和求解過程見文獻[10]。

1.2 壓縮機和膨脹機

壓縮機和膨脹機的響應時間與換熱器相比要小很多,按穩(wěn)態(tài)元件處理。帶預冷的氮膨脹液化流程中,壓縮機為預冷循環(huán)和氮氣制冷循環(huán)提供動力,膨脹機是主要的制冷元件。壓縮機與膨脹機的熱力學過程相反,原理相似,熱力學模型形式相同。壓縮機特性方程[11]如下:

多變壓縮功為

出口溫度為

多變效率為

式中,Hpol為多變壓頭,J·kg-1;ηpol為多變效率;m為多變指數(shù);R為氣體常數(shù),J·mol-1·K-1;cp為質量定壓熱容;下角標1代表入口狀態(tài),2代表出口狀態(tài)。

1.3 節(jié)流閥

節(jié)流閥的節(jié)流過程近似等焓過程,響應時間短,同樣按穩(wěn)態(tài)元件處理。節(jié)流閥安裝在冷箱的原料氣出口,用于節(jié)流LNG,降低其壓力,便于儲存。節(jié)流閥特性方程[12-13]為

式中,W為質量流量,kg·s-1;Af為閥門過流面積,m2;Kv為流量系數(shù);ζ為閥門開度;C0～C4為常數(shù)。

2 冷箱降溫動態(tài)模擬的并行方法

高效并行算法的研究必須與具體的機器特征、物理過程特點相結合。本文中動態(tài)模擬計算均在DELL Precision Work Station T7500 Tower上進行,處理器為Intel Xeon X5675@3.07 GHz(六核);天然氣(原料氣)降溫主要通過冷箱內的6個串聯(lián)板翅式換熱器完成,流程如圖1所示。高溫原料氣通過6個換熱器逐級降溫至液化溫度,而氮氣制冷劑通過膨脹機對外做功為原料氣液化提供冷量,制冷劑和原料氣通過換熱器進行熱交換。

圖1 冷箱系統(tǒng)工藝流程及計算單元劃分Fig.1 Process flow diagram of cold box and division of computing units

2.1 板翅式換熱器的結構參數(shù)與網(wǎng)格劃分

6個板翅式換熱器的翅片參數(shù)相同:翅片厚度、高度、間距和有效寬度分別為0.20、6.40、1.40 mm 和0.50 m。采用簡化后的一維模型進行板翅式換熱器的動態(tài)傳熱計算。圖2以單元2為例給出了板翅式換熱器的流道分布與求解區(qū)域劃分,在邊界網(wǎng)格上進行信息的接收與發(fā)送。冷箱中板翅式換熱器參數(shù)見表1。

圖2 板翅式換熱器(單元2)流道分布及求解區(qū)域劃分Fig.2 Channel arrangement of plate-fin heat exchanger(unit 2)and grid division

表1 冷箱中板翅式換熱器參數(shù)Table 1 Parameters of plate-fin heat exchangers in cold box

由于板翅式換熱器通道多,結構參數(shù)多,即便采用一維網(wǎng)格進行動態(tài)計算,計算量也相當巨大。如果將6個換熱器串聯(lián)的冷箱降溫過程動態(tài)模擬計算集中在一個程序中進行,計算任務只被分配在1個處理器(計算核心)上,其余5個處理器閑置,多核CPU并不能發(fā)揮其性能,無法提高計算速度。因此,實現(xiàn)冷箱中降溫過程動態(tài)模擬的多核并行計算是提高計算速度的有效途徑。

2.2 計算單元劃分

根據(jù)液化流程的特點將冷箱降溫動態(tài)計算模型劃分為6個計算單元,各計算單元通過接收端和發(fā)送端進行通信。圖1給出了計算單元的劃分方法和物流參數(shù)的通信方式,通過主進程的控制將6個計算單元的計算任務分配到6個處理器上,并進行數(shù)值計算的同步控制和單元間的通信,從而實現(xiàn)流程模擬的分單元計算和邊界耦合。

2.3 并行計算模式

通過計算單元的劃分和計算任務的分配使動態(tài)模擬的并行計算方法與多核CPU相匹配,形成如圖3中多指令流多數(shù)據(jù)流(multiple instruction multiple data,MIMD)的體系結構。每個計算單元執(zhí)行的子進程對應一個處理器結點,各進程具有獨立的堆棧和代碼段,可獨立執(zhí)行。利用同步路障對每個積分步進行同步控制,使進程間的通信同步,進程之間的通信通過共享內存中顯式的函數(shù)調用來完成。

圖3 多核并行計算的體系結構Fig.3 Architecture of multi-core parallel computation

圖4給出了動態(tài)模擬的多核并行計算流程圖。其中,參數(shù)設置:在并行程序的用戶界面上輸入動態(tài)計算的自定義參數(shù)(總步數(shù)N、步長s、模型保存的間隔時間Si、動態(tài)操作參數(shù)等)。

圖4 多核并行計算流程圖Fig.4 Flowchart of multi-core parallel computation

程序運行步驟:

(1)讀取自定義參數(shù);對計算單元Unit1～Unit6賦初值;初始化已計算步數(shù)i和動態(tài)操作編號j。

(2)同步路障計數(shù)器初始化。

(3)判斷已完成計算步數(shù)是否達到用戶設置;若是,則計算結束;若否,則各計算單元同時開始第i步計算。

(4)利用集中式同步路障計數(shù)器對各計算單元的進程進行同步控制,已完成第i步計算的進程阻塞,直至所有進程完成第i步計算。

(5)保存計算單元Unit1～Unit6的第i步計算結果;計算單元Unit1～Unit6的各進程間通信;用已計算步數(shù)i計數(shù)。

(6)判斷是否到達模型保存時間;若是,保存模型;若否,不保存。

(7)判斷是否到達動態(tài)操作時刻;若是,則在計算單元內進行操作參數(shù)設置;若否,則不進行參數(shù)設置。

(8)同步路障計數(shù)器初始化;循環(huán)進入下一步,直至完成所有超步(同步的周期)。

3 計算結果分析

利用小型撬裝天然氣液化裝置進行冷箱降溫過程的試驗研究,試驗系統(tǒng)詳見文獻[14]?？紤]到試驗的安全性,采用氮氣作為原料氣代替天然氣進行液化。試驗過程中的主要操作為氮氣壓縮機的啟動,動態(tài)模擬過程中通過設置冷箱入口的氮氣壓力來模擬。氮氣壓縮機出口壓力變化情況:時刻分別為0、0.25、0.50、0.78、1.39、1.75、2.25、3.50 h時對應的壓力為191.09、230.16、573.87、772.28、843.1、935.29、1 004.88、1 110.5 kPa。

圖5為冷箱降溫過程中原料氣、氮氣制冷劑系統(tǒng)中試驗溫度與模擬值對比。由圖5可見,模擬結果與試驗值在大多數(shù)節(jié)點上吻合較好,驗證了并行計算中應用的數(shù)學模型的合理性,板翅式換熱器等設備的特性參數(shù)設置準確,物性、換熱系數(shù)等參數(shù)計算準確,單元劃分、路障同步等并行方法合理。

R22預冷后的原料氣(圖5(b)中M1)溫度先上升后下降。這是由于初始時刻整個冷箱溫度約為0℃,啟動氮氣壓縮機后,其排氣溫度較高,經(jīng)水冷后溫度下降到20℃,進入冷箱后對換熱器頂部流體有加熱效果,隨著預冷機組和膨脹機的啟動,R22(二氟一氯甲烷)和低溫氮氣又使冷箱頂部的溫度降低,因此M1的溫度呈現(xiàn)先上升后下降的趨勢。

節(jié)流后的原料氣(圖5(b)中M5)試驗溫度與模擬值相異。模擬過程中,根據(jù)壓差計算流量,節(jié)流閥后的壓力設為定值,2.75 h后節(jié)流閥出液,因此,模擬結果中原料氣節(jié)流后溫度從此保持不變。而試驗中由于原料氣壓力不穩(wěn)定,節(jié)流閥后的分離器壓力是變化的,從而造成試驗結果與模擬值的差異。

圖6為通過模擬計算得到的膨脹機多變效率變化曲線。綜合分析可知:冷箱降溫過程中隨著膨脹機溫度降低和入口壓力升高,其多變效率不斷增加;在高、低溫膨脹機的進出口壓力基本一致的情況下,高溫膨脹機(膨脹機1)進出口溫差不斷增加,低溫膨脹機(膨脹機2)溫差基本不變的情況(圖5)與圖6中高溫膨脹機多變效率不斷增加而低溫膨脹機多變效率基本不變的情況相符。上述分析表明:膨脹機的多變效率決定了冷箱降溫過程中物流的溫降速率,而多股流換熱過程的影響較為次要。因此,可以通過提高膨脹機(特別是高溫膨脹機)的多變效率或多變效率的增加速率來優(yōu)化冷箱的降溫速度,但為了防止換熱器的熱應力破壞,應保證降溫速率在30℃/h以內。

圖5 并行計算的氮氣制冷劑及原料氣溫度與試驗值對比Fig.5 Comparison of temperature curves of parallel computation and experiment results of nitrogen refrigerant and natural gas

圖6 模擬過程中膨脹機多變效率變化曲線Fig.6 Variations of polytropic efficiency of expanders in simulation

4 并行計算性能測試

采用時間比率f(n,M)、加速比S(n,M)和結點效率因子e(n,M)對并行計算性能進行評價,定義分別為

式中,tphy為物理過程的實際時間;tsim為模擬計算用時;n為計算結點(核心)數(shù)目;M為計算規(guī)模,用網(wǎng)格數(shù)表示;f(n,M)為n個核心并行計算的時間比率;f(1,M)為最優(yōu)單核串行計算的時間比率。

時間比率是實際物理現(xiàn)象時間與模擬計算時間的比值,反映計算速率的快慢;加速比表示并行計算與串行計算時間比率的比值,反映并行計算的加速倍率;結點效率因子e(n,M)表示并行計算加速比與計算結點數(shù)的比值,反映并行計算中計算結點效率提高的倍數(shù)。

在冷箱降溫動態(tài)模擬計算過程中采用單核串行計算對5 h的降溫過程進行模擬,最快需要7 d;而采用六核并行計算7.25 h內即可完成。并行計算加速比為23,計算核心的平均效率提高為原來的3.83倍。表2為單核串行計算和多核并行計算的性能測試結果對比。

表2 串、并行計算性能對比Table 2 Performance comparison of parallel and serial computation

5 結論及應用前景

(1)模擬結果與試驗數(shù)據(jù)吻合程度較好,驗證了動態(tài)模型的準確性和并行方法的合理性。

(2)LNG冷箱中工質的降溫速率主要取決于膨脹機的多變效率,多股流換熱過程的影響較為次要。

(3)并行計算使模擬進程加速23倍,計算結點效率是單核計算的3.83倍,并行計算性能優(yōu)越。

(4)分單元計算、邊界耦合的并行方法是解決大數(shù)據(jù)量動態(tài)模擬計算的有效途徑,而高效并行計算技術的研發(fā)必須與并行機性能、應用程序執(zhí)行特點和實際物理模型特征相結合。

本文的并行方法不拘泥于單機多核環(huán)境,可以擴展到計算機集群環(huán)境下進行并行計算,但應注意異構環(huán)境中的負載平衡和通信開銷問題。

[1] 顧安忠.液化天然氣技術手冊[M].北京:機械工業(yè)出版社,2010:474-484.

[2] 王坤,徐風雨,李紅艷,等.小型MRC天然氣液化裝置中板翅式換熱器動態(tài)特性仿真研究[J].低溫工程, 2007,157(3):44-49.

WANG Kun,XU Fengyu,LI Hongyan,et al.Dynamic performance simulation of plate-fin heat exchangers in small scale MRC-LNG plant[J].Cryogenics,2007,157 (3):44-49.

[3] 王松漢.板翅式換熱器的電子計算機計算[J].石化技術,1980(3):367-380.

WANG Songhan.The electronic computer calculation of plate-fin heat exchanger[J].Petrochemical Industry Technology,1980(3):367-380.

[4] 周偉明.多核計算與程序設計[M].武漢:華中科技大學出版社,2009.

[5] 胡峰,胡保生.并行計算技術與并行算法綜述[J].電腦與信息技術,1999(5):47-59.

HU Feng,HU Baosheng.Parallel computing technology and parallel algorithms review[J].Computer and Information Technology,1999(5):47-59.

[6] 李津,李忠澤,朱自強,等.分區(qū)并行跨聲速流的計算[J].空氣動力學學報,1998,16(3):374-378.

LI Jin,LI Zhongze,ZHU Ziqiang,et al.The parallel computation of transonic flow[J].Acta Aerodynamica Sinica,1998,16(3):374-378.

[7] 朱自強,麻曉波,付鴻雁,等.流場分析和設計的并行計算[J].科學技術與工程,2003,3(6):582-587.

ZHU Ziqiang,MA Xiaobo,FU Hongyan,et al.Parallel computation of flow field analysis and design[J].Science Technology and Engineering,2003,3(6):582-587.

[8] 周乃春,鄧有奇,楊其德.戰(zhàn)術導彈繞流流場并行計算[J].空氣動力學學報,2002,20(sup l):64-69.

ZHOU Naichun,DENG Youqi,YANG Qide.Parallel computation of tactical missile[J].Acta Aerodynamica Sinica,2002,20(sup l):64-69.

[9] 溫建華,朱自強,吳宗成,等.基于N-S方程串并行計算的機翼優(yōu)化設計[J].北京航空航天大學學報, 2008,34(2):127-130.

WEN Jianhua,ZHU Ziqiang,WU Zongcheng,et al. Wing's optimization design using serial and parallel computations based on N-S equations[J].Journal of Beijing University of Aeronautics and Astronautics,2008,34(2): 127-130.

[10] 徐益峰,蔡祖恢.平行流多流體板翅式換熱器的動態(tài)數(shù)學模型[J].化工學報,1998,49(6):721-728.

XU Yifeng,CAI Zuhui.Research on dynamic model of multi-fluid plate-fin heat exchangers[J].Journal of Chemical Industry and Engineering,1998,49(6):721-728.

[11] 錢錫俊,陳弘.泵和壓縮機[M].東營:中國石油大學出版社,2007:127-131.

[12] 柯明純.液壓馬達和電液比例節(jié)流閥性能分析與測試的研究[D].杭州:浙江大學機械與能源工程學院,2006.

KE Mingchun.Study on performance analysis and testing of hydraulic motor and electro-hydraulic proportional throttle[D].Hangzhou:The Institute of Mechatronic Control Engineering of Zhejiang University,2006.

[13] 丁延鵬.天然氣管輸系統(tǒng)快瞬變流動特性研究[D].青島:中國石油大學儲運與建筑工程學院,2011.

DING Yanpeng.Study on the characteristics of fast transients in natural gas transmission system[D].Qingdao: College of Pipeline and Civil Engineering in China University of Petroleum,2011.

[14] 朱建魯,李玉星,王武昌,等.晃蕩條件下氮膨脹液化過程冷箱運行可靠性試驗[J].化工學報,2013,64 (4):1183-1190.

ZHU Jianlu,LI Yuxing,WANG Wuchang,et al.Reliability experiments in a cold box with nitrogen expansion liquefaction process running under sloshing conditions [J].CIESC Journal,2013,64(4):1183-1190.

(編輯 沈玉英)

Parallel computation of dynamic simulation of cool down process in LNG cold box

WANG Lin,LI Yuxing,ZHU Jianlu,WANG Yating,SHENG Huanhuan,WANG Wuchang
(College of Pipeline and Civil Engineering in China University of Petroleum,Qingdao 266580,China)

As the single-core can't afford the high load of dynamic simulation of the cool down process in LNG cold box,a symmetric multi-core parallel computing method under shared memory was proposed to implement the dynamic simulation of the cool down process in cold box of nitrogen expansion liquefaction process with pre-cooling.In dynamic mathematical models,one-dimensional model is used in the simulation of plate-fin heat exchanger,and compressor,expander and throttle were treated as steady-state operations.In parallel methods,compute units were divided and coupled in the boundary based on the characteristics of parallel machine and liquefaction process.The simulation processes of parallel units were controlled using synchronization barrier,and the data were communicated through explicit function calls.The comparison of simulation results with experimental data shows that the dynamic models and parallel methods were reasonable.The simulation results show that the cooling rate of the working fluid depends mainly on the polytropic efficiency of the expanders.The result of performance test of parallel computation shows that the parallel method accelerates the simulation process 23 times,and the efficiency of computing nodes is 3.83 times the efficiency of the single-core computing.

natural gas liquefaction;cool down process in cold box;dynamic simulation;multiple instruction multiple data (MIMD);multi-core;parallel computation

TE 646

A

1673-5005(2014)04-0148-06

10.3969/j.issn.1673-5005.2014.04.022

2013-11-23

國家重大科技專項(2011ZX05026-006-07);中央高?；究蒲袠I(yè)務費專項(13CX06064A)

王琳(1986-),男,博士研究生,主要從事多相管流及油氣田集輸技術研究。E-mail:lincw_wang@qq.com。

王琳,李玉星,朱建魯,等.LNG冷箱中降溫過程的動態(tài)模擬并行計算[J].中國石油大學學報:自然科學版,2014,38(4):148-153.

WANG Lin,LI Yuxing,ZHU Jianlu,et al.Parallel computation of dynamic simulation of cool down process in LNG cold box [J].Journal of China University of Petroleum(Edition of Natural Science),2014,38(4):148-153.